• 대한전기학회
Mobile QR Code QR CODE : The Transactions of the Korean Institute of Electrical Engineers
  • COPE
  • kcse
  • 한국과학기술단체총연합회
  • 한국학술지인용색인
  • Scopus
  • crossref
  • orcid

  1. (School of Mechanical and Control Engineering, Handong Global University, South Korea.)
  2. (LIG Nex1, South Korea.)



Long-term trajectory prediction, close-in-weapon system (CIWS), predicted impact point, anti-ship missiles, terminal maneuver

1. 서 론

중동전과 포클랜드분쟁을 거치면서 세계 각국은 대함유도탄위협을 효과적으로 방어할 수 있는 함정방어체계 구축에 심혈을 기울여왔다. 이에 따라 현대 전투함정들은 대부분 사거리 별 방어 수단을 순차적으로 가동하는 소위 다층방어개념을 도입하고 있다(1). 함정용 다층방공망의 한 축을 담당하고 있는 근접방어체계(CIWS: closed-in weapon system)는 지역방어(area defense) 망을 돌파한 표적에 대해 근거리에서 중・소구경 탄을 연사하여 요격하는 직접파괴(hard-kill) 방식의 요격체계이다. 다층방어체계가 상대적으로 취약한 우리나라의 상황을 고려할 때, 함정 전투력 제고 측면에서 CIWS가 차지하는 비중과 의미는 매우 크다고 할 수 있다.

일반적으로 CIWS는 탐색/추적 레이더, 사격통제컴퓨터와 함포로 구성되며, 교전가능 시간이 수십 초 내외로 매우 짧아 표적 탐색-추적-사격에 이르는 일련의 사격통제(fire control) 과정이 완전자동화 되어 있다(2). 이때, CIWS가 대함유도탄(ASM: anti-ship missile)을 성공적으로 요격하기 위해서는 표적 궤적변화와 요격탄의 탄착시간을 고려하여 예상 명중점(PIP: predicted impact point)을 산출해야 한다. 함포의 대응시간에 비해 요격 대상인 ASM이 고속으로 비행하므로 산출된 PIP의 정확도에 따라 CIWS의 전체 성능이 좌우되는 경향이 있다(3). 불행하게도 대부분의 최신 ASM은 CIWS 레이더에 의한 피탐 및 요격 가능성을 최소화하기 위해 다양한 종말 회피기동을 수행한다(4-8). 해면밀착 비행(sea-skimming) 직후 수평면 급선회(weaving)를 하거나, 함정에 접근하여 수직면으로 급상승(pop-up) 혹은 고입사각 급하강(dive) 하는 것이 그 대표적 사례이다. 이에 따라 호밍궤적에 영향을 끼치는 속력, 고도, 받음각, 추력 등 비행조건이 종말 단계에서 크게 변화하여 레이더 측정치를 이용한 장시간 표적궤적 예측이 매우 어렵다(9). 따라서, ASM의 비행동특성에 대한 고려 없이 위치-속도-가속도로 단순하게 ASM의 운동을 모델링하는 전통적인 표적추적필터 설계 방식으로는 만족할만한 ASM 방어성능을 확보할 수 없다(10).

이에 대한 대안으로 종말 유도법칙 추론에 기초한 표적상태추정 방법의 도입이 검토되었다(11,12). 이 방법은 표적이 채택할 수 있는 종말 유도법칙에 대한 가설을 설정하고 다중모델 필터링 혹은 EM(expectation maximization) 기법 등을 활용하여 적절한 가설을 취사선택하는 방식을 취하고 있다. 이를 통해, 기존 표적추적필터에 비해 향상된 궤적 예측 성능을 제공하는 것으로 알려져 있다. 하지만, 이들 방법은 수평면 혹은 수직면 상에서의 표적 회피기동만을 고려하고 있어 유도법칙이 복합적으로 사용되는 경우에 대응이 용이하지 않다. 더욱이 ASM의 속력이 일정하게 유지된다는 가정 하에 설계되어 표적 급기동 시 항력의 영향을 적절히 반영할 수 없다. 따라서, CIWS와 같이 수 초 가량의 장시간 표적궤적 예측이 불가피한 상황에서는 우수한 성능을 담보하기 어렵다.

이러한 한계를 극복하기 위해 본 논문에서는 ASM 표적의 동특성변화와 다양한 회피기동에 대응할 수 있는 새로운 형태의 장시간 궤적 예측 알고리듬을 제안한다. 이를 위해 일반적인 ASM 표적의 6자유도 운동모델이 사용된다. 궤적 예측 필터 설계 시 복잡도를 완화하기 위해, ASM의 추력제어 및 공력학적 특성을 반영하여 최소한의 미지 변수를 이용하여 ASM 비행동특성을 근사한다. 또한 수평면-수직면 상에서의 회피기동을 아우를 수 있는 일반화된 유도기법을 고려하여 표적 가속도를 시선각, 시선변화율, 편향항 및 중력의 함수로 기술한다. ASM 표적의 비행동특성 파라미터는 매우 느리게 변화하는 특성을 가지며 유도기법 구성에 필요한 이득들은 회피기동 전략에 따라 그 값이 달라지는 구간상수(piecewise constant)로 간주할 수 있으므로, 랜덤워크 형태의 불규칙 프로세스로 모델링할 수 있다. 이렇게 얻어진 시스템 모델의 비선형성을 감안하여 무향칼만필터(UKF: unscented Kalman filter) 이론을 적용하여 표적 궤적 예측 필터가 설계된다. 제안된 기법은 표적 비행동특성 변화 혹은 다양한 형태의 회피기동이 존재하는 상황에 능동적으로 대처할 수 있어, 장시간에 걸친 고기동 표적의 궤적 예측 성능, 더 나아가 CIWS의 ASM 표적 요격확률(kill probability)를 획기적으로 향상시킬 수 있다. 제안된 방법의 유용성을 확인하기 위해 전형적인 ASM 표적 조우 시나리오에 대한 시뮬레이션을 수행하고, 기존 기법과의 성능 비교분석 결과를 제시한다.

2. 고기동 대함유도탄의 운동모델

2.1 좌표계 및 기호 정의

다양한 형태의 종말기동을 수행하는 대함유도탄 표적의 운동을 모델링하기 위해 그림 1의 상대기하를 고려하자. 그림에 사용된 주요 좌표계는 다음과 같다.

2.1.1 관성좌표계 (I-frame)

함정에 장착된 CIWS용 탐색・추적 레이더를 원점으로 하는 NED(north-east-downward) 좌표계로, 표적의 관성운동을 기술하기 위한 기준 좌표계이다.

2.1.2 시선좌표계 (L-frame)

CIWS 레이더로부터 표적을 잇는 시선(LOS: line-of-sight) 벡터를 $X_{L}$축으로 하는 오른손 좌표계이다. 수평/수직면 시선각을 $\lambda_{h},\:\lambda_{v}$로 쓰면 I-frame으로부터 L-frame까지의 좌표변환행렬 및 회전각속도는 각각 다음과 같이 정의된다 (그림 1(a) 참조).

(1)
$C_{I}^{L}=R_{y}(\lambda_{v})R_{z}(\lambda_{h}),\: bold\omega_{IL}^{L}=\begin{bmatrix}-\dot\lambda_{h}s\lambda_{v}&\dot\lambda_{v} &\dot\lambda_{h}c\lambda_{v}\end{bmatrix}^{T}$

여기서 $R_{\chi}(\epsilon)$은 $\chi$축을 중심으로 $\epsilon$만큼 회전시키는 회전변환행렬(rotation matrix)을 의미하며, $c\eta =\cos(\eta),\:s\eta =\sin(\eta)$이다.

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

표적 무게중심을 원점으로 하고, 표적의 속도벡터 방향을 $X_{V}$축, 수평면 우측방향을 $Y_{V}$축으로 하는 오른손 좌표계이다. 그림 1(a)에 도시한 바와 같이 표적의 수평면, 수직면 비행경로각(flight path angle)이 각각 $\gamma_{h},\:\gamma_{v}$라면, I-frame에서 V-frame 으로의 좌표변환행렬 및 회전각속도는 다음과 같이 정의된다.

(2)
$C_{I}^{V}= R_{y}(\gamma_{v})R_{z}(\gamma_{h}),\: bold\omega_{IV}^{V}=\begin{bmatrix}-\dot\gamma_{h}s\gamma_{v}&\dot\gamma_{v}&\dot\gamma_{h}c\gamma_{v}\end{bmatrix}^{T}$

그림 1(b)에 도시한 바와 같이 받음각(AOA: angle of attack) $\alpha$는 ASM의 동체축($X_{B}$)과 속도벡터 간의 사잇각, 뱅크각 $\phi_{a}$는 $-Z_{V}$축과 가속도 벡터 $bold a_{L}$이 이루는 각도를 나타낸다.

그림. 1. 좌표계 정의

Fig. 1. Definitions of coordinate systems

../../Resources/kiee/KIEE.2021.70.2.408/fig1.png

2.2 대함유도탄의 6자유도 운동모델 근사

전술한 바와 같이 장시간 표적궤적 예측 필터의 정확도를 향상시키기 위해서는 ASM 표적동체에 작용하는 외력을 고려해야 한다. 실제 ASM의 6자유도 운동모델은 매우 복잡한 비선형 시변 미분방정식으로 표현되므로, 이를 직접 표적궤적 예측에 활용하는 것은 CIWS 교전 소프트웨어의 실시간 구현 측면에서 바람직하지 않다. 따라서, 본 논문에서는 표적 궤적 예측 정확도를 크게 저해하지 않는 선에서 6자유도 운동모델을 간략화하고, 이를 필터 설계에 활용하는 방법을 택한다.

ASM에 작용하는 외력은 양력(lift), 항력(drag), 추력(thrust), 중력(gravity)으로 구분할 수 있다. 이들 외력에 의한 가속도를 V-frame 성분으로 다시 쓰면 다음과 같다.

(3)
$bold a_{L}^{V}=\begin{bmatrix}0\\a_{y}\\a_{z}\end{bmatrix},\: bold a_{D}^{V}=\begin{bmatrix}-a_{D}\\0\\0\end{bmatrix},\: bold a_{T}^{V}=a_{T}\begin{bmatrix}c\alpha \\s\alpha s\phi_{a}\\-s\alpha c\phi_{a}\end{bmatrix},\: bold\mathrm{g}^{V}=\mathrm{g}\begin{bmatrix}-s\gamma_{v}\\0\\c\gamma_{v}\end{bmatrix}$

그림 1(b)에서 양력에 의한 가속도 $bold a_{L}$은 ASM의 속도에 수직한 방향으로 인가되는 유도명령 $(a_{y},\: a_{z})$에 의한 것이며, 항력에 의한 가속도 $bold a_{D}$는 속도벡터의 반대 방향으로 작용한다. ASM 표적의 특성 상 대부분의 운용 구간에서 받음각이 비교적 작게 유지되므로 다음 근사식을 적용해도 무방하다(13).

(4)
$a_{L}= DL\in E bold a_{L}DL\in E\simeq q\dfrac{S_{ref}}{m}C_{L_{\alpha}}\alpha ,\:$ $a_{D}= DL\in E bold a_{D}DL\in E\simeq q\dfrac{S_{ref}}{m}(C_{D_{0}}+C_{L_{\alpha}}\alpha^{2}),\:$

여기서 $q$는 동압(dynamic pressure), $S_{ref}$와 $m$은 각각 ASM의 동체 단면적과 질량, $C_{L_{\alpha}}$ 및 $C_{D_{0}}$는 공력미계수를 의미한다.

ASM 표적의 비행동특성은 다음과 같이 기술된다.

(5)
$\left .\dfrac{d bold V}{dt}\right |_{I}^{V}=\left .\dfrac{d bold V}{dt}\right |_{V}^{V}+ bold\omega_{IV}^{V}\times bold V^{V}$

여기서 V-frame 표적 속도벡터는 $bold V^{V}=\begin{bmatrix}V & 0 & 0\end{bmatrix}^{T}$이다.

식(3)으로부터 식(5)의 좌변을 계산할 수 있다.

(6)
$\left .\dfrac{d bold V}{dt}\right |_{I}^{V}= bold a_{L}^{V}+ bold a_{D}^{V}+ bold a_{T}^{V}+ bold g^{V}=\begin{bmatrix}a_{T}c\alpha - a_{D}- gs\gamma_{v}\\ a_{y}+ a_{T}s\alpha s\phi_{a}\\a_{z}- a_{T}s\alpha c\phi_{a}+ g c\gamma_{v}\end{bmatrix}$

따라서, 식(2)식(6)을 이용하여 ASM 표적의 동특성 방정식(5)를 다음 미분방정식으로 재정리할 수 있다.

(7)
$\begin{bmatrix}\dot V \\\dot\gamma_{h}\\\dot\gamma_{v}\end{bmatrix}=\begin{bmatrix}a_{T}c\alpha - a_{D}- gs\gamma_{v}\\\dfrac{1}{Vc\gamma_{v}}\left(a_{y}+ a_{T}s\alpha s\phi_{a}\right)\\ -\dfrac{1}{V}\left(a_{z}- a_{T}s\alpha c\phi_{a}+ gc\gamma_{v}\right)\end{bmatrix}$

표적 상대기하 변화를 고려하기 위해 앞서와 마찬가지로 코리올리 방정식을 적용하면 다음 관계식을 얻는다.

(8)
$\left .\dfrac{d bold R}{dt}\right |_{I}^{L}=\left .\dfrac{d bold R}{dt}\right |_{L}^{L}+ bold\omega_{IL}^{L}\times bold R^{L}$

여기서 L-frame 표적 위치벡터는 $bold R^{L}=\begin{bmatrix}R&0&0\end{bmatrix}^{T}$ 이다. 식(8)의 좌변을 다시 쓰면,

(9)
$\left .\dfrac{d bold R}{dt}\right |_{I}^{L}=C_{I}^{L}C_{V}^{I}bold V^{V}=V\begin{bmatrix}c\lambda_{v}c\gamma_{v}c(\gamma_{h}-\lambda_{h})+s\lambda_{v}s\gamma_{v}\\c\gamma_{v}s(\gamma_{h}-\lambda_{h})\\s\lambda_{v}c\gamma_{v}c(\gamma_{h}-\lambda_{h})-c\lambda_{v}s\gamma_{v}\end{bmatrix}$

식(1)을 이용해 식(8)의 우변을 계산한 후, 식(9)를 대입하면 ASM 표적의 상대기하 변화를 기술할 수 있다.

(10)
$\begin{bmatrix}\dot R \\\dot\lambda_{h}\\\dot\lambda_{v}\end{bmatrix}= \begin{bmatrix}V\left(c\lambda_{v}c\gamma_{v}c(\gamma_{h}-\lambda_{h})+ s\lambda_{v}s\gamma_{v}\right)\\\dfrac{V}{R c\lambda_{v}}c\gamma_{v}s(\gamma_{h}-\lambda_{h})\\ -\dfrac{V}{R}\left(s\lambda_{v}c\gamma_{v}c(\gamma_{h}-\lambda_{h})- c\lambda_{v}s\gamma_{v}\right)\end{bmatrix}$

식(7)식(10)으로부터 ASM 표적의 동특성과 상대운동이 6차 비선형 미분방정식으로 간략화 됨을 알 수 있다.

3. 종말기동을 고려한 장시간 표적 궤적 예측 필터링 기법

앞절에서 유도된 ASM의 6자유도 운동모델은 유도명령 $(a_{y},\: a_{z})$, 공력, 추력과 관련한 미지변수를 포함하고 있다. 따라서, ASM 표적의 장시간 궤적 예측을 위해 레이더에서 제공하는 상대거리 및 시선각 측정치를 이용하여 이들 미지 변수를 추정해야 한다.

3.1 대함유도탄의 종말 회피기동 고려방안

ASM의 종말 회피기동은 수직면과 수평면 상에서 서로 다른 패턴을 갖는다. CIWS가 다양한 형태의 회피기동에 대응하기 위해서는 이들 기동패턴을 아우를 수 있는 일반화된 ASM 유도기법을 고려하여 궤적 예측 필터를 구성할 필요가 있다.

∙ 수평면 회피기동

ASM의 수평면 회피기동은 급선회 유도기법에 의해 모델링 가능하다. 이 유도기법을 채택하는 경우 ASM은 해면밀착 비행과 동시에 주기적인 횡방향 급기동을 통해 선회궤적을 갖게 된다. 궤적의 곡률이 달라지는 영역에서 CIWS에 의한 궤적 예측오차가 크게 증가하는 경향이 있으므로 효과적인 종말 회피 기동 중 하나라 할 수 있다. 현존하는 기법들을 아우를 수 있는 급선회 유도명령은 다음과 같이 기술된다(7,8).

(11)
$a_{y}\simeq\eta_{7}V\dot\lambda_{h}+\eta_{8}R\lambda_{h}$

여기서 $\eta_{7}$과 $\eta_{8}$은 각각 수평면 유도오차를 0으로 만들기 위한 비례항법상수 및 선회기동에 사용된 제어이득을 나타낸다.

∙ 수직면 회피기동

ASM 표적은 레이더에 의한 피탐 가능성과 CIWS 함포의 대응사격 시간을 최소화하기 위해 원거리에서는 해면밀착비행 혹은 고고도 비행을 하다가 근거리에서 급상승(pop-up) 혹은 급하강(dive) 형태의 회피기동을 수행한다. 이러한 형태의 수직면 종말회피기동을 구현하기 위해서는 랑데부(Rendezvous), 잔여비행시간 다항식(time-to-go polynomial), 궤적변조(trajec- tory modulation), 입사각제어(impact angle control) 등 다양한 유도기법이 활용된다(4-6). 이들 기법에서 사용되는 유도명령을 일반화하면 다음과 같다.

(12)
$a_{z}\simeq\eta_{9}V\dot\lambda_{v}+\eta_{10}\dfrac{V^{2}}{R}\lambda_{v}+\eta_{11}\dfrac{V^{2}}{R}-gc\gamma_{v}$

여기서 $\eta_{9}\sim\eta_{11}$은 수직면 급기동 및 입사각 제어를 위한 이득을 의미한다. 우변 마지막 항은 중력보상항이다.

ASM 표적의 유도법칙에 사용된 제어이득 $\eta_{7}\sim\eta_{11}$은 종말 유도단계에 따라 그 값이 달라지는 구간상수(piecewise constant)로 간주할 수 있다. 따라서, 이들 변수는 랜덤워크(random walk) 프로세스로 모델링 가능하다.

(13)
$\dot\eta_{i}=w_{i},\: w_{i}\sim N(0,\: q_{i}),\: i=7,\:\cdots ,\: 11$

여기서 $N(m,\:\sigma^{2})$을 평균이 $m$, 분산이 $\sigma^{2}$인 정규분포를 의미한다.

3.2 고기동 표적 궤적 예측 필터 설계

종말회피기동을 수행하는 ASM의 궤적 예측 필터 설계에 앞서 시스템 모델을 유도해보자. 만일, 유도명령 $(a_{y},\: a_{z})$이 주어져 있다면, 그림 1(b)에서 $a_{L}=\sqrt{a_{y}^{2}+ a_{z}^{2}}$이므로, 식(4)로부터 ASM의 받음각을 계산할 수 있다.

(14)
$\alpha =\dfrac{a_{L}}{\eta_{12}},\:\eta_{12}=q\dfrac{S_{ref}}{m}C_{L_{\alpha}}$

식(14)를 동특성 방정식(7)에 대입하면 표적 속력변화에 관한 식을 다시 쓸 수 있다.

(15)
$\dot V =a_{T}c\alpha - a_{D}- gs\gamma_{v}=\eta_{13}-\dfrac{a_{y}^{2}+a_{z}^{2}}{\eta_{12}}- gs\gamma_{v}$,

여기서 $\eta_{13}=a_{T}c\alpha - q\dfrac{S_{ref}}{m}C_{D_{0}}$.

ASM의 받음각이 비교적 작은 값을 유지하므로 식(7)의 표적 비행경로각 변화를 다음과 같이 근사할 수 있다.

(16)
$\begin{bmatrix}\dot\gamma_{h}\\\dot\gamma_{v}\end{bmatrix}\approx\begin{bmatrix}\dfrac{1}{Vc\gamma_{v}}a_{y}\\-\dfrac{1}{V}\left(a_{z}+ gc\gamma_{v}\right)\end{bmatrix}$,

필터 설계를 위해 ASM의 동특성을 받음각, 뱅크각 등 빠르게 변화하는 변수 대신 $\eta_{12},\:\eta_{13}$을 이용하여 근사적으로 기술하는 것이 보다 유리하다. 식(15)에 사용된 변수 $\eta_{12}$는 동압과 공력미계수의 함수로 표현되므로 ASM의 고도와 속력에 따라 그 값이 달라질 수 있다. 하지만, ASM의 운용조건을 고려할 때 실제 상황에서 $\eta_{12}$의 변화폭은 그다지 크지 않다. 또한 ASM은 속도제어루프를 포함하고 있어 엔진 추력과 항력에 의한 가속도가 거의 상쇄되는 특성이 있다. 이러한 관찰 결과에 따라 미지수 $\eta_{12}$과 $\eta_{13}$를 매우 느리게 변화하는 불규칙변수로 모델링해도 무방하다.

(17)
$\dot\eta_{i}=w_{i},\: w_{i}\sim N(0,\: q_{i}),\: i=12,\: 13$

식(10),(13), (15)~(17)을 종합한 후, 이를 필터 계산주기 $T$로 이산화하면 장시간 표적궤적 예측필터의 시스템 모델을 얻을 수 있다.

(18)
$\left\{\begin{matrix}bold x_{k+1}& = bold f(bold x_{k})+ bold u_{k},\: \\ bold y_{k}&= H bold x_{k}+ bold v_{k}\end{matrix}\right .$

여기서 $N_{x}=13,\: N_{y}=3$라 하면 상태변수 $bold x\in ℝ^{N_{x}\times 1}$, 비선형함수 $bold f(\bullet)\in ℝ^{N_{x}\times 1}$, 레이더에서 제공되는 측정치 $bold y\in ℝ^{N_{y}\times 1}$, 측정행렬 $H\in ℝ^{N_{y}\times N_{x}}$와 상호비상관되어 있는 공정잡음 $bold u\in ℝ^{N_{x}\times 1}$ 및 측정잡음 $bold v\in ℝ^{N_{y}\times 1}$의 정의는 다음과 같다.

$bold x =\begin{bmatrix}R\\\lambda_{h}\\\lambda_{v}\\V\\\gamma_{h}\\\gamma_{v}\\\eta_{7}\\\vdots \\\eta_{13}\end{bmatrix},\: bold f(\bullet)= bold x + T\begin{bmatrix}V\left(c\lambda_{v}c\gamma_{v}c(\gamma_{h}-\lambda_{h})+s\lambda_{v}s\gamma_{v}\right)\\\dfrac{V}{Rc\lambda_{v}}c\gamma_{v}s(\gamma_{h}-\lambda_{h})\\-\dfrac{V}{R}\left(s\lambda_{v}c\gamma_{v}c(\gamma_{h}-\lambda_{h})-c\lambda_{v}s\gamma_{v}\right)\\\eta_{13}-\dfrac{a_{y}^{2}+a_{z}^{2}}{\eta_{12}}- gs\gamma_{v}\\\dfrac{1}{Vc\gamma_{v}}a_{y}\\-\dfrac{1}{V}\left(a_{z}+ gc\gamma_{v}\right)\\0^{7\times 1}\end{bmatrix},\:$ $bold y =\begin{bmatrix}\widetilde R \\\widetilde\lambda_{h}\\\widetilde\lambda_{v}\end{bmatrix},\:H=\begin{bmatrix}I^{3\times 3}&0^{3\times 10}\end{bmatrix},\: bold u\sim N(0,\: Q),\: bold v\sim N(0,\: R)$

이상의 모델링 결과로부터 장시간 표적 궤적 예측 문제가 불규칙 비선형 상태공간 방정식(18)에 대한 비선형 상태추정 문제로 귀결됨을 알 수 있다. 불행하게도 시스템 방정식에 사용된 함수 $bold f(bold x_{k})$로 인해 상태변수에 관한 확률분포의 비선형성이 야기된다. 따라서, 표적 상태변수에 관한 예측성능을 제고하기 위해서는 $bold f(bold x_{k})$에 의한 확률분포의 비선형성을 적절히 다룰 수 있어야 한다. 이를 위해 본 논문에서는 UT(unscented transform)를 통해 확률분포의 4차 적률까지 근사 가능한 무향칼만필터 이론을 활용하여 궤적 예측 필터를 설계한다. 잘 알려져 있듯이 무향칼만필터는 전통적인 확장칼만필터에 비해 상태변수에 관한 확률분포의 비선형성을 보다 정확히 반영할 수 있어 상대적으로 우수한 추정성능을 제공할 수 있을 뿐만 아니라, 필터 구현과정에서 자코비안 행렬 계산이 불필요하다는 장점이 있다(14).

이제, 불규칙 비선형 상태공간 방정식(18)에 무향칼만필터링 기법을 적용하면 다음과 같이 미지의 종말 회피기동과 표적동특성을 고려한 장시간 궤적 예측 필터를 설계할 수 있다.

∙ 사전추정치 확률분포의 시그마 포인트 산출

(19)
$bold x_{k|k-1}^{0}=\hat bold x_{k|k-1},\:$ $bold x_{k|k-1}^{i}=\hat bold x_{k|k-1}+\sqrt{N_{x}+\lambda}[\sqrt{P_{k|k-1}}]_{i},\:$ $bold x_{k|k-1}^{i+N_{x}}=\hat bold x_{k|k-1}-\sqrt{N_{x}+\lambda}[\sqrt{P_{k|k-1}}]_{i},\: i=1,\:\cdots ,\:N_{x}$

여기서 $\hat bold x_{k|k-1}\in ℝ^{N_{x}\times 1}$과 $P_{k|k-1}\in ℝ^{N_{x}\times N_{x}}$는 각각 사전추정치 및 추정오차 공분산을 의미한다. 상태변수의 확률분포 근사를 위해 필요한 시그마 포인트는 추정치 $\hat bold x$ 주위에서 추정오차공분산 $P$의 제곱근에 비례하도록 추출되며, 그 개수는 $2N_{x}+1$개 이다(15). 또한, $\sqrt{A}$는 QR 혹은 Cholesky 분해 등을 이용하여 계산되는 행렬 $A$의 제곱근, $[A]_{i}$는 행렬 $A$의 $i$번째 열, $\lambda =\alpha_{x}^{2}(N_{x}+\kappa)- N_{x}$는 시그마 포인트 선정에 사용된 비례 계수를 나타낸다. 참고로, 변수 $\alpha_{x}$는 시그마 포인트가 평균 주변에 퍼진 정도를 나타내며, $\kappa$는 근삿값의 고차 적률 일치과정에서 추가적인 자유도를 확보하기 위해 도입된 것이다.

∙ 측정치 갱신

(20)
$\hat bold y_{k|k-1}=\sum_{j=0}^{2N_{x}}W_{j}^{(m)}bold y_{k|k-1}^{j},\:$ $bold y_{k|k-1}^{i}=H bold x_{k|k-1}^{j},\:$ $P_{k|k-1}^{y}=\sum_{j=0}^{2N_{x}}W_{j}^{(c)}(bold y_{k|k-1}^{j}-\hat bold y_{k|k-1})(bold y_{k|k-1}^{j}-\hat bold y_{k|k-1})^{T}+R_{k},\:$ $P_{k|k-1}^{xy}=\sum_{j=0}^{2N_{x}}W_{j}^{(c)}(bold x_{k|k-1}^{j}-\hat bold x_{k|k-1})(bold y_{k|k-1}^{j}-\hat bold y_{k|k-1})^{T},\:$ $\hat bold x_{k|k}=\hat bold x_{k|k-1}+K_{k}(bold y_{k}-\hat bold y_{k|k-1}),\:$ $K_{k}=P_{k|k-1}^{xy}(P_{k|k-1}^{y})^{-1},\:$ $P_{k|k}=P_{k|k-1}-K_{k}P_{k|k-1}^{y}K_{k}^{T}$

여기서 $R\in ℝ^{N_{y}\times N_{y}}$는 식(18)에 정의된 측정잡음 $bold v$의 분산, $P_{k|k-1}^{y}\in ℝ^{N_{y}\times N_{y}}$는 잔차 분산, $P_{k|k-1}^{xy}\in ℝ^{N_{x}\times N_{y}}$는 측정치 $bold y_{k}$와 사전추정치 $\hat bold x_{k|k-1}$ 간의 공분산이다. 사전추정치 $\hat bold x_{k|k-1}$ 및 추정오차 공분산 $P_{k|k-1}$은 시그마 포인트 $bold x_{k|k-1}^{j}(j=0\sim 2N_{x})$의 가중 합으로 표현되며, 이때 사용된 가중치 $W_{j}^{(m)} $,$W_{j}^{(c)}$는 다음과 같이 설정된다.

(21)
$W_{0}^{(m)}=\lambda /(N_{x}+\lambda),\: W_{0}^{(c)}=\lambda /(N_{x}+\lambda)+(1 -\alpha_{x}^{2}+\beta_{x})$ $W_{j}^{(m)}=W_{j}^{(c)}= 1 /(2(N_{x}+\lambda)),\: j = 1,\:\cdots ,\: 2N_{x}$

위의 식에서 $\beta_{x}$는 상태변수 $bold x$의 확률 분포에 대한 4차 적률을 일치시키기 위해 도입된 계수이다.

∙ 사후추정치 확률분포의 시그마 포인트 산출

(22)
$bold x_{k|k}^{0}=\hat bold x_{k|k},\:$ $bold x_{k|k}^{i}=\hat bold x_{k|k}+\sqrt{N_{x}+\lambda}[\sqrt{P_{k|k}}]_{i},\:$ $bold x_{k|k}^{i+Nx}=\hat bold x_{k|k}-\sqrt{N_{x}+\lambda}[\sqrt{P_{k|k}}]_{i},\:i=1,\:\cdots ,\:N_{x}$

여기서 $\hat bold x_{k|k}\in ℝ^{N_{x}\times 1}$과 $P_{k|k}\in ℝ^{N_{x}\times N_{x}}$는 각각 사후추정치 및 추정오차 공분산을 의미하며, 시그마 포인트 $bold x_{k|k}^{j}(j=0\sim 2N_{x})$의 선정 방법은 앞서 식(19)에 사용된 방식과 동일하다.

∙ 시스템 전파

(23)
$bold x_{k+1|k}^{j}=f(bold x_{k|k}^{j}),\:j=0,\:\cdots ,\:2N_{x},\:$ $\hat bold x_{k+1|k}=\sum_{j=0}^{2N_{x}}W_{j}^{(m)}bold x_{k+1|k}^{j}$, $P_{k+1|k}=\sum_{j=0}^{2N_{x}}W_{j}^{(c)}(bold x_{k+1|k}^{j}-\hat bold x_{k+1|k})(bold x_{k+1|k}^{j}-\hat bold x_{k+1|k})^{T}+Q_{k}$

여기서 $Q\in ℝ^{N_{x}\times N_{x}}$는 앞서 정의된 공정잡음 $bold u$의 공분산이다.

CIWS 교전통제 소프트웨어의 PIP 산출에 필수적인 ASM 표적의 궤적 예측치는 최근 레이더 측정치를 이용한 사후추정치 산출결과 (20)을 기준으로 식(23)을 주어진 예측시간만큼 반복함으로써 손쉽게 계산 가능하다.

4. 성능 분석

제안한 기법의 유용성을 확인하기 위해 전형적인 ASM 교전시나리오에 대한 모의실험을 수행하였다. 다양한 회피기동 시나리오에 대한 대응 능력을 확인하기 위해 CIWS 성능평가에서 가장 중요하게 고려되는 수직면 급상승(pop-up) 및 수평면 급선회(weaving) 시나리오 각각에 대한 성능 분석을 수행한다. 시뮬레이션에 사용된 주요 파라미터는 표 1에 정리하였다. 궤적 예측오차 분석을 위해 몬테칼로 100회 반복실험을 수행한다. 실제 상황에서 빈번하게 사용되는 CT(constant turn) 모델에 기반하여 설계된 표적 궤적 예측필터를 비교대상으로 삼는다.

표 1. 시뮬레이션 조건

Table 1. Simulation conditions

item

parameters

radar

R}= 4.572[\mathrm{m}]$, $\sigma_{\lambda_{h}}=0.0573^{\circ}$, $\sigma_{\lambda_{v}}=0.0573^{\circ}$

target

pop-up

$bold R^{I}(0)=[12.97,\: 0,\: 0]^{T}[k\mathrm{m}],\:$

$bold V^{I}(0)=[-273.97,\: -93.70,\:0]^{T}[\mathrm{m}/\mathrm{s}]$

weaving

$bold R^{I}(0)=[9.26,\: 0,\: 0]^{T}[k\mathrm{m}],\:$

$bold V^{I}(0)=[-241.74,\: -82.68,\:0]^{T}[\mathrm{m}/\mathrm{s}]$

trajectory predictor

\begin{align*} P_{0|0}=diag(5\mathrm{m}2^{\circ}2^{\circ}10[\mathrm{m}/\mathrm{s}]2^{\circ}2^{\circ}\\ 1.00.52.02.01.50.1252.5)^{2} \end{align*} $Q = diag\left(0.00.00.00.0(1.0^{\circ}\cdot T)^{2}(1.0^{\circ}\cdot T)^{2}\right.$ $0.0(0.1^{2}\cdot T)(0.1^{2}\cdot T)(0.1^{2}\cdot T)$ $\left.(0.1^{2}\cdot T)(1.0^{2}\cdot T)(0.01\cdot T)^{2}\right)$, $R=diag(\sigma_{R}^{2},\:\sigma_{\lambda_{h}}^{2},\:\sigma_{\lambda_{v}}^{2})$, $T=20[m\mathrm{s}]$ $\alpha_{x}=0.5,\:\beta_{x}=2,\:\kappa = N_{x}- 3 = 10$

그림 2에 도시한 바와 같이, 급상승 시나리오는 초기 속력 $291.6[\mathrm{m}/\mathrm{s}]$로 해면밀착 비행을 하다 약 $31.3[\sec]$에 급기동하여 최대 고도 약 $0.33[k\mathrm{m}]$까지 상승한 후 최종 입사각이 약 $30^{\circ}$로 하강하는 시나리오이다. 이를 통해 조우직전 급격한 표적 비행동특성 변화가 궤적 예측성능에 미치는 영향을 분석할 수 있다.

그림 3은 고도 $4.5[\mathrm{m}]$로 해면밀착 비행을 하며 약 $0.05[\mathrm{Hz}]$ 주기로 수평면 급선회 기동을 수행하는 시나리오이다. 급선회 시 궤적곡률이 크게 변화하는 구간이 지속적으로 발생하므로 상대운동과 관련된 변수들만을 추정하는 기존 기법의 한계와 제안한 방법의 효과를 동시에 확인할 수 있다.

그림. 2. 표적 궤적: 급상승 기동

Fig. 2. Target trajectory: pop-up maneuver

../../Resources/kiee/KIEE.2021.70.2.408/fig2.png

CIWS 레이더의 최대 탐지거리는 $30[k\mathrm{m}]$, 요격가능 최대거리는 $3[k\mathrm{m}]$로 가정하였다. 현존하는 CIWS에 포함되어 있는 함포의 제원을 참고할 때 포탄의 발사관 이탈속도(muzzle speed)는 $900\sim 1100[\mathrm{m}/\mathrm{s}]$이므로, $2[\sec]$ 후의 궤적 예측치를 기준으로 성능 평가를 수행한다. CIWS에 의한 요격성능을 수치적으로 정량화하기 위해 PIP 예측오차의 거리 RMSE를 기준으로 요격확률이 산출된다. 실제 CIWS 함포는 초당 수천 발의 연사(simultaneous shots) 능력를 지니고 있어 요격확률 산출이 매우 복잡하지만, 본 논문에서는 다음과 같이 개산(槪算) 식을 사용하였다(16).

$P_{k}=1-(.5)^{(r^{2}/CEP^{2})}$, $CEP=\sigma\sqrt{2\ln 2}$,

여기서 $r$은 치명거리(lethal radius), $\sigma$는 궤적 예측치의 거리오차 표준편차를 의미한다. 참고로, 위와 같이 간략화된 식을 사용하더라도 요격확률의 상대적 비교에는 무리가 없다.

수직면 회피기동 시나리오에 대한 표적 동특성 및 유도명령 추정 결과는 그림 4에 도시된 바와 같다. 약 $31.28[\sec]$ 이전에는 ASM 표적이 속도와 고도를 일정하게 유지하는 순항비행 구간이므로 유도명령과 비행동특성 파라미터의 변화가 상대적으로 작다 (그림 2(b) 참조). 이로 인해 해당 구간에서 표적 동특성 추정오차 역시 작은 수준으로 유지된다. 급상승에 따라 표적 비행동특성 및 유도명령이 급격히 변화하는 경우 추정오차 및 표준편차도 순간적으로 증가하지만 이내 추정오차가 0으로 수렴하는 것을 확인할 수 있다. 이는 제안 방법이 상대운동과 관련된 변수뿐만 아니라 표적의 비행동특성 변화에 능동적으로 대처할 수 있음을 보여주는 결과이다.

그림. 3. 표적 궤적: 급선회 기동

Fig. 3. Target trajectory: weave maneuver

../../Resources/kiee/KIEE.2021.70.2.408/fig3.png

그림. 4. 표적 동특성 추정오차: 급상승 기동

Fig. 4. Target dynamics estimation: pop-up maneuver

../../Resources/kiee/KIEE.2021.70.2.408/fig4_1.png

../../Resources/kiee/KIEE.2021.70.2.408/fig4_2.png

그림. 5. 급상승 궤적 예측 성능

Fig. 5. Pop-up trajectroy prediction performance

../../Resources/kiee/KIEE.2021.70.2.408/fig5.png

수직면 급상승 궤적 예측 성능은 그림 5와 같다. 샘플링 주기 $20[m\mathrm{s}]$에 비해 예측시간이 $2[\sec]$로 상대적으로 매우 길게 설정되었음에도 불구하고 기동 시점을 제외한 전 구간에서 제안한 기법의 궤적 예측오차가 영평균에 가까울 정도로 만족할 만한 추정성능을 보인다. 반면, 기존 기법의 경우 종말 회피기동과 관련된 표적 동특성 및 유도명령을 고려하지 못해 급상승 기동이 지속되는 $31.28[\sec]$ 이후에 예측오차가 크게 증가한다. 기존 방법은 필터 구동 초기를 제외하면 순항구간에서는 궤적 예측거리 RMSE가 $10\sim 20[\mathrm{m}]$사이를 유지하다 급상승 기동구간에서 $30\sim 80[\mathrm{m}]$수준으로 급격히 증가한다.

그림. 6. 표적 동특성 추정오차: 급선회 기동

Fig. 6. Target dynamics estimation: weaving maneuver

../../Resources/kiee/KIEE.2021.70.2.408/fig6.png

제안 기법에 비해 기존 기법의 성능이 떨어지는 것은 궤적 예측 필터 설계 시 순항구간의 경우 추력, 항력 등 비행동특성 파라미터, 기동 구간에서는 표적 유도기법에 관한 적극적 고려가 부족했기 때문으로 판단된다. 이와 대조적으로, 제안한 기법은 급상승 기동 직후 필터 수렴에 필요한 시간(약 $1.5[\sec]$)를 제외하면 전 비행구간에서 궤적 예측거리 RMSE가 $5[\mathrm{m}]$내외로 억제된다. 그림 5(b)는 수직면 회피기동 상황에서의 요격확률을 도시한 것이다. 제안한 기법은 우수한 궤적 예측 성능을 바탕으로 표적 기동과 상관없이 표적이 근접함에 따라 요격확률이 약 $30\%$에서 $90\%$가량으로 지속적으로 증가한다. 이와 달리, 기존 방법은 교전가능 최대거리 $3[k\mathrm{m}]$에서 약 $15\%$수준에서 급상승 기동 직전에는 약 $50\%$까지 요격확률이 점차 증가하지만, 급기동을 궤적 예측에 적절히 반영하지 못해 표적이 $1.5[k\mathrm{m}]$이내로 근접했음에도 불구하고 오히려 요격확률이 $30\%$대로 감소하는 경향을 보인다. 이러한 결과는 비행동특성 및 유도기법의 적응 추정이 궤적 예측 정확도와 요격확률 향상에 결정적 역할을 함을 암시하는 것이다.

그림 6~7은 수평면 급선회 기동에 대한 성능분석 결과를 도시한 것이다. 제안한 기법을 적용할 경우 앞서와 마찬가지로 비행동특성 파라미터($\eta_{12},\:\eta_{13}$)와 유도명령을 성공적으로 추정하여 추정오차가 0으로 점근 수렴한다.

그림. 7. 급선회 궤적 예측 성능

Fig. 7. Weaving trajectroy prediction performance

../../Resources/kiee/KIEE.2021.70.2.408/fig7.png

그림 7은 $2[\sec]$ 후의 궤적 예측오차를 나타낸 것이다. 제안 기법과 기존 기법 모두 표적의 현재 상대기하를 추정하는 데는 큰 문제가 없지만, 궤적 예측 성능에는 매우 큰 차이를 보인다. 예상대로 기존 기법의 궤적 예측거리 RMSE는 급선회 유도명령에 동조하여(그림 3(a) 참조) 주기적으로 증감하는 것을 관찰할 수 있으며, 최대 오차가 약 $80[\mathrm{m}]$에 이른다. 반면, 제안된 궤적 예측 필터의 경우 전 비행구간에서 급선회 기동 여부와 상관없이 궤적 예측오차 RMSE가 $3[\mathrm{m}]$내외로 유지된다. 앞서 추정한 표적 궤적 예측오차 특성을 활용한 요격확률은 그림 7(b)와 같다. 기존 기법의 경우 ASM의 수평면 회피 기동에 따른 궤적 변화가 매우 큰 장거리에서 요격확률이 $25\% $내외로 매우 낮은 수준을 보인다. 이후, 표적이 함정에 접근하면서 급선회 기동의 크기가 줄어듦에 따라 요격확률이 약 $50\%$수준을 회복한다. 하지만, 이러한 저조한 궤적 예측 성능으로는 동일 표적에 대해 CIWS가 수차례 반복적으로 요격을 시도해야 하므로 다표적 교전 상황에서 CIWS가 충분한 교전시간을 확보하는데 어려움을 겪을 수밖에 없다. 반면 제안된 필터의 경우 표적 비행동특성 및 급선회 유도명령에 대한 우수한 추정결과를 바탕으로 유효 사거리 내에서 CIWS의 요격확률을 향상시키는데 기여한다. 제안된 기법을 활용할 경우 요격확률은 수평면 급선회 기동 여부와 상관없이 약 $65\% $수준으로 유지되어 기존 방법에 비해 장거리에서 동일표적 요격에 필요한 반복 교전회수와 교전시간을 효과적으로 감소시킬 수 있을 것으로 판단된다.

5. 결 론

CIWS의 요격확률 제고를 위한 고기동 ASM 표적의 장시간 궤적 예측 필터 설계방법을 제안하였다. 궤적 예측 성능이 ASM 표적의 비행동특성에 따라 크게 달라진다는 점에 착안하여 간략화된 6자유도 운동모델을 필터 설계에 활용하였다. 이때, 다양한 형태의 종말 회피기동에 효과적으로 대응하기 위해 ASM 표적의 일반화된 유도명령이 운동모델에 반영되었다. 비선형 상태추정 필터 설계를 통해 미지의 비행동특성 파라미터와 유도법칙에 사용된 제어이득을 추정이 가능함을 확인하였다. 제안한 기법을 사용하는 경우 예측시간 $2[\sec] $를 기준으로 ASM 표적의 급상승 및 급선회 기동 상황에서 기존 기법대비 요격확률이 최대 2배가량 증가한다. 따라서 제안한 장시간 궤적 예측 기법은 현재 우리 해군이 개발을 추진 중인 차세대 CIWS의 방공능력을 획기적으로 향상시킬 수 있는 기술적 해법 중 하나가 될 것으로 기대된다.

Acknowledgements

본 연구는 고속/고기동 해상 소형 표적 추적 기술(LIGNEX1-2020- 0741(00))의 일환으로 수행되었습니다. 연구지원에 감사드립니다.

References

1 
R. J. Prengaman, E. C. Wetzlar, 2001, Integrated ship defense, Johns Hopkins APL technical digest, pp. 523-535Google Search
2 
W. D. Blair, 1988, Some Alternates in Fire Control Processing for the PHALANX CIWS, NSWC TR88-143Google Search
3 
M. E. Pittelkau, 1984, Tracking and Track Prediction of Maneuvering Targets for Gun Engagement, NSWC TR84-393, Naval Surface Warfare CenterGoogle Search
4 
C. Song, C. K. Ryoo, M. J. Tahk, 2004, Optimal guidance law for impact angle control, IFAC Proceedings, pp. 653-658Google Search
5 
H. S. Shin, J. I. Lee, H. J. Cho, 2011, Trajec- troy modulation guidance law for anti-ship missiles, AIAA Guidance, Navigation and Control Conf.DOI
6 
T. H. Kim, C. H. Lee, 2013, Time-to-go poly- nomial guidance with trajectory modulation for observ- ability enhancement, in IEEE Trans. Aero. Electr. Sys., pp. 55-73DOI
7 
C. H. Lee, T. H. Kim, 2015, Biased PNG for target observability enhancement against non- maneuvering targets, IEEE Trans. Aero. Electr. Sys., pp. 2-17Google Search
8 
Joon Heo, Han Ho Shin, Seong-Jeub Jeon, 2017, Output Characteristics of Wireless Power Transfer System with Semi-random Magnetic Flux, KIEE Summer Conference 2017, pp. 12-14Google Search
9 
B. L. Clark, 1976, The development of an adaptive Kalman target tracking filter, AIAA Guidance and Control Conf., pp. 365-382DOI
10 
G. A. Watson, W. D. Blair, 1991, Constant speed prediction for maneuvering targets using a three dimensional turning rate, IEEE Southeastern Symposium on System Theory, pp. 239-243DOI
11 
J. Yun, C. K. Ryoo, K. Choi, 2010, Anti-ship missile’s guidance law estimation filter using IMM estimator, AIAA Guidance, Navigation, and Control Conf.DOI
12 
G. Xue, X. Wang, W. Zhan, Y. Liang, C. Zhao, 2016, Anti-ship missile target tracking based on guidance law estimation, IEEE Conf. Information FusionGoogle Search
13 
D. McLean, 1990, Automatic Flight Control Systems, PrenticeHallGoogle Search
14 
S. J. Julier, J. K. Uhlmann, 2004, Unscented filtering and nonlinear estimation, Proceedings of the IEEE, pp. 401-422DOI
15 
Simon J. Julier, Jeffrey K. Uhlmann, Hugh F. Durrant Whyte, 1995, A new approach for filtering nonlinear systems., American Control Conf.DOI
16 
A. R. Washburn, 2002, Notes on Firing Theory, Naval Postgraduate SchoolGoogle Search

저자소개

함다혜(Dahye Ham)
../../Resources/kiee/KIEE.2021.70.2.408/au1.png

2019년 한동대학교 기계제어공학부(공학사),

2019년~현재 동 대학원 기계제어공학과 석사과정.

관심분야는 레이다 표적추적, 라이다 확장표적추적, 차량용 센서융합 등

조형찬(Hyung-Chan Cho)
../../Resources/kiee/KIEE.2021.70.2.408/au2.png

2020년 한동대학교 기계제어공학부(공학사),

2020년~현재 동 대학원 기계제어공학과 석사과정. 관심분야는 이종센서 정보융합,

자율이동체 유도제어 등

나원상(Won-Sang Ra)
../../Resources/kiee/KIEE.2021.70.2.408/au3.png

1998년 연세대학교 전기공학과(공학사),

2000년, 2009년 동 대학원 전기컴퓨터공학과(공학석사), 전기전자공학과(공학박사).

2000년~2009년 국방과학연구소 유도조종부 선임연구원.

2009년~현재 한동대학교 기계제어공학부 교수.

2015년 Cran- field University 방문교수.

관심분야는 상태추정 및 정보융합 이론, 표적추적필터, 자율이동체 유도조종기법 등.

안지훈(Ji-Hoon An)
../../Resources/kiee/KIEE.2021.70.2.408/au4.png

2014년 한양대학교 전자시스템공학(공학사),

2016년 동 대학원 전자시스템공학(공학석사).

2016년~현재 LIG넥스원(주) 레이다연구소 선임 연구원.

관심분야는 레이더 제어, 표적추적 알고리즘 SW 개발 등

김상현(Sang-Hyun Kim)
../../Resources/kiee/KIEE.2021.70.2.408/au5.png

1998년 한국항공대학교 항공우주공학과(공학사),

2015년 고려대학교 국방기술경영학과(공학석사).

2007년~현재 LIG넥스원 CIWS-II 사업단 수석연구원.

관심분야는 체계설계, 사격통제체계, 구동장치, 진동/충격/소음 등

박준현(June-Hyune Park)
../../Resources/kiee/KIEE.2021.70.2.408/au6.png

1996년 홍익대학교 전자공학과(공학 석사).

1996년 3월~현재 LIG넥스원(주) 레이다연구소 수석연구원.

관심분야는 능동 위상 배열 레이다, 항공기 레이다, 레이다 신호처리 등