Evaluation of the DCT-PLS Method for Spatial Gap Filling of Gridded Data

격자자료 결측복원을 위한 DCT-PLS 기법의 활용성 평가

Abstract

Long time-series gridded data is crucial for the analyses of Earth environmental changes. Climate reanalysis and satellite images are now used as global-scale periodical and quantitative information for the atmosphere and land surface. This paper examines the feasibility of DCT-PLS (penalized least square regression based on discrete cosine transform) for the spatial gap filling of gridded data through the experiments for multiple variables. Because gap-free data is required for an objective comparison of original with gap-filled data, we used LDAPS (Local Data Assimilation and Prediction System) daily data and MODIS (Moderate Resolution Imaging Spectroradiometer) monthly products. In the experiments for relative humidity, wind speed, LST (land surface temperature), and NDVI (normalized difference vegetation index), we made sure that randomly generated gaps were retrieved very similar to the original data. The correlation coefficients were over 0.95 for the four variables. Because the DCT-PLS method does not require ancillary data and can refer to both spatial and temporal information with a fast computation, it can be applied to operative systems for satellite data processing.

지구환경 변화를 파악하는 데 있어서는 장기 시계열의 격자자료가 필수적이며, 기후 재분석장과 위성자료는 대기 및 지표면 변수에 대하여 전 지구 규모에서 주기적이고 정량적인 정보로 활용되고 있다. 본 연구에서는 위성자료의 결측 문제를 해결하기 위한 방안으로 DCT-PLS (penalized least square regression based on discrete cosine transform) 기반의 결측복원 기법을 서로 다른 특성을 가진 복수의 격자자료에 적용하고, 정량적인 검증을 통하여 그 활용성을 평가하였다. 원본 자료와의 객관적인 비교를 위하여 결측이 없는 LDAPS (Local Data Assimilation and Prediction System) 모델로부터 상대습도, 풍속 일자료를 추출하고, MODIS (Moderate Resolution Imaging Spectroradiometer)의 월간 합성 LST (land surface temperature), NDVI (normalized difference vegetation index) 영상을 사용하여, 임의로 생성된 결측 블록이 원본에 매우 가깝게 복원됨을 확인하였고, 4가지 변수 모두에서 상관계수 0.95 이상의 일치도를 나타내었다. DCT-PLS 기반 결측복원 기법은 별도의 보조자료를 필요로 하지 않고, 필요시 시간 및 공간 정보를 모두 활용할 수 있으며, 처리속도가 비교적 빠르기 때문에 현업시스템에 사용될 수 있을 것으로 사료된다.

요약

지구환경 변화를 파악하는 데 있어서는 장기 시계열의 격자자료가 필수적이며, 기후 재분석장과 위성 자료는 대기 및 지표면 변수에 대하여 전 지구 규모에서 주기적이고 정량적인 정보로 활용되고 있다. 본 연구에서는 위성자료의 결측 문제를 해결하기 위한 방안으로 DCT-PLS (penalized least square regression based on discrete cosine transform) 기반의 결측복원 기법을 서로 다른 특성을 가진 복수의 격자자료에 적용하고, 정량적인 검증을 통하여 그 활용성을 평가하였다. 원본 자료와의 객관적인 비교를 위하여 결측이 없는 LDAPS (Local Data Assimilation and Prediction System) 모델로부터 상대습도, 풍속 일자료를 추출하고, MODIS (Moderate Resolution Imaging Spectroradiometer)의 월간 합성 LST (land surface temperature), NDVI (normalized difference vegetation index) 영상을 사용하여, 임의로 생성된 결측 블록이 원본에 매우 가깝게 복원됨을 확인하였고, 4가지 변수 모두에서 상관계수 0.95 이상의 일치도를 나타내었다. DCT-PLS 기반 결측복원 기법은 별도의 보조자료를 필요로 하지 않고, 필요시 시간 및 공간 정보를 모두 활용할 수 있으며, 처리속도가 비교적 빠르기 때문에 현업시스템에 사용될 수 있을 것으로 사료된다.

1. 서론

최근 기후변화로 인해 한반도를 포함하여 전 세계 기상 및 해양현상의 특성이 빠르게 변화하고 있다 (Lee, 2018). 이러한 지구환경 변화를 파악하는 데 있어서 장기 시계열 자료가 필수적이며(Kim and Cho, 2019), 기후재분석장과 위성자료는 대기 및 지표면 변수에 대하여 전 지구 규모에서 주기적이고 정량적인 정보로 활용되고 있다(Tralli et al., 2005; Cho et al., 2013). 또한 가용한 자료의 종류가 늘어나고 시공간 해상도가 향상되고 있기 때문에 다양한 규모에서의 시공간 변화를 분석하기에 적합하다. 그러나 광학센서 기반의 위성자료는 구름 등에 따른 결측을 피할 수 없는데, 일정 수준 이상의 결측은 자료의 불확실성을 크게 하므로 시공간적으로 연속적인 지구환경 분석에 어려움이 따른다. 따라서 위성자료의 결측 복원을 통한 시공간적 연속성의 확보는 매우 중요한 연구 대상이다.

위성자료를 비롯한 시공간 데이터 셋의 결측 복원을 위하여 여러 기법의 연구가 진행되어 왔다. 경험직교함수(Empirical Orthogonal Function, EOF)로부터 도출한 공간 변동성 정보를 이용하여 결측 복원을 수행하거나 (Beckers and Rixen, 2003; Alvera-Azcárate et al., 2007), 크리깅(Kriging)과 같은 공간통계법을 이용한 결측 복원 (Kostopoulou, 2020), 그리고 시공간 변동성에 대해 단순 스펙트럼 분석(Singular Spectrum Analysis, SSA)을 이용하기도 했다 (Schoellhamer, 2001; Kondrashov and Ghil, 2006; Hocke and Kämpfer, 2009). 또 다른 접근으로, 타깃 변수와 설명인자 사이의 지구물리학적 관계를 이용하는 회귀크리깅(Regression Kriging) 또는 그 관계의 비선형성을 고려하는 인공지능 기반의 연구도 진행 중인데 (Ren et al., 2019; Llamas et al., 2020), 이러한 고급 통계법에는 보다 더 높은 컴퓨팅 파워가 요구되며(Katzfuss and Cressie, 2012; Eidsvik et al., 2014), 설명인자를 확보하지 못하면 사용할 수 없다는 단점도 존재한다.

이처럼 지구환경정보를 포함하는 격자자료 결측 복원의 중요성과 여러 노력에도 불구하고, 특히 우리나라 지역을 대상으로 위성자료나 기상격자 자료의 결측 복원에 대한 객관적인 평가나 정량적인 검증이 제시된 사례는 거의 없다. 이에, 본 연구에서는 동일한 시공간 해상도의 설명 변수가 존재하지 않는 상황에서 시계열 격 자자료의 시간적 및 공간적 연관성 정보를 활용하는 DCT-PLS(penalized least square regression based on discrete cosine transform) (Garcia, 2010; Wang, 2012) 기법을 사용하여 격자 자료의 결측 복원을 수행하고, 그 정확도를 정량적으로 평가하고자 한다. 다양한 경우에 대한 정량적인 평가를 위해서는 결측에 해당하는 시공간에 참값이 존재해야 하기 때문에, 실제로 결측이 거의 없는 격자자료를 대상으로 임의의 시공간에 결측을 생성하고 그 복원 값에 대한 검증을 수행하는 것이 효과적이다. 이 점에 착안하여, 본 연구에서는 기상청 국지예보모델 LDAPS(Local Data Assimilation and Prediction System) 자료와 MODIS(Moderate Resolution Imaging Spectroradiometer) 월 단위 산출물을 선택하여 결측 복원 실험을 실시하였다. LDAPS는 모델 자료이므로 기상장에 결측이 없고, MODIS의 LST(Land Surface Temperature), NDVI(Normalized Difference Vegetation Index) 월 단위 산출물은 NASA의 후처리에 의해 결측이 거의 없는 상태로 가공되므로, DCT-PLS 기반의 결측 복원 성능 평가에 적합할 것으로 판단하였다. 기법을 적용하여 원본 자료와의 정량적인 대조를 통해 성능 평가를 수행하였다. LDAPS는 2019년 3월 1일부터 11월 30일까지 275일에 대하여 매일 정오의 상대습도, 풍속격자자료를 사용하였고, MODIS 월간 합성LST는 MYD11C3v006, 월간 합성NDVI는 MYD13C2v006 산출물에서 2003년부터 2019년까지 각 204개월 영상을 사용하였으며, 매장면 마다 임의의 결 측 블록을 생성하고 이를 복원한 후, 검증정확도를 산출하는 방식으로 실험을 수행하였다.

2. 자료와 방법

1) LDAPS 기상자료

우리나라 기상청 수치예보모델 중에서 단기예보 모델은 영국 기상청에서 개발한 통합수치모델(Unified Model, UM)을 사용하며, 예보 지역에 따라 전 지구 예보 시스템(Global Data Assimilation and Prediction System, GDAPS), 지역 예보시스템(Regional Data Assimilation and Prediction System, RDAPS), 국지 예보시스템(LDPAS)으로 구분한다. LDAPS는 1.5 km 격자로 연직 40 km 높이까지 70개의 층으로 구성되며, 수치 예보 모델의 특성상 초기 및 경계조건은 GDAPS와 동일하게 사용한다. 예보는 1일 8회(00, 03, 06, 09, 12, 15, 18, 21 UTC) 3시간 간격으로 수행되는데, 이때 예보 오차를 실시간으로 계산하여 다음 예보에 반영하는 3차원 변분자료 동화기법 (3D variational data assimilation, 3DVAR)을 사용하여 오차를 최소화한다(https://data.kma.go.kr/cmmn/main.do). 단일면(지표면) 변수는 78종이며, 그 중에서 2019년 3월 1일부터 11월 30일까지 275일에 대하여 매일 정오의 상대습도, 풍속 변수를 사용하였다.

2) MODIS 위성자료

Aqua 위성은 고도 705 km에서 1일 2회 한반도 상공을 지나며, 낮 시간대의 경우 남에서 북쪽으로 13:30분경 통과한다. Aqua 위성의 LST 산출물 MYD11C3v006, NDVI 산출물 MYD13C2v006은 0.05° 해상도의 월간 합성자료로서, LP DAAC(Land Processes Distributed Active Archive Center)로부터 2003년부터 2019년까지 각 204개월 영상을 취득하여 사용하였다.

3) 결측복원 기법

DCT-PLS 결측 복원 기법은 다차원 자료의 평활화(smoothing)를 위한 것으로서 (Garcia, 2010; Wang, 2012), 이산코사인변환(Discrete Cosine Transform, DCT)을 이용한 벌점최소제곱법회귀(Penalized Least SquareRegression, PLS)를 기반으로 한다. PLS의 기본 원리는 원자료와 평활화 자료 사이의 균형을 위하여, 원자료와 평활화 자료 사이의 편차를 나타내는 잔차제곱합(residual sumof-squares, RSS)과 평활화 자료 간의 거칠기를 반영하는 페널티항(P)으로 구성된 식(1)을 최소화하는 것이다 (Whittaker, 1923; Wahba, 1990; Eilers, 2003).

\(F(\hat{y})=R S S+s P(\hat{y})=\|\hat{y}-y\|^{2}+s P(\hat{y})\)       (1)

최근 연구에서 PLS이 DCT에 의해 공식화될 수 있다는 것이 입증되었는데 (Garcia, 2010; Wang, 2012), DCT 는 다차원 적용이 가능하므로, DCT-PLS를 시공간으로 확장하여 격자자료의 결측 복원을 적용할 수 있다

\(F(\hat{X})=\left\|\mathrm{W}^{1 / 2^{\circ}}(\hat{X}-X)\right\|^{2}+\mathrm{s}\left\|\nabla^{2} \hat{X}\right\|^{2}\)      (2)

여기에서, X는 결측이 존재하는 원자료이며 W는 X와 동일한 크기의 가중치 행렬이다. 가중치는 결측이 없는 픽셀은 1, 결측치는 0을 주어 이진 행렬을 구성한다. 또한 || · ||는 유클리드 놈(Euclidean norm), ∇2는 라플라스 연산자(Laplace operator)를 나타내며 °는 슈어(Schur)의 곱을 나타낸다. s는 \(\hat X\)의 평활화 정도에 관여하는 매개변수로 이를 적절히 조정함으로써 과대 또는 과소평 활화(over- or under-smoothing)를 피할 수 있다.

\(\begin{array}{l} \hat{X}=\operatorname{IDCT}\left(\Gamma^{\circ} D C T\left(W^{\circ}(\hat{X}-X)+X\right)\right), \\ \Gamma_{i_{1}, i_{2}, i_{3}}=\left(1+\mathrm{s}\left(\sum_{j=1}^{3}\left(2-\cos \frac{\left(i_{j}-1\right) \pi}{n_{j}}\right)\right)^{2}\right)^{-1} \end{array}\)        (3)

자료가 3차원인 경우 식(2)는 DCT를 이용한 식 (3)과 같이 나타낼 수 있으며, i_j 는 j차원의 i번째 성분, n_j 는 3차원 자료 X행열의 크기를 나타낸다. 이때 3차원은 XYZ 공간뿐만 아니라 XYT의 시공간에도 적용할 수 있다.

4) 실험과정

DCT-PLS 기반 결측 복원 기법의 객관적인 성능 평가를 위하여, 결측이 거의 없는 LDAPS 자료와 MODIS 월간 합성자료에 10×10픽셀의 결측블록 2개를 랜덤하게 생성하고, 이를 복원한 후, 원자료와 픽셀별로 비교하여 검증통계량을 계산함으로써 결측 복원 성능을 평가하였다. 연구영역은 북위 33~38.5°, 동경 125.8~129.8°의 남한 일대를 대상으로 하였으며, LDAPS 상대 습도, 풍속 자료와 MODIS LST, NDVI 자료의 해상도는 0.05°로 통일하여 110×80픽셀을 구성하였으며, LST와 NDVI는 결측 블록이 육상에만 생성되도록 설정하였다. 상대 습도, 풍속, LST, NDVI에 대하여 2차원 DCT-PLS 결측 복원을 수행하고, 복원된 모든 픽셀을 원자료와 비교하여 평균편의오차(mean bias error, MBE), 평균 절대오차(mean absolute error, MAE), 평균제곱근오차(root mean square error, RMSE) 그리고 상관계수(correlation coefficient, CC)를 계산하였다. 또한, 매우 전형적인 코싸인 곡선 형태의 점진적인 시계열 변화를 가지는 NDVI의 경우에는 XYT의 시공간 3차원 DCT-PLS를 추가적으로 수행하여 정확도 향상의 효과가 있는지 검토하였다. 이때 결측 복원하려는 시점을 t라고 할 때, t-1, t+1 시점의 자료는 2차원 DCT-PLS 처리를 한 뒤, t-1, t, t+1 장면에 대하여 3차원 DCT-PLS를 수행하는 방식이며, 이동창(moving window)의 크기는 3으로 하였다 (Fig. 1).

Fig. 1. 3D DCT-PLS gap filling with a moving window.

3. 결과 및 고찰

각 275장의 상대 습도, 풍속격자 격자자료 그리고 각 204장의 LST, NDVI 격자자료에 대하여, 각 장면마다 임의의 생성한 10×10 의결 측 블록 2개, 즉 200개 픽셀에 대하여 DCT-PLS 결측 복원을 수행하였다. LDAPS는 55,000개(275×200), MODIS는 40,800개(204×200)의 결측화소가 생성되지만, 1장의 영상 내에서 결측블록의 교집합이 일부 발생할 수 있고, LST, NDVI와 같은 지표면 변수의 경우 해양에 속하는 픽셀은 제거하였기 때문에, 검증에 사용된 픽셀 수는 상대 습도 54,522개, 풍속 54,502개, LST 37,690개, NDVI 36,831개이다. 이에 대한 정확도 통계량은 Table 1~4에 나타냈고, 또한 Table 5는 NDVI에 대한 3차원 DCT-PLS 결과를 보여준다. 상관계수의 경우, 상대습도 0.967, 풍속 0.958, LST 0.986, NDVI 0.924로서, 복원된 픽셀과 원자료와의 일치도가 상당히 높은 것으로 보인다. MAE는 상대 습도 2.997%, 풍속 0.584 m/s, LST 1.087 K, NDVI 0.054로서 복원 오차가 상당히 적어, 실제 활용이 가능할 것으로 사료된다. 또한 일반적으로 매우 점진적인 시계열 변화를 보이는 NDVI의 경우, 시공간 3차원 DCT-PLS로 결측 복원을 수행한 결과, 상관계 수 0.966, MAE 0.039로서, 2차원 DCT-PLS에 비해 다소 향상된 정확도가 산출되었다. 이처럼 공간 변동성이 비교적 일정하고 시간에 따른 변화가 점진적인 자료의 경우에는, 전후 시간대를 고려한 3차원 DCT-PLS를 추가적으로 사용함으로써 보다 더 개선된 결과를 얻을 수도 있다.

Table 1. Accuracy statistics of gap-filled LDAPS relative humidity using DCT-PLS (Mar-Nov, 2019)

Table 2. Accuracy statistics of gap-filled LDAPS wind speed using DCT-PLS (Mar-Nov, 2019)

Table 3. Accuracy statistics of gap-filled MODIS LST using DCT-PLS (2003-2019)

Table 4. Accuracy statistics of gap-filled MODIS NDVI using DCT-PLS (2003-2019)

Table 5. Accuracy statistics of gap-filled MODIS NDVI using 3D DCT-PLS (2003-2019)

Fig. 2~5는 상대 습도, 풍속, LST, NDVI의 2차원 결측 복원 결과를 산점도로 나타낸 것이고, Fig. 6은 NDVI의 3차원 결측 복원 결과이다. 4가지 기상변수 모두, 복원 픽셀과 원자료의 관계가 1:1선에 가깝게 분포하는 것을 확인할 수 있다.

Fig. 2. Observed vs. predicted LDAPS relative humidity using DCT-PLS (Mar-Nov, 2019).

Fig. 3. Observed vs. predicted LDAPS wind speed using DCT-PLS (Mar-Nov, 2019).

Fig. 4. Observed vs. predicted MODIS LST using DCT-PLS (2003-2019).

Fig. 5. Observed vs. predicted MODIS NDVI using DCT-PLS (2003-2019).

Fig. 6. Observed vs. predicted MODIS LST using 3D DCT-PLS (2003-2019).

Fig. 7~10은 각각 상대 습도, 풍속, LST, NDVI에 대한 원본 자료, 임의 생성된 결측블록, 그리고 복원 사례를 나타내는 지도이다. 대부분의 경우에 10×10의 결측 블록 2개가 원자료와 상당히 유사한 분포로 복원된 것을 확인할 수 있다. 상대 습도는 해양에서 전반적으로 높고, 육지가 해양보다는 낮은 값의 대조적인 패턴이 유지되기 때문에, 해양과 육지가 만나는 해안선 부근에서의 결측 복원이 쉽지 않음에도 불구하고, Fig. 7의 3월 20일 사례에서처럼 해안선 부근의 결측 블록도 상당히 원자료에 가깝게 복원되는 것을 확인하였다. 또한 110×80의 8,800픽셀로 구성된 1장의 영상을 처리하는 데 걸리는 시간이 i5 또는 i7 PC에서 0.3초 미만으로 소요되므로 처리 속도가 매우 빠르다.

Fig. 7. Maps of LDAPS relative humidity retrieved by DCT-PLS.

Fig. 8. Maps of LDAPS wind speed retrieved by DCT-PLS.

Fig. 9. Maps of MODIS LST retrieved by DCT-PLS.

Fig. 10. Maps of MODIS NDVI retrieved by 3D DCT-PLS.

4. 결론

본 연구에서는 위성자료의 결측 문제를 해결하기 위한 방안으로 DCT-PLS 기반의 결측 복원 기법을 서로 다른 특성을 가진 복수의 격자자료에 적용하고, 정량적인 검증을 통하여 그 활용성을 평가하였다. 원본 자료와의 객관적인 비교를 위하여 결측이 없는 LDAPS 모델로부터 상대 습도, 풍속 변수를 추출하고, MODIS의 월간 합성 LST, NDVI 영상을 사용하여, 임의로 생성된 결측 블록이 원본에 매우 가깝게 복원됨을 확인하였고, 이로써 DCT-PLS 기반의 결측 복원 기법이 실제로 온도, 습도, 식생 등의 위성 산출물에 적용 가능함이 입증되었다. 특히 시계열 변화 패턴이 매우 점진적인 식생지수의 경우에는 전후 시점의 영상을 참조하는 3D DCT-PLS를 통하여 보다 더 복원 정확도가 향상될 수 있을 것으로 보인다. DCT-PLS 기반 결측 복원 기법은 별도의 보조자료를 필요로 하지 않고, 필요시 시간 및 공간정보를 모두 활용할 수 있으며, 처리 속도가 비교적 빠르기 때문에 현업시스템에 사용될 수 있을 것으로 사료된다. 향후에는 고해상도 위성영상을 포함하여 다양한 위성영상을 대상으로 적용 및 평가가 필요할 것이다.