Terrain Shadow Detection in Satellite Images of the Korean Peninsula Using a Hill-Shade Algorithm

음영기복 알고리즘을 활용한 한반도 촬영 위성영상에서의 지형그림자 탐지

Abstract

In recent years, the number of users has been increasing with the rapid development of earth observation satellites. In response, the Committee on Earth Observation Satellites (CEOS) has been striving to provide user-friendly satellite images by introducing the concept of Analysis Ready Data (ARD) and defining its requirements as CEOS ARD for Land (CARD4L). In ARD, a mask called an Unusable Data Mask (UDM), identifying unnecessary pixels for land analysis, should be provided with a satellite image. UDMs include clouds, cloud shadows, terrain shadows, etc. Terrain shadows are generated in mountainous terrain with large terrain relief, and these areas cause errors in analysis due to their low radiation intensity. previous research on terrain shadow detection focused on detecting terrain shadow pixels to correct terrain shadows. However, this should be replaced by the terrain correction method. Therefore, there is a need to expand the purpose of terrain shadow detection. In this study, to utilize CAS500-4 for forest and agriculture analysis, we extended the scope of the terrain shadow detection to shaded areas. This paper aims to analyze the potential for terrain shadow detection to make a terrain shadow mask for South and North Korea. To detect terrain shadows, we used a Hill-shade algorithm that utilizes the position of the sun and a surface's derivatives, such as slope and aspect. Using RapidEye images with a spatial resolution of 5 meters and Sentinel-2 images with a spatial resolution of 10 meters over the Korean Peninsula, the optimal threshold for shadow determination was confirmed by comparing them with the ground truth. The optimal threshold was used to perform terrain shadow detection, and the results were analyzed. As a qualitative result, it was confirmed that the shape was similar to the ground truth as a whole. In addition, it was confirmed that most of the F1 scores were between 0.8 and 0.94 for all images tested. Based on the results of this study, it was confirmed that automatic terrain shadow detection was well performed throughout the Korean Peninsula.

최근 지구관측 위성이 급격히 발전함에 따라 사용자의 수가 증가하고 있다. 이에 따라 지구관측위성위원회(Committee on Earth Observation Satellites, CEOS)에서는 분석준비자료(Analysis Ready Data, ARD)라는 개념을 제안하고 분석준비자료의 요구 조건을 CEOS ARD for Land (CARD4L)로 정의하여 사용자 친화적인 위성영상을 제공하기 위해 노력하고 있다. 분석준비자료에는 육상분석에 불필요한 픽셀이 식별된 마스크(Unusable Data Mask, UDM)가 영상과 함께 제공되어야 한다. UDM의 종류는 구름, 구름 그림자, 지형그림자 등이 있다. 지형그림자는 지형기복이 큰 산악지형에서 발생되며 지형그림자가 생긴 지역은 복사조도가 낮기 때문에 분석 결과에 오류를 야기시킨다. 기존 지형그림자 탐지연구는 지형그림자 보정을 위해 지형그림자 픽셀을 탐지하는데 목적을 두었지만, 이것은 지형보정 기법으로 대체 가능하다. 따라서 지형그림자 탐지 목적을 확장할 필요가 있다. 산림과 농업분석을 목적으로 한 차세대중형위성 4호(CAS500-4)의 활용을 위해 본 연구에서는 지형그림자 탐지 범위를 태양의 영향을 적게 받는 지역까지 확장하였다. 본 논문은 남북한을 대상으로 지형그림자 마스크 생성을 위해 지형그림자 탐지 가능성을 분석하는데 목적이 있다. 지형그림자 탐지를 위해서 태양의 위치, 지표면의 경사와 경사방향을 이용한 음영기복 알고리즘을 사용하였다. 한반도를 촬영한 5 m급 공간해상도의 RapidEye 영상과 10 m급 공간해상도의 Sentinel-2 영상들을 대상으로 참값과 비교하며 최적의 음영기복 임계값을 결정하였다. 결정된 임계값을 사용하여 지형 그림자 탐지를 수행하고 결과를 분석하였다. 정성적 결과로는 전체적으로 참값과의 형상이 유사함을 확인하였다. 정량적 실험결과는 F1 score가 대부분 0.8에서 0.94 사이인 것을 확인하였다. 본 연구 결과를 바탕으로 남북한을 대상으로 자동적인 지형그림자 탐지가 잘 수행됨을 확인하였다.

1. 서론

위성영상은 데이터 전처리에 많은 시간이 소비되기 때문에 사용자가 위성영상을 쉽게 접근하는데 한계가 있다. 이를 위해, 지구관측위성위원회(Committee on Earth Observation Satellites, CEOS)는 사용자 친화적인 위성영상을 제공하기 위해 분석준비자료(Analysis Ready Data, ARD)라는 개념과 분석준비자료의 요구사양을 CEOS ARD for Land (CARD4L) 로 명명하여 제시하였다(Committee on Earth Observation Satellites, 2021).

ARD 제품은 위성영상과 분석에 필요 없는 데이터 마스크(Unusable Data Mask, UDM)를 함께 제공해야 한다. UDM은 구름, 지형그림자 등 지상관측에 필요 없는 픽셀을 식별한 마스크 데이터를 의미한다. 그 중 지형그림자란 산림지에서 지형의 기복에 의해 발생하는 음영지역이다. 지형그림자는 지표면의 복사조도 손실을 야기시키기 때문에 지상분석에 부정확한 영향을 미친다(Li et al., 2016). 따라서 이를 식별한 마스크 데이터가 필요하다.

기존 지형그림자 탐지 연구 Giles(2001)는 복사조도가 손실된 픽셀을 보정하기 위해서 직달일사가 지형에 의해 차단된 그림자 픽셀을 탐지하는 연구를 수행하였다. 하지만 Tan et al. (2013)은 모든 픽셀을 지형의 기하 상태를 고려하여 모든 픽셀에 대해 지형보정을 수행하면 변화탐지 정확도를 향상시킨다고 하였다. 그러므로 지형그림자 보정을 위한 지형그림자 탐지는 지형보정 기법으로 대체 가능하다. 따라서 지형그림자 탐지 목적을 확장할 필요가 있다.

Sternberg and Shoshany (2001)는 태양의 영향을 많이 받는 남 사면(South-facing slope)과 태양의 영향을 적게 받는 북 사면(North-facing slope)을 기준으로 산림의 바이오매스(Biomass) 차이를 분석하였다. Yang et al. (2020)도 남 사면과 북 사면을 기준으로 식생 종류와 생태다양성이 차이가 나는 것을 분석하였다. 산림과 농업분석을 목적으로 한 차세대중형위성 4호(이하 농림위성)는 산림 생태와 산림현황 모니터링 등의 활용계획이 있다(Kim et al., 2019). 따라서 농림위성 활용을 위한 지형그림자 탐지 연구는 산림의 특성 등을 파악하기 위해서 태양의 영향을 적게 받는 경사진 지역까지 포함하여 탐지할 필요가 있다.

본 연구에서는 지표면의 경사로 인해 반사도가 감소된 지역을 반그림자라고 정의하였다. 연구의 탐지 범위는 지형그림자를 빛이 지표면에 도달하지 않아서 발생한 그림자 지역과 반그림자 영역까지 확장하였다. 연구의 목적은 남북한을 대상으로 2025년에 발사될 5 m 공간해상도를 가지는 농림위성의 지형그림자 마스크 생성을 위해서 지형그림자 자동 탐지 가능성을 분석하는 것이다. 지형그림자 탐지 기법은 모델기반 기법으로 불리는 광선 추적(ray tracing) 기법, 음영기복(hill-shade) 기법 등으로 구분된다(Mostafa, 2017). 본 논문은 광선 추적 기법 대비 처리시간이 짧은 음영기복 기법을 적용하였다. 실험 데이터는 국토지리정보원(National Geospatial Information Institute)에서 제작한 10 m 공간해상도의 digital elevation model (DEM) 데이터, 5 m급 공간해상도의 RapidEye 위성영상 및 10 m급 공간해상도의 Sentinel2 위성영상을 활용하였다. 영상의 전처리를 위해서 기하 및 정사보정이 수행되었으며 DEM 데이터는 리샘플링(resampling)하여 위성영상과 공간해상도를 일치시켰다. 한반도를 촬영한 위성영상을 대상으로 태양의 고도각을 이용하여 반그림자 임계값을 분석하였다. 탐지정확도 분석을 위해 본 연구에서는 지형그림자 탐지결과를 참값 데이터와 비교하였다. 결과적으로, 탐지 결과의 F1 score가 2개의 실험지역을 제외하고 0.8–0.94 사이임을 확인하였다. 본 논문의 결과가 향후 농림위성 영상의 지형그림자 마스크 생성에 활용될 수 있을 것이라고 기대한다.

2. 연구자료 및 방법

2.1. 연구 방법

본 연구에서 적용한 연구방법은 Fig. 1과 같다. 첫 번째로 위성영상과 DEM 데이터에 전처리 과정이 수행된다. 두 번째로 DEM을 이용하여 픽셀마다 X축과 Y축방향의 고도 변화량을 산출하고 이를 이용하여 지형의 경사각(slope)과 경사방향(aspect)을 계산한다. 다음으로 지형경사각과 경사방향 그리고 태양 고도각과 방위각을 이용하여 음영기복도를 생성한다. 최종적으로 음영기복 임계값을 적용하여 지형그림자를 탐지한다.

Fig. 1. A workflow of DEM-based terrain shadow detection algorithm.

본 논문에서 탐지하고자 하는 지형그림자는 태양광이 지표면에 도달되지 않아서 발생하는 순그림자 영역과 태양광이 지표면에 도달하지만 지표면의 경사로 인하여 반사도가 감소되는 반그림자 영역을 모두 포함한다. 순그림자는 태양과 지표면의 기하학적 관계를 통해 음영기복도의 밝기값이 0 이하인 픽셀이다. 그러나 순그림자는 그 영역이 매우 적어 지형기복에 따른 음영지역을 모두 표현하기에 한계가 있다. 따라서 본 연구에서는 순그림자 영역과 반그림자 영역을 모두 포함할 수 있는 최적의 음영기복 임계값을 선정하고 이 값을 이용하여 지형그림자 마스크를 생성하였다. 정확도 확인은 탐지결과를 참값 데이터와 비교하여 정성 및 정량적 분석을 하였다. 정성적 분석은 참값 데이터와 육안 판독을 통해 판단하였으며 정량적 분석은 F1 score, 재현율(Recall), 정밀도(Precision) 성능지표를 이용하여 분석하였다.

2.1.1. 데이터 전처리

위성영상의 기하학적인 왜곡을 보정하기 위해서 위성영상은 전처리과정이필요하다.왜곡을 보정하기 위해서 위성영상은 지상기준점 칩(Ground Control Points Chip)과 DEM을 이용하여 기하 및 정사보정이 수행된다. 본 연구에서는 기하 및 정사보정이 수행되지 않은 RapidEye 영상을대상으로선행된연구(Yoon, 2019; Park et al., 2020; Son et al., 2021; Choi and Kim, 2022)에 따라 기 구축된 지상기준점 칩을 활용하였다. 남한지역의 지상기준점 칩은 25 cm급 항공정사영상과 국토지리정보원에서 제작하는 통합기준점, 삼각점 및 항공정사영상 제작에서 사용된 사진기준점의 지상좌표를 활용하였으며, 북한지역의 지상기준점 칩은 50 cm급의 위성 정사영상과 정사영상의 평면 좌표 및 DEM의 타원체고 높이값을 활용하였다. 사용한 지상기준점 칩의 개수는 영상마다 다르며 최소 765점, 최대 1,594점이 사용되었다. Fig. 2는 사용한 지상기준점 칩의 예시이다. DEM 데이터는 클리핑(clipping)과정을 통해 DEM 데이터를 위성영상의 영역에 맞게 일치시키고 공간해상도 일치를 위해 공일차 보간법(Bilinear Interpolation)이 수행된다.

Fig. 2. Example of ground control points (GCP) chips.

2.1.2. 고도 변화량(Elevation Gradient) 추정

고도 변화량 추정 기법에는 여러 기법들(Horn, 1981; Ritter, 1987; Zevenbergen and Thorne, 1987)이 존재한다. 그 중에서 Horn (1981)의 기법은 DEM의 높이 값에 오차가 있더라도 결과에 민감하게 영향을 주지 않아 변화량을 가장 잘 표현할 수 있는 장점이 있다(Horn, 1981). Jones (1998)는 변화량 계산 기법들의 평균제곱근(Root Mean Square Error, RMSE)를 비교하는 연구를 수행하였으며, Horn (1981)이 제안한 기법이 가장 좋은 결과를 보여주었다. 따라서 본 논문에서는 Horn (1981)의 기법을 사용하였다. Horn (1981)이 제안한 고도 변화량 수식은 식(1), (2)와 같다. 식(1)은 X축 방향의 고도 변화량을 의미하며, 식(2)는 Y축 방향의 고도 변화량을 의미한다. Δx, Δy는 픽셀 사이즈를 의미한다. A부터 I는 Fig. 3과 같이 인접 픽셀의 고도값을 의미한다.

\(\begin{aligned}\frac{d z}{d x}=\frac{(A+2 D+G)-(C+2 F+I)}{8 \times \Delta x}\end{aligned}\)       (1)

\(\begin{aligned}\frac{d z}{d y}=\frac{(A+2 B+C)-(G+2 H+I)}{8 \times \Delta y}\end{aligned}\)       (2)

Fig. 3. Horn’sgradient estimation methodof x-axisdirection.

Fig. 3은 Horn (1981)의 기법에서 X축의 고도 변화량 과정을 나타낸다. 3×3 마스크를 이용하여 중앙 픽셀(E)을 기준으로 인접 픽셀들의 고도값을 차분하여 중앙 픽셀의 높이 변화량을 계산한다.

2.1.3. 경사각 및 경사방향 계산

본 논문에서 경사각은 지표면의 가파른 정도를 측정한 수치이다. 경사방향은 지표면이 북쪽을 기준으로 시계방향으로 0°에서 360°까지 표현되며 지표면이 향하는 방향을 나타낸다(Doumit, 2017). 식(3)과 (4)는 각각 경사각(γ)과 경사방향(δ) 계산식이다. \(\begin{aligned}\frac{d z}{d x}\end{aligned}\)와 \(\begin{aligned}\frac{d z}{d y}\end{aligned}\)는 각각 X축, Y축의 고도 변화량이다.

\(\begin{aligned}\gamma=\operatorname{Tan}^{-1}\left(\sqrt{\left(\frac{d z}{d x}\right)^{2}+\left(\frac{d z}{d y}\right)^{2}}\right)\end{aligned}\)       (3)

\(\begin{aligned}\delta=\operatorname{Tan}^{-1}\left(-\left(\frac{d z}{d y}\right) /\left(\frac{d z}{d x}\right)\right)\end{aligned}\)       (4)

2.1.4. 음영기복도 생성

음영기복도를 생성하기 위해서 태양의 방향벡터와 지표면의 법선벡터가 필요하다. 태양의 방향벡터(\(\begin{aligned}\vec{S}\end{aligned}\))와 지표면의 법선벡터(\(\begin{aligned}\vec{N}\end{aligned}\))는 각각 식(5), (6)과 같이 계산된다. \(\begin{aligned}\vec{S}\end{aligned}\)는 Fig. 4(a)와 같이 태양의 고도각(α)과 태양의 방위각(β)을 이용하여 표현할 수 있다. \(\begin{aligned}\vec{N}\end{aligned}\)는 Fig. 4(b)와 같이 지표면의 경사각(γ)과 지표면의 경사방향(δ)을 이용하여 표현 가능하다.

\(\begin{aligned}\vec{S}=\left[\begin{array}{c}\cos (\alpha) \times \cos (\beta) \\ \cos (\alpha) \times \sin (\beta) \\ \sin (\alpha)\end{array}\right]\end{aligned}\)       (5)

\(\begin{aligned}\vec{N}=\left[\begin{array}{c}\cos \left(90^{\circ}-\gamma\right) \times \cos (\delta) \\ \cos \left(90^{\circ}-\gamma\right) \times \sin (\delta) \\ \sin \left(90^{\circ}-\gamma\right)\end{array}\right]\end{aligned}\)       (6)

Fig. 4. Description of vectors required for hill-shade method: (a) the sun’s directional vector (\(\begin{aligned}\vec{S}\end{aligned}\)) and (b) a surface’s normal vector (\(\begin{aligned}\vec{N}\end{aligned}\)) on 3-dimensional space.

음영기복 알고리즘은 램버시안 반사(Lambertian reflectance)와 램버시안 코사인 법칙(Lambertian Cosine’s Law)을 적용하여 태양의 방향벡터(\(\begin{aligned}\vec{S}\end{aligned}\))와 지표면의 법선 벡터(\(\begin{aligned}\vec{N}\end{aligned}\))의 관계를 설명한 알고리즘이다(Yoeli, 1967; Farmakis-Serebryakova and Hurni, 2020). 램버시안 반사는 지표면에 의해 반사된 태양의 빛을 모든 방향에서 동일하게 반사한다고 가정한다. 램버시안 코사인 법칙에 따르면 반사된 빛의 양은 태양의 방향벡터(\(\begin{aligned}\vec{S}\end{aligned}\))와 지표면의 법선벡터(\(\begin{aligned}\vec{N}\end{aligned}\))의 사이각(입사각)이 코사인 값과 비례하며 두 벡터를 내적(dot product)하여 구할 수 있다. 따라서 지표면에 의해 반사된 빛의 양은 식(7)과 같이 구할 수 있다. 식(7)에서 Cos(i)의 i는 입사각(incidence angle)을 의미한다.

\(\begin{aligned} \vec{S} \cdot \vec{N}= & \operatorname{Cos}(i)=\left(\sin (\alpha) \times \sin \left(90^{\circ}-\gamma\right)+\right. \\ & \left.\cos (\alpha) \times \cos \left(90^{\circ}-\gamma\right) \times \cos (\beta-\delta)\right)\end{aligned}\)       (7)

2.1.5. 지형그림자 탐지

식(7)에 따르면 음영기복도의 밝기값은 두 벡터의 입사각(i)이 0°이면 코사인 값이 1이고 180°이면 –1로 구분된다. 이러한 이론을 바탕으로 지표면은 두 벡터의 입사각이 작아질 수록 빛을 많이 받는 것을 의미하며 입사각이 커질수록 빛을 적게 받는다(Konnelly, 2002). Konnelly(2002)는 이 관계를 통해 두 벡터의 입사각이 90° 이상일때 그림자로 정의하였다. 본 연구에서는 선행연구 Konnelly (2002)에 따라 식(8)과 같이 입사각이 90° 이상일때를 태양광이 지표면에 도달하지 않아서 발생하는 순그림자라고 판단하였다.

\(\begin{aligned} \text {Terrain Shadow Detection} =\\\left\{\begin{array}{lrlrl}\text { True Shadow, } & -1 & \leq \operatorname{Cos}(i) & \leq 0 \\ \text { others, } & 0 & <\operatorname{Cos}(i) & \leq 1\end{array}\right.\end{aligned}\)       (8)

그러나 본 논문의 2.1.절에서 제시된 것처럼 순그림자 만으로는 영상에서 관측되는 기복변이에 따른 음영지역을 모두 표현하지 못한다. 따라서 본 연구에서는 지형그림자를 Fig. 5와 같이 순그림자(True shadow)와 반그림자(Half shadow)로 구분하였다. 그림에서 Sunlit은 태양복사에너지가 온전한 비그림자영역을 의미한다. 순그림자는 태양광이 지표면에 도달되지 않는 조건인 입사각이 90°이상일 때 해당한다(Li et al., 2012). 그러나 반그림자는 태양광이 지표면에 도달하지만 지표면의 경사도로 인해 반사도가 낮아지기 경향이 나타나며 그 경계를 확인하기 어렵다.

Fig. 5. Types of terrain shadow.

본 연구에서는 반그림자의 경계까지 모두 포함하여 탐지하고자 하였다. 비그림자 영역(Sunlit area)과 반그림자의 경계를 찾기 위하여 Fig. 6과 같이 가정을 세웠다. Fig. 6(a)는 지형 경사각이 거의 없는 평지를 나타내며 비그림자 영역에 해당한다. 평지에서는 태양의 천정각이 입사각이 된다. Fig. 6(b)는 태양의 반대편으로 경사가 있는 지표면의 입사각을 나타낸다. 이 때의 입사각은 태양의 천정각에서 지표면의 경사각(γ)을 더한 값이 된다. 본 논문에서 반그림자 영역은 태양 반대 사면의 경사각이 특정값보다 커질 경우로 가정하여 식(9)와 같이 정의하였다. 즉, 입사각이 태양의 천정각 대비 특정 각도 이상인 경우, 지형기복에 따른 음영지역으로 가정하였다.

Halfshadow = incidence angle ≥ Sun’szenith + γ       (9)

Fig. 6. A supposition of setting a threshold of half terrain shadow: (a) a case of flat terrain and (b) a case of the inclined terrain facing away from the sun in mountainous terrain.

Richter et al. (2009)은 본 논문의 가정과 유사하게 태양의 천정각이 55° 이상일 때 Cos(태양의 천정각+10°), 45°에서 55° 사이일때 Cos (태양의 천정각+15°), 45°보다 작을 때 Cos (태양의 천정각+20°)의 값을 경우에 따라 구분하여 지형보정을 수행하였다. 본 논문에서는 영상에서 추출한 참값 지형그림자를 이용하여 반그림자 영역추출을 위한 최적의 경사각(γ)을 추출하고 이 결과를 Richter의 연구결과와 유사한지 비교하였다.

2.1.6. 정확도검증

정확도검증은 참값 데이터와 육안 판독을 통해 정성적 평가를 하였으며 오차행렬(Confusion Matrix)을 이용하여 정량적 평가를 하였다. 오차행렬은 정탐지(True Positive, TP), 정무탐지(True Negative, TN), 오탐지(False Positive, FP), 미탐지(False Negative, FN)로 표현되며 4가지 지표를 이용하여 성능측정지표를 계산할 수 있다. 성능측정지표는 재현율(Recall), 정밀도(Precision), F1 score를 사용하였다. 정밀도는 예측결과가 True라고 탐지한 것 중에서 실제 True인 것의 비율로 식(10)과 같이 정의되며 오탐지에 대한 정보를 제공한다. 재현율은 참값이 True인 것 중에서 예측 결과가 True라고 예측한 것의 비율로 식(11)과 같이 정의되며 미탐지에 대한 정보를 제공한다. F1 score는 정밀도와 재현율의 조화평균을 의미하며 식(12)와 같다.

\(\begin{aligned}\text {Precision}=\frac{T P}{T P+F P}\end{aligned}\)       (10)

\(\begin{aligned}\text {Recall}=\frac{T P}{T P+F N}\end{aligned}\)       (11)

\(\begin{aligned}\text {F1 score}=2{\times} \frac{Precision {\times} Recall}{Precision + Recall}\end{aligned}\)       (12)

2.2. 실험 데이터

본 연구에서는 Table 1과 같이 공간해상도가 농림위성과 유사한 5 m 공간해상도의 RapidEye 위성영상과 10 m 공간해상도의 Sentinel-2 위성영상을 사용하였다. 한반도의 산악지형이 복잡한 지형구조를 가지고 있기 때문에 안정적인 지형그림자 산출물을 얻기 위하여 한반도에 고루 분포된 영상으로 실험을 할 필요가 있다. 본 연구에서는 Fig. 7과 같이 남북한을 대상으로 널리 분포된 4장의 RapidEye 영상과 16장의 Sentinel-2 영상을 사용하였다. Table 2에 실험에 사용한 영상의 종류, 촬영일, 태양 방위각 및 천정각을 나열하였다. 태양의 위치가 계절별로 달라지기 때문에 실험에 사용된 영상의 시기는 4월부터 11월까지 정하였다. 한국천문연구원의 자료에 따르면, 오전 11시의 서울을 기준으로 연간 태양의 천정각은 최소 24.4°, 최대 65.6°이며 방위각은 최소 118.1°, 최대 158.8도이다. 실험의 사용된 태양의 천정각이 최소 19.7°부터 최대 54.6°, 방위각이 최소 132.9°부터 167.1°까지 사용되기 때문에 실험 데이터가 계절별 지형그림자 탐지 분석에 적합하다고 판단하였다. Fig. 8은 왼쪽부터 실험에 사용된 위성영상과 해당 위치의 DEM, 참값 데이터를 의미한다. DEM 자료의 검은색 픽셀은 바다나 평지대와 같은 저고도 지역을 의미하고 흰색 픽셀은 산악지형과 같은 고고도 지역을 의미한다. 참값데이터는 Quantum Geographic Information System (QGIS)에서 제공하는 음영기복도 기법을 이용하였으며, 전처리된 DEM과 각 영상마다 Table 2와 같이 동일한 태양의 위치를 사용하여 제작하였다. QGIS에서 제공하는 음영기복도의 값은 0–255 사이의 값으로 표현되어 있으며 음영이 많이 진 곳은 0, 음영이 없는 곳은 255를 의미한다. 참값데이터는 QGIS에서 제공하는 음영기복도의 값을 수동으로 조절하며 육안으로 식별되는 지형그림자 영역과 가장 잘 맞는 임계치를 영상별로 결정한 뒤, 이후 원본 위성 영상에서 보이는 구름 및 구름 그림자를 제외하고 육안으로 미세 조정하였다. QGIS로 제작한 참값의 임계값은 Table 3에 정리하였다. Scene-1~5의 참값데이터에서 보이는 빨간색 영역은 DEM의 고도값이 불확실성으로 해당 위치를 참값에서 수동으로 제거하였다. 본 논문에서의 고도값 불확실성의 기준은 비접근지역인 비무장 지대(Demilitarized Zone) 인근 지역으로 정하였다.

Fig. 7. Study area used in the experiment.

Table 1. Specification of satellite image used in the experiment and (proposal) CAS500-4 (Kwon et al., 2021)

Table 2. Metadata of satellite image used in the experiment

Fig. 8. Satellite image used in the experiment. From the left: satellite image, DEM, ground truth.

Table 3. Hillshade threshold of ground true data used in the experiment

3. 연구결과 및 토의

3.1. 음영기복도 생성 및 순그림자 탐지 결과

Fig. 9는 음영기복도 및 중간산출물 생성결과를 생성과정을 나타낸다. DEM을 이용하여 X와 Y축의 고도 변화량 맵(a), (b)를 생성하였다. 이후 (a)와 (b)를 이용하여 경사도 맵(c)와 경사방향 맵(d)를 생성하였다. (c)와 (d) 그리고 태양의 위치를 이용하여 최종적으로 음영기복도 맵(e)가 생성된다. 이러한 절차대로 위성영상 20장의 음영기복도를 제작하였다.

Fig. 9. The process of generating a hillshade map: (a) an elevation gradient map of X-axis, (b) an elevation gradient map of Y-axis, (c) a slope map, (d) an aspect map, and (e) a hillshade map.​​​​​​​

Fig. 10은 3개의 영상에 대한 음영기복도 생성 및 순그림자 탐지 결과이다. (a)는 Scene-1, (b)는 Scene-2, (c)는 Scene-3 영상을 의미한다. 왼쪽부터 순서대로 위성영상, 위성영상의 확대 영상, 음영기복도, 입사각이 90.0° 이상일 때 순그림자 탐지 결과이다. (a), (b), (c)의 순그림자 탐지 결과를 보면 탐지결과가 매우 적은 것을 확인할 수 있다. 이 결과를 바탕으로 반 그림자를 포함한 지형그림자 임계값은 입사각이 90.0° 이상일 때가 아니며 최적의 입사각을 확인해볼 필요가 있음을 알 수 있다.

Fig. 10. Results of hillshade map and true terrain shadow detection in each image : (a) Scene-1, (b) Scene-2, and (c) Scene-3 (From the left: original image, subset image, hillshade map of subset image, result of true shadow detection overlaid with subset image).​​​​​​​

3.2. 반그림자 탐지를 위한 임계각 분석결과

반그림자 탐지를 위해서 각 영상별로 제작된 음영기복도에서 입사각 임계값을 1°씩 조절하여 지형그림자를 탐지하고, 각각의 탐지결과를 참값 지형그림자와 비교하여 탐지 정확도를 비교하였다. Fig. 11은 20개 데이터 중 일부 실험데이터의 정확도 분석을 한 그래프를 나타낸다. Y축은 성능지표의 값을 의미하며 X축은 입사각 임계값을 의미한다. 즉 입사각이 그래프의 X축의 값보다 큰 경우는 지형그림자로 판단하였다. 각 그래프에서 파란색 선은 정밀도이며 회색 선은 재현율, 빨간색 점선은 F1 score를 의미한다. 각 그래프 마다 빨간색 화살표는 F1 score가 가장 높은 지점을 가리키며 초록색 화살표은 입사각이 태양의 천정각일 때를 가리킨다. (a)부터 (j)는 각 영상마다 정확도 결과를 나타낸다.

Fig. 11. A graph of accuracy analysis results according to incidence angle, (a–j) mean the accuracy of each scene (red arrow indicated when F1 score is highest, green arrow indicated when the incidence angle is the sun’s zenith angle).​​​​​​​

Table 4는 각 영상마다 F1 score가 가장 높을 때의 입사각과 정밀도, 재현율을 나타낸다. Table 4와 Fig. 11을 보면 각 영상마다 입사각의 임계값이 다른 것을 확인할 수 있다. 동일지역인 Scene 13과 14는 임계값이 같지만 동일지역인 Scene 8과 9는 임계값이 다른 것을 확인하였다. 이러한 이유는 동일지역일지라도 태양의 위치가 달라지면 임계값이 크게 달라지는 것이라고 판단하였다. 또한 모든 그래프의 초록색 화살표를 보면 정밀도와 F1 score가 매우 낮아지며 그 뒤로 값이 수렴하는 것을 확인할 수 있다. 이러한 결과로 지형그림자의 입사각 임계값 결정은 태양의 천정각보다 커야 된다고 판단하였다.

Table 4. Results of accuracy analysis when F1 score for each satellite image is highest​​​​​​​

마지막으로 영상마다 입사각과 태양의 천정각의 관계를 파악하였다. 식(9)에 따라서 반그림자를 결정하기 위한 지표면 고도각을 결정하기 위해서 Table 5와 같이 각 영상마다 태양 천정각, 임계값, 임계값-천정각으로 구분하였다.

Table 5. Results of determining the optimal terrain shadow threshold

SD: standard deviation

임계값-천정각은 지형그림자 임계값을 결정하는 고도각과 같다. 각 영상마다 고도각은 편차가 존재하였으나 편차에 특별한 경향성은 관측되지 않았다. 한반도를 촬영한 영상을 대상으로 자동적인 지형그림자 탐지결과를 얻기 위해서 고도각은 평균과 표준편차를 계산하였다. 계산결과는 Table 5와 같이 평균이 8.7°이고 표준편차가 ±2.9°로 산출되었다.

본 논문의 최적 임계값이 Richter et al. (2009)의 결과와 일치되지 않았으나 각 임계값이 태양의 천정각 +10°와 유사한 결과를 보여주었다. Balthazar et al. (2012)의 결과에 따르면, 이러한 결과가 실험지역의 특성에 따라 다르기 때문이라고 하였다. 따라서 본 연구는 실험에 사용된 영상에서 한정적으로 적용되는 전제조건을 가진다. 향후에는 실험에 사용된 영상 외의 다른 시기의 영상을 이용하여 재검증하는 과정이 필요하다. 결론적으로, 본 연구에서는 언급한 전제조건에 한하여 한반도를 대상으로 하는 최적 임계값이 태양의 천정각 +8.7°로 결정하였다.

3.3. 최적 임계값을 적용한 지형그림자 탐지 정확도 결과

최적 임계값을 적용한 지형그림자 탐지 정확도 분석은 정량적, 정성적 분석 방법을 통해 수행하였다. 정량적 분석 결과는 Table 6에 정리하였다. 앞서 1°씩 임계각을 조절하며 참값과 비교하여 영상 별로 찾은 최적의 정확도 결과에서 20개 영상의 F1 score는 평균 0.91였다. 반면, 앞선 분석을 통해서 얻어진 최적 임계값을 20개 영상에 공통 적용하여 얻은 F1 score는 평균 0.86였다. F1 score가 낮아지는 영상들이 존재하였으나 대부분 F1 score가 0.8 이상임을 확인하였다. F1 score가 많이 낮아지는 영상들은 Scene-2, 3, 8 영상이었다. 그 이유는 최적 임계값이 앞서 분석한 최고 성능의 입사각 임계값보다 최대 6.3°로 차이 나기 때문이며 이러한 차이의 원인은 추후 연구가 필요하다.

Table 6. Accuracy results of total terrain shadow detection with optimal threshold​​​​​​​

Fig. 12는 정성적 결과를 보여준다. 위쪽부터 위성영상, 확대된 위성영상, 참값데이터, 탐지결과를 나타낸다. (a)는 Scene-3 영상, (b)는 Scene-4의 영상, (c)는 Scene-12 영상, (d)는 Scene-19 영상이다.

Fig. 12. Terrain shadow detection result from RapidEye image and Sentinel-2 image (From the upper: original image, subset image, ground truth in subset image, proposed method results of subset image): (a) Scene-3, (b) Scene-4, (c) Scene-12, and (d) Scene-19 (Gray is non-shadow, and black is terrain shadow).​​​​​​​

정성적 분석 결과, 전체적으로 참값데이터와 유사하나 부분적으로 과탐지나 미탐지가 발생한 것을 확인하였다. 정량적 분석 결과로는 F1 score 결과가 대부분 0.8–0.94를 유지하였다. 따라서 남북한 지역을 대상으로 자동적으로 지형그림자를 탐지하기위해서 지형그림자 탐지 최적 임계값이 천정각 +8.7°일때 지형그림자 탐지가 잘 수행된 것을 확인하였다.

4. 결론

본 논문은 향후 발사될 농림위성의 지형그림자 마스크 생성을 위해 농림위성의 공간해상도와 유사한 Rapideye 위성영상과 Sentinel-2 위성영상을 이용하여 지형그림자 탐지 방법을 제안하였다. 지형그림자 탐지를 위해 본 논문은 음영기복 알고리즘을 적용하였다. 보정을 위한 지형그림자 탐지 연구는 지형보정 기법으로 대체 가능하기 때문에 지형그림자 탐지 목적을 확장시켰다. 이를 위해 본 논문에서는 지형그림자를 태양광이 지표면에 도달되지 않는 순그림자 영역과 태양광이 지표면에 도달하지만 지표면의 경사로 인해 반사도가 감소되는 반그림자 영역으로 구분하였고 두가지 영역 모두를 포함하여 지형그림자라고 정의하였다. 이후 반그림자 영역에 대한 임계값을 분석하였다. 최종적으로 실험에 사용된 한반도를 촬영한 영상을 대상으로 가장 적합한 임계값을 찾았으며 대부분 F1 score가 0.8에서 0.94 사이임을 확인하였다. 정성적 분석으로는 지형그림자의 형상이 참값과 매우 유사함을 확인하였다. 따라서 본 연구결과를 통해 한반도를 대상으로 자동적으로 균일한 품질의 지형그림자 마스크가 제공될 것이라 기대한다. 향후 연구는 Richter et al. (2009)의 연구결과처럼 태양 천정각에 따른 임계값 비교연구와 최적의 임계값을 자동으로 찾을 수 있는 기법에 대한 연구가 필요하다.