감시정찰 (Surveillance and Reconnaissance) 영역 계산

한국군 최초의 정찰위성 1호가 2023년 12월 2일 발사에 성공하였다. 2024년에는 정찰위성 2호부터 4호까지 차례로 발사될 예정이라고 한다. 최근에는 전통적인 군사 영역은 물론이고 민간 영역에서도 우주 자산을 이용한 감시 및 정찰 활동이 빠르게 증가하고 있다.

참고로 감시(surveillance)와 정찰(reconnaissance)은 모두 정보 수집을 위한 활동이지만 시간과 임무의 구체성에서 차이가 있다. 감시는 주로 장기적으로 특정 지역이나 대상을 모니터링하여 정보를 수집하는 과정이다. 예를 들어 군사적인 목적으로 특정 지역을 지속적으로 관찰하는 것이 여기에 해당한다. 정찰은 일반적으로 특정 목적을 위한 한시적이고 전략적인 정보 수집에 중점을 둔다. 예를 들어 새로운 또는 알려지지 않은 지역의 상황을 파악하거나 적의 위치 등을 알아내는 것등이 이에 포함된다.

 

 

그렇다면 인공위성이 수행할 수 있는 감시정찰 영역은 어떻게 계산할 수 있을까.

먼저 인공위성의 고도가 큰 영향을 미칠 것이다. 고도가 높을 수록 인공위성이 수평선까지 '볼 수 있는' 시야 범위가 넓어지기 때문이다. 두번째는 인공위성 센서의 시야각(FOV, field of view)이다. 센서의 시야각이 넓을 수록 더 넓은 영역이 시야에 들어올 것이다. 실제 계산에는 특정 센서의 해상도, 대기 굴절, 타원형 지구의 곡률 등 세부적인 내용도 고려해야겠지만, 여기서는 원형 지구 가정하에 센서의 시야각과 인공위성의 고도를 통해 지표면에 닿는 센서의 범위를 이용하여 감시정찰 영역을 계산해보고자 한다.

 

 

위 그림에서 지구에 표시된 파란색 선은 궤도의 지상궤적(ground track)이다. A점은 인공위성의 위치를 지상에 투사시킨 점으로서 지표점(nadir)이다. 센서의 관측 중심선이 지표점을 향한다면 '센서가 볼 수 있는 영역'은 A점을 중심으로 한 노란색으로 빗금친 영역이다. 위성이 지구를 공전하고 있으므로 센서가 볼 수 있는 영역은 이동하면서 지상궤적을 중심으로 파란색 점선으로 표시된 일정 밴드폭을 형성할 것이다. 만약 센서의 관측 중심선이 약간 기울어져 C점을 향한다면 그 영역은 지상궤적에서 이격된 선을 중심으로 일정 밴드를 형성할 것이다.

이제 특정 싯점에서 형성되는 '센서가 볼 수 있는 영역', 즉 감시정찰 가능 영역을 계산해 보자.

 

 

위 그림에서 지구 표면에 두 개의 위치가 표시되어 있다. A점은 위성의 지표점(nadir)이고 B점은 위성 센서의 조준점이다. \(\lambda_{lon,A}, \lambda_{lon,B}\) 는 각각 A점과 B점의 경도, \(\lambda_{lat,A}, \lambda_{lat,B}\) 는 A점과 B점의 위도이다. 여기서 지구는 원형으로 가정한다. \(\Lambda\) 는 A점과 B점의 지상 거리(ground range)다. \(\Lambda\) 는 각도의 단위일 수도 있고 거리의 단위일 수도 있다. 지구를 원형으로 가정했으므로 둘 사이의 변환은 간단하다. 방위각 \(\beta\) 는 A점에서 본 B점의 방향을 나타낸다. 또는 위성 센서가 보고 있는 방향을 나타낸다고 볼 수 있다. 북쪽으로 기준으로 시계방향으로 측정한다.

여기서 목표는 위성의 고도와 센서의 시야각(FOV), 그리고 지표점(nadir)의 위치 좌표가 주어졌을 때, 지상 거리인 \(\Lambda\) 를 계산하는 것이다.

A점과 B점을 연결한 선을 포함하는 원을 기준으로 지구의 절단면을 그리면 아래 그림과 같다.

 

 

그림에서 \(R_{er}\) 은 지구의 반지름, \(h_{sat}\) 는 위성의 고도, \(r_{sat}\) 는 지구 중심에서 위성까지의 거리, \(\eta_{FOV}\) 는 센서의 시야각(FOV), H점은 지구 지평선의 위치, A점은 위성의 지표점(nadir), B점은 위성 센서의 조준점 또는 관심지점, \(\eta\) 는 위성의 고도를 기준으로한 관심지점의 조준각(boresight angle), \(\eta_{hor}\) 는 지구 수평선의 조준각, \(\Lambda_{hor}\) 는 A점에서 H점까지의 지상 거리 (각도), \(\Lambda_{FOV}\) 는 시야각 범위의 지상 거리 (각도), \(\Lambda\) 는 A점에서 B점까지의 지상 거리(각도)다.

\(\rho_{hor}\) 는 위성에서 지평선 지점까지의 거리, \(\rho\) 는 위성에서 관심지점까지의 거리고, \(\gamma\) 는 계산상의 편의를 위하여 도입한 각이다.

지구 수평선의 조준각 \(\eta_{hor}\) 와 A점에서 H점까지의 지상 거리 \(\Lambda_{hor}\) 는 위성의 고도 \(h_{sat}\) 와 지구의 반지름 \(R_{er}\) 을 이용하여 쉽게 계산할 수 있다.

 

\[ \begin{align} & \sin \eta_{hor} = \frac{R_{er}}{R_{er}+h_{sat}} \tag{1} \\ \\ & \cos \Lambda_{hor} = \frac{R_{er}}{R_{er}+h_{sat}} \end{align} \]

 

\(\eta_{hor}\) 는 정의상 \(\pm 90^0\) 를 넘지않는다. 지상 거리 \(\Lambda_{hor}\) 는 센서가 관측할 수 있는 최대 지상 거리를 의미한다. \(\rho_{hor}\) 는 \(\eta_{hor}\) 을 이용하여 다음과 같이 계산할 수 있다.

 

\[ \begin{align} \rho_{hor}=r_{sat} \cos \eta_{hor} =(R_{er}+h_{sat} ) \cos \eta_{hor} \tag{2} \end{align} \]

 

피타고라스 정리를 이용하여 \(\rho_{hor}\) 를 계산할 수도 있다.

 

\[ \begin{align} \rho_{hor} &= \sqrt{r_{sat}^2-R_{er}^2 } = \sqrt{(R_{er}+h_{sat} )^2-R_{er}^2 } \tag{3} \\ \\ &= \sqrt{2R_{er} h_{sat}+h_{sat}^2 } \end{align} \]

 

이제 여기서부터는 평면(plane)에서의 삼각함수 공식이 필요하다. 다음과 같은 삼각형이 있을 때, 세 각 \(A, B, C\) 와 마주보는 변의 길이 \(a, b, c\) 에 관한 사인법칙은 다음과 같다.

 

\[ \begin{align} \frac{a}{\sin⁡ A} = \frac{b}{ \sin B} = \frac{c}{ \sin C} \tag{4} \end{align} \]

 

 

제1코사인법칙은 다음과 같다.

 

\[ \begin{align} & a= b \cos C+ c \cos B \tag{5} \\ \\ & b=a \cos C+c \cos A \\ \\ & c=a \cos B+b \cos A \end{align} \]

 

\(\gamma\) 는 사인법칙 (4)를 이용하면 다음과 같이 계산할 수 있다.

 

\[ \begin{align} \sin \gamma = \frac{r_{sat} \sin \eta}{R_{er} } \tag{6} \end{align} \]

 

여기서 \(\eta= \frac{\eta_{FOV}}{2}\) 이다. 각도를 계산할 때 사분면을 고려하는 것이 중요한데, \(\gamma\) 는 정의상 \(90^0\) 이상의 값을 갖는다.

제1코사인법칙 (5)를 이용하면 \(\rho\) 를 계산할 수 있다.

 

\[ \begin{align} \rho =R_{er} \cos \gamma+r_{sat} \cos \eta \tag{7} \end{align} \]

 

\(\rho\) 를 구했으므로 마지막으로 사인법칙 (4)를 이용하여 \(\Lambda\) 를 계산할 수 있다.

 

\[ \begin{align} \sin \Lambda = \frac{ \rho \sin \eta}{R_{er} } \tag{8} \end{align} \]

 

그림을 보면 정의상 \(\Lambda\) 는 \(90^0\) 보다 작아야 한다.

지상 거리인 \(\Lambda\) 까지 구했으므로 지표점(nadir)에서 가능한 모든 방위각 (\(0 \le \beta \le 360^0\) )을 조사하면 위성 센서가 볼 수 있는 B점의 좌표를 계산할 수 있다.

 

 

B점의 좌표를 계산하기 위해서는 구면(sphere)에서의 삼각함수 공식이 필요하다. 다음과 같은 삼각형이 구면에 있을 때, 세 각 \(A, B, C\) 와 마주보는 각 \(a, b, c\) 에 관한 사인법칙은 다음과 같다.

 

\[ \begin{align} \frac{ \sin a}{\sin A} = \frac{ \sin b}{\sin B} = \frac{ \sin c}{ \sin C} \tag{9} \end{align} \]

 

 

코사인법칙은 다음과 같다.

 

\[ \begin{align} & \cos a= \cos b \cos c+ \sin b \sin c \cos A \tag{10} \\ \\ & \cos b= \cos c \cos a+ \sin c \sin a \cos B \\ \\ & \cos c= \cos a \cos b+ \sin a \sin b \cos C \end{align} \]

 

지상 거리 \(\Lambda\) 와 방위각 \(\beta\), 지표점(nadir)의 좌표, B점의 좌표의 기하학적인 관계는 다음 그림과 같다.

 

 

B점의 위도 \(\lambda_{lat,B}\) 는 코사인법칙 (10)에서 \(a \to (90^0-\lambda_{lat,B})\), \(b \to \Lambda\), \(c \to (90^0-\lambda_{lat,A})\), \(A \to \beta\) 로 놓으면 다음과 같이 계산할 수 있다.

 

\[ \begin{align} & \cos (90^0-\lambda_{lat,B} ) = \cos \Lambda \cos (90^0-\lambda_{lat,A} ) \\ & \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ + \sin \Lambda \sin (90^0-\lambda_{lat,A} ) \cos \beta \\ \\ \to \ & \sin \lambda_{lat,B} = \cos \Lambda \sin \lambda_{lat,A} + \sin \Lambda \cos \lambda_{lat,A} \cos \beta \tag{11} \end{align} \]

 

B점의 경도 \(\lambda_{lon,B}\) 를 알기 위해서는 먼저 \(\Delta \lambda_{lon}\) 를 구해야 한다. 사인법칙 (9)를 이용하면 \(\Delta \lambda_{lon}\) 를 다음과 같이 계산할 수 있다.

 

\[ \begin{align} & \frac{\sin \Lambda }{ \sin \Delta \lambda_{lon}}= \frac{\sin⁡ (90^0-\lambda_{lat,B}) }{ \sin \beta} \\ \\ \to \ & \sin \Delta \lambda_{lon} = \frac{ \sin \beta \sin \Lambda}{\cos \lambda_{lat,B} } \tag{12} \end{align} \]

 

또는 코사인법칙 (10) 에서 \(a \to \Lambda\), \(b \to (90^0-\lambda_{lat,B})\), \(c \to (90^0-\lambda_{lat,A})\), \(A \to \Delta \lambda_{lon}\) 로 놓으면 다음과 같이 계산할 수도 있다.

 

\[ \begin{align} & \cos \Lambda = \cos (90^0-\lambda_{lat,B} ) \cos (90^0-\lambda_{lat,A} ) \\ & \ \ \ \ \ \ \ \ \ \ \ \ \ \ + \sin (90^0-\lambda_{lat,B} ) \sin (90^0-\lambda_{lat,A} ) \cos \Delta \lambda_{lon} \\ \\ \to \ & \cos \Lambda = \sin \lambda_{lat,B} \sin \lambda_{lat,A} + \cos \lambda_{lat,B} \cos \lambda_{lat,A} \cos⁡ \Delta \lambda_{lon} \\ \\ \to \ & \cos \Delta \lambda_{lon} = \frac{ \cos \Lambda- \sin \lambda_{lat,B} \sin \lambda_{lat,A}}{ \cos \lambda_{lat,B} \cos \lambda_{lat,A} } \tag{13} \end{align} \]

 

식 (12) 또는 (13)을 통하여 \(\Delta \lambda_{lon}\) 를 계산한 후, B점의 경도 \(\lambda_{lon,B}\) 는 다음식으로 계산할 수 있다.

 

\[ \begin{align} \lambda_{lon,B} = \lambda_{lon,A}+ \Delta \lambda_{lon} \tag{14} \end{align} \]

 

만약, 센서의 관측 중심선이 지표점(nadir)을 향하지 않고 기울어져 있다면 전체적인 기하학적 형상이 비대칭이 되어서 약간 더 복잡한 계산이 필요하다.

 

 

관측 중심선의 관측 중심각(center boresight angle) \( \eta_{ctr}\) 이 주어지면 max 조준각과 min 조준각을 계산할 수 있다.

 

\[ \begin{align} \eta_{max} &= \eta_{ctr}+ \frac{\eta_{FOV}}{2} \tag{15} \\ \\ \eta_{min} &= \eta_{ctr} - \frac{\eta_{FOV}}{2} \end{align} \]

 

그런 다음, 식 (6), (7), (8)을 이용하여 차례로 \( (\gamma_{max}, \gamma_{min})\), \( ( \eta_{max}, \eta_{min})\), \( (\Lambda_{max}, \Lambda_{min})\) 을 계산하고 최대 및 최소 지상 거리 (각도)의 차이인 \(\Lambda =\frac{1}{2} (\Lambda_{max}-\Lambda_{min})\) 를 구하면 된다.

관측 중심선이 지면과 닿는 위치 C의 좌표 \((\lambda_{lat,ctr}, \lambda_{lon,ctr})\) 는 센서가 바라보는 방향각 \(\beta_{sens}\) 가 주어지면 식 (11), (12), (14)로 계산할 수 있다.

 

\[ \begin{align} & \sin \lambda_{lat,ctr} = \cos \Lambda_{ctr} \sin \lambda_{lat,A} + \sin \Lambda_{ctr} \cos \lambda_{lat,A} \cos \beta_{sens} \tag{16} \\ \\ & \sin \Delta \lambda_{lon,ctr} = \frac{ \sin \beta_{sens} \sin \Lambda_{ctr} }{ \cos \lambda_{lat,ctr} } \\ \\ & \lambda_{lon,ctr} = \lambda_{lon,A} + \Delta \lambda_{lon,ctr} \end{align} \]

 

여기서 \(\Lambda_{ctr}\) 은 지표점(nadir)에서 C점까지의 지상 거리로서 식 (8)에 의해서

 

\[ \begin{align} \sin \Lambda_{ctr} = \frac{ \rho_{ctr} \sin \eta_{ctr}}{R_{er}} \tag{17} \end{align} \]

 

이다. 여기서 \(\rho_{ctr}\) 은 식 (7)에 의해서

 

\[ \begin{align} \rho_{ctr}=R_{er} \cos \gamma_{ctr} +r_{sat} \cos \eta_{ctr} \tag{18} \end{align} \]

 

이고, \(\gamma_{ctr}\) 은 식 (6)에 의헤서

 

\[ \begin{align} \sin \gamma_{ctr} = \frac{ r_{sat} \sin \eta_{ctr} }{R_{er}} \tag{19} \end{align} \]

 

이다.

지금까지의 논의는 지구가 원형 구체라는 가정하에서 진행됐지만 실제로 지구는 적도면이 부풀어 오른 타원형 구체이므로 감시정찰 영역 계산에 있어서 약간의 오차가 내재되어 있다는 것을 알아야 한다.

정리하면 다음과 같다.

위성의 고도, 센서의 시야각이 주어졌을 때, 지상 거리는 다음과 같은 순서로 구할 수 있다.

   [1] 식 (1)로 \(\eta_{hor}\) 와 \(\Lambda_{hor}\) 을 계산한다.
   [2] 식 (2)로 \(\rho_{hor}\) 를 계산한다.
   [3] 식 (6)으로 \(\gamma\) 를 계산한다.
   [4] 식 (7)로 \(\rho\) 를 계산한다.
   [5] 식 (8)로 \(\Lambda\) 를 계산한다.
   [6] 가능한 방위각 (\(0 \le \beta \le 360^0\)) 범위안에서 식 (11), (12), (14)로 감시정찰 영역을 계산한다.

이번에는 위성의 지표점(nadir)의 좌표, 고도, 센서의 시야각, 센서의 관측 중심각, 센서가 바라보는 방향각이 주어졌을 때, 지상 거리와 관측 중심의 지상 좌표는 다음과 같은 순서로 구할 수 있다.

   [1] 식 (1)로 \(\eta_{hor}\) 와 \(\Lambda_{hor}\) 을 계산한다.
   [2] 식 (2)로 \(\rho_{hor}\) 를 계산한다.
   [3] 식 (15)로 \((\eta_{max}, \eta_{min})\) 을 계산한다.
   [4] 식 (6), (7), (8)로 \((\gamma_{max}, \gamma_{min})\), \((\eta_{max}, \eta_{min})\), \((\Lambda_{max}, \Lambda_{min})\) 을 계산한다.
   [5] \((\Lambda_{max}, \Lambda_{min})\) 로 \( \Lambda\) 를 계산한다.
   [6] 식 (19)로 \(\gamma_{ctr}\) 를 계산한다.
   [7] 식 (18)로 \(\rho_{ctr}\) 를 계산한다.
   [8] 식 (16)으로 관측 중심의 지상 좌표를 계산한다.
   [9] 가능한 방위각 \((0 \le \beta \le 360^0)\) 범위안에서 식 (11), (12), (14)로 감시정찰 영역을 계산한다. 이 때 위성의 지표점(nadir)의 좌표 대신에 관측 중심의 지상 좌표를 사용한다.

 

 

다음은 위 알고리즘을 매트랩 코드로 구현한 것이다.

 

surveillance.m

function [Lambda_hor, rho_hor, Lambda, lat_ctr, lon_ctr] = ...
    surveillance(latA, lonA, h, fov, eta_ctr, beta_sens)

% surveillance ground range computation
% input:
%       latA, lonA, h: latitude-longitude-height of nadir in (deg, km)
%       fov: sensor field of view in deg
%       eta_ctr: center boresight angle in deg
%       beta_sens: sensor azimuth look-angle in deg
% output:
%       Lambda_hor: maximum ground range in deg
%       rho_hor: maximum slant range in km 
%       Lambda: ground range in deg
%       lat_ctr, lon_ctr: latitude-longitude of center boresight in deg
%
% coded by st.watermelon

% constants
R_er = 6378.1363; % km

if nargin == 4
    eta_ctr = 0; 
    beta_sens = 0;
elseif nargin == 5
    beta_sens = 0;
end


% deg to rad
fov = fov*pi/180;
eta_ctr = eta_ctr*pi/180;

latA = latA*pi/180;
lonA = lonA*pi/180;

beta_sens = beta_sens*pi/180;


% horizon-boresight angle, eta_horizon
r_sat = R_er + h;
eta_hor = asin(R_er/r_sat);

% horizon-ground-range, Lambda_hor
Lambda_hor = acos(R_er/r_sat) * 180/pi;

% horizon slant range, rho_hor
rho_hor = r_sat*cos(eta_hor);


if eta_ctr == 0 % if sensor points nadir

    % intermediate angle, gamma (>= pi/2)
    eta = fov/2;
    gamma = asin(r_sat*sin(eta)/R_er);
    if gamma < pi/2
        gamma = pi-gamma;
    end

    % slant range, rho
    rho = R_er*cos(gamma) + r_sat*cos(eta);

    % ground-range, Lambda (< pi/2)
    Lambda = asin(rho*sin(eta)/R_er); 
    if Lambda >= pi/2
        Lambda = pi-Lambda;
    end 
    lat_ctr = latA*180/pi;
    lon_ctr = lonA*180/pi;
    Lambda = Lambda*180/pi;

else % if not

    % max and min boresight angle, eta_max, eta_min
    eta_max = eta_ctr + fov/2;
    eta_min = eta_ctr - fov/2;

    % intermediate angle, gamma (>= pi/2)
    gamma_max = asin(r_sat*sin(eta_max)/R_er);
    if gamma_max < pi/2
        gamma_max = pi-gamma_max;
    end

    gamma_min = asin(r_sat*sin(eta_min)/R_er);
    if gamma_min < pi/2
        gamma_min = pi-gamma_min;
    end

    % slant range, rho_max, rho_min
    rho_max = R_er*cos(gamma_max) + r_sat*cos(eta_max);
    rho_min = R_er*cos(gamma_min) + r_sat*cos(eta_min);

    % ground-range, Lambda_max, Lambda_min (< pi/2)
    Lambda_max = asin(rho_max*sin(eta_max)/R_er); 
    if Lambda_max >= pi/2
        Lambda_max = pi-Lambda_max;
    end 

    Lambda_min = asin(rho_min*sin(eta_min)/R_er); 
    if Lambda_min >= pi/2
        Lambda_min = pi-Lambda_min;
    end 

    % total ground-range, Lambda_fov 
    Lambda_fov = Lambda_max - Lambda_min;
    Lambda = Lambda_fov/2;

    % ceteral line position
    gamma_ctr = asin(r_sat*sin(eta_ctr)/R_er);
    if gamma_ctr < pi/2
        gamma_ctr = pi-gamma_ctr;
    end
    rho_ctr = R_er*cos(gamma_ctr) + r_sat*cos(eta_ctr);
    Lambda_ctr = asin(rho_ctr*sin(eta_ctr)/R_er); 
    if Lambda_ctr >= pi/2
        Lambda_ctr = pi-Lambda_ctr;
    end 
    
    lat_ctr = asin(cos(Lambda_ctr)*sin(latA) ...
                +sin(Lambda_ctr)*cos(latA)*cos(beta_sens));

    del_lon_ctr = asin(sin(beta_sens)*sin(Lambda_ctr)/cos(lat_ctr));
    lon_ctr = lonA+del_lon_ctr;

    lat_ctr = lat_ctr*180/pi;
    lon_ctr = lon_ctr*180/pi;
    Lambda = Lambda*180/pi;

end