1. 서론

백두산은 900년대부터 꾸준히 화산활동을 기록한 것으로 알려져 있다. 백두산은 1903년 발생한 소규모의 분화를 마지막으로 화산활동에 대한 징후가 발생하지 않은 것으로 보고되었다(Horn and Schmincke, 2000). 그러나 2000년대 들어서부터 백두산의 분화 전조현상의 징후가 나타나기 시작했다. 보고된 화산 전조현상들은 화산성 지진활동, 지표의 변위, 온천수 온도 증가 및 화학성분 변화 등이 있다(Xu et al., 2012; Wei et al., 2013; Yun and Lee, 2012). 백두산의 불안정한 상태로 인한 분화가능성이 증가하면서 백두산의 분화 전조현상 연구가 지속적으로 진행되고 있다. 특히, Xu et al. (2012)는 천지 아래에 존재하는 마그마방으로 인해 밑바닥에서 뜨거운 온천과 기체가 분출되고 있으며 화산활동에 의해 마그마방이 변화하여 온천수의 온도가 증가하며, 화학성분이 변화하는 것을 분석한 바 있다. 이로 인해 백두산 백운봉 아래 천지호 지역은 다른 천지호 지역에 비해 늦은 시기에 얼며, 지표온도 또한 높게 나타나는 이상지역이 존재한다(KMA, 2021). 따라서, 해당 이상지역이 온천 분출 지역일 가능성이 크며, 천지가 결빙 되어 있는 겨울시기에 이상지역의 얼음변화도를 관측 및 분석한다면 화산활동의 정도를 유추할 수 있다.

백두산의 칼데라호인 천지는 북한과 중국의 경계에 위치하며, 해발 2,450 m의 고산지대에 존재하며, 천지의 넓은 지역의 정보를 얻기에 현장 조사를 통해서는 정치적 지리적 공간적 한계점이 존재한다(Yun, 2013). 따라서 현장관측의 한계점을 극복하고, 화산활동에 대한 정보를 신속하게 전달할 수 있는 위성영상자료를 활용한 원격탐사 기술을 적용하였다. 원격탐사 기술은 위성영상자료를 활용하기 때문에 현장관측의 한계뿐 아니라 주기적으로 생성이 되기 때문에 시계열 정보생성이 가능하다. 일반적으로 광학위성과 레이더 위성을 활용되는데, 광학영상의 경우 기상조건에 따라 많은 영향을 받는 한계점이 있지만, 레이더 영상에 비해 자료의 접근성이 높으며, 수십년 간 축적된 자료로부터 시계열 분석이 가능하며 매개체로부터 획득한 다양한 분광정보를 활용할 수 있다(Han and Lee, 2018).

광학위성 영상에서의 가시광선 대역 및 근적외선 대역의 수체, 짙은 그림자 및 얼음 분광반사도를 비교해 보면 수체, 짙은 그림자, 얼음 순으로 분광반사도가 낮으며 얼음의 두께에 따라서도 분광반사도가 다르다(Han et al., 2019). 특히, 수체는 주변보다 분광반사도가 낮기 때문에 얼음이 많이 융해되거나 얼음의 두께가 얇아 물에 가까워질수록 어둡게 나타나며, 완전히 결빙된 상태는 밝게 나타난다. 그러나 시계열영상의 경우 이상 지역이 존재하지 않은 경우도 있기 때문에, 수체와 분광반사도가 비슷한 짙은 그림자를 활용하였다. 따라서 가시광선 및 근적외선 대역에서 짙은 그림자와 얼음의 상대적인 분광반사도 차이를 활용하면 얼음의 융해정도를 정량적으로 나타내는 얼음변화도 값을 관측하는데 좋은 지표로 활용할 수 있다. 그러나 모든 시계열 영상 및 추후에 획득하는 광학위성영상 내의 가장 두꺼운 얼음값, 짙은 그림자값 그리고 수체값을 각각 측정하여 얼음변화도 값을 관측하는 데는 어려움이 존재한다.

최근 이미지 관련 분야에서 합성곱 신경망을 활용한 딥러닝 기반 기술은 다양한 활용분야에서 기존 탐지 기법의 성능을 상회하는 것으로 나타났다. 특히. 저차원 및 고차원 기능을 통합할 수 있는 U자형 대칭 네트워크(U-Net)를 제안하여 다양한 이미지 시멘틱 분야에서 높은 정확도를 나타냈다(Wu et al., 2018). 그러나 단일 인코더를 가지는 기본 U-net의 경우 다양한 분광정보를 가지고 있는 Landsat 광학위상영상의 입력데이터를 하나의 인코더에 입력하면서 합성곱 필터에 의해 각 채널이 가지고 있는 분광정보가 혼합된다. 이에 따라 각 입력 데이터가 가지고 있는 독립적인 분광정보는 소실된다(Hou et al., 2019). 이러한 독립적인 정보 소실 문제를 해결하기 위하여 Hou et al. (2019)는 두 개의 인코더를 가지는 딥러닝 모델을 제안하였다. 입력 노드를 분리하고 서로 다른 인코더에 배치함으로써, 각 입력데이터가 가지고 있는 독립적인 정보를 최대한 유지하여 분류 성능을 향상시켰다.

따라서 본 연구에서는 기존의 U-Net 모델의 구조적 한계를 극복한 Modified U-Net을 활용하여 Landsat 광학위성영상의 VNIR 대역으로부터 얼음변화도를 관측 및 모델 성능을 확인했다. 1988년부터 2020년까지 천지가 결빙된 총 83장의 영상을 수집하여 전처리 과정을 거쳐 Blue, Green, Red 그리고 Near Infrared (NIR) 대역 정보를 가진 학습데이터를 생성하였다. 이에 대한 정답 값으로는 짙은 그림자와 얼음의 상대적인 분광반사도 차이를 통해 천지 내 얼음변화도 값을 측정하여 활용했다. 설계한 모델에 대해서 천지 내 얼음변화도 값을 추정하는 회귀 모델을 학습하였으며, 최종 성능은 모델 학습에 활용되지 않은 시기의 영상에서의 실제 얼음 및 짙은 그림자 값을 관측해 획득한 얼음변화도 값과 회귀모델을 통해 예측한 얼음변화도 값의 차이를 RMSE 값과 상관계수로 나타냈다.

2. 연구지역 및 데이터

1) 연구지역

천지는 백두산의 과거 화산활동에 의해 백두산 화구에 생긴 칼데라 호수로 북한의 삼지연군과 중국 동북지방의 국경지대에 위치한다(Fig. 1). 백두산 천지는 수심이 373 m, 둘레는 13.11 km 면적은 9.82 km² 그리고 수면의 해발고도는 2189.1 m이며, 높은 위도와 고도로 인해 기온은 연 평균 -6~8°C, 최고 온도는 18~20°C이다(Yun and Oh, 2014). 천지지역 기후는 우리나라와 비슷하게 여름에는 덥고 겨울에는 춥기 때문에 천지는 계절에 따라 해빙과 결빙을 반복하는 시계열 패턴을 보여준다(KMA, 2017). 천지호 물의 약 61.5%는 천지아래의 마그마에 의해 증발된 지하수이며 나머지 대부분은 천지에 내리는 강수에 의해 결정된다(Suh et al., 2013).

Fig. 1. Location map of the Cheonji at Mt. Baekdu.

2) 데이터

본 연구에서는 과거 1985년부터 현재까지 백두산을 지속적으로 촬영했으며, 다양한 분광대역 센서를 가지고 있는 Landsat -5/-7/-8 광학위성영상을 활용하였다. 사용된 모든 Landsat -5/-7/-8 광학위성영상은 단순한 지형보정과 기하보정이 수행된 Level 1 산출물로써, 각각의 다중 분광 센서인 Thematic Mapper (TM)/Enhanced Thematic Mapper Plus (ETM+)/Operational Land Imager (OLI)의 다중 분광 센서를 통해 획득한 가시광선 및 근적외선 대역을 활용하였다. 이는 가시광선 및 근적외선 대역에서 물과 얼음의 분광반사도 차이가 커서 천지의 이상지역을 뚜렷하게 관측할 수 있기 때문이다(Han et al., 2019). TM/ETN+/OIL의 3가지 다중 분광 센서 모두 30 m 공간해상도를 가지며, 자세한 밴드 정보는 Table 1과 같다.

Table 1. Specifications of spectral bands (Blue, Green, Red and NIR) of Landsat -5/-7/-8

연구를 위해 1985년부터 2020년까지 운량 10% 미만의 영상 중 12~3월 영상만을 수집했다. 천지는 겨울철 낮은 온도로 의해 천지 표면이 10월부터 얼기 시작하여 4월까지 얼음이 유지되는데 그 중 12~3월 영상만이 천지표면이 완전히 결빙되어 온천수의 영향으로 이상지역이 뚜렷하게 관측되기 때문이다. 최종적으로 수집된 영상 중 육안분석을 통해 천지에 구름의 영향이 미치는 영상을 제외하여 총 83장을 선정하였다(Table 2). 선정된 모든 시계열 영상들은 천지를 모두 포함하면서 자세히 관측할 수 있도록 동일한 좌상단 좌표를 통해 200 × 200 크기로 잘라서 구축하였다(Fig. 2).

Table 2. List of image information used in this study

Fig. 2. Landsat RGB images cropped 200 × 200 size: (a) 20140223, (b) 20180218, (c) 20210226, (d) 20150314, (e) 20160128, and (f) 20180316.

3. 연구방법

백두산 천지의 이상지역 내의 얼음변화도를 Modified U-Net을 활용하여 관측하기 위해 사용된 데이터는 Landsat -5/-7/-8 Level 1 광학위성영상이다. 수집한 영상을 딥러닝 기법에 적용하여 얼음변화도 값을 관측하기 위해 입력데이터를 200 × 200 크기로 자른 뒤, 전처리 과정을 거쳐 Blue, Green, Red 그리고 NIR 대역 정보를 가진 입력데이터를 생성하였다. 먼저 다중 시기 광학위성영상자료에서 필수적인 전처리 과정인 대기보정을 수행하였다. 그 후, 관측하고자 하는 얼음변화도 값을 측정하기 위해서 대기보정이 완료된 영상의 수체, 짙은 그림자 그리고 얼음의 상대적인 분광반사도를 활용하였다. 측정된 얼음변화도 값은 모델의 ground truth로 활용된다. 제작된 입력데이터와 ground truth는 같은 시계열 순서로 통해 모델 훈련에 필요한 학습데이터와 모델을 검증하기 위한 테스트 데이터를 8:2 비율로 각각 66장과 17장으로 나누었다(Table 2). 딥러닝 기반 모델에 학습 시, 학습용 데이터셋이 부족한 경우 모델의 강건성이 떨어지며 오버피팅(overfitting)과 같은 다양한 문제가 발생하므로 데이터 증강을 통해 학습데이터셋을 19배 증강시켜 1,254장의 데이터 셋을 구축했다. 이 과정을 거친 학습데이터셋을 Modified U-Net에 적용하였고 테스트 데이터셋을 통해 모델을 정량적으로 검증하였다. 연구는 다음의 순서를 통하여 진행되었다(Fig. 3).

Fig. 3. Data flow used in this study.

1) 대기보정

시계열 데이터를 통해 분석을 수행하기 위해서는 대기보정은 필수적이다. Landsat 시리즈 광학위성영상 Level-1의 경우 센서에서 생성된 영상의 화소 값(Digital Number, DN)을 제공한다. 센서에서 감지된 영상의 DN값은 대기 입자에 의해 발생한 산란, 흡수 및 굴절에 의해서 많은 영향을 받아 데이터를 왜곡시킨다(Chavez, 1996). 특히, 천지와 같은 호수에서 반사되는 복사량은 매우 낮으므로, 대기보정은 필수적이다. 따라서 전체 시계열 영상에 동일한 방법으로 대기에 의한 영향을 제거해야 한다(Park et al., 2012).

대기의 영향을 보정하는 방법은 크게 대기관측 자료를 통해 보정하는 대기모델을 이용하는 방법과 영상을 이용하는 방법이 있다. 대기모델을 이용하는 방법은 Second Simulation of the Satellite Signal in the Solar (6S) 등과 같은 모델을 통해 보정하는 방법이지만 영상 촬영 지역의 대기관측자료가 필요하며, 많은 시간, 비용 그리고 인력 등이 소요된다(Park et al., 2012). 본 연구지역은 지형적으로 접근이 어려운 해발 2,450 m의 산악지형에 있으므로 영상을 이용한 대기에 의한 영향을 보정하였다. 이 경우는 직접적인 대기 관측자료가 필요하지 않은 장점이 있다. 본 연구에서 활용된 대기보정 방법은 US Geological Survey (USGS)에서 Landsat 위성의 종류에 따라 제공한 모델을 활용하였다. 위성의 센서에 기록된 DN을 대기 상층복사도로 변환해줌으로써 대기 보정이 가능하며 광학위성센서에 따라 다른 식을 적용해야 한다(Roy et al., 2016). Landsat -5/-7 위성의 경우, 식(1)을 통해 변환해주었다.

\(\begin{aligned}L_{\lambda}=\frac{L M A X_{\lambda}-L M I N_{\lambda}}{D N M A X-D N M I N}(D N-D N M I N)+L M I N_{\lambda}\\\end{aligned}\)     (1)

다음의 식에서 L_λ는 복사휘도(Radiance)이며 DN은 센서의 화소 값이다. DNMAX와 DNMIN는 DN의 최소값, 최대값이다. LMIN_λ과 LMAX_λ는 각각 DN 최소값, 최대값에 대하여 spectral radiance scales이다. 최종적으로 L_λ(Radiance)는 식(2)를 통해 대기 상층복사도로 변환하였다.

\(\begin{aligned}\rho_{\lambda}=\frac{\pi L_{\lambda} d^{2}}{\operatorname{ESUN}_{\lambda} \operatorname{Cos}\left(\theta_{s}\right)}\\\end{aligned}\)     (2)

다음의 식에서 ρ_λ는 대기 상층복사도이고 d는 지구와 태양까지의 거리 AU (천문 단위)이다. ESUN_λ는 대양으로부터 복사된 복사조도(Solar exoatmospheric irradiance)이다. θ_s는태양의천정각(Solar zenith angle)이다(D’Allestro and Parente, 2015).

Landsat -8 위성의 경우, 식(3)을 통해 태양각도 보정을 하지 않은 대기 상층복사도를 구한 뒤 최종적으로 식(4)를 통해 태양각도를 고려한 대기 상층복사도를 구한다.

ρ′_λ = M_pDN + A_ρ (3)

다음의 식에서 ρ′_λ는 태양각도 보정을 하지 않은 대기 상층복사도이며, M_p는 Band-specific multiplicative rescaling factor로 본 연구에서 활용된 모든 영상은 0.00002로 동일하였다. DN은 센서의 화소값이고 A_ρ는 Band-specific additive rescaling factor로, 모든 영상이 –0.1로 동일하였다. 이후에 식(4)를 통해서 태양의 각도를 고려하여 최종 대기 상층복사도를 구한다.

\(\begin{aligned}R_{\lambda}=\frac{\rho_{\lambda}^{\prime}}{\operatorname{Sin}\left(\theta_{S E}\right)}\\\end{aligned}\)     (4)

다음의 식에서 R_λ는 태양의 각도를 고려한 최종 대기 상층복사도이며, (Sin(θ_SE)은 태양의 천정각이다.

2) 상대적인 분광반사도를 활용한 얼음변화도 값 생성

대기보정이 완료된 영상의 분광반사도 활용하여 정량적으로 얼음변화도 값을 측정하였다. Red, Green, Blue 그리고 NIR대역에서는 수체에서 매우 낮은 분광 반사도를 나타내며, 얼음에서는 비교적 높은 분광반사도를 나타낸다. 또한 천지내에서 짙게 진 그림자 부분은 수체보다 평균적으로 0.01–0.02 정도 높고, 얼음보다는 0.4–0.5의 낮은 분광반사도를 나타내므로 시계열 영상 중 얼음이 융해되지 않은 시기의 영상의 경우 수체 대신 그림자가 짙게 진 부분의 값을 대체할 수 있다. 이에 따라, 물 또는 짙은 그림자와 얼음의 상대적인 분광 반사도 차이를 이용하면 결빙된 천지내에서 융해된 지역의 경계를 뚜렷하게 구분할 수 있으며, 이상지역에서의 얼음이 상대적으로 융해된 정도를 정량적으로 측정할 수 있다. 먼저 영상의 입력데이터 Blue, Green, Red 그리고 NIR 밴드의 반사도 값을 모두 더하고, 천지호 내에서 수체 부분의 반사도 값을 1로, 호수의 얼음 반사도 값을 0으로 가정한다. 이 때 천지의 시계열 영상 중에는 얼음이 전혀 녹지 않은 시기의 영상의 경우에는 수체 대신 그림자가 짙게 진 부분의 값을 활용하였다. 그 후 선형보정식을 사용하여 얼음과 수체 사이의 값을 추정하고, 0보다 작은 값은 0으로, 1보다 큰 값은 1로 설정하였다. 최종적으로 1–선형보정값을 통해 최종 얼음변화도 값을 측정하였다. 따라서 얼음변화도 값은 0–1 사이의 값을 가지며, 1에 가까운 값을 가질수록 얼음이 더 많이 융해 되었음을 의미하고 0에 가까운 값일수록 융해되지 않은 상태임을 의미한다. 그러나, 강설 등으로 인하여 피복의 변화가 큰 육지 영역이 패치에 포함되면 딥러닝 훈련 시 오차를 야기할 수 있다. 이를 저감하기 위하여, ground truth 및 입력데이터에서 육지 영역을 제외하였다(Fig. 4).

Fig. 4. Ice gradient using the relative spectral reflectance: (a) 19850122, (b) 20030124, (c) 20070212, (d) 20090217, (e) 20110130, and (f) 20170106.

Table 3는 전체 시계열 영상에 대해 각각 그림자가 가장 짙게 진 부분의 값과 호수 얼음의 값을 나타낸다.

Table 3. Times-series water and ice value

3) 데이터 분할 및 데이터 확장

딥러닝 모델을 훈련 및 검증하기 위해서 전처리된 입력자료와 제작한 얼음변화도 값(ground truth)을 같은 순서로 랜덤으로 섞은 뒤 훈련데이터 66장과 테스트데이터 17장으로 나누었다. 훈련데이터를 통해 딥러닝 모델을 훈련하고, 테스트 데이터를 통해 학습이 완료된 모델의 최종 성능을 평가하였다.

본 연구에서 학습 데이터 제작을 위해 사용된 광학위성영상은 구름 등 기상 상황에 민감하여 동일 지역에 대한 많은 데이터를 확보하기 어렵다. 이에 따라, 활용할 수 있는 훈련데이터수가 66장으로 학습된 모델이 훈련데이터에 대해 과적합되어 테스트셋의 데이터를 잘 설명하지 못하는 문제점이 발생하여 데이터 확장 기법이 필수적이다(Wang and Perez, 2017). 데이터 확장 기법은 원본 학습데이터에 임의의 기하학적 변환을 통하여 데이터셋의 양을 증가시킴으로써, 과적합을 방지하고 모델의 강건성을 향상시켜준다. 따라서 본 연구에서는 먼저 –45°~45°까지 5°단위로 무작위로 회전하여 19배 증강하여 1,254장의 훈련데이터 셋을 구축하였다. 시계열적으로 특별한 화산활동이 없는 경우 이상지역의 모양이 비슷한데 우선적으로 다양한 회전변환을 해줌으로써 이상지역의 모양을 변환해줌으로써 다양한 사례에 대해 정보를 추출할 수 있는 능력을 확보하였다. 그 후 증강된 훈련데이터셋에 추가적으로 기하 왜곡(affinetransformation), 임의 회전, 임의 영역 삭제(cut off)의 데이터 확장 기법을 적용하였다. 특히, 학습시에 증강하기 전 원본데이터는 사용하지 않고 증강된 데이터로만 학습함으로써 모델의 강건성을 더욱 향상시켰다.

4) 딥러닝 기법 적용

본 연구에서는 얼음변화도 값을 측정할 때, 모든 시계열 영상을 한 장씩 확인하여 수체지역의 분광반사도 픽셀 값. 가장 짙게 진 그림자 분상반사도 픽셀 그리고 얼음의 분광반사도 픽셀을 확인하였다. 위성영상 특성상 추후에도 지속적으로 데이터를 제공하며, 제공된 데이터에 대해 이러한 후처리 작업을 통해 얼음변화도 값을 측정하는 것은 많은 시간과 노력이 필요하다. 이러한 한계점을 극복하기 위해 Modified U-Net을 활용하여 추후의 시계열 광학위성영상에 대해 자동적으로 얼음 변화도 값을 측정하는 모델을 학습하고자 한다.

최근 딥러닝은 원격 감시를 비롯한 다양한 영상 처리 분야에서 널리 사용되어 왔으며 기존의 비딥러닝 알고리즘 성능을 능가한다(Yu et al., 2018; Baek, 2022). 또한, 원격 감시 이미지의 복잡한 기능을 더 잘 처리할 수 있으며 강력한 자율 학습 능력을 가지고 있다(Yao et al., 2018). 특히, 딥러닝 기법 중 영상분야에서 가장 활발하게 활용되는 기술은 합성곱 신경망이다. 이미지 인식 및 대상 감지 외에도 합성곱 신경망은 픽셀 단위의 문제에서도 성공적으로 적용되었다(Schmidt et al., 2019). 따라서 본 연구에서는 얼음변화도 값을 측정하는 모델을 구축하기 위해 합성곱 신경망 모델을 활용한 구조 중 하나인 U-Net의 구조를 사용하였다. U-Net은 원격 감시 이미지 분할을 포함하여 픽셀 단위 분류 문제를 해결하는데 매우 효과적이며 뿐만 아니라, 픽셀 단위 회귀 문제에서도 좋은 결과를 보여주었다(Yao et al., 2018). 그러나 기존의 U-Net 구조는 한 개의 인코더로 이루어져 있어 합성곱 필터과정에서 각 채널이 가지고 있는 정보가 혼합된다(Hou et al., 2019). 이에 따라 각 입력데이터가 가지고 있는 독립적인 정보가 소실된다. 이러한 독립적인 정보 소실 문제를 해결하기 위하여 본 연구에서는 두개의 인코더를 가지는 Modified U-Net 모델을 활용하였다(Hou et al., 2019).

본 연구에서는 천지 이상지역에서 각 분광대역에서 가지고 있는 정보를 최대한 유지하기 위하여 상관도가 높은 가시광선 대역들은 하나의 인코더로 입력하고 가시광선대역과 상관도가 낮은 NIR 대역은 다른 인코더에 배치하였다. 본 연구에서 활용된 모델의 구조는 Fig. 5와 같다.

Fig. 5. Modified U-Net architecture in this study.

본 연구에서 제안한 모델의 구조는 기존의 U-Net 구조와 비슷하게 U자형 구조를 가진다. 전체 모델의 구조는 인코더와 디코더로 구성되며 두 구조는 서로 대칭적이다. 각 인코더와 디코더는 합성곱 계층, 전치 합성곱계층, 풀링 계층으로 구성되어 있다. 본 연구에서는 비선형 문제, 과적합 우려 및 계산 효율성을 해결하고 복잡한 특징을 학습하기 위한 고차원 피처맵을 만들기 위해 3 × 3 커널을 사용하여 3개의 합성곱 계층을 연속으로 적용하였고 연산이 끝남과 동시에 1개의 풀링 계층을 지속적으로 계산하였다(Simonyan and Zisserman, 2014). 또한, 매 인코딩 레이어에서 연산을 통해 생성된 특성맵은 오른쪽의 디코딩 레이어에서 합성곱 필터 연산전에 결합함으로써 해상도 정보를 유지하면서 새로운 정보도 함께 활용할 수 있도록 하였다. 매번 맥스풀링(Max-pooling)을 적용할 때는 이전 계층보다 특성맵의 크기를 반으로 줄였으나 특성맵의 수는 두 배로 늘렸다. 가장 위의 계층에서는 특성맵의 크기는 200 × 200이며, 특성맵의 수는 64개 생성되도록 구성하였으며, 이를 반복적으로 네번 적용하여 가장 아래 계층에서 특성맵의 크기는 12 × 12이며, 특성맵의 수는 1,024이다. 디코더 과정에서는 이를 반대로 적용하였다. 최종 출력 값은 얼음의 변화값이 0–1 사이의 값을 가지므로, 동일하게 0–1 사이의 값으로 예측하도록 회귀모델 형태로 설계하였다. 최종 출력층의 활성화 함수는 ReLU, 손실함수는 평균절대값오차, 최적화함수는 Nadam, 학습률은 0.00001, 미니배치 크기는 5, 학습 반복수는 250으로 설정하였다(Table 4).

Table 4. Hyperparameter used for Modified U-Net

5) 모델 성능 검증

훈련데이터를 통해 만들어진 딥러닝 모델에 대한 신중하게 평가 및 분석이 필요하다. 이에 따라, 학습에 활용되지 않은 테스트 데이터셋을 활용하여 모델의 성능을 정량적으로 평가하였다. 본 연구에서 활용한 모델은 회귀 모델로써 평균 제곱근 오차(RMSE)와 피어슨 상관 계수(Pearson correlation coefficient analysis)를 통해 정량적인 평가를 실시했다. 평균 제곱근 오차는 실제 값과 예측값의 편차를 제곱하여 평균한인 평균 제곱 오차에 제곱근한 값으로, 이상치에 대해 민감한 평균 제곱 오차보다 조금 더 정확한 평가 지표이다. y를 실제값, \(\begin{aligned}\hat{y}\\\end{aligned}\)를 모델이 예측한 값, n이 데이터 개수라고 할 때, 평균 제곱근 오차는 식(5)를 통해 계산된다(Chai and Draxler, 2014).

\(\begin{aligned}\operatorname{RMSE}=\sqrt{\frac{(y-\hat{y})^{2}}{n}}\\\end{aligned}\)     (5)

또한, 실제값과 모델이 예측한 값의 상관성 분석을 위해 두 변수 사이의 상관계수 γ를 계산했다. 상관계수란 두 변수 X, Y 간의 선형 상관 관계를 계량화한 수치로써, 본 연구에서 활용한 피어슨 상관계수는 코시-슈바르츠 부등식에 의해 +1과 –1 사이의 값을 가지며 +1은 완벽한 양의 선형 상관 관계, 0은 선형 상관 관계없음, –1은 완벽한 음의 선형 상관관계를 의미한다. x, y는 각각 변수의 값을 나타내며, \(\bar{x}\), \(\begin{aligned}\bar{y}\\\end{aligned}\)는 각 변수집합의 표본 평균이라고 할 때. 상관 계수 γ는 식(6)을 통해 계산된다(Egghe and Leydesdorff, 2009). 본 연구에서는 평균 제곱근 오차 및 상관계수를 구할 때, 영상 한 장 내의 픽셀값을 모두 더해서 비교하였다.

\(\begin{aligned}\gamma=\frac{\Sigma(x-\bar{x})(y-\bar{y})}{\sqrt{\Sigma(x-\bar{x})^{2}(y-\bar{y})^{2}}}\\\end{aligned}\)     (6)

4. 결과

Fig. 6는 천지를 면밀하게 관찰하기 위해 공간해상도가 30 m 크기인 Landsat -5/-7/-8 영상을 천지가 모두 포함하도록 200×200으로 영상들을 자른 뒤 대기 보정 및 육지지형을 제외한 True color map (Red, Green, Blue composite band)이며, False color map (NIR, Blue, Green composite band) 그리고 상대적인 분광반사도를 활용하여 얼음변화도 값을 측정한 ground truth map이다. Ground truth 값은 모두 0–1값을 보여준다. Fig. 6의 흰색 box는 온천수 방출 예상지역으로 얼음이 녹아 있는 것을 뚜렷하게 확인할 수 있는 지역이다. 특히, Fig. 6(c1, c2, c4)에서 가장 많이 녹은 부분이 두 군데로 나뉘는데 이 지역은 온천수 방출되는 지역으로 예상된다. 이상지역의 모양을 확인해보면 일반적으로 온천수 방출 예상지역에서 뜨거운 온천수로 인해 열이 사방으로 천천히 전도되면서 얼음이 융해되는 모양이 Fig. 6(1–4)와 같다. 그러나 Fig. 6(1–4)와 다르게 Fig. 6(5–6)을 확인해보면 이상지역이 가늘고 길게 뻗어 나타나 것을 알 수 있는데 이는 온천수가 상당한 압력을 지니고 방출되고 있는 것으로 판단된다. 이는 온천수의 양도 늘고, 마그마의 영향으로 온도가 더 높아져서 나타나는 현상으로 해석이 가능하다. 또한, 시계열 영상 중인 얼음변화도 값이 가장 큰 값인 시기를 확인해보면 2002년 12월~2003년 1월 시기의 영상 4장이 유사한 패턴으로 이상지역이 넓고, 얼음변화도 값이 높게 나타난다. 이는 이 당시 온천수가 매우 뜨겁고 많은 양을 방출하고 있었다는 것을 의미한다(Fig. 7). 다른 화산 전조현상결과를 확인해본 결과 2002년에서 2006년 사이에 백두산화산 활동이 활발했음을 나타났는데, 이 시기와 정확히 일치한다(Xu et al., 2012). 따라서 천지호가 얼었을 때 이상지역의 모양, 크기, 얼음변화도 값의 변화는 백두산 화산활동을 분석하는데 매우 중요한 요소로 판단된다.

Fig. 6. Ice gradient input data for Modified U-Net regression model: (a) True color map, (b) False color map (NIR, Blue, Green), (c) Ice gradient (GT), (1) 19871230, (2) 19970123, (3) 20010323, (4) 20080114, (5) 19980110, and (6) 19980126.

Fig. 7. Ice gradient in 2002–2003: (a) 20021231, (b) 20030108, (c) 20030124, and (d) 20030201.

다음으로, 앞서 제작한 입력자료를 활용하여 본 연구에서 제안한 Modified U-Net에 학습하였고 모델의 결과를 RMSE와 상관 계수를 통해 평가하였다. RMSE는 예측 모델에 대한 테스트 데이터 세트를 기반으로 계산되었으며 140 값으로 한 영상 내의 픽셀 개수가 40,000개인 것을 감안하면 상당히 낮은 값이다(Fig. 8(a)). 상관계수를 확인해봐도 약 0.9968로 실제 값과 거의 일치하는 사실을 알 수 있다(Fig. 8(b)). 또한 추가적으로 정성적인 평가를 통해 확인해본 결과 이상지역 내의 값의 분포는 거의 일치하는 것을 확인할 수 있었다(Fig. 9). 따라서 본 연구에서 활용된 모델이 연구 영역의 얼음변화도 값을 관측하는데 상당히 정확하다는 것을 알 수 있다. 이는 Modified U-Net의 특징인 두 개의 인코더에서 특성맵을 추출함으로써, 분광대역에서의 특성을 최대한 활용하여 이상지역 내에서 융해된 부분을 명확하게 측정할 수 있었다. 이를 통해 추후에 획득하는 Landsat 광학위성영상을 학습되어진 Modified U-Net을 활용하여 이상지역의 얼음변화도 값을 정량적으로 관측할 수 있으며, 백두산 화산활동 모니터링에 좋은 지표로 활용될 수 있다.

Fig. 8. Model performance evaluation results: (a) Root Mean Square Error (RMSE) and (b) Pearson’s correlation coefficient.

Fig. 9. Qualitative analysis for test data: (a) true color map, (b) prediction map using the modified U-Net, (c) ground truth, (1) 20160112, (2) 20170106, (3) 20180117, and (4) 20200131.

4. 결론 및 토의

백두산의 화산활동으로 인해 천지내 방출되는 온천수의 온도와 수압 등이 변하게 된다. 이에 따라 천지의 백운봉 아래 천지지역에는 다른 천지 부분보다 얼음이 빨리 융해되고 늦게 얼며 지표온도가 높은 이상지역이 존재한다. 따라서 본 연구에서는 화산활동으로 의해 이상지역의 얼음변화도를 모니터링하기 위해 Landsat 광학위성영상을 활용하여 얼음변화도 값을 관측했으며, 이를 딥러닝 모델인 Modified U-Net에 적용하여 추후의 영상에서도 효율적이고 정확한 얼음변화도 값을 관측하고자 한다. 이를 위하여 천지가 얼어 있는 1985년부터 2020년까지 12~3월 영상을 활용하였다. 먼저 광학위성영상의 수체, 그림자 그리고 얼음의 분광반사도 차이를 통해 얼음변화도 값을 정량적으로 측정한 후 분석하였다. 또한, 전처리한 데이터와 제작한 얼음변화도 값을 Modified U-Net 회귀 모델에 적용하여 회귀결과를 정량적·정성적으로 검증하였으며 다음과 같은 결론을 얻을 수 있었다.

첫째, 상대적인 분광반사도를 통해 얼음변화도 값을 정량적으로 측정하고 시계열적으로 분석하였다. 값이 1에 가까울수록 얼음이 많이 융해되었으며, 0에 가까울수록 얼음에 가까운 것을 확인한 결과, 시계열적으로 이상지역 내에서 얼음이 많이 녹는 지역 두 군데를 확인할 수 있었으며 이는 온천수가 방출되는 지역으로 예상된다. 또한, 시계열적으로 이상지역의 모양이 일반적이지 않은 일자형의 영상을 확인할 수 있는데 이 시기에 온천수의 수압이 평소보다 높았던 것으로 해석할 수 있다. 특히, 2002~2003년 시기에는 얼음변화도 값이 높았는데 이 시기는 백두산의 화산활동이 활발했던 시기이다. 이에 따라, 얼음변화도 값이 화산활동을 모니터링할수 있는 좋은 지표임을 알 수 있다.

다음으로 전처리한 데이터와 제작한 얼음변화도 값을 두 개의 인코더를 가진 Modified U-Net에 적용하였으며 훈련에 사용되지 않은 테스트 데이터만으로 딥러닝 모델의 성능을 정성적·정량적으로 평가했다. 모델 성능의 정량적 평가 지표는 RMSE와 상관계수를 활용했으며 각각 140, 0.9968로 높은 회귀 정확도를 보여줬다. 또한, 정성적으로 평가한 경우에도 실제 값과 예측된 모델값이 거의 동일한 분포를 나타냈다. 따라서 Modified U-Net을 활용하면 높은 정확도로 얼음변화도를 관측할 수 있으며 이는 백두산의 화산활동을 모니터링하는 방법 중 하나로 사용될 수 있다. 그러나 본 연구에는 몇 가지 한계를 지니고 있다. 첫째, 광학위성영상 특성상 기상의 영향을 많이 받아 연구지역의 정기적인 관측에 한계가 있으며 적시성이 저하된다. 둘째, 천지호 얼음의 경우 단시간에 어는 것이 아니라 장시간동안 얼음으로 변하므로 장시간의 온도, 습도, 강우량 및 적설량 등에 따른 해석의 오류를 최소화할 필요가 있다. 따라서 향후 연구에서는 복사해상도, 시간해상도 및 공간해상도가 높은 Sentinel-2, KOMPSAT-2/-3/-3A, 및 차세대 중형위성 등의 광학위성영상을 함께 활용할 필요가 있다. 또한 주변 지역별상세관측자료(Automatic Weather Station)로부터 얻은 장기간의 데이터를 활용하여 데이터에 대한 보정을 할 필요가 있다.

사사

이 연구는 기상청 「화산활동 감시·예측 및 한반도 지각활동 진단 기술개발」(과제번호: 1365003682)의 지원을 받아 수행되었습니다.

이 연구는 국토교통부의 『스마트시티 혁신인재육성 사업(’19–’23)』으로 지원되었습니다.

이 연구는 2020년 대한민국 정부의 재원으로 한국연구재단(NRF-2020R1A2C1006593)의 지원을 받아 수행되었습니다.