본 연구의 대상 지역은 한반도를 포함하는 북위 33°~38.5°, 동경 126°~130°에 해당하는 지역이며, 2020년 한 해 동안의 GK2A 위성관측 자료를 이용하였다 (그림 1). 한반도 지역은 산악 지형과 바다가 둘러싸여 있어 해양성 기후와 산지성 기후가 복잡하게 영향을 미치며, 계절적 기후 경향이 나타나므로 온도 변화의 다양성 및 주기성이 함께 나타나는 지역이다.

Fig. 1. 
Map of the study area. Rectangular boundaries for East Asian region (24.0°N~49.0°N, 112.0°E~148.0°E) and south Korean region (33.0°N~38.5°N, 126.0°E~130.0°E) are shown.

본 연구에서는 표 1과 같이 GK2A 위성자료와 GK2A 위성자료 처리를 위하여 필요한 보조 입력자료로서 MODIS 지표정보 자료와 재분석 모델자료를 사용하였다. GK2A 기반의 온도 산출 알고리즘의 실행을 위하여 주요 입력자료인 GK2A 위성 관측 자료와 보조 입력자료로서 MODIS emissivity, MODIS landcover, 그리고 지표 물리특성 정보를 포함한 재분석 모델자료인 Global Land Data Assimilation System (GLDAS) 자료가 사용되었다. GK2A 기반의 온도 산출물의 비교검증을 위하여 MODIS Land Surface Temperature (LST) 자료도 사용되었다. 개별 자료에 대한 상세설명은 아래와 같다.

한반도를 관측하는 정지궤도 위성 중 GK2A는 기상 및 환경 모니터링을 위한 목적으로 운영되고 있으며, 한반도 및 주변 지역에 대하여 실시간으로 관측하고 모니터링하는 데 사용되고 있다. GK2A 위성은 지구에서 약 45000 km 떨어진 정지궤도관측을 통해 대기오염물질의 현황과 장거리 이동을 관측할 수 있다. GK2A의 탑재 센서인 Advanced Meteorological Imager (AMI)는 수동형 원격탐사 기기로서 가시광선부터 적외선 영역의 복사에너지를 관측한다. AMI의 가시광선 영역의 채널은 지구로부터의 반사된 태양 복사를 관측하며, 태양광이 없는 야간이나 대상의 열적 변화 정보를 구별하는 능력은 없다. 반면에 적외선 영역의 관측 채널은 지구 복사를 측정하여 대상의 방사 프로파일 정보와 열적 변화 정보를 제공할 수 있으나 대기중의 장파장 에너지에 민감하지 못한 대상은 탐지가 어렵다.

본 연구에서 사용된 인공위성 자료는 국산 정지궤도위성인 GK2A 위성의 AMI 원시 관측 자료를 복사 보정 (radiometric calibration)한 자료인 Level-1B (L1B) 자료를 사용하였다. GK2A 위성의 관측영역은 확장된 동아시아 (Extended East Asia) 모드이며, 본 연구의 대상 지역인 한반도 주변 영역만 추출하여 사용하였다. GK2A AMI 센서는 16개 채널 (중심 파장 0.47 μm, 0.51 μm, 0.64 μm, 0.86 μm, 1.38 μm, 1.6 μm, 3.9 μm, 6.2 μm, 6.9 μm, 7.3 μm, 8.6 μm, 9.6 μm, 10.4 μm, 11.2 μm, 12.4 μm, 13.3 μm)에서 지구 관측을 수행하며, 약 10분 간격으로 전구 영역을 관측한다. GK2A AMI의 공간 해상도는 0.470 μm~0.856 μm 파장 범위의 가시 채널에서 1 km 공간 분해능을 가지고 있고, 0.640 μm 채널은 500 m의 공간 분해능, 적외 채널에서는 2 km 공간 분해능을 가지고 있다. 그리고 최대 2.5분에 한 번씩 관측하고 있으며 동아시아 영역에서는 10분 간경으로 관측자료를 제공한다. 기존의 정지궤도 위성이 약 2 km~4 km 정도의 공간 해상도를 가지는 3~5개의 채널 정보를 제공하였으나, GK2A AMI는 상대적으로 고해상도의 다중채널 영상정보를 제공하여 보다 상세한 수준의 지구환경 관측정보를 획득할 수 있게 되었다.

MODIS 지표 방사율 자료 (코드명: MOD11C3, 버전 6.1)는 적외선 채널 관측을 기반으로 만들어졌으며, 지구의 표면 온도와 방사율을 산출하기 위해 다양한 알고리즘과 물리적 모델링을 사용한다 (Wan, 1999). 이 자료는 전 지구를 격자화하여 한 격자당 0.05° (적도 기준 5600 m)를 가지는 월별 지표 방사율값을 제공한다. 지표의 방사율 산출을 위하여 흑체 (black body)로 가정된 지표면에 대한 온도와 열복사 에너지와의 관계를 해석하여 지표의 방사율을 계산한다 (Wolfe et al., 1998).

MODIS 토지피복 자료 (코드명: MCD12C1)는 지구의 지표를 여러 유형 (예: 숲, 초원, 도시, 농경지 등)으로 분류하여 유사한 피복들을 총 7개의 클래스로 나누어 결정한다. 전 세계에 대해 0.05° (약 5600 m) 공간 해상도를 가지는 각 토지 피복 정보를 포함하고 있다. MCD12C1 데이터는 1년에 한 번 측정된 지표 유형 데이터로서 제공되므로 연도별 지표 유형 변화를 파악할 수 있다 (Strahler et al., 1999).

전 지구 토지 데이터 동화 시스템 (Global Land Data Assimilation System; GLDAS)의 목표는 첨단 지표면 모델링과 데이터 동화 기법을 사용하여 위성 및 지상 기반 관측 데이터 제품을 수집하여 최적의 지표면 상태와 플럭스 필드를 생성하는 것이다 (Rodell et al., 2004). 토지 정보 시스템 (Land Information Sys-tem; LIS)과 여러 개의 지표면 모델을 구동하고 방대한 양의 관측 기반 데이터를 통합하며, 전 세계 규모의 격자화된 데이터를 생성한다. GLDAS 자료는 하루 동안 3시간 간격으로 8회 (00, 03, 06, 09, 12, 16, 18, 21 UTC) 생산되고 있으며, 공간 해상도는 경위도 0.25°의 등간격으로 총 경도별 1440개, 위도별 721개의 격자 단위로 구성된 자료를 제공한다.

MODIS LST 자료 (코드명 MOD11)는 지표 온도 산출물이 저장되어 있으며, 공간 해상도 1 km를 가진다 (Wan, 1999). MODIS LST 산출물은 적외선 채널에서 방사율을 활용하여 흑체 방사율과 비교를 통해 지구 표면의 온도를 추정하는 수치해석적 계산 과정을 거쳐서 지구 표면의 온도를 산출한다. 이 자료는 GK2A 지표 온도 산출 결과물과의 비교를 통하여 결과물의 오차범위 수준을 분석하였다.

본 연구에서 제안한 GK2A 온도 산출 알고리즘은 1) 입력자료 준비 및 처리 과정, 2) 청천 화소 탐지, 3) 지표 온도 산출 및 4) 대기 온도 산출 과정으로 구성된다. 그림 2는 지표 및 대기 온도 산출을 위한 자료처리 흐름도를 나타낸다. 입력자료 준비 및 처리 과정에서는 GK2A AMI L1b 자료를 주요 입력자료로 읽은 후, 사용자가 지정한 위도와 경도에 해당하는 영역에 대하여 잘라내어 자료처리를 위한 자원을 줄인다. 이때 GK2A 위성 자료뿐 아니라, 표 1에서 나열된 보조 입력자료에 대해서도 동일한 영역으로 추출된 자료를 재구성하여 사용하였다. 그리고 이렇게 준비된 입력자로는 GK2A AMI L1 파일명으로부터 읽은 관측 시간과 각 화소별 기하정보 자료를 이용하여 태양 천정각 (Sun zenith angle)과 태양 방위각 (sun azimuth angle)을 계산하였다. 태양의 위치각도는 주야간 경계값 설정, 반사도값 보정, 태양 거울 반사각 (sunglint angle) 계산에 사용되므로 정확한 수치가 필요하다.

Fig. 2. 
Flow chart of GK2A land surface and air temperature retrieval process.

정지궤도 위성인 GK2A의 청천 화소 탐지는 여러 가지의 단위 테스트를 적용하여 해당 화소가 구름 영역인지 여부를 판단한다. 구름탐지 테스트는 먼저 각 화소별 가시채널 반사도와 적외채널 방사율을 사용하여 일련의 임계값 테스트 (threshold test)를 수행한다. 구름탐지 테스트에서 사용되는 임계값은 각 화소의 주요 지역 분류 또는 지표면 유형에 따라 달라진다. 이러한 GK2A 구름탐지의 목적은 구름이나 광학적으로 두꺼운 에어로솔이 주어진 화소 위치에서 지구 표면을 가리는지 여부를 결정하는 것이다. 구름 탐지 테스트는 기존의 AVHRR (Stowe et al., 1995) 및 MODIS (Ackerman et al., 1997)에서 사용되어 신뢰성이 검증된 테스트에 대하여 GK2A에 적용하였다. 그리고 개별 구름탐지 테스트 결과 (Gi)를 모두 곱한 값은 청천 화소의 신뢰도 값을 결정하는 데 사용되며 (식 1), 구름이 존재하는 화소에 대한 레벨 플래그로 정의된다. 식 (1)에서 사용된 구름탐지 마스크(Cloud Mask; CM)에 대한 접근방식은 해당 화소에 적용된 여러 탐지 테스트 결과를 통해 구름에 대한 신뢰도값을 정량적으로 결정하는 것이다.

GK2A/AMI 센서의 16개 채널 중 10.4 μm 채널은 대기의 창 영역으로서 대기투과도에 덜 민감하며, 청천 화소 조건에서 지표의 온도 감지에 유용하다. 청천 화소 조건에서는 대기 중 산란 효과가 매우 적고 지표면은 흑체로 가정할 수 있다. 이 파장대의 적외선 채널에서 복사 전달 방정식 (Radiative Transfer Equation; RTE)의 지표와 대기 중의 미량기체에 의한 방출량이 주요 요소값이다. 따라서, 위성에 도달하는 적외 영역의 복사 휘도는 지표와 대기 방출에 관한 합으로 다음과 같은 식 (2)로 표현될 수 있다.

위 식에서는 위성 방향과 지표면 방향으로 향하는 복사를 각각 설명하는 두 개의 항 (L_u 및 L_d)이 포함된다. 여기서 지표면 수준에서 복사량은 하향 복사와 지표면의 반사율을 곱한 것과 (지표의 투과도는 0으로 가정하면, 지표반사도: r_g=1-ε_g) 지표에서 방출되는 복사량의 합이다. 그리고 지표에서 방출되어 지구 대기로 향하는 복사량은 대기투과도 (t_a)에 의해 줄어든다. 청천 화소 조건에서 대기투과도와 상향 복사량은 무시할 만한 수준이므로, 지표 수준에서의 밝기온도값 (BT_g)은 방출률의 함수로 다음과 같은 식 (3)과 같이 정리된다.

지표 근처의 대기 온도 산출 이론은 에너지 균형 모델에 의하여 설명될 수 있다. 즉, 지표면에서 흡수된 순복사 (Rn)는 다음의 식 (4)와 같이 지열 (G), 현열 (H), 잠열 (LE)의 유출 플럭스의 균형을 유지해야 한다.

순복사는 지표로 입사하는 태양 복사와 방출되는 지구 장파 복사의 합계이며, 대기로부터 지표면으로 전달되는 복사에너지량과 지표면으로부터 다시 대기로 방출되는 양의 차이이다. 그리고 지하층의 온도가 일반적으로 지표 온도보다 낮기 때문에 낮 동안 지표에 가열된 열량은 지하로 전달된다. 이렇게 낮 동안 지하에 저장된 열량은 밤에 지표가 식으면 지표와 대기로 전달된다. 잠열의 경우, 주로 식물의 증산과 지표면의 증발로 인하여 생성된 물에 의하여 생성된다. 그리고 이 과정에서 남은 에너지량 (=Rn-G-LE)은 지표면 온도를 높이고 지표면에서 대기 중의 현열으로 전달된다. 현열 플럭스와 잠열 플럭스는 LST와 지표 기온의 함수로 간주할 수 있으며, 토양 열 플럭스는 순 복사량으로 표현할 수 있다. 따라서, 각 항목에 대한 물리량이 확보된 후, 지표 부근의 대기 온도에 관한 값을 역으로 산출한다.

청천 화소 결정을 위하여 사용되는 구름탐지 기법은 주로 주간 동안의 구름의 영향이 최소화된 화소를 찾기 위하여 구름 화소나 구름과 유사한 영역에 대한 강한 제거 방법이 사용되나, 이러한 방법론은 주야간이 바뀌는 여명기 (또는 황혼기)나 야간에 적용이 어렵게 된다. 따라서, 본 연구에서 사용된 구름탐지 방법은 주야간 경계에서의 구름탐지 불연속성 최소화와, 육상/해상 경계에 대한 구름탐지 불연속성 최소화, 그리고 AMI 관측 조건에 최적화된 연속 구름탐지 수준을 제시함으로써, 사용자에게 선택가능한 청천 화소 수준을 제공할 수 있도록 설계되었다.

그림 3은 에어로솔 산출을 위하여 사용된 구름탐지 결과와 본 연구 결과에 대해서 주간 (2020/04/05 00:00 UTC)과 야간 (2020/04/05 09:00 UTC) 사례를 비교한 사례이다. 먼저 청천 화소 탐지 산출물은 주로 가시채널 반사도와 적외채널 밝기온도 경계값 테스트가 주야간 경계조건 없이 사용되며, 공간균일도 테스트와 같은 통계적 유사도 테스트가 사용되어 매우 가혹한 조건하에서 구름 화소로 의심되는 모든 화소를 제거한다. 반면에 본 연구에서 사용된 다양한 구름테스트는 주간과 야간, 육지와 해양 등의 처리 프로세스가 나뉘어져 있으므로 연속탐지가 가능하다. 이러한 차이점은 가시채널 반사도 테스트 적용 시 태양고도각이 낮은 부분에서 과대탐지를 하는 것과 주야간 경계면에서 잡음신호가 발생하는 것, 그리고 상층운 탐지를 위하여 사용되는 1.38 μm 채널의 경계값 테스트 경우, 기존의 알려진 경계값을 주야간 경계 없이 사용할 경우 오류가 발생함을 확인하였다. 특히 여명기가 발생하는 영역에서 큰 오류가 나타난다. 이것은 아주 낮은 태양고도각에서 반사도값이 크게 증가하므로 나타나는 현상이므로 여명기에 대한 경계값 조정을 통하여 보완이 가능하였다.

Fig. 3. 
Comparison of cloud detection score results generated by applying the cloud detection test. Cloud mask values used for aerosol calculation (left) and used in this study (right) for daytime at 2020/04/05 00:00 UTC (top) and nighttime at 09:00 UTC (bottom).

지표에서 방사되는 복사량은 최초 플랑크 함수로 표현될 수 있으며, 대기투과도가 1에 가까운 경우 (투명한 대기), 대기 상단 (Top of Atmosphere; TOA)과 지표에서의 밝기온도값은 거의 동일한 값을 가지게 된다 (식 2). 그러나 대기 조성물질에 따른 대기투과도에 따라 복사량이 달라지게 되므로, TOA에서의 복사량에 대해 대기투과도에 의한 보정을 실시하면 지표 수준에서 방사되는 복사량으로 간주할 수 있다. 그리고 이렇게 대기투과도가 보정된 복사량은 다시 지표 방사율을 보정하여 지표면에서의 열량을 가정한 지표면 온도값으로 환산할 수 있다.

GK2A 지표 온도 산출 알고리즘을 이용하여 산출된 지표 온도 산출 결과 사례는 그림 5와 같다. 이 결과에서는 GK2A 기반의 주간 TOA와 지표 수준의 밝기온도 산출 결과 사례를 나타내며, TOA에서의 밝기온도는 지표에서 방출된 적외 복사량이 대기 중에서 감쇄된 양과 다시 대기 조성물질에 의하여 방사된 양이 더해져 센서가 감지한 결과이다. 따라서 지구표면과 태양의 기여도가 다소 포함될 수 있다. 지표 온도는 대기에 의하여 감쇄된 양과 지표 방출률을 고려하여 보정된 값으로 TOA 밝기온도보다는 약간의 증가된 온도값을 확인할 수 있다. 그림 4에서 주간 동안 TOA에서의 평균 밝기온도값은 289±5.0 K이며, 지표의 평균 밝기온도는 290±5.7 K로서 약 1 K의 차이가 발생한다.

Fig. 4. 
GK2A observation based BT (unit=Kelvin) at TOA (left) and surface (right) levels during daytime.

그림 5에서는 GK2A 기반의 야간 TOA와 지표 수준의 밝기온도 산출 결과 사례를 나타내며, 태양의 간섭이 없는 상태이므로 주간에 비하여 온도가 낮다. TOA와 지표 모두 지표에서 방출된 적외복사량이 주요 에너지원이며, 대기투과도와 방출률의 영향은 TOA 온도에서만 영향을 받는다. 주간과 마찬가지로 야간에서도 TOA에서의 평균 밝기온도값은 277±5.3 K이며, 지표의 평균 온도는 279±5.0 K로서 약 2 K의 차이가 발생한다.

Fig. 5. 
GK2A observation based BT (unit=Kelvin) at TOA (left) and surface (right) levels during nighttime.

그림 6은 GK2A 지표 온도 산출 결과와 NASA MODIS L2 LST 자료와 비교한 산점도이다. MODIS L2 LST 자료는 코드명 MOD11 또는 MYD11으로 산출된 버전 6.1의 최신 업데이트된 산출알고리즘이 적용된 자료이다. MODIS LST는 토지피복분류, 지표 방출율, 식생, 기상 모델 등 다양한 외부 입력자료를 이용하여 생산된 지표 온도자료이다. GK2A 기반의 지표 온도값과 MODIS LST와의 비교 결과는 선형회귀 직선의 기울기 1.15, MSE=2.30 K, BIAS=2.08 K, RMSE=2.88 K 범위 수준의 오차범위를 나타내었다.

Fig. 6. 
Comparison of GK2A-based surface temperature with MODIS L2 LST.

현열 플럭스 및 잠열 플럭스는 벌크 저항 접근법을 사용하여 구할 수 있으며, 현열 플럭스 H는 다음과 같이 계산된다.

여기서 ρ는 공기의 밀도 (kg/m³), C_p는 일정한 압력에서 공기의 비열 (J/kg K), T_s는 지표면 온도 (K), r_a는 공기 역학적 저항 (s/m)이다. r_a는 다음 식을 사용하여 계산한다 (Brutsaert, 1982):

여기서 z_u 및 z_t는 풍속 u 및 대기 온도가 측정되는 높이이고, d₀은 변위 높이이며, z_0M 및 z_0H는 각각 운동량과 열 수송에 대한 거칠기 길이이다. 모든 측정값은 미터 (m) 단위이며, 또한 Ψ_M과 Ψ_H은 운동량과 열에 대한 안정성 보정 함수로서, Monin-Obukhove length (Brutsaert, 1982)에 따라 달라지고, k는 칼만상수 (=0.4)이다. 거칠기 길이는 현장에서 관측값을 사용하는 것이 가장 이상적이나, 위성의 공간 해상도에 따른 불확성이 존재하므로, 각 토지 피복별로 분류된 문헌값이 사용되었다 (Kondo, 2000, 1994; Hansen, 1993; Brutsaert, 1982). 위성 영상에서 각 화소별 거칠기 길이는 다음과 같이 로그 평균화 절차 (Noilhan and Lacarrère, 1995)에 의해 조정된다.

여기서 z_0,i는 지표면 토지피복 유형에 대한 거칠기 길이이고, n은 유효 화소 수 (3 by 3 화소에 해당하는 9개)이다. 여기서 사용되는 거칠기 길이는 일반적으로 사용되는 토지 피복별 대표값을 사용하였다. 여기서 z_u와 z_t 및 복잡한 지형 영역의 d₀은 DEM 데이터를 사용하며, 관심 지역 내 지형 높이의 합과 지역 최대 높이와 지역 평균 높이의 차이 평균 높이의 차이는 z_u와 z_t로 사용된다. d₀의 값은 Macdonald et al. (1988)이 제안한 모델을 사용한다.

여기서 z_h는 장애물의 높이, A는 상수 (=4.43), λP는 장애물의 평면 면적 밀도, 즉 전체 면적에 대한 지형 장애물 면적의 비율이다. 이 방정식을 DEM 데이터에 적용하기 위해 z_h는 최대 높이와 최소 높이의 차이로 근사화한다. 또한, 총 픽셀 수에 대한 지역 평균 높이보다 높은 고도를 가진 픽셀 수의 비율 (즉, x /9)이 λP 대신 사용된다. 이상의 모델식을 따르면, 지표 온도 T_s와 현열 플럭스 H를 알면 T_a를 역으로 계산할 수 있다.

따라서, GK2A 대기 온도 산출 알고리즘은 식 (9)를 풀기 위하여 필요한 관측값이나 모델값을 이용하여 대기 온도를 산출하기 위하여 설계되었다. 여기서 필요한 관측값은 지표 온도로서 3.2절에서 설명된 GK2A 지표 온도 산출물을 사용하였다. 모델값은 현열플럭스 H이며, 시공간적으로 규격화된 정보를 얻기 위하여 NASA에서 제공하고 있는 재분석 자료인 GLDAS로부터 위성관측 지점과 동일한 위치에 해당하는 현열플럭스 값을 획득하였다. 그리고 지표 거칠기와 공기역학적 저항값은 각 화소별 지리적 형태값을 이용하여 계산하였다.

그림 7은 2020년 4월 15일 00UTC에 GK2A 대기 온도 산출 알고리즘을 이용하여 산출된 결과를 나타낸다. 칼라합성 영상에서는 동아시아 영역에서 육지와 해양에 위치한 구름에 관한 위치를 가시적으로 확인할 수 있으며, GK2A 연속구름탐지 결과를 통하여 구름이 있는 위치에 대한 마스크를 생성하였다. 그리고 지표 온도 산출을 위하여 청천 화소 정보가 사용되며, 산출된 지표 온도는 다시 대기 온도 산출을 위하여 사용되었다. 이상의 과정을 통하여 구름탐지, 지표 온도, 대기 온도 산출과정은 유기적으로 연결되어 있으며, 지역별 온도변화 분포정보를 위성관측 시간별로 생성할 수 있다. 그러나, GK2A 대기 온도 산출 결과에서는 입력자료로 사용된 GK2A 위성자료의 지표 온도 해상도 (2 km), GLDAS 현열플럭스의 공간 해상도 (0.25°)로서 상이하고 화소와 격자에 대한 위치정보가 정확히 일치하지 않기 때문에 해안선 근처에서 과대평가된 대기 온도가 산출됨을 알 수 있다. 이렇게 과대평가된 지점을 대기 온도 산출에서 제외시키기 위하여 지표 온도와 GLDAS 현열플럭스가 공간적으로 일치하여 존재하는 지점값만을 사용하여 대기 온도를 산출하였다.

Fig. 7. 
Example of (a) RGB color composite image, (b) cloud mask, (c) land surface termperature, and (d) air temperature images from the GK2A data measured on 2020/04/15/0000UTC.

그림 8에서는 GK2A와 GLDAS내부에 포함된 재분석 자료 중 대기 온도를 비교한 결과이다. 두 자료의 비교 결과는 BIAS=3.33 K, RMSE=5.3 K 정도 수준의 편차를 가지는 것으로 나타났다. 지표 온도의 비교사례에 비하여 오차범위가 다소 높은 수준이나, 제한된 리소스를 이용하여 대기 온도가 산출가능함을 확인하였으며 정확도 평가를 위한 방법론과 결과 제시라는 의의를 가진다. 또한, 이러한 결과는 정지궤도위성 관측값을 이용하여 지표 근처의 대기 온도 산출을 위한 최초의 결과이며, 정확도 개선을 위하여 대기 온도 산출을 위하여 사용된 모델 계산식의 입력값들의 개선이 요구될 것이다.

Fig. 8. 
Scatter plot of GK2A air temperature versus GLDAS air temperature values.