The Journal of Korean Institute of Electromagnetic Engineering and Science

The Korean Institute of Electromagnetic Engineering and Science

J. Korean Inst. Electromagn. Eng. Sci. 2021; 32(1):66-74

pISSN: 1226-3133, eISSN: 2288-226X

DOI: https://doi.org/10.5515/KJKIEES.2021.32.1.66

논문/REGULAR PAPERS

발사체 분류를 위한 탄도계수 추정 방법 비교

오유근¹^,^†

, 전주환¹

, 정광용^*

, 김남문^*

Ballistic Coefficient Estimation for Classification of Projectiles

You Geun Oh¹^,^†

, Joohwan Chun¹

, Kwang Yong Jung^*

, Nam Moon Kim^*

¹한국과학기술원 전기전자공학부

^*한화시스템

¹Electrical Engineering, KAIST

^*Hanwha System

^†Corresponding Author: You Geun Oh (e-mail: yuk0320@kaist.ac.kr)

© Copyright 2021 The Korean Institute of Electromagnetic Engineering and Science. This is an Open-Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/4.0/) which permits unrestricted non-commercial use, distribution, and reproduction in any medium, provided the original work is properly cited.

Received: Aug 10, 2020; Revised: Aug 24, 2020; Accepted: Nov 27, 2020

Published Online: Jan 31, 2021

요약

탄도계수는 발사체가 공기 저항을 뚫고 얼마나 잘 나아갈 수 있는지를 수치화한 것이다. 발사체의 운동방정식은 탄도계수에 대한 함수이기 때문에 발사체의 탄도계수를 추정하는 방법은 발사체의 궤적을 추적할 수 있는 등의 이유로 많은 연구가 이루어졌다. 본 논문에서는 발사체를 분류하기 위해 탄도계수를 추정하는 방법을 제안했다. 사용된 주요 기법은 확장 칼만 필터와 상미분방정식 파라미터 추정이다. 그리고, 실제 발사체의 궤적을 시뮬레이션을 통해 만들고, 그 궤적을 사용하여 제안한 방법들로 발사체의 탄도계수를 추정하고, 그 성능을 비교 분석했다.

Abstract

The ballistic coefficient is a measure of a projectile’s ability to overcome air resistance. Because the motion equations of the projectile are functions of the ballistic coefficient, several studies estimating the projectile’s ballistic coefficient have been conducted to track its trajectory. In this paper, two methods are proposed to estimate the ballistic coefficient for classification of the projectile. The primary techniques used are the extended Kalman filter and ordinary differential equation parameter estimation. The true projectile trajectory is created via simulation, and the proposed method is used to estimate the ballistic coefficient. Furthermore their performances are compared and analyzed.

Keywords: Ballistic Coefficient; Extended Kalman Filter; Ordinary Differential Equation Parameter; Motion Equation; Trajectory

Ⅰ. 서 론

탄도계수는 발사체가 공기 저항을 뚫고 얼마나 잘 나아갈 수 있는지를 수치화한 것이다. 발사체는 종류에 따라서 다양한 탄도계수를 가지고 있다. 높은 탄도계수를 가지고 있는 발사체는 공기 중을 날아갈 때, 공기 저항을 이겨내는 힘이 강하므로 속도 감소율이 낮다고 할 수 있다. 탄도계수는 발사체의 질량에 비례하고, 단면적과 항력 계수에 반비례하다.

이전까지 각종 발사체에서 탄도계수를 추정하려는 연구들이 있었다. 발사체의 운동방정식은 탄도계수와 밀접한 관련이 있다. 따라서, 발사체의 탄도계수를 추정하면 발사체의 궤적을 잘 추적할 수 있기 때문에 중요성이 크다. 참고문헌 [1]은 high accuracy satellite drag model (HASDM)을 세우기 위해 탄도계수를 추정했다. 참고문헌 [1]에서 HASDM은 정확한 기온과 공기 밀도를 측정하기 위해 사용되었다. 참고문헌 [2]는 과거의 two line elements(TLE) 정보를 통해 저고도 파편의 탄도계수를 추정했다. 참고문헌 [2]가 탄도계수를 추정한 이유는 파편의 궤적을 추적하기 위해서이다. 참고문헌 [3]은 2D 환경에서 여러 종류의 칼만 필터를 사용하여 발사체의 탄도계수를 추정하고, 그 성능을 비교분석하였다. 참고문헌 [4]와 참고문헌 [5]는 재진입 단계에서 발사체의 탄도계수를 추정했다. 발사체의 종류에 따라서 발사체는 각기 다른 궤적을 가진다. 또한, 다른 종류의 발사체는 각각 사용성이 다르므로 발사체의 종류를 알아내면 그 목적에 맞게 대응할 수 있을 것이다. 본 논문에서는 발사체를 분류하기 위한 목적으로 발사체의 탄도계수를 추정한다. 참고문헌 [1]은 발사체가 아닌 공기의 특성을 알아내기 위한 연구이고, 참고문헌 [2]는 추력이 꺼진 후, 더 이상 발사체의 모습을 가지지 않는 파편에 대해서 탄도계수를 추정하였다. 참고문헌 [3]은 2D 상황을 가정했고, 참고문헌 [4]와 참고문헌 [5]는 재진입 단계라는 특정한 상황을 가정하였다. 본 논문은 3D 공간의 일반적인 상황에서 발사체의 분류를 위해 탄도계수를 추정한 새로운 시도를 했다.

본 논문의 2장에서는 레이다와 발사체의 개념도를 제시하고, 탄도계수를 추정하기 위한 발사체의 궤적을 만들기 위해 운동방정식을 제안했다. 3장에서는 널리 쓰이는 방법인 확장 칼만 필터 방법을 소개했다. 4장에서는 확장 칼만 필터 방법을 응용한 non-state 방법을 새로 제시했다. 이 방법은 확장 칼만 필터로 추정된 속도를 이용해서 탄도계수를 추정하는 방식이다. 5장에서는 상미분방정식 파라미터 추정 방법을 새로 제시했다. 이 방법은 비용함수를 설정하고, 비용함수를 최소화시키는 최적화 문제를 푸는 방법이다. 6장에서는 시뮬레이션 결과를 통해 얻은 논문의 결론과 탄도계수 추정의 추후 연구 방향을 제시했다.

Ⅱ. 발사체의 운동방정식

탄도계수를 추정하기 위해서는 시간에 따른 발사체의 궤적이 필요하다. 이번 장에서는 발사체의 운동방정식을 세우고, 그 운동방정식을 토대로 발사체의 궤적을 시뮬레이션했다. 발사체에 가해지는 힘은 항력, 중력으로 제한했고, 발사체의 추력이 꺼진 후의 상황으로 가정했다. 이 시뮬레이션 결과는 3장, 4장, 5장에서 탄도계수를 추정할 때 실제 궤적으로 여겨진다.

그림 1에 지구 표면에 위치한 레이다와 상공에 있는 발사체를 나타냈다. 그림 1에는 두 개의 좌표계가 있는데, 하나는 지구 중심을 원점으로 하는 ECEF 좌표계 (earth-centered earth-fixed coordinates)이고, 다른 하나는 레이다의 위치를 원점으로 하는 SEZ 좌표계(topocentric-horizon coordinates)이다. 레이다의 위도와 경도는 각각 λ, L이고, 지구는 반지름 r_e인 완전한 구라고 가정했다.

그림 1. | Fig. 1. 레이다와 발사체의 개념도 | Geometry of the radar and the target.

Download Original Figure

발사체의 운동방정식은 다음 뉴튼의 운동방정식으로부터 얻을 수 있다^[6]. 발사체는 하나의 질점이라고 가정되었다.

d v e f d t | e f = d v e f d t | i − ω × v e f = f − 2 ω × v e f + g

(1)

여기서 v_ef는 ECEF 좌표계에 대한 발사체의 속도, $d v e f d t | e f$ 는 ECEF 좌표계에 대한 v_ef의 시간에 대한 미분, $d v e f d t | i$ 는 ECI 좌표계에 대한 v_ef의 시간에 대한 미분, ω는 지구 자전 각속도를 의미한다. g, f는 각각 중력과 중력을 제외한 힘을 의미한다. ω×v_ef는 코리올리 효과를 의미한다.

발사체의 운동방정식은 SEZ 좌표계에서 나타냈으며, 3개의 위치방정식과 3개의 속도방정식으로 이루어져 있다. 위치의 시간에 대한 미분은 속도로 나타냈고, 속도의 시간에 대한 미분은 식 (1)을 활용하여 나타냈다. 아래 첨자 s, e, z는 SEZ 좌표계의 축을 의미하고, p와 v는 각각 발사체의 위치와 속도를 의미한다.

[p ˙ s p ˙ e p ˙ z] = [v s v e v z] [v ˙ s v ˙ e v ˙ z] = [f s f e f z] − 2 ω [− v e sin L v s sin L + v z cos L − v e cos L] + [g s g e g z]

(2)

s축과 e축에 가해지는 중력은 너무 작으므로 무시했다. z축으로 가해지는 중력은 지구 표면에서의 높이에 따라 달라지는 함수이며, 지구 자전으로 인해 생기는 원심력은 무시했다. 중력에 대한 식은 식 (3)과 같다. h는 지구 표면에서의 높이를 의미한다.

g s = 0, g e = 0 g z = − g = − 9.8 ∗ (1 − 2 h r e) (m / s 2)

(3)

f에 해당하는 힘은 항력만 남게 되고, 그 크기는 다음의 식과 같이 얻을 수 있다. D, ρ, C_D, S는 각각 항력의 크기, 공기 밀도, 항력계수, 발사체의 단면적을 의미한다. α, β는 발사체의 속도 방향을 알려주는 각도로 그림 2에 그 의미를 더 자세히 나타냈다. M은 발사체의 질량을 의미한다.

v 2 = v s 2 + v e 2 + v z 2, D = 12 ρ v 2 C D S

(4)

cos α = v s 2 + v e 2 v s 2 + v e 2 + v z 2, cos β = v s v s 2 + v e 2

(5)

[f s f e f z] = − 1 M [D cos α cos β D cos α sin β D sin α]

(6)

그림 2. | Fig. 2. 발사체의 속도 방향 | Direction of the velocity.

Download Original Figure

공기 밀도는 지구 표면에서의 높이가 25,000 m 이상일 때 다음과 같은 식으로 나타낼 수 있다^[7]. T, P는 각각 온도와 압력을 나타낸다.

T = − 131.21 + 0.00299 h [∘ C]

(7)

P = 2.488 (T + 273.1 216.6) − 11.388 [kg / m 2]

(8)

ρ = P 0.2869 (T + 273.1) [kg / m 3]

(9)

탄도계수는 식 (10)과 같이 정의된다^[1]. 발사체의 질량과 단면적이 일정하게 유지될 때, 항력계수가 변화함에 따라서 탄도계수도 변화하게 된다. 항력계수는 물체의 속도에 따라서 크게 변하기 때문에 물체가 일정한 속도를 가지고 있을 때, 탄도계수가 일정하다고 가정할 수 있다. 식 (4)와 식 (10)을 결합하면 $D M = 1 2 b ρ v 2$ 을 얻을 수 있다. 식 (2), 식 (3), 식 (5)와 논문의 가정들을 결합하면 다음과 같은 7개의 미분방정식으로 발사체의 운동을 표현할 수 있다. 여기서 ρ, υ, α, β, g는 모두 x의 함수이다.

b = M C D S

(10)

x = [x 1, x 2, x 3, x 4, x 5, x 6, x 7] T = [p s, p e, p z, v s, v e, v z, b] T x ˙ 1 = x 4, x ˙ 2 = x 5, x ˙ 3 = x 6 x ˙ 4 = − 1 2 x 7 ρ v 2 cos α cos β − 2 ω x 5 sin L x ˙ 5 = − 1 2 x 7 ρ v 2 cos α cos β − 2 ω (x 4 sin L + x 6 cos L) x ˙ 6 = − 1 2 x 7 ρ v 2 sin α + 2 ω x 5 cos L − g x ˙ 7 = 0

(11)

높은 상공에서는 공기의 밀도가 낮아지기 때문에 발사체의 운동이 달라지고, 지구의 자전효과도 지표면에서보다 훨씬 큰 영향을 가지게 된다. 이번 장에서 생성된 궤도는 이 영향들을 모두 고려했다. 또한, 식 (1)에서 계산이 되는 중력의 원심력 요소와 같이 지표면의 운동방정식에서 흔히 무시되는 요소들을 전부 고려해 주었기 때문에 높은 상공에서도 실제 발사체의 궤적을 잘 표현하고 있다. 만약, 하드웨어적인 문제로 인해 실제 미사일의 궤도가 달라지더라도 그에 맞춰 탄도계수도 함께 달라지기 때문에 탄도계수를 추정하는 것에는 문제가 되지 않는다.

이렇게 만들어진 궤도는 공기의 밀도가 작아지는 높은 상공에서의 달라지는 발사체의 운동, 지구의 본 논문에서는 2가지 탄도계수를 사용했다. 그 값은 각각 1,800 kg/m²와 3,200 kg/m²이다. 식 (11)을 이용해서 만든 발사체의 궤적은 그림 3과 같다. 그림 3에 사용된 탄도계수는 3,200 kg/m²이다. 시뮬레이션에 사용한 나머지 파라미터들은 표 1에 보였다.

그림 3. | Fig. 3. 발사체의 궤적 | True trajectory of the projectile.

Download Original Figure

표 1. | Table 1. 시뮬레이션 파라미터 | Simulation parameters.

Sampling interval	0.1(s)
Number of samples	4,500
ω	7.2921159×10⁻⁵(rad/s)
r_e	6,378,137 (m)
L	36.6424°
x _(0,1:3)	[10,000 10,000 43,000](m)
x _(0,4:6)	[400 400 350](m/s)

Download Excel Table

Ⅲ. 확장 칼만 필터-State 방법

이번 장에서는 탄도계수를 추정하는데 널리 쓰이는 방법인 확장 칼만 필터 state 방법을 소개하고, 시뮬레이션 결과를 보였다. State 방법은 탄도계수를 상태 벡터에 포함하여 확장 칼만 필터를 진행하고, 그 결과로 얻은 추정 탄도계수 값을 그대로 추정값으로 한다.

확장 칼만 필터의 모델과 과정은 다음과 같다^[8]. k는 샘플링 시간, $x^k ∣ k$ 는 관측치 업데이트를 통해 추정된 시간 k에서의 보정 추정값, $x^k ∣ k − 1$ 은 시간 업데이트를 통해 추정된 시간 k에서의 예측 추정값이다. P_k|k는 시간 k에서의 보정 오차 공분산 행렬이고, P_k|k-1는 시간 k에서의 예측 오차 공분산 행렬이다. e_k는 시간 k에서의 이노베이션 벡터, z_k는 시간 k에서의 관측치를 의미한다. w_k는 공분산이 Q인 프로세스 노이즈, v_k는 공분산이 R인 관측 노이즈를 의미한다. 프로세스 노이즈와 관측 노이즈는 모두 가우시안 노이즈이다.

x k = f ˜ (x k) + w k (상태 모델) z k = h (x k) + v k (관측 모델) x^0 | 0, P 0 | 0 지정 for k = 1,2,3, … x^k | k − 1 = f ˜ (x^k − 1 | k − 1) P k | k − 1 = F k − 1 P k − 1 | k − 1 F k − 1 T + Q R e, k = R + H k P k | k − 1 H k T e k = z k − h (x^k | k − 1) x^k | k = x^k | k − 1 + P k | k − 1 H k T R e, k − 1 e k P k | k = P k | k − 1 − P k | k − 1 H k T R e, k − 1 H k P k | k − 1 end

(12)

F k = ∂ f ˜ ∂ x T | x = x^k | k = [∂ f ˜ 1 ∂ x 1 ... ∂ f ˜ 1 ∂ x 7 ... ∂ f ˜ 7 ∂ x 1 ... ∂ f ˜ 7 ∂ x 7] | | x^k | k

(13)

H k = ∂ h ∂ x T | x = x^k | k − 1 = [∂ h r ∂ x 1 ... ∂ h r ∂ x 7 ∂ h θ ∂ x 1 ... ∂ h θ ∂ x 7 ∂ h ϕ ∂ x 1 ... ∂ h ϕ ∂ x 7] | | x^k | k − 1

(14)

함수 f ̃는 발사체의 운동방정식 식 (11)을 나타내는 함수로 다음과 같이 정의된다. △t는 샘플링 간격을 의미한다.

f ˜ (x) = x + f (x) Δ t

(15)

f (x) = [f 1 (x), f 2 (x), f 3 (x), f 4 (x), f 5 (x), f 6 (x), f 7 (x)] f 1 (x) = x 4, f 2 (x) = x 5, f 3 (x) = x 6 f 4 (x) = − 1 2 x 7 ρ v 2 cos α cos β + 2 ω x 5 sin L f 5 (x) = − 1 2 x 7 ρ v 2 cos α sin β − 2 ω (x 4 sin L + x 6 cos L) f 6 (x) = − 1 2 x 7 ρ v 2 sin α + 2 ω x 5 cos L − g f 7 (x) = 0

(16)

함수 h는 레이다로 들어오는 관측치를 나타내는 함수이다. 레이다는 레이다에서 표적까지의 거리 r, 고각 θ, 방위각 ϕ를 관측치로 받는다. 함수 h는 다음과 같은 식으로 나타낼 수 있다.

[r, θ, ϕ] T = h (x)

(17)

r = h r (x) = x 12 + x 22 + x 32

(18)

θ = h θ (x) = sin − 1 (x 3 r)

(19)

ϕ = h ϕ (x) = tan − 1 (x 2 x 1)

(20)

F_k와 H_k는 각각 식 (13)과 식 (14)와 같이 정의된 함수 $f ˜$ 와 h의 자코비안 행렬이다.

시뮬레이션에서 사용한 값은 모두 표 1과 동일한 값을 사용했다. 확장 칼만 필터에서 초기값을 설정해 주어야 하는 $x^0 | 0$ 와 P_0|0는 각각 다음과 같은 값을 주었다.

x^0 | 0 = x 0 + n, n ~ N (0, P 0 | 0) P 0 | 0 = d i a g (1000, 1000, 1000, 100, 100, 100, 1000)

(21)

프로세스 노이즈와 측정 노이즈는 각각 다음과 같다. 측정치의 단위는 (m, rad, rad)이다.

Q = 0, R = d i a g (10, 0.01, 0.01)

(22)

그림 4에 b=3,200일 때, 시간에 따른 탄도계수 추정을 보 였다. 시간이 지날수록 칼만 필터의 추정값이 실제값으로 잘 추정이 되고 있음을 볼 수 있다. 100번의 몬테카를로 방법으로 시뮬레이션한 결과, 추정 탄도계수의 평균과 평균 제곱근 편차는 각각 3200.1, 2.01555였다. b=1,800일 때의 두 값은 1,799.8, 0.8484였다. 추정이 비교적 잘 되었다고 볼 수 있는 기준인 평균 제곱근 편차가 5 이하로 떨어지는 데 필요한 샘플링 개수는 b=1,800일 때 1,100개였다.

그림 4. | Fig. 4. b=3,200일 때 탄도계수 추정 (3장) | Ballistic coefficient estimation when b=3,200.

Download Original Figure

확장 칼만 필터의 상태 벡터에는 탄도계수뿐만 아니라, 발사체의 위치와 속도도 들어가기 때문에 이 항들도 역시 추정된다. 시뮬레이션 결과, 20초 이후에는 속도의 오차가 0.5 m/s 이내로 추정됨을 확인했다. 속도의 추정이 비교적 안정적이기 때문에 이 사실을 이용해 4장에서 새로운 방법을 제시한다.

Ⅳ. 확장 칼만 필터-Non-State 방법

이번 장에서는 탄도계수를 확장 칼만 필터 non-state 방법을 새로 제시하고, 시뮬레이션 결과를 보였다. 3장의 시뮬레이션에서 탄도계수보다는 속도의 추정이 더 안정적임을 확인할 수 있었다. Non-state 방법도 state 방법과 마찬가지로 탄도계수를 상태 벡터에 포함하여 확장 칼만 필터를 진행한다. 그러나 그 결과로 얻은 추정 탄도계수 값을 추정값으로 사용하지 않는다. 그 대신, 추정 속도를 이용해서 일정한 계산 과정을 통해 탄도계수 추정값을 얻는다. 상태 벡터에 들어있는 탄도계수는 확장 칼만 필터의 과정에서 위치, 속도 추정에 영향을 주기 때문에 여전히 의미가 있다. 식 (4)와 식 (10)을 결합하면 다음과 같이 쓸 수 있다.

b = M C D S = ρ v 2 M 2 D = ρ v 2 2 a D

(23)

여기서 a_D는 항력에 의한 가속도를 의미한다. 칼만 필터를 통해 추정한 속도를 통해 발사체의 추정 가속도를 계산할 수 있다. 추정 가속도에서 항력에 의한 가속도만 구하는 방법은 다음과 같다. a_aero를 발사체의 가속도에서 중력가속도와 지구의 자전으로 인한 가속도를 뺀 나머지 가속도라고 하자. 그렇다면 a_aero는 오차가 섞여 있기 때문에 발사체의 속도 방향과 완전한 반대 방향을 보이지는 않을 것이다. a_aero를 속도에 대한 방향으로의 크기를 추출하면 a_D가 된다. 그림 5에 발사체의 가속도들을 나타냈다. 이를 식으로 나타내면 다음과 같다.

a a e r o = a − [00 − g] − [2 ω v e sin L − 2 ω (v s sin L + v z cos L) 2 ω v e cos L]

(24)

a D = a a e r o · v / | | v | |

(25)

그림 5. | Fig. 5. 발사체의 가속도 | Acceleration of the projectile.

Download Original Figure

이렇게 얻은 $a^D$ 를 식 (23)에 대입하면 non-state 방법으로 추정한 탄도계수를 얻을 수 있다.

시뮬레이션에서 사용한 값은 모두 3장의 시뮬레이션과 동일하다. 그림 6에 b=3,200일 때, 시간에 따른 탄도계수 추정을 보였다. Non-state 방법 역시, 그래프의 파동이 크지만 시간이 지남에 따라 실제값으로 잘 추정됨을 알 수 있다. 100번의 몬테카를로 방법으로 시뮬레이션한 결과, b=3,200일 때, 추정 탄도계수의 평균과 평균 제곱근 편차는 각각 3,200.0, 1.5037였고 b=1,800일 때, 추정 탄도계수의 평균과 평균 제곱근 편차는 각각 1,800.0, 0.3498였다. b=1,800일 때, 평균 제곱근 편차가 5 이하로 낮아지는 샘플링 개수는 1,350개였다.

그림 6. | Fig. 6. b=3,200일 때 탄도계수 추정 (4장) | Ballistic coefficient estimation when b=3,200.

Download Original Figure

평균 제곱근 편차로 보았을 때 non-state 방법은 state 방법과 비교해서 더 좋은 추정을 했다고 할 수 있다. 하지만, non-state 방법의 경우, 수렴하는 과정에서 큰 파동이 생김을 볼 수 있다. 그 이유는 SEZ 좌표계의 각 축에서의 시간에 따른 속도의 추정오차가 매우 낮았다고 하더라도 추정 속도의 방향은 실제 속도의 방향과 크게 달라질 수 있다. 이는 결과적으로 a_aero를 제대로 추정할 수 없게 만들고 큰 파동을 만든다. 수렴하는 과정에서도 같은 이유로 파동이 생기고, 결국 정확한 추정을 위해서는 많은 샘플링 개수, 즉 많은 시간이 필요하다. b=1,800일 때, 평균 제곱근 편차가 5 이하로 낮아지는 샘플링 개수도 non-state 방법이 state 방법보다 250개 정도 더 많았다. 이는 올바른 추정을 위해서 non-state 방법이 더 많은 시간이 필요하다는 것을 알려준다.

Ⅴ. 상미분방정식 파라미터 추정 방법

이번 장에서는 탄도계수를 추정하는 방법인 상미분방정식 파라미터 추정 방법을 새로 제안하고, 시뮬레이션 결과를 보였다. 발사체의 운동방정식 식 (11)은 상미분방정식이다. 최대우도추정을 사용하여 비용함수를 설정하고, 비용함수를 최소화하는 파라미터 x₀를 구했다. 탄도계수는 일정하다고 가정했기 때문에 x₀의 7번째 항인 x_0,7이 탄도계수 추정값이 된다.

발사체의 운동방정식 식 (26)은 $x ˙ (t) = f (x (t), t), x (t 0) = x 0$ 형식의 상미분방정식이다. 발사체의 실제 궤적을 만든 운동방정식은 랜덤함수가 아니므로 궤적은 x₀에 의존하게 된다. 다시 말해, 동일한 초기 상태 벡터 x₀로 만든 궤적은 모두 같다.

x = [x 1, x 2, x 3, x 4, x 5, x 6, x 7] T = [p s, p e, p z, v s, v e, v z, b] T x ˙ 1 (t) = x 4 (t), x ˙ 2 (t) = x 5 (t), x ˙ 3 (t) = x 6 (t) x ˙ 4 (t) = − 1 2 x 7 (t) ρ v 2 cos α cos β + 2 ω x 5 (t) sin L x ˙ 5 (t) = − 1 2 x 7 (t) ρ v 2 cos α sin β − 2 ω (x 4 (t) sin L + x 6 (t) cos L) x ˙ 6 (t) = − 1 2 x 7 (t) ρ v 2 sin α + 2 ω x 5 (t) cos L − g x ˙ 7 (t) = 0

(26)

x(t;x₀)을 초기 x가 x₀일 때 시간 t에서의 x라고 하자. 매 시간마다 레이다에는 발사체에 대한 3개의 측정치 z_t가 들어온다. 따라서, 시간에 따라 h(x(t;x₀))가 z_t와 가장 비슷할 때, 탄도계수를 잘 추정했다고 말할 수 있다. 이를 표현하는 비용함수는 최대우도추정 방법을 사용했고, 다음의 식 (27)과 같다^[9]. 함수 h의 식은 식 (17)–식 (20)과 같다. 이 비용함수를 최소화 시키는 최적화 문제는 식 (28)과 같다. 비용함수를 달라지게 하는 파라미터는 x₀이다. r_j는 관측치 노이즈의 분산을 의미한다. 이 최적화 문제의 해의 7번째 항 $x^0, 7$ 이 추정값이 된다.

L (x 0) = ∑ t ∑ j = r, θ, ϕ (z t, j − h j (x (t; x 0))) 2 2 r j

(27)

min L x 0 (x 0) = min x 0 ∑ t ∑ j = r, θ, ϕ (z t, j − h j (x (t; x 0))) 2 2 r j

(28)

시뮬레이션에 쓰인 실제 x는 x₀=(10000, 10000, 43000, 400, 400, 350, 3200 or 1800)으로 3장, 4장의 시뮬레이션과 동일하게 주었다. x₀의 상한과 하한은 확장 칼만 필터의 시뮬레이션에서 사용된 P_0|0의 각 대각항의 제곱근에 3배 값을 x₀에 더하고 뺀 값으로 했다. 이는 3장과 4장의 시뮬레이션 결과와 비교하기 위해 최대한 비슷한 시뮬레이션 환경을 주기 위함이다. 표준편차의 3배를 실제 값에 더하고 뺀 값 사이 구간에서 3장과 4장의 시뮬레이션 초기 추정값이 있을 확률이 99% 이상이기 때문이다. 관측치 노이즈는 다음과 같이 확장 칼만 필터 시뮬레이션과 동일한 값을 주었다.

r r = 10, r θ = 0.01, r ϕ = 0.01

(29)

6장의 시뮬레이션에서 푼 최적화 문제는 다음의 식과 같이 쓸 수 있다.

min x 0 L (x 0) = min x 0 ∑ t ∑ j = r, θ, ϕ (z t, j − h j (x (t; x 0))) 2 2 r j l b ≦ x 0 ≦ u b l b = (9905, 9905, 42905, 370, 370, 320, 2900 or 1500) u b = (10095, 10095, 43095, 430, 430, 380, 3500 or 2100)

(30)

그림 7은 b=3,200일 때, 상미분방정식 파라미터 추정 방법을 사용하여 추정된 탄도계수의 예시를 보여준다. 100번의 몬테카를로 방법으로 시뮬레이션 결과, b=3,200일 때, 추정 탄도계수의 평균과 평균 제곱근 편차는 각각 3,200.5, 1.1680였다. b=1,800일 때, 두 값은 1,800.3, 0.3344였다. 3장과 4장의 시뮬레이션 결과보다 평균 제곱근 편차가 훨씬 낮음을 확인할 수 있다. 그러나, 이 방법의 경우, 실시간 추적이 불가능하고, 샘플링한 시간만큼 모든 관측치를 받고 나서 추정이 가능하다는 단점이 있다. 그러나 확장 칼만 필터 방법도 정확한 추정까지 시간이 걸린다. b=1,800일 때, 평균 제곱근 편차가 5 이하로 낮아지는 샘플링 개수는 1,100개로 앞선 방법과 비교했을 때 작은 값이므로 이 경우에는 문제가 되지 않는다.

그림 7. | Fig. 7. b=3,200일 때 탄도계수 추정 예시 (5장) | Ballistic coefficient estimation when b=3,200.

Download Original Figure

Ⅵ. 결 론

본 논문에서는 발사체 분류를 위한 탄도계수 추정 방법을 제안하고, 그 방법들을 비교해 보았다. 제안한 방법은 확장 칼만 필터 non-state 방법, 상미분방정식 파라미터 추정 방법이다.

State 방법과 non-state 방법의 경우, 모두 확장 칼만 필터를 사용하여 탄도계수를 추정했다. 시뮬레이션 결과, 충분한 시간이 있을 때, non-state 방법이 state 방법보다 추정이 더 정확하다는 것을 확인할 수 있었다. 그러나, state 방법이 non-state 방법보다 추정 그래프가 안정적이었고, 올바른 추정까지 걸리는 시간도 더 적게 걸렸다. 상미분방정식 파라미터 추정의 경우, 성능이 앞선 두 확장 칼만 필터 방법들보다는 좋았지만, 실시간 추정이 불가능하다는 단점이 있었다. 하지만, 올바른 추정까지 걸리는 시간이 확장 칼만 필터 방법들에 비해서 느리지 않았다.

본 논문에서는 발사체를 분류하기 위해 탄도계수를 추정했지만, 만약 발사체의 파편에 대한 운동방정식을 올바르게 모델링한다면, 발사체와 파편의 탄도계수를 같은 방법을 통해 추정하고 분류할 수 있을 것이다. 앞으로는 속도에 따라 달라지는 탄도계수 등을 고려하여 운동방정식을 더 정밀하게 세워 탄도계수를 추정하는 연구가 진행 되어야 할 것이다.

Acknowledgements

본 연구는 한화시스템 협동과제로 수행되었음.

References

[1].

B. R. Bowman, “True satellite ballistic coefficient determination for HASDM,” in AIAA/AAS Astrodynamics Specialist Conference and Exhibit, Monterey, CA, Aug. 2002.

[2].

J. Sang, J. C. Bennett, and C. H. Smith, “Estimation of ballistic coefficients of low altitude debris objects from historical two line elements,” Advances in Space Research, vol. 52, no. 1, pp. 117-124, Jul. 2013.

[3].

R. Jesionowski, P. Zarchan, Comparison of Filtering Options for Ballistic Coefficient Estimation. Arlington, VA, Sparta, 1998.

[4].

P. Dodin, P. Minvielle, and J. P. Le Cadre, “Estimating the ballistic coefficient of a re-entry vehicle,” IET Radar, Sonar & Navigation, vol. 1, no. 3, pp. 173-183, Jun. 2007.

[5].

T. H. Kim, K. R. Moon, and T. L. Song, “Variable- structured interacting multiple model algorithm for the ballistic coefficient estimation of a re-entry ballistic target,” International Journal of Control, Automation and Systems, vol. 11, no. 6, pp. 1204-1213, 2013.

[6].

C. Jekeli, Inertial Navigation Systems with Geodetic Applications, Berlin, Germany, Walter de Gruyter, 2000.

[7].

J. J. Martin, Atmospheric Reentry an Introduction to Its Science and Engineering, Englewood Cliffs, NY, Prentice- Hall, 1966.

[8].

A. H. Jaswinski, Stochastic Processes and Filtering Theory, New York, NY, Academic Press, 1970.

[9].

M. Peifer, J. Timmer, "Parameter estimation in ordinary differential equations for biochemical processes using the method of multiple shooting," IET Systems Biology, vol. 1, no. 2, pp. 78-88, Mar. 2007.