1. 서 론

탄도탄 기술 고도화에 따라 점증하는 안보위협에 효과적으로 대처하기 위해, 지난 수십 년 간 탄도탄 방어체계(BMDS: ballistic missile defense system)에 관한 연구가 활발히 수행되어 왔다. BMDS를 개발, 운용하기 위해서는 잠재적 위험성을 지닌 표적을 조기에 탐지/식별한 후, 요격탄에 예상 궤적을 제공하는 고성능 추적시스템의 개발이 전제되어야 한다. 불행하게도 탄도탄 표적은 항공기 표적에 비해 매우 빠른 속력으로 비행할 뿐만 아니라, 대기권 재진입 시 항력(drag force)에 의한 기동특성이 크게 변화하므로 궤적추정 및 예측의 정확도를 담보하기가 쉽지 않다. 더욱이, 의도적인 기동을 수행하여 BMDS를 무력화하는 신형 탄도탄의 등장은 표적추적 필터링 문제의 복잡도를 가중시키는 요인으로 작용하고 있다. 이러한 도전과제를 효과적으로 해결하기 위해서는 표적추적 필터 설계에 앞서 탄도탄 표적에 적합한 기동모델링 방법에 관한 기술적 검토가 이루어져야 하지만, 아직까지 이와 관련한 연구는 매우 미진한 상황이다.

기존 탄도탄 표적추적 기법에 관한 연구들은 주로 시스템 모델의 비선형성으로 인한 성능저하를 억제하는데 초점을 맞추고 있다. 대표적인 초창기 연구결과 중 하나인 RPCC (radar principal Cartesian coordinate) 기법은 추적시스템을 구성하고 있는 능동위상배열 레이다의 신호처리 특성에 대한 이해를 바탕으로, 기존 직교좌표계 위치 대신 상대거리 및 방향 코사인을 상태변수로 하는 추적필터 설계 방법이다. 이 경우, 레이다 측정모델이 선형 측정방정식으로 정의되므로, 레이다 측정잡음의 상호 상관성에 의한 EKF(extended Kalman filter)의 성능저하 문제를 해결할 수 있는 것으로 알려져 있다. 하지만, 변침기동 탄도탄의 경우와 같이 표적이 복잡한 기동패턴을 보일 경우에는 도리어 운동방정식의 비선형성이 크게 강화되어, RPCC 기법 적용에 따른 효과가 퇴색될 수 있다 ^(1,2).

탄도탄 표적추적 필터 설계를 위한 시스템 모델의 비선형성에 의한 성능저하를 완화하기 위해 다양한 형태의 비선형 상태추정 기법이 제안되었다 ^(3-5). IEKF(iterative extended Kalman filter), IKF(interval Kalman filter) 기법은 탄도탄 표적의 비선형 운동방정식을 선형화하는 과정에서 발생하는 오차를 감소시키기 위해 고안되었다 ⁽³⁾. 하지만, EKF가 지닌 구조적 문제, 즉 측정치와 필터이득 간의 상관성으로 인한 추정치 편향 문제를 여전히 안고 있다. 그 대안으로 비선형 확률밀도 함수를 직접 근사하는 UKF(unscented Kalman filter), CKF (cubature Kalman filter), PF(particle filter) 등의 적용이 시도되었다. 이 방법들은 모두 탄도탄 표적의 미지 항력특성, 물성치 및 공기밀도 등의 함수로 기술되는 탄도계수(ballistic coefficient)를 추가 상태변수로 사용하여 재진입 구간에서 표적의 감속기동에 대응하고자 하였다. 하지만, 탄도계수 불확실성에 의해 추적필터의 성능이 저하되거나, 탄도탄이 변침 기동하는 최악의 상황에서는 궤적추정치가 발산하는 등, 실제 적용에 한계를 노출한 바 있다 ^(4,5).

이러한 이유로, 일각에서 추적필터의 잔류오차를 이용하여 탄도탄 표적의 기동 상황에 대응하는 적응 표적추적 필터링 기법이 연구되었다 ^(6,7). 이 접근법은 표적 기동 보상 방법에 따라 공정잡음 적응형 필터링 기법과 입력 추정기법으로 나뉜다. 첫 번째 방법의 핵심 아이디어는 탄도탄 표적의 감속기동 모델에 포함된 공기밀도의 불확실성을 공정잡음으로 반영하자는 것이다. 이 방법을 사용할 경우, 표적의 감속기동으로 인해 잔류오차가 커지면 공정잡음 분산을 적절히 증가시켜 추적필터의 불안정성을 손쉽게 억제할 수는 있지만, 궤적 추정치의 분산이 커지는 단점이 있다. 이와 달리, 두 번째 방법은 일정한 시구간의 잔류오차 데이터를 분석하여 표적 기동을 탐지함과 동시에 표적기동에 따른 추정오차를 보상한다. 하지만, 이 방법은 탄도탄 표적의 급기동에 빠르게 대처하기 어렵고 시스템 모델의 비선형성으로 인해 추정오차 보상에 필요한 기동시점과 크기의 해석적으로 산출하는 과정이 상당히 복잡해 실제 적용이 용이하지 않다.

앞서 소개한 선행 연구결과들은 탄도탄 표적추적 문제의 본질이 필터링 이론 보다는 표적 감속/변침을 모두 아우를 수 있는 일반화된 탄도탄 기동모델의 수립에 있음을 암시한다. 이러한 통찰로부터 소위 MaRV(maneuvering re-entry vehicle) 필터가 고안되었다 ^(8,9). 이 기법은 기존 기법들과 달리, 탄도탄에 작용하는 종/횡 가속도를 각각 별도의 탄도계수로 모델링하여 탄도탄 표적 운동방정식에 반영하였다. MaRV 필터는 탄도탄 표적이 임의의 방향으로 기동하더라도 비교적 안정적인 트랙 유지 성능을 제공한다. 하지만, 탄도탄 표적의 공력특성이 탄도계수에 명시적으로 반영되지 않아 급기동 시점에서 추정치 편차가 커지거나 필터의 수렴속도가 늦어지는 단점이 있다. 또한, 추정된 탄도계수가 추적필터의 발산을 억제할 뿐 물리적 의미를 갖지 못하는 경우가 많아, 이를 활용한 체계적인 표적 식별/분류가 어렵다 ⁽⁹⁾.

본 논문에서는 탄도탄 표적의 공력특성을 사전지식으로 활용하는 명시적 탄도탄 기동모델에 기반한 추적필터링 방법을 현실적 대안으로 제시한다. 사전에 표적 분류를 통해 해당 탄급의 개략적인 공력특성을 활용할 수 있다는 가정 하에 양력 및 항력과 관련된 탄도계수들 간의 구속조건을 명시적으로 도출하고, 이를 탄도탄 표적 기동 모델링에 반영한다. 이 경우, 탄도탄 표적의 운동방정식은 3축 관성위치/속도와 양력에 관한 2개의 탄도계수를 포함한 총 8개의 상태변수로 기술된다. 제안한 모델링 기법을 적용할 경우, 기존 MaRV 기법에 비해 필터 차수를 낮출 수 있을 뿐만 아니라 운동방정식 구성 시

그림. 1. 레이다-표적 상대기하

Fig. 1. Relative target geometry with respect to the radar

기동모델의 불확실성이 줄어들어 기동시점에서 표적추적 필터의 수렴속도와 안정성을 동시에 개선할 수 있다. 실제 상황에서 탄도탄 표적의 공력특성에 관한 사전지식이 불완전하므로, 사전지식에 대한 민감도 분석을 수행한다. 이를 통해, 명시적 기동모델이 완벽하지 않더라도 추적성능에 긍정적 영향을 끼침을 확인한다. 제안한 기법의 타당성을 입증하기 위해 재진입 변침기동 궤적에 대한 모의실험이 수행된다.

2. 기동 탄도탄 추적을 위한 사전지식 모델링

본 장에서는 기동 탄도탄의 운동특성을 모델링한다. 재진입 구간에서 탄도탄 표적에 작용하는 구심력 및 코리올리 힘에 의한 가속도는 항력 및 양력에 비해 매우 작으므로 평면 지구를 가정하고, 지구 자전 효과는 무시한다.

2.1 좌표계 및 운동특성

레이다에 대한 탄도탄의 상대기하는 그림 1에 도시한 바와 같다. 그림에서 사용된 변수들의 정의는 다음과 같다.

$[x,\:y,\:z]^{T}$ : 탄도탄 표적 위치(I-frame)

$[\dot x ,\:\dot y ,\:\dot z]^{T}$ : 탄도탄 표적의 속도(I-frame)

$r,\:\lambda_{h},\:\lambda_{v}$ : 레이다-표적 간 상대거리 및 방위각/고각(I-frame)

탄도탄 상대운동 모델링을 위해 도입된 좌표계의 정의는 다음과 같다.

· 관성좌표계 (I-frame)

탄도탄의 관성운동을 기술하기 위한 기준 좌표계로, 편의상 레이다를 원점으로 하고 $X_{I}$축이 정북방향, $Y_{I}$축이 정동방향, 그리고 $Z_{I}$축이 지구 타원체에 접하는 평면의 수직 하 방향으로 정의되는 NED(North-East-Down) 좌표계로 정의한다.

· 시선좌표계 (L-frame)

그림. 2. 받음각 정의

Fig. 2. Definition of the angle of attack

원점이 레이다 위치와 일치하고, $X_{L}$ 축은 레이다와 표적을 잇는 시선벡터, $Y_{L}$ 축은 I-frame 수평면 상에 위치하는 오른손 좌표계이다.

· 표적속도좌표계 (V-frame)

탄도탄 무게 중심을 원점으로 하고, $X_{V}$ 축은 탄도탄 속도방향, $Y_{V}$ 축은 I-frame 수평면 상에 위치하는 오른손 좌표계이다. 탄도탄 표적 속력이 $v=\sqrt{\dot x^{2}+\dot y^{2}+\dot z^{2}}$, 지면속력이 $v_{g}=\sqrt{\dot x^{2}+\dot y^{2}}$라면 V-frame에서 I-frame으로의 좌표변환 행렬 $C_{V}^{I}$은 다음과 같이 정의된다.

$C_{V}^{I}=\left[\begin{matrix}\dot x /v&-\dot y /v_{g}&-\dot x\dot z /(vv_{g})\\\dot y /v&\dot x /v_{g}&-\dot y\dot z /(vv_{g})\\\dot z /v&0&v_{g}/v\end{matrix}\right]$

· 표적동체좌표계 (B-frame)

원점이 탄도탄 중심과 일치하고, $X_{B}$축은 탄도탄 코 방향, $Y_{B}$축은 탄도탄 우측 날개방향과 일치하는 오른손좌표계이다.

그림 1 및 그림 2의 좌표계 및 기호정의에 따라 공력 및 중력에 의한 탄도탄 표적의 가속도를 다음과 같이 쓸 수 있다 ⁽⁸⁾.

(1)

\begin{align*} f_{a}(v)=\dfrac{d}{dt}\left[\begin{aligned}\begin{aligned}\dot x \\\dot y\end{aligned}\\\dot z\end{aligned}\right]=Q(h,\:v)C_{V}^{I}\left[\begin{aligned}\begin{aligned}-\alpha_{D}\\-\alpha_{T}\end{aligned}\\-\alpha_{C}\end{aligned}\right]+\left[\begin{aligned}\begin{aligned}0\\0\end{aligned}\\g(h)\end{aligned}\right] \end{align*}

여기서 $v =[\dot x \dot y \dot z]^{T}$이며, $\rho$와 $g$는 각각 고도 $h$에 따른 공기 밀도와 중력, $Q(h,\:v)=0.5\rho(h)v^{2}$는 동압(dynamic pressure), $(\alpha_{D},\:\alpha_{T},\:\alpha_{C})$은 표적의 변침기동을 모델링하기 위해 도입된 일반화된 탄도계수(generalized ballistic coefficients)를 의미한다. 참고로 탄도계수에 동압을 곱한 값이 가속도 물리량에 대응한다.

식 (1)에서 탄도계수 $(\alpha_{D,\:}\alpha_{T,\:}\alpha_{C})$는 탄도탄 표적의 무차원 공력계수 $(C_{D,\:}C_{L})$과 다음 관계를 갖는다.

$\begin{bmatrix}\alpha_{D}\\\alpha_{T}\\\alpha_{C}\end{bmatrix}=\dfrac{S_{ref}}{m}\begin{bmatrix}C_{D}(h,\: v,\:\alpha_{t})\\C_{L}(h,\: v,\:\beta)\\C_{L}(h,\: v,\:\alpha)\end{bmatrix}$

여기서 $m$과 $S_{ref}$는 표적 질량 및 단면적, $\alpha$와 $\beta$는 수직/수평 받음각, $\alpha_{t}\approx\sqrt{\alpha^{2}+\beta^{2}}$는 총 받음각, $C_{D}$와 $C_{L}$은 표적 형상에 의해 결정되는 항력(drag) 및 양력(lift) 계수이다. 참고로, 탄도계수의 부호에 따른 기동 방향은 다음과 같다.

부호	$\alpha_{D}$	$\alpha_{T}$	$\alpha_{C}$
(+)	감속	우선회	상승
(-)	가속	좌선회	하강

2.2 공력 사전지식 모델링

앞서 살펴본 바와 같이 탄도탄 표적의 기동은 표적 형상에 의해 결정되는 공력계수로 기술된다. 따라서, 공력특성 사전지식의 활용정도에 따라 탄도탄 표적 기동 모델링 방법이 달라진다. 본 논문에서는 공력특성 사전지식을 거의 사용하지 않는 방법을 암시적 기동모델로, 사전지식을 적극적으로 사용하는 방법을 명시적 기동모델로 명명한다.

· 암시적 탄도탄 기동모델

변침기동 탄도탄 기동 모델링은 유도항력(lift induced drag)의 반영 방법을 찾는 문제로 귀결된다. 횡기동에 의한 양력은 탄도탄 표적의 항력을 증가시키므로, 재진입 단계에서 탄도궤적을 그리는 통상적인 탄도탄 표적과는 매우 다른 형태의 감속기동을 야기한다. 만일 이러한 상황에서, 전통적인 방법들에서와 같이 항력에 의한 탄도계수만을 고려하여 추적필터를 구성하면 표적추적 성능 저하, 더 나아가 최악의 경우 표적추적 실패에 이를 수 있다 ^(7,8).

대부분의 탄도탄 표적은 십자 대칭형상(cruciform)을 가지므로, 유도항력을 다음과 같이 기술할 수 있다 ⁽¹⁰⁾.

(2)

$C_{D}=C_{D_{0}}+KC_{L}^{n}$

여기서 $C_{D_{0}}$는 양력이 없을 경우의 항력(zero-lift drag), $K$는 유도항력 계수(induced drag parameter), $n$은 표적 형상에 따라 달라지는 양항곡선의 지수 항으로 초음속 탄도탄의 경우 $n\approx 2$를 만족한다 ⁽⁹⁾. 이 경우, 양항비 $LD$는 다음 식과 같이 표현된다.

(3)

$LD=C_{L}/C_{D}=C_{L}/(C_{D_{0}}+KC_{L}^{2})$

암시적 기동모델에서는 탄도탄 표적이 항상 최대 양항비로 기동한다고 가정한다. 임계 양력계수(critical lift coefficient) $C_{L}^{*}$에 대한 $K$, 즉 식 (3)을 $C_{L}$에 대하여 미분한 값이 0이 되도록 하는 $K$는 다음과 같다.

(4)

$\dfrac{\partial LD}{\partial C_{L}}=0\to C_{L}^{*},\: K^{*}=\dfrac{C_{D_{0}}}{(n-1)(C_{L}^{*})^{n}}\approx\dfrac{C_{D_{0}}}{(C_{L}^{*})^{2}}$

식 (4)를 식 (3)에 대입하면, 최대 양항비 $LD_{\max}$를 얻는다.

(5)

$LD_{\max}=\dfrac{1}{2}\left(\dfrac{C_{L}^{\ast}}{C_{D_{0}}}\right)$

위의 결과를 이용하면, 양력 및 항력계수를 $\lambda =(C_{L}/C_{L}^{*})$로 다시 쓸 수 있다.

(6)

$C_{L}=2C_{D_{0}}LD_{\max}\lambda$, $C_{D}=C_{D_{0}}(1+\lambda^{2})$

식 (6)으로부터 항력 및 양력에 의한 표적의 3차원 기동을 탄도계수 간의 상관관계로 기술할 수 있다.

(7)

$\alpha_{D}=\delta_{D}(1+\lambda^{2})$, $\alpha_{L}=\lambda\delta_{L}$

여기서 탄도계수 $\delta_{D}$와 $\delta_{L}$의 정의는 다음과 같다.

$\delta_{D}=\dfrac{S_{ref}}{m}C_{D_{0}}$, $\delta_{L}=\dfrac{S_{ref}}{m}2C_{D_{0}}LD_{\max}$

양력에 의한 가속도는 속도벡터의 법평면에 존재하므로, $\alpha_{L}$을 수평방향 성분 $\alpha_{C}$와 수직방향 성분 $\alpha_{T}$로 분기할 수 있다.

(8)

$\alpha_{T}=\lambda\delta_{T},\:\alpha_{C}=\lambda\delta_{C}$

여기서 $\alpha_{L}=\sqrt[]{\alpha_{T}^{2}+\alpha_{C}^{2}}$ 혹은 $\delta_{L}=\sqrt[]{\delta_{T}^{2}+\delta_{C}^{2}}$이 만족된다.

Remark 2.1. 식 (7) 및 식 (8)로부터 탄도탄 표적의 기동 가속도를 모델링하기 위한 3개의 탄도계수 $(\alpha_{D},\:\alpha_{T},\:\alpha_{C})$는 4개의 미지 공력특성 파라미터 $(\lambda$,$\delta_{D},\:$$\delta_{T}$,$\delta_{C})$의 함수임을 알 수 있다. 만일, 공력특성 파라미터를 직접 추정할 수 있다면, 탄도탄 표적 궤적추정 뿐만 아니라 서로 다른 공력특성, 즉 형상이 다른 표적의 분류(classification)/식별(recognition)도 시도해 볼 수 있다 ⁽⁹⁾. 하지만, 이들 간의 매핑 관계를 유일하게 결정할 방도가 없으므로, 실제로는 공력특성 파라미터에 대한 추적필터의

그림. 3. 양항곡선 사전지식 예

Fig. 3. An example of a model of the drag polar a priori information

가관측성이 확보되지 않는다. 다시 말해, 공력특성 파라미터 $(\lambda$,$\delta_{D},\:$$\delta_{T}$,$\delta_{C})$를 상태변수로 설정하는 경우 물리적 의미를 갖는 추정치를 얻을 수 없으며, 이에 따라 표적 분류/식별 가능성도 사라진다.

이상의 관찰로부터 암시적 탄도탄 기동모델로는 추적필터 설계 시 양항 구속조건 (3)을 제대로 고려하기 어렵다는 결론을 얻는다. 차선책 중 하나는 가관측성이 결여된 공력특성 파라미터 $(\lambda$,$\delta_{D},\:$$\delta_{T}$,$\delta_{C})$ 대신, 식 (9)의 탄도계수 $(\alpha_{D},\:\alpha_{T},\:\alpha_{C})$를 탄도탄 표적의 상태변수로 간주하여 추정하는 것이다. 다만, 이 경우에는 항력-양력이 독립적으로 표적 동체에 작용한다고 가정하고 탄도탄 운동방정식을 구성해야 하므로, 다소 간의 성능저하를 감수해야 한다. 따라서, 암시적 기동모델에 기반한 추적필터의 상태변수는 $9$차원 벡터로 정의된다.

(9)

$x_{k}=\left[x,\: y,\: z,\:\dot x ,\:\dot y ,\:\dot z ,\:\alpha_{D},\:\alpha_{T},\:\alpha_{L}\right]^{T}$

탄도계수의 변화 특성은 위치/속력에 비해 매우 느리게 변화하므로, 탄도탄 운동모델은 다음과 같이 기술된다.

(10)

\begin{align*} \dfrac{d}{dt}x = f(x)=\left[\begin{aligned}\begin{aligned}\begin{aligned}\begin{aligned}\begin{aligned}\begin{aligned}\dot x \\\dot y\end{aligned}\\\dot z\end{aligned}\\Q(h,\:v)\left(-\dfrac{\dot x}{v}\alpha_{D}+\dfrac{\dot y}{v_{g}}\alpha_{T}+\dfrac{\dot x\dot z}{vv_{g}}\alpha_{C}\right)\end{aligned}\\Q(h,\:v)\left(-\dfrac{\dot y}{v}\alpha_{D}-\dfrac{\dot x}{v_{g}}\alpha_{T}+\dfrac{\dot y\dot z}{vv_{g}}\alpha_{C}\right)\end{aligned}\\Q(h,\:v)\left(-\dfrac{\dot z}{v}\alpha_{D}-\dfrac{v_{g}}{v}\alpha_{C}\right)+g(h)\end{aligned}\\\begin{aligned}\begin{aligned}0\\0\end{aligned}\\0\end{aligned}\end{aligned}\right] \end{align*}

· 명시적 탄도탄 기동모델

식 (2)에 기술된 바와 같이 탄도탄 표적의 양력과 항력에 의한 기동 특성은 $C_{D_{0}}$와 $K$에 의해 결정된다. 이때, $C_{D_{0}}$는 초음속 충격파에 의해 발생하는 조파항력(wave drag), 탄도탄 코와 꼬리의 기압차에 의해 발생하는 기저항력(base drag), 표면 점성마찰에 의해 결정되는 마찰항력(friction drag) 등으로 구성된다. 이들 성분들은 탄두부 형상, 탄 길이 및 지름 등, 개략적인 탄도탄 표적 형상정보로 근사 가능하다 ⁽¹¹⁾. 대칭 형상 탄도탄의 $K$ 역시 마찬가지로 표적 형상의 함수로 기술된다 ⁽¹⁰⁾. 따라서, 추적대상 탄도탄 표적 군(群)이 주어져 있다면, 알려진 형상정보를 이용하여 개략적으로 산출된 $C_{D_{0}}$및 $K$를 사전지식화 할 수 있다. 일례로 그림 3은 임의의 탄도탄 형상에 대한 양항곡선을 식 (2)의 형태로 근사한 것이다. 이를 고도-속력를 변화시키며 사전지식화 한 결과는 그림 4와 5에 도시한 바와 같다. 그림에서 볼 수 있듯, $C_{D_{0}}$의 경우 고도 $30[km]$ 미만, 속력 $1,\:000[m/s]$미만에서 유의미한 차이가 발생하지만, $K$의 경우 고도 및 속력에 따른 변화 폭이 크지 않아 상수 $\bar{K}$로 근사해도 무방하다. 이 경우, 식 (2)는 다음과 같이 간략화 된다.

(11)

$C_{D}=C_{D_{0}}(h,\:v)+\bar{K}C_{L}^{2}$

식 (11)을 탄도계수의 정의를 이용해 다시 쓰면

(12)

$\alpha_{D}=\alpha_{D_{0}}+\bar{k}\left(\alpha_{T}^{2}+\alpha_{C}^{2}\right)$

여기서 $\alpha_{D_{0}}=\dfrac{S_{ref}}{m}C_{D_{0}}$, $\bar{k}=\dfrac{m}{S_{ref}}\bar{K}$를 의미한다. 식 (12)는 상태변수로 사용되는 탄도계수들 사이의 구속조건(state equality constraint)이다. 따라서, 이를 변침기동 탄도탄 표적추적 필터 설계에 활용할 경우 상태추정 성능 향상을 도모할 수 있다. 만일, 상수 $\bar{k}$와 더불어 $\alpha_{D_{0}}$까지 사전지식화 할 수 있다면 추정해야 할 탄도계수의 개수는 $(\alpha_{D},\:\alpha_{T},\:\alpha_{C})$ 3개에서 $(\alpha_{T},\:\alpha_{C})$ 2개로 감소하며, 이에 따라 추적필터의 상태변수는 다음과 같이 8차원의 벡터가 된다.

(13)

$x =\left[x,\: y,\: z,\:\dot x ,\:\dot y ,\:\dot z ,\:\alpha_{T},\:\alpha_{C}\right]^{T}$

(14)

\begin{align*} \dfrac{d}{dt}x = f(x)=\left[\begin{aligned}\begin{aligned}\begin{aligned}\begin{aligned}\begin{aligned}\begin{aligned}\dot x \\\dot y\end{aligned}\\\dot z\end{aligned}\\Q(h,\:v)\left(-\dfrac{\dot x}{v}\alpha_{D}+\dfrac{\dot y}{v_{g}}\alpha_{T}+\dfrac{\dot x\dot z}{vv_{g}}\alpha_{C}\right)\end{aligned}\\Q(h,\:v)\left(-\dfrac{\dot y}{v}\alpha_{D}-\dfrac{\dot x}{v_{g}}\alpha_{T}+\dfrac{\dot y\dot z}{vv_{g}}\alpha_{C}\right)\end{aligned}\\Q(h,\:v)\left(-\dfrac{\dot z}{v}\alpha_{D}-\dfrac{v_{g}}{v}\alpha_{C}\right)+g(h)\end{aligned}\\\begin{aligned}0\\0\end{aligned}\end{aligned}\right] \end{align*}

그림. 4. 고도 및 속력에 따른 $C_{D_{0}}$

Fig. 4. $C_{D_{0}}$ with respect to altitude and absolute velocity

그림. 5. 고도 및 속력에 따른 $K$

Fig. 5. $K$ with respect to altitude and absolute velocity

Remark 2.2. 해외 BMDS는 탄도탄 표적의 상승단계 혹은 외기권을 비행하는 중기단계에서 위성센서 혹은 조기경보레이다 에서 획득된 신호정보를 이용한 탄급(彈級) 분류(classification)를 시도하는 것으로 알려져 있다 ⁽¹²⁾. 탄도탄 표적의 클래스 정보가 제공되면 해당 클래스의 사전지식 $\alpha_{D_{0}},\:\bar{k}$로 정의된 명시적 기동모델 (12) 혹은 운동방정식 (14)를 사용하여 추적필터를 설계할 수 있다.

Remark 2.3. 명시적 기동모델의 사용은 암시적 기동모델에 비해 표적 탄종 식별(identification) 측면에서 여러 장점을 갖는다. 암시적 기동모델과 달리, 명시적 기동모델을 표적추적에 사용하면 필터의 잔류오차(residue) 분석을 통해 표적이 어떤 공력특성을 지니고 있는지 판단할 수 있으므로, 이를 표적식별에 활용할 수 있다. 실제 상황에서 상존하는 사전지식의 불완전성으로 인해 표적식별 성능이 다소 저하될 수 있으나, 교전 레이다에서 제공되는 RCS(radar cross-section), HRRP(high resolution range profile) 등 표적 신호정보와의 정보융합을 통해 이를 보완할 수 있다 ⁽¹³⁾.

3. 기동 탄도탄 표적추적 필터 설계

식 (1)의 가속도 운동 모델 $f_{a}(·)$에 기반하여, 기동 탄도탄 추적을 위한 이산 시스템 모델은 다음과 같이 정의된다.

여기서 $x_{k}$는 상태변수, $T$는 필터 샘플링 주기, $w_{k}$는 모델링 오차를 반영하기 위한 공정잡음으로 분산이 $Q_{k}$인 영평균 백색잡음으로 가정한다. $N_{\alpha}$를 탄도계수 상태변수의 차원수라 하면, 위의 식에 사용된 행렬의 정의는 다음과 같다.

$F_{0}=\begin{bmatrix}I^{3\times 3}&T· I^{3\times 3}&0^{3\times N_{\alpha}}\\0^{3\times 3}&I^{3\times 3}&0^{3\times N_{\alpha}}\\0^{3\times 3}&0^{3\times 3}&I^{3\times N_{\alpha}}\end{bmatrix}$,\begin{align*} G_{a}=\left[\begin{aligned}\begin{aligned}0^{3\times 3}\\T· I^{3\times 3}\end{aligned}\\0^{3\times N_{\alpha}}\end{aligned}\right] \end{align*},$G_{w}=\left[\begin{aligned}0^{3\times(3+N_{\alpha})}\\I^{6\times(3+N_{\alpha})}\end{aligned}\right]$

탄도계수의 변화 특성은 위치, 속도 변화에 비해 무시할 만하므로 white-jerk로 간주된다.

탄도탄 표적추적과 관련한 기존 연구결과들에 따르면, 재진입 구간에서 항력 및 양력에 의한 가속도의 크기가 증가할수록 비선형 시스템 방정식 (15)의 자코비안의 행렬놈(matrix norm)이 커져 $w_{k}$의 전파가 매우 크게 일어난다 ⁽⁷⁾. 이러한 문제는 표적의 수평속력이 줄어들수록 더욱 두드러진다. 따라서, 수평방향 속력 및 기동 가속도의 크기에 따라 공정잡음 분산 $Q_{k}$의 체계적 선정이 필요하다. 이를 위해, 본 논문에서는 선형화 된 시스템 모델 $F_{k}$를 활용하여 매 시점 공정잡음 분산 $Q_{k}$를 근사적으로 계산하는 방법을 사용한다.

여기서 $Q$는 표적가속도 및 탄도계수 변화율의 연속시간 공정잡음 공분산을 의미하며, 어떤 기동모델을 필터에 사용하는지에 따라 다음과 같이 정의된다.

$Q$ $=$ $diag([q_{a_{x}}^{2},\:q_{a_{y}}^{2},\:q_{a_{z}}^{2},\:q_{\dot\alpha_{D}}^{2},\:q_{\dot\alpha_{T}}^{2},\:q_{\dot\alpha_{C}}^{2}]^{T})$ $\cdots$ 암시적 기동모델

$Q$ $=$ $diag([q_{a_{x}}^{2},\:q_{a_{y}}^{2},\:q_{a_{z}}^{2},\:q_{\dot\alpha_{T}}^{2},\:q_{\dot\alpha_{C}}^{2}]^{T})$ $\cdots$ 명시적 기동모델

위 식에서 $(q_{a_{x}},\:q_{a_{y}},\:q_{a_{z}})$는 3축 가속도의 공정잡음 표준편차를, $(q_{\dot\alpha_{D}},\:q_{\dot\alpha_{T}},\:q_{\dot\alpha_{C}})$는 탄도계수 변화율에 해당하는 공정잡음 표준편차를 의미한다. 식 (16)에서 추적필터의 $k$시점 사후추정치를 $\hat x_{k|k}$라 하면 선형화 된 시스템행렬 $F_{k}$는 다음과 같이 정의된다.

$F_{k}=\left .\dfrac{\partial f(·)}{\partial x}\right |_{x =\hat x_{k|k}}$

한편, 재진입 단계 표적추적을 위한 지상레이다에서는 표적까지의 상대거리 및 수평/수직 시선각이 측정되므로, 측정방정식을 다음과 같이 쓸 수 있다.

여기서 측정치 $y_{k}$와 비선형함수 $h(·)$의 정의는 다음과 같다.

$y_{k}=[\widetilde r ,\:\widetilde\lambda_{h},\:\widetilde\lambda_{v}]^{T}$, \begin{align*} h(x)=\left[\begin{aligned}\begin{aligned}\sqrt{x^{2}+y^{2}+z^{2}}\\\tan^{-1}(y/x)\end{aligned}\\-\tan^{-1}(z/\sqrt{x^{2}+y^{2}})\end{aligned}\right] \end{align*}

또한, 레이다 측정잡음 $n_{k}$는 분산이 $R_{k}$인 영평균 백색잡음으로 가정해도 그 일반성을 잃지 않는다.

따라서, 탄도탄 표적추적 필터는 레이다 측정치 (17)과 운동방정식 (15)를 만족하는 상태변수 $x_{k}$에 대한 비선형 상태추정 문제로 귀결된다. 시스템 모델의 비선형성에 대처하기 위해 상태천이 및 측정치 갱신 시 UT(unscented transform)가 적용된 $\sigma -$포인트 필터를 활용한다. 이 경우 사전지식으로 주어진 $C_{D_{0}}$ 및 $\bar{k}$에 대한 미분 연산 없이도 시스템 전파 및 측정치 갱신이 가능하다. 추적필터 순환 식은 다음과 같다.

· 시스템 전파

$k-1$ 시점에서 식 (15)에 의해 전파된 $\sigma$-포인트 $\chi_{i}^{x}$는 다음과 같이 계산된다.

여기서 $L$은 상태변수의 차원 수이며, $\kappa$는 $\sigma$-포인트의 분포를 결정짓는 설계 파라미터이며, $\left(\sqrt[]{(L+\kappa)P_{x}}\right)_{i}$는 $(L+\kappa)P_{x}$의 제곱근 행렬의 $i$번째 열을 의미한다. 식 (18) 에서 $\chi_{i}^{x}$ 산출 시, 사전지식 (12)을 활용하여 예측 항력을 계산할 수 있다.

위의 식에서 $\hat h_{k-1}$, $\hat v_{k-1}$는 $k-1$시점 사후추정치로 계산된 고도 및 속력을 의미한다. 이 경우, $k-1$ 시점의 사전추정치 산출결과는 다음과 같다.

· 측정치 갱신

$\chi_{i}^{y}$을 식 (17)에 의해 측정치 공간으로 변환된 $\sigma$-포인트라 하면, $k$ 시점에서의 필터 이득은 다음과 같이 계산된다.

이 때, $k$ 시점의 사후추정치 및 사후추정오차 공분산은 다음과 같다.

4. 모의실험

제안기법의 유용성을 확인하기 위해, 그림 6과 그림 7의 변침기동 탄도탄 재진입 궤적에 대한 모의실험을 수행한다. 탄도탄 표적은 37.9초부터 우선회 및 급상승(pull-up) 한 후, 175.1초에는 좌선회, 182.0초부터는 급하강(pull-down) 기동으로 전환된다. 기동 가속도 크기는 최대 $4.04[g]$이다. 탄도탄 표적 기동모델이 추적성능에 미치는 영향을 분석하기 위해 다음 4가지 추적필터의 성능을 비교하였다.

방법 #1) 항력 탄도계수 $\alpha_{D}$만을 고려하는 전통적인 BRV(ballistic reentry vehicle) 필터

방법 #2) 3축 가속도 운동을 반영하기 위해 Singer 기동모델을 사용한 mBRV(modified BRV) 필터

방법 #1) 암시적 기동모델을 사용하여 설계된 MaRV 필터

방법 #2) 명시적 기동모델을 활용하는 제안기법

시나리오에 따른 사후 추정오차 공분산의 크기가 유사하도록 공정잡음 수준을 표 1과 같이 선정하였다. 필터의 샘플링 주기는 $0.1[s]$이다.

그림. 6. 변침기동 탄도탄 표적의 수평-수직 궤적

Fig. 6. Maneuvering ballistic missile horizontal-vertical trajectory

그림. 7. 기동 궤적 시나리오 속력 및 가속도 프로파일

Fig. 7. Absolute velocity and acceleration profile of maneuvering trajectory scenario

· 탄도탄 표적 기동모델의 적합성 평가

주어진 시나리오에 대한 각 추적필터의 단일시행(single run) 시뮬레이션 결과는 그림 8 및 그림 9와 같다. 전통적인 BRV 필터의 경우, 탄도탄 표적의 횡기동을 고려하지 못하므로 변침기동 구간에서 추정 오차가 급격하게 증가하여 거리오차는 약 70m, 속력오차는 약 –200m/s에 이른다. 이러한 오차는 필터에서 계산된 이론적 추정오차 표준편차의 2.2배와 10배에 이르는 매우 큰 값이다. 특히, 변침기동 구간에서 잔류치가 편향되는 특성이 있다. 실제 레이다 추적시스템에서는 필터에서 산출된 사전 추정치와 이론적 추정오차 공분산을 이용하여 게이팅(gating) 과정을 수행 한 후, 게이팅 된 유효측정치 만을 이용해 사후추정치를 산출한다. 따라서, 전통적인 BRV 필터 (방법 #1)는 사실상 탄도탄 표적추적에 실패할 가능성이 높다. 탄도탄 표적의 공력학적 특성을 전혀 반영하지 않은 Singer 기동모델을 이용한 mBRV 필터(방법 #2) 역시 추정오차의 첨두치가 미소하게 줄어들지만 BRV와 마찬가지로 추적실패 가능성에서 자유롭지 못하다.

그림. 8. 양항 특성이 미고려된 방법 #1,#2의 추적성능

Fig. 8. Tracking performance of method #1, #2 (without drag-to-lift ratio)

그림. 9. 양항 특성이 고려된 방법 #3,#4의 추적성능

Fig. 9. Tracking performance of method #3, #4 (with drag-to-lift ratio)

반면 암시적 기동모델에 기반한 MaRV 필터(방법 #3)는 최대 거리오차 25m, 속력오차 30m/s 미만의 추정성능을 보여준다. 각 추정치에 대한 잔류치 역시 거의 0으로 유지되어 변침기동에도 불구하고 전반적으로 안정적으로 표적추적이 이루어지고 있음을 알 수 있다. 다만 기동구간에서 속력오차가 크게 요동치며 증가하므로 탄도탄 표적의 궤적 예측 정확도가 다소 떨어질 가능성이 있다. 이와 달리, 명시적 기동모델을 사용한 제안기법(방법 #4)은 거리 추정오차는 암시적 기동모델을 사용한 추적필터와 대등하거나 다소 향상된 수준을 유지하지만, 속력오차가 크게 개선되는 양상을 보인다. 이는 명시적 기동모델이 양항 구속조건 식 (12)를 이용하여 탄도탄 표적 기동 모델링 불확실성을 줄일 수 있기 때문이다. 이상의 결과는 탄도탄 표적추적 필터 설계 시, 공력학적 특성을 반영한 기동 모델의 중요성을 다시 한 번 보여주는 것이다.

· 명시적 기동모델에 의한 추적성능 개선 효과

기동모델의 영향을 상세히 살펴보기 위해 200회 반복시행을 통해 추적성능을 분석하였다. 그림 10과 그림 11은 방법 #3과 #4의 추정성능을 도시한 것이다. 실제 상황에서는 표적양항 특성에 대한 사전지식이 불완전하므로, 그 영향을 확인하기 위해 식 (12)의 ($\alpha_{D_{0}},\:\bar{k}$)에 의도적으로 오차 $(\Delta\alpha_{D_{o}},\:\Delta\bar{k})$를 더해준 값을 이용해 제안한 추적필터를 설계하고 그 성능을 도출하였다. 이때, $\Delta\alpha_{D_{o}}$,$\Delta\bar{k}$의 크기는 참값의 10%, 20%로 제한하고, 해당 범위 내에서 균등분포를 따르도록 설정하였다.

그림에서 확인할 수 있듯이, 두 추적필터의 위치 및 속도 오차 평균은 재진입 구간에서 모두 0으로 수렴하는 것을 볼 수 있다. 흥미로운 점은 공력 불확실성이 5% 미만으로 작은 경우, 횡기동에 의한 유도항력이 존재하는 37.9초~175.1초 구간에서 제안기법의 수평면 위치/속도 오차 표준편차가 1.9~3.2배 개선된다는 점이다. 횡기동 종료 후, 항력만 작용하는 탄도비행 구간(175.1초 이후)에서는 두 방법이 유사한 오차수준을 보인다. 다만, 해당구간에서 명시적 기동모델을 사용하는 경우 $z$ 방향 위치/속도오차의 표준편차가 암시적 기동모델을 사용한 경우에 비해 1.8~2.0배 줄어든다. 이렇듯 추정오차 표준편차가 줄어드는 이유는 앞서 언급한 것과 마찬가지로 명시적 기동모델이 유도항력 구속조건 (12)를 최대한 활용하기 때문이다. 구속조건을 미 고려하는 암시적 기동모델을 사용한 추적필터는 구조적으로 ($\alpha_{D}$,$\alpha_{T}$,$\alpha_{C}$)의 불확실성이 운동방정식 (10)에 의해 3축 속도로 빠르게 전파될 수밖에 없다. 반면, 명시적 기동모델을 추적필터 설계에

그림. 10. 사전지식 불완전성에 따른 위치 추정성능

Fig. 10. The position estimation performance with respect to the uncertainty of a priori information

사용하면, 사전지식을 통해 $\alpha_{D}$가 구속조건 (12)로 대체되므로 기동모델의 불확실성이 전파되는 현상을 현저히 감소시킬 수 있다. 그림 12 및 그림 13의 탄도계수 추정결과에서 두 방법의 수직/수평양력 탄도계수 $\alpha_{T}$, $\alpha_{C}$ 추정성능은 서로 유사하지만, $\alpha_{D}$ 추정오차는 큰 차이를 보이는 이유도 동일한 맥락에서 설명 가능하다. 참고로, 속력 추정오차 수준은 탄도탄 표적에 대한 교전 우선순위 결정과 요격에 필수적인 표적궤적 예측성능에 지대한 영향을 끼친다. 따라서, 명시적 기동모델의 적용은 탄도탄 표적의 탄착점과CEP(circular error probability), 예상조우점(PIP: predicted impactpoint) 산출의 정밀도를 향상시켜, 탄도탄

그림. 11. 사전지식 불완전성에 따른 속도 추정성능

Fig. 11. The velocity estimation performance with respect to the uncertainty of a priori information

방어체계 전반의 성능 개선에 도움을 줄 수 있다. 한편, 그림 10 및 그림 11에서 볼 수 있듯이, 약 15%의 사전지식의 불완전성이 존재하더라도 명시적 기동모델이 암시적 기동모델에 비해 더 좋은 추정성능을 제공함을 알 수 있다. 사전지식의 불완전성이 커질수록, 명시적 기동모델에 의한 추정오차 평균은 암시적 기동모델에 비해 커지는 경향이 있다. 이와 달리, 명시적 기동모델을 적용한 추적필터의 추정오차 표준편차는 사전지식의 불완전성이 커지더라도 크게 증가하지 않는다. 이는 사전지식의 불완전성과 상관없이 구속조건 (12)가 표적 추적의 안정성 확보에 큰 도움이 됨을 암시하는 것이다. 특히, 사전지식 불

그림. 12. 탄도계수 참값 및 추정치 평균 $(\Delta\alpha_{D_{0}}=\Delta\bar{k}=0)$

Fig. 12. The true and average estimated values of the ballistic coefficients

그림. 13. 탄도계수 추정오차 표준편차 $(\Delta\alpha_{D_{0}}=\Delta\bar{k}=0)$

Fig. 13. Standard deviation of the estimation error of the ballistic coefficients

완전성에 의해 위치 추정오차가 소폭 증가 하더라도, 속도 추정오차 표준편차를 작게 유지 할 수 있으므로, 기동 탄도탄 표적 트랙유지와 궤적예측 측면에서 명시적 기동모델이 지닌 장점은 여전히 유효하다.

· 사전지식을 활용한 탄종 분류/식별 가능성

앞서의 시뮬레이션 결과를 토대로 명시적 기동모델을 사용하여 거둘 수 있는 부수적 효과를 생각해 볼 수 있다. 사전지식의 불일치는 추정오차의 편향, 즉 추적필터의 잔류치를 편향시킨다. 만일 탄도탄 표적의 탄종간 공력특성이 큰 차이를 보인다면,

그림. 14. 표적 형상 및 양항 특성 사전지식

Fig. 14. The target shapes and a priori informations of the drag-to-lift ratio

그림. 15. 다중모델 필터의 모델 확률

Fig. 15. Model probabilities of the multiple model filter

이러한 성질을 역으로 탄도탄 표적의 탄종 분류/식별에 활용할 수 있다. 그 가능성을 타진하기 위해, 그림 14의 가상 탄도탄 표적 2종 각각에 대해 방법 #4를 이용하여 다중모델 필터 뱅크를 구성하였다 ⁽¹⁴⁾. 탄종 별 양항 특성은 주어진 형상에 대해 공력 경험식을 적용하여 산출되었다 ⁽¹¹⁾. 궤적 생성에는 참표적의 공력특성이 사용되었다. 그림 15는 다중모델 필터 뱅크의 모델 확률을 나타낸 것이다. 대기밀도가 낮은 약 80초 이하 구간(고도 30km 이상)에서는 탄도계수에 대한 가관측성이 존재하지 않으므로, 모델 확률이 0.5를 중심으로 거의 동일하게 유지된다. 하지만, 탄도탄 표적이 대기권으로 진입하여 항력 및 양력이 작용하는 80초 이후 구간에서는 모델 확률의 차이가 발생하여 탄종 분류/식별이 가능함을 확인할 수 있다. 참고로, 기존 연구에 따르면 MaRV 필터(방법 #3)의 경우 기동 구간 발생하는 가속도를 추정할 수는 있지만, 이를 이용한 표적 분류/식별은 사실상 불가능한 것으로 알려져 있다. 따라서, 본 논문에서 제시한 명시적 기동모델에 근거한 표적추적 필터 설계 방법론을 채택할 경우, 변침기동 탄도탄 표적에 대한 안정적 트랙유지 및 궤적예측 정확도 향상, 그리고 표적 조기 식별에 따른 추적/요격 자원의 효율적 활용 등 BMDS의 성능 전반을 개선할 수 있다.

5. 결 론

기동 동특성을 공력 사전지식을 활용하여 명시적으로 기술하는 새로운 변침기동 탄도탄 추적필터를 제안하였다. 기동 탄도탄 표적추적의 성능을 좌우하는 핵심 요소가 일반화된 탄도탄 기동모델의 선정에 있다는 사실에 착안하여 탄도탄의 공력학적 특성에 의해 표적 기동 모델링에 사용되는 탄도계수 간의 명시적 상관관계를 추적필터의 구속조건으로 사용하였다. 이를 통해 변침기동 탄도탄에 대해 우수한 궤적추정 성능을 얻을 수 있었다. 특히, 명시적 기동모델의 사용을 통해 기동시점에서 수렴속도가 개선되어, 기존 기법에 비해 추정오차를 크게 감소시킬 수 있다. 공력 사전지식 모델이 실제 공력과 불일치 할 경우, 편향 오차가 유발될 수 있지만, 그럼에도 불구하고 안정적인 트랙유지가 가능함을 확인하였다. 제안한 기법은 탄도탄 표

적의 기동에 대한 대응능력과 함께 잔류오차를 활용한 표적 분류/식별의 가능성을 제공하므로, 탄도탄 방어체계 구축에 크게 기여할 수 있을 것으로 기대된다.