1. 서론

산림생태계 내의 총일차생산량은 바이오매스(목재, 임산물 등) 자원 생산량과 직결되고, 산림생태계의 건강성, 산림식물계절 및 생태계 서비스의 중요한 지표가 된다(Ruidisch et al., 2014; Robinson et al., 2014). 총일차생산량은 일정 기간 동안 식물의 광합성을 통해 발생된 탄소의 총량으로 생물권과 대기권을 연결하는 탄소순환과정을 평가할 수 있는 요소이다(Beer et al., 2010; Wu et al., 2010; Balzarolo et al., 2015). 전지구적 및 지역적 차원의 총일차생산량의 측정과 관측이 필요한 이유는 육상생태계의 이산화탄소 고정에 대한 보다 정확한 관측정보가 다양한 생태계 기작들을 평가할 수 있는 기초자료가 되기 때문이다(Coops et al., 1998; Running et al., 1999).

우리나라 산림생태계에서도 기후변화 영향 평가 보고 법적 의무화(매 5년)에 따른 전국 단위 산림자원 및 산림생태계 변화 조사를 매년 실시하게 됨에 따라 산림생산성도 산림생태계서비스 즉, 생물량 기반의 에너지와 영양물질 교환, 생물다양성 보존, 온실가스 감축 등에 연계된 지표로서 전반적인 현황 파악이 중요하게 되었다.

산림의 총일차생산량을 관측하는 방법에는 에디공분산 타워 같은 대기 중 이산화탄소 흐름을 측정하여 계산하거나 인공위성영상자료를 사용한다(Kang et al., 2014). 에디공분산 측정 방법은 생태계와 대기 간의 이산화탄소를 비롯한 물질과 에너지의 흐름을 관측하여 정량화 할 수 있는 검증된 방법이다(Aubinet et al., 2000) 그러나 에디공분산타워는 다른 실측 방법과 마찬가지로 비용적인 측면에서 그 수가 제한적일 수 밖에 없고 타워 주변의 관측 범위(footprint)만 대표한다는 공간적한계가 있다(Wang et al., 2017; Ogutu et al., 2013). 인공위성 자료의 경우, 넓은 지역과 직접 측정이 불가능한 곳을 관측할 수 있으나 인공위성자료에서 산출된 총일차생산량 역시 자료의 신뢰성에 지속적인 논란이 되고 있고 특히 복잡지형에서는 대부분의 연구자들은 연구 지역에 맞게 수정 보완하여 사용하고 있다(Cramer et al., 1999; Reeves et al., 2007; Zhang et al., 2007). 특히 우리나라의 경우, 산림이 국토의 약 63%를 차지하고 있고 그 대부분이 복잡지형에 위치하고 있어 우리나라 산림에서의 인공위성 총일차생산량 산출물은 불확실성이 높다. 인공위성 자료는 1950년대 이후로 지속적으로 발달하여 축적된 자료의 양이 많아 빅데이터로서의 가치를 가진다. 지구 관측에 관한 위성자료로부터 새로운 정보와 지식을 얻고자 하는 움직임이 미국과 유럽 등 주요 국가들과 구글과 같은 기업에서 발빠르게 이루어 지고 있다(Liu et al., 2018). 위성 자료의 활용은 생태계의 분류와 변화 예측을 보다 정확하게 해주고 기존의 관측 데이터가 가지는 공간적인 한계도 대부분 극복 가능할 수 있게 한다. 특히 산림 생태계와 같이 복잡하고 다양한기작들이 얽혀 있고 공간적으로 직접적인 관측 한계가 있는 분야에서는 그 활용도가 더욱 기대된다(Melaas et al., 2013). 따라서 산림분야에서의 인공위성 자료를 활용하기에 앞서 위성자료로부터 추출된 자료와 지상에서 관측된 정밀도 높은 결과(에디공분산 등)를 비교하여 위성자료가 관측결과를 얼마나 정확하게 예측 할 수 있는지부터 연구되어야 할 것이다.

본 연구에서는 인공위성 자료와 기계학습 알고리즘을 활용하여 우리나라의 산림 유역의 총일차생산량을 추정해 보고 그 활용가능성을 고찰해 보고자 한다. 우리나라의 산림 지역에서 미항공우주국(NASA)의 MODIS (Moderate Resolution Imaging Spectroradiometer) 자료를 입력 자료로 사용하여 기계학습 알고리즘 Support Vector Machine (SVM) 모델을 훈련하고 모델의 예측 결과를 지상 측정 자료(에디공분산 플럭스 타워)의 총일차생산량 및 MODIS 총일차생산량과 비교해 보았다. 모델의 입력자료는 MODIS의 지상 산출물들을 이용하여 새로운 산출물로 계산한 MODIS 처리 입력자료 (Processed MODIS)와 NASA 데이터센터에서 제공되는 MODIS 지상 산출물 자료를 재계산 없이 그대로 사용한 입력자료(Unprocessed MODIS)를 사용하였다. 또한 구축된 모델을 다른 지역들로 적용하였을 때 기계학습 알고리즘 모델의 예측력을 분석하고 향후 전국 산림으로의 확장 가능성을 평가해 보고자 한다.

이 연구에서 사용된 기계학습 알고리즘은 Support Vector Machine (SVM)으 로 Artificial Neural Network (ANN), Random Forest (RF)와 함께 산림생태계 연구에서 가장 많이 사용되고 있는 알고리즘이다. SVM 알고리즘은 자료 중심 모델에서의 그 능력을 이미 입증한 기계학습 알고리즘 중 하나로 수문(Jajarmizadeh et al., 2014; Granata et al., 2016; Langhammer and Cesak, 2016), 산림과 비산림(Ichii et al., 2017; Dou et al., 2018), 생태계 잠열 (Wang et al., 2017) 등 다양한 생태학적 연구에서 사용되었다. 기존 연구에 의해 세 가지 알고리즘 모두 예측력은 이미 증명되었고 (Lee et al., 2019) 이 연구에 앞서 세 가지 알고리즘을 예비 비교를 해보았는데 본 연구의 주제에는 SVM이 가장 높은 예측력을 가진 것으로 나타나 SVM 알고리즘 모델을 좀 더 자세하게 다루었다.

2. 방법

1) 연구지역 및 에디공분산 자료

연구지역은 대기의 난류를 직접 측정하여 총일차생산량을 추정하는 에디공분산 타워가 있는 산림유역의 6개 지점, 경기도 포천시 소흘읍 광릉숲(GDK), 경기도 파주시 적성면 설마천유역(SMC), 강원도 홍천군(HC) 매화산, 강원도 삼척시(SC) 검봉산, 전라남도 완도군 완도수목원(WD), 제주도 서귀포시 난대시험림(Jeju)을 선택하였다. 본 연구에서 기계학습 알고리즘 모델의 변수의 정답 즉, 목표 변수와 결과의 비교 자료로 사용된 에디공분산 플럭스 자료는 각 지점에 설치된 에디공분산플럭스 타워에서 관측된 것으로 자세한 정보는 Table 1에 제시되어 있다. 광릉숲유역(GDK)은 수령 80년 이상의 고령림으로 주요 수종은 졸참나무(Quercus serrata)와 서어나무(Carpinus laxiflora)이고 설마천유역(SMC)은 수령 20년~50년 사이로 굴참나무(Quercus variaabillis), 신갈나무(Quercus mongolica), 졸참나무(Quercus serrata)가 주 종이다(Kwon et al., 2009). 홍천(HC) 매화산은 수령 40~50년 사이의 혼효림으로 일본잎갈나무(Larix kaempfere), 쪽동백나무(Styrax obassia), 층층나무(Cornus controversa) 등이 우점종이다. 삼척(SC) 검봉산은 수령이 1~10년으로 영급이 낮은 혼효림으로 소나무(Pinus densiflora), 굴참나무 (Quercus variaabillis), 신갈나무(Quercus mongolica), 완도수목원(WD)와 제주의 난대시험림(Jeju) 지역은 활엽수종과 상록 활엽수종이 우점종으로 WD 지역은 붉가시나무 (Quercus acuta), 구실잣밤나무(Castanopsis sieboldii), 후박나무(Machilus thunbergii), 참식나무(Neolitsea sericea), Jeju 지역은 붉가시나무(Quercus acuta), 동백나무(Camellia japonica), 서어나무(Carpinus laxiflora) 등이 혼효되어 있다(Yun and Chun, 2018). 에디공분산 플럭스 타워는 해발고도 260 – 1190 m 사이에 분포하며 전국으로 비교적 넓게 분포 되어 있다. 관측항목은 풍향, 풍속, 3차원 풍속, 온습도, 대기압, 강우량, 복사량, 퀀텀량, CO₂ 농도, H₂O 농도 등이 있고 타워의 높이는 20~40 m이다. 해당 지역의 30년 (1981~2010) 연평균기온은 10.3°C(영월(평창))~16.6°C (제주), 강수량은 1224~1923 mm 이다. 30년 평균 기온과 강수량에 대한 정보는 기상청 홈페이지(http://data. kma.go.kr)를 참조로 하였으며 30년 평년자료가 존재 하지 않는 지역은 가까운 측정지점의 정보를 제시하였다.

Table 1. Site information of eddy covariance flux sites at Gwangneung, Seolmacheon, Hongcheon, Samcheok, Wando, and Jeju, South Korea. Ta is mean air temperature, P is total rainfall

각 연구지역 플럭스 타워 정상의 에디공분산 시스템은 3차원 초음파 풍향풍속 온도계(SAT, Model CSAT3, Campbell Scientific, Inc., USA)와 고속반응 적외선 기체분석기(IRGA, Models LI-7500 (LI-COR, Inc., USA), EC150, and EC155 (Campbell Scientific, Inc.))로 구성되었다. 에디공분산 시스템의 샘플링 속도는 10Hz, 평균시간은 30분, 10Hz 관측과 30분 평균 자료는 자료 집록기(Model CR3000, Campbell Scientific Inc.)에 저장된다. 저장된 자료는 표준화된 KoFlux 자료처리 방법에 따라 에디 공분산 자료처리 응용 프로그램 EddyPro® (LI-COR, Inc.)와MATLAB® (The MathWorks Inc., USA)으로 작성된 KoFlux 자료처리 프로그램으로 처리하였다(Hong et al., 2009; Kang et al., 2018). 표준화된 KoFlux 자료처리 방법은 추세 제거(30분 블록 평균), 시간 지연 보상(추정된 시간 지연 범위 내에서 공분산이 최대가 되는 시점 탐지), 좌표 변환(평면 맞추기 회전(Wilczak et al., 2001; Yuan et al., 2007), 풍향 별로 8개 평면의 경사를 한 달 간격으로 계산하여 적용), SAT 온도의 수증기 보정, 주파수 반응 보정(추세 제거 시 발생하는 고주파수 필터링 효과 보정, SAT와 IRGA 간 분리 효과 보정, 기기 경로 길이, 튜브 감쇄 등의 이유로 발생한 저주파수 필터링 효과 보정), 공기 밀도 보정(Webb et al., 1980), 품질 관리(정상 상태 및 난류 발달 상태 점검), 프로파일 관측자료를 이용한 저류항 계산 및 튀는 자료 제거(Papale et al., 2006), 결측자료 메우기(주변 분포 표본 추출법)와 야간 CO₂ 플럭스 자료 보정 및 배분(마찰속도 필터링 방법, Reichsteinet al., 2005; Kang et al., 2019)을 포함한다. 총일차생산량은 야간 CO₂ 플럭스 자료 보정과정에서 얻어진 기온과 호흡(야간 CO₂ 플럭스) 간의 관계가 주간에도 변화가 없다는 가정 하에 추정한 주간의 호흡량과 주간에 관측된 CO₂ 플럭스의 차이로부터 계산한다.

2) MODIS 입력자료

본 연구에서 사용된 MODIS 자료의 경우 공간해상도(500~1000 m)가 에디공분산 타워가 측정할 수 있는 범위와 다르므로 MODIS의 한 픽셀 값이 얼마나 대표성을 갖는지 확인할 필요가 있다. 에디공분산 플럭스타워가 위치한 곳을 임상도와 함께 확인해 본 결과, 해당 지역의 MODIS 픽셀은 제주(Jeju) 지역을 제외하고 나머지 모두 90% 이상이 산림으로 구성된 픽셀이었다. 제주(Jeju) 지역도 픽셀의 80%가 산림으로 해당 지역들의 MODIS 값들은 산림지역의 값으로 대표 될 수 있다고 가정하였다. 산림의 종류는 완도(WD)의 경우 94%가 활엽수림, 삼척(SC)은 약 92%가 침엽수림으로 구성되어 있고 나머지 지역들은 약 40~68% 활엽수림, 약 11~48 % 침엽수림, 약 9~27 %가 혼효림으로 분류되어 있다.

이 연구에서 사용된 MODIS자료는 표 2에 자세히 나와 있다. MODIS는 NASA에서 개발된 극궤도 위성 Terra와 Aqua에 탑재된 다중분광센서로 전지구적 규모의 대기와 육상 및 해양 모니터링 목적으로 개발되었다. Terra는 10:30, Aqua는 1:30분 하루에 두 번 우리나라 상공을 통과하면서 영상을 수집한다. MODIS 영상은 250 m – 5600 m 의 공간해상도를 가지며 Terra는 1999년부터 Aqua는 2002년부터 다양한 지구 모니터링 자료를 축적하였다.

(1) Processed MODIS

MODIS 자료는 대기보정 및 방사보정을 하여 제공되는 자료를 사용하였고 총일차생산량(Gross Primary Production, GPP) 추정을 위해 필요한 변수들은 MODIS 총일차생산량 알고리즘을 기준으로(Running and Zhao, 2015) 선정하였다.

\(G P P=\left[\varepsilon_{\max } f\left(T_{\min }\right) f(V P D)\right]\left[F P A R\left(0.45 R_{s}\right)\right]\)

총일차생산량 알고리즘에 사용되는 입력자료는 다음의 설명과 같이 몇 단계의 산출 과정을 거쳐 생산된다. 이 자료들은 기계학습 알고리즘의 설명 변수들로 사용되며 MODIS 처리 입력자료(Processed MODIS)로 명명하였다. ε_max는 최대광이용효율(KgC m-² day^-1 MJ^-1)이고 Running et al. (2000)과 White et al. (2000)에 의해 제안된 식생형에 따른 생물리특성 색인표(Biome Property Look Up Table, BPLUT)를 따랐다. Tmin은 일최저기온으로 Aqua위성에 탑재된 MODIS의 대기프로파일(MYD07_L2)의 지표와 가장 가까운 지점에서 산출된 기온을 사용하였다. Aqua 위성은 새벽 1시 30분경 우리나라 상공을 지나므로 MYD07_L2 자료는 일 최저 기온으로 볼 수 있다(Lee et al., 2011). VPD(Vapor Pressure Deficit)는 포화수증기압차, FPAR(Fraction of Photosynthetically Active Radiation)는 광합성유효복사량 흡수비율(%)을, Rs(Solar radiation)는 일사량(MJ Day^-1)을 나타낸다. VPD는 대기중 수증기와 관련된 기상인자로 이 인자의 생산을 위해 MODIS의 대기프로파일자료인 MYD07_L2로부터 산출한 지표 부근의 기온과 이슬점온도를 사용하였다. 두 가지 기온 자료를 바탕으로 Dingman (2008)의 방법에 따라 현재수증기압(Vapor Pressure, VP), 포화수증기압(Saturated Vapor Pressure, SVP), VPD를 산출하였다. 일 단위 일사량 추정에는 MODIS 알베도(MCD43)와 AMSR-E(Advanced Microwave Scanning Radiometer)를 사용하였고, Bird and Hulstrom (1981)에 의해 제안되고 Houborg and Soegaard (2004)와 Ryu et al. (2008) 등에 의해 검증된 simple clear sky model (Bird model)을 적용하였다. 광합성유효복사량(Photosynthetically Active Radiation, PAR)은 지표에 입사하는 일사량을 총량에서 식물이 흡수하는 가시광선의 평균 비율인 45%(0.45)를 일사량에 곱하여 산출할 수 있다(Larcher, 2003). 자세한 방법은 Jang et al. (2009)에 제시되어 있으므로 참조할 수 있다.

(2) Unprocessed MODIS 입력자료

포화수증기압차(VPD)나 일사량 등은 비교적 복잡한 산출 과정을 거쳐야 한다. 따라서 이 알고리즘에 대해 익숙하지 못하면 입력 자료를 생산하는데 어려움이 있다. 이 연구에서는 MODIS 생산 자료들 간의 복잡한처리 과정을 거치지 않고 MODIS 데이터 센터에서 제공되는 자료를 별도의 계산 과정 없이 그대로 사용 가능한 입력자료들로 선정해 기계학습 알고리즘 모델의 설명변수로 사용하였고 MODIS 무처리 입력자료 (Unprocessed MODIS)로 명명하였다. 설명변수들은 Processed MODIS 자료와 마찬가지로 MODIS 총일차 생산량 알고리즘을 기준으로 별도의 계산 과정 없이 사용 가능한 것과 대체 가능한 인자들을 선택하였고 지표면 온도(Land surface temperature, LST), FPAR, 잠재증발산(Potential Evapotranspiration, PET), 개량식생지수 (Enhanced Vegetation Index, EVI)를 선택하였다. FPAR의 경우, Processed MODIS를 만들기 위해 사용된 전천 후 기상입력자료가 Unprocessed MODIS에서는 사용되지 않으므로 10:30분경에 상공을 통과하는 위성인 Terra의 자료를 사용하였다(MOD15). 사용된 MODIS 자료에 대한 정보는 Table 2에 제시되어 있다. 잠재증발산 (PET)은 표면 알베도, 공기 중 수분, 에너지 등의 적절한 조건이 갖추어지면 발생하는 증발산으로 잠재증발산 (PET)이 증가하면 포화수증기압차(VPD)가 증가한다 는 두 요인 간의 상관관계는 많은 연구들에 의해 증명되었기 때문에 잠재증발산(PET)을 포화수증기압차(VPD)의 대체 변수로 선택하였다(Granger and Gray, 1989; Novick et al., 2015; Massmann et al., 2018). 개량식생 지수(EVI)는 식생 커버를 민감하게 표현해 주는 지수로 특히 우리나라 산림과 같이 식생이 밀접하게 모여 있는 곳에서 식생 커버와 높은 연관성을 보이고 특히 산림에서 총일차생산량과 높은 상관성을 가지고 있다(Sims et al., 2006; Wardlow and Egbert, 2010; Lee et al., 2017).

Table 2. Summary of MODIS Products used in the input variables from 2006 to 2017

* MYD11 used for both input data sets

3) 기계학습 알고리즘

이 연구에서 사용된 기계학습 알고리즘은 Support Vector Machine (SVM)으로 Random Forest (RF)와 함께 산림생태계 연구에서 가장 많이 사용되고 있는 알고리즘이다. SVM 알고리즘은 단순회귀식과 사례기반의 최근접이웃 학습 방법을 조합한 알고리즘으로 비교적 복잡한 관계도 모사 가능하다. SVM 알고리즘은 다차원에 그려진 자료들 간의 설명 가능한 특징을 가지는 최적화된 평면 경계(초평면, hyperplane)를 그린다고 가정 할 수 있다. 예를 들어 Fig. 1과 같은 자료가 있다고 가정할 때 이 자료들을 두 그룹으로 나누는 A, B, C와 같은 기준선을 정해야 한다면 SVM 알고리즘은 A, B, C 중 어떤 기준선이 최적의 기준선을 판단해야 한다. 이 때 ‘초평면’으로 명명된 이 경계는 데이터 사이의 마진(margin)이 가장 큰 기준선을 선택하고 그 기준선으로부터 최적화된 분류 기준으로 그룹화 되도록 분류한다. 이 때 서포트벡터(support vector)라고 명명된 포인트가 초평면에 가장 가까운 데이터로 최대 마진을 설명할 수 있는 포인트가 된다(Fig. 1(b)). 따라서 SVM 알고리즘의 목표는 서포트벡터(support vectors)의 최대 마진이 되는 포인트를 구하기 위해 두 평면 사이의 거리(유클리디안거리, Euclidean norm, 즉 원점부터 벡터 w까지의 거리)를 최소화하는 값을 구하는 방향으로 학습한다. 선형적으로 자료들이 나누어지는 관계에서는 모든 입력 자료에 대해 유클리디안거리가 최소화하는 값을 구하면 되지만 대부분의 자료들은 정확히 선형적으로 나누어 지지 않는다. 비선형의 관계가 많은 생태계 시스템에서는 비선형 관계의 독립 변수들을 고차원 공간으로 배치하여 선형관계로 변환시키기 위해 커널 기능(kernel function)을 사용하게 된다. 커널 기능에는 시그모이드, 가우시안, 다항식, 선형 커널 등이 있고 자료와 그 목적에 맞게 커널을 선택할 수 있다. 본 연구에서는 기본 커널 radial 커널을 사용하였다(Cortes and Vapnik, 1995). 본 연구에서는 프로그램 R (version 3.3.3), SVM알고리즘을 위한 패키지 “e1071”을 사용하였고(Dimitriadou et al., 2008; Smola and Schölkopf, 2002), 자료에 맞는 최적의 SVM 알고리즘 모델을 만들기 위해 10-fold cross validation 실행하여 best parameters를 찾아서 적용하였다.

Fig. 1. (a) There are two groups of points, shown in closed circle and opened squares. Three possible separating hyperplanes are labeled A, B, and C. (b) Example of maximum margin with support vectors in gray circle and squares that lie on the dashed lines.

3. 결과 및 토의

1) 입력자료에 따른 총일차생산량 추정

위성자료 산출물의 추가적인 가공 정도 즉 복잡한 계산과정을 거쳐 나오는 산출물(Processed MODIS)과 기본적인 보정과 품질검증만 거친 위성자료를 추가적인 계산을 하지 않고를 그대로 사용하였을 때(Unprocessed MODIS) 산림 총일차생산량의 예측력 차이를 알아보기 위해 단일 지역에서 먼저 SVM 알고리즘으로 총일차생산량 추정 모델을 구축하였다. 비교적 장기 에디공분산 플럭스 자료가 있는 광릉숲(GDK)를 대상으로 MODIS 입력자료를 설명변수로 하는 – Processed MODIS, Unprocessed MODIS – 두 가지 세트로 구성하여 모델을 구축하고 모델의 예측결과를 검증하기 위하여 에디공분산 플럭스 타워에서 관측 및 추정된 총일차 생산량을 사용하였다.

모델의 훈련자료는 8일 단위 2006년부터 2013년까지 (n=468),검증 및 비교자료는 2014년과 2015년으로(n=92)설정하였다. Processed MODIS와 Unprocessed MODIS로 각각 훈련된 모델로 2014년과 2015년 산림 총일차생산량을 예측해 보았을 때 에디공분산 총일차생산량과 비교하여 피어슨 결정계수가 0.97로 동일하고 기울기 (Slope)는 Processed MODIS 입력자료로 훈련된 모델은 0.94와 Unprocessed MODIS 입력자료로 훈련된 모델은 0.992로 후자가 근소한 차이로 위성영상으로부터의 추정치와 지상 측정 자료 값의 변화의 차이가 더 작은 것으로 나타났다(Fig. 2). 이는 위성영상에서 추출된 자료들을 복잡한 계산과정을 거친 2차 산출물로 가공하는 단계를 생략하여도 결과값에 매우 큰 영향을 주지 않는 것으로 보인다. 따라서 이후 연구 결과는 Unprocessed MODIS 입력자료를 구축된 모델에 사용하였다.

Fig. 2. The relationship between GPP from SVM model using processed MODIS input data set and GPP from eddy covariance measurement (EC GPP, right panel), and GPP from SVM model using unprocessed MODIS input data set and EC GPP (left panel) for the test years indicated in Table 3 (2014-2015).

NASA에서 MODIS 산출물로 총일차생산량(MOD
17A2 & MYD17A2)이 제공된다. 그러나 MODIS 총일차 생산량의 불확실성은 지속적으로 제기되었는데 이는 전지구적으로 같은 알고리즘을 적용하고 BPLUT에 의존하여 모수화하므로 각 지역마다 다른 기상 조건이나 토지 이용 종류를 반영하기 보다는 전지구적 산출물이기 때문이거나 알고리즘이 식생의 비선형적인 현상을 너무 단순화하여 모사하므로 발생하는 오차, 측정 자료의 오차 등도 있는 것으로 보인다(Zhao et al., 2005; Kimbell et al., 2017). MODIS 총일차생산량의 정확도를 보기 위해 Unprocessed MODIS 입력자료를 사용한 SVM모델 (SVM GPP)이나 에디공분산 타워의 총일차생산량(Eddy GPP)을 비교해 보았다. MODIS 총일차생산량은 비성장기와 성장기의 시작하는 지점에서 과소평가 하는 경향을 보였다. 또한 성장기의 최고점인 여름 동안에는 갑작스러운 떨어짐을 보였다(Fig. 3, left). MODIS 총일차생산량과 에디공분산 총일차생산량 피어슨 결정계수가 각각 0.74로 에디공분산 총일차생산량–SVM 총일차생산량 관계(R² = 0.94)보다 약한 상관관계를 보였다(Fig. 3, right). MODIS 총일차생산량은 SVM 모델의 결과 보다 에디공분산 총일차생산량과의 상관성이 낮았고 이는 기계학습 알고리즘 모델이 현장 예측에 사용되는 것이 적합하다고 판단할 수 있다.

Fig. 3. GPP from eddy covariance measurement (Eddy GPP), SVM model using unprocessed MODIS input (SVM GPP), and MODIS product for the test years, 2014-2015 (right panel). The relationship between MODIS GPP and EC GPP (open diamond, left panel) and SVM GPP (gray circles, left panel). The coefficients, R² of the regression are given.

2) SVM 모델의 확장 적용

(1) 단일 지역 훈련 모델

기계학습 알고리즘 SVM으로 구축된 총일차생산량 예측 모델의 공간적 확장 적용을 위해 두 가지 방법으로 모델을 구축하였다. 첫 번째 모델은 단일지역에서 훈련된 모델을 다른 지역으로 적용해 보는 것이고 두 번째 모델은 다중 지역의 자료들을 모아 훈련시키고 예측력을 테스트 해보는 것이었다. 먼저 단일 지점에서 훈련된 모델을 다른 지역으로 확장하여 적용하는 방법은 해당 지역들에 자료가 충분치 않거나 자료의 수집이 어려울 경우에 적용해 볼 수 있고 예측하고자 하는 지역범위가 매우 클 경우 모든 지역을 다 훈련 할 수 없을 때 적용 가능성을 볼 수 있다. 또한 지상 측정 자료의 부재나 측정 기간이 짧아 모델의 입력자료나 검증자료로 사용 하지 못하는 지역에 확장 적용이 가능한지 고찰해 볼 수 있도록 생태적 구성이 비슷하고 식물 생리 현상이 비슷한 곳에 적용이 가능한지를 고찰해 볼 수 있다.

에디공분산 플럭스 타워의 자료가 장기간(2006~현재) 측정되어 온 지점인 광릉 KoFlux 타워 활엽수림 지점(GDK)에서 구축된 SVM 총일차생산량 추정 모델을 다른 에디공분산 타워 지점들(GDK, SMC, HC, SC, WD, Jeju)에 적용해 보았다. 8일 단위 Unprocessed MODIS 입력자료로 2006년부터 2015년까지 훈련하고(n=460) 에디공분산 자료의 유무에 따라 테스트 연도는 지역별로 다르게 하여 총일차생산량을 예측하였다(Table 3). 모델의 총일차생산량과 에디공분산 총일차생산량 간의 상관관계는 피어슨 상관계수 0.75 – 0.95, RMSE 7.4 – 13.3 사이로 시계열 패턴은 완도(WD) 지역이 가장 맞지 않았고 홍천(HC) 지역이 가장 잘 모사하였다(Fig. 4). 삼척 (SC) 지역을 제외하고 SVM 모델의 총일차생산량은 계 절 변화 패턴은 잘 모사하지만 최고점 예측은 모두 과소평가 하는 것으로 나타났다. 광릉숲(GDK) 지역의 경우, 봄철과 초여름에 총일차생산량이 최고점을 보이는데 SVM 모델의 결과는 패턴은 모사하고 있으나 정량적으로는 과소평가하는 경향을 보인다. 이는 SVM 모델을 사용하였을 경우 지속적으로 나타나는 문제로 커널 기능을 사용하는 알고리즘에서 커널 파라미터와 커널 선택에 따라 모델의 파라미터가 결정되는데 이 과정에서 종종 모델의 과적합을 발생시키고 과적합으로 인해 모델의 결과는 과대/과소평가 된다(Cawley and Talbot, 2010). 이러한 결과는 유출량 예측(Granata et al., 2016)이나 수문 이벤트 예측(Langhammer and Cesak, 2016), 산림탄소예측(Dou et al., 2018) 등에서도 나타난다. 또한 지역 생태의 특성상 총일차생산량이 높다거나 에디공분산 플럭스 자료 자체의 과대평가 역시 배제할 수 없는 이유가 될 수 있다. 광릉숲(GDK)의 에디공분산 타워에서 측정된 8일 단위 최고 총일차생산량은 2017년 (85.2 gC m^-2 8-day^-1)을 제외하고 58.9 gC m^-2 8-day^-1~74.9gC m^-2 8-day^-1 사이에 분포한다. 2017년을 포함한 평균 69. 7 gC m^-2 8-day^-1로 2017년의 최고 총일차생산량은 다른 해에 비해 약 18% 높은 값을 보인다. 2016년 8일 단위 최고 총일차생산량 역시 74.3 gC m^-2 8-day^-1으로 평균보다 높은 값을 보인다. 이 두 해는 다른 해보다 갑자기 더 높은 최고 총일차생산량값을 보이는데, 광릉숲(GDK)지역의 2014년과 2015년에는 SVM 모델의 총일차생산량 예측력이 높은 결과(Fig. 3)와 비교해 갑작스러운 상승은 SVM모델이 정확하게 예측하지 못하는 것으로 보인다. Dou et al. (2018)의 연구에서도 노령림과 활엽수림, 혼효림에서 결과 값들이 영급이 낮거나 침엽수림인 곳에 비해 산발적인 관계를 보였다. 삼척(SC)의 경우 2000년 4월에 있었던 동해안 산불로 인해 파괴된 지역으로 생태계 복원이 이루어 지고 있는 곳이다. 국립산림과학원의 흉고직경 2 cm 이상의 진계생장목 개체수 조사에 의하면 2014부터 지속적으로 진계생장목의 생물량이 증가하고 있으며 흉고직경을 기준으로 2018년에는 급격한 증가가 발생하였다(국립산림과학원 내부자료). 이는 2016년에 비해 2017년에 높은 탄소고정량을 발생시킬 수 있는 근거로, 다른 해보다 높은 총일차생산량을 보이는 원인 중 하나일 가능성이 있음으로 판단된다. 그러나 SVM 모델은 이 이벤트를 모사하지 못하였는데 한 가지 이유는 훈련에 사용된 자료에는 산불이나 기타 다른 이유로 인해 생태계 교란이 발생하였음을 인지시켜주는 입력자료는 사용되지 않아서 인 것으로 판단된다.

Fig. 4. Relationship between eddy covariance GPP and SVM GPP utilizing data from unprocessed MODIS input for the test years indicated in Table 3. Pearson coefficient (R) is given. GDK = Gwangneung, SMC = Seolmacheon, HC = Hongcheon, SC = Samcheok, WD = Wando, and Jeju = Jeju.

Fig. 5. Seasonal variation of Gross Primary Production (GPP) obtained with eddy covariance measurement (black circle) and with Support Vector Machine (SVM) algorithm (gray circle). GDK = Gwangneung, SMC = Seolmacheon, HC = Hongcheon, SC = Samcheok, WD = Wando, and Jeju = Jeju.

(2) 다중 지역 훈련 모델

단일 지역의 자료로 훈련하였을 때 각 지역별로 특별히 발생하는 이벤트를 모사하지 못하는 한계점을 극복하기 위해 두 번째로 전체 지역의 자료를 훈련자료로 하여 모델을 구축하고 에디공분산 총일차생산량과 비교해 보았다(Fig. 6과 Fig. 7). 각 지역마다 목표 변수가 될 수 있는 에디공분산 플럭스 자료가 충분치 않아 광릉숲 (GDK) 지역을 제외하고 나머지 지역에서는 각각 1년씩만(2016년) 훈련자료로 사용하였다. 그 결과 피어슨 상관계수를 기준으로 설마천(SMC), 완도(WD), 제주 (Jeju) 지역에서는 모델의 결과가 개선되었고 삼척의 경우 약간 감소하였지만 RMSE나 평균제곱오차(MSE)를 기준으로 개선된 결과를 보였다(Table 3). 설마천 지역은 겨울 기간과 성장 시즌이 시작되는 시기에 단일 지역 훈련자료로 훈련된 SVM모델의 결과가 과소평가된 총일차생산량을 보였으나 전체 지역의 훈련자료로 구축된 SVM모델의 결과는 약간 개선된 것을 볼 수 있다. 다른 지역에 비해 총일차생산량 값이 높게 나타난 제주 (Jeju) 지역에서는 모델의 훈련 자료에 제주 지역의 훈련 자료가 포함된 이후 총일차생산량이 증가하는 결과를 보였다. 삼척(SC) 지역의 경우, 2016년의 매우 낮은 시계열 총일차생산량의 영향으로 2017년에 단일 지역모델의 결과보다 전체적으로 낮은 총일차생산량 값을 보였다. 완도(WD) 지역에서는 단일 지역 모델과 다중지역 모델 모두에서 가장 낮은 예측력을 보이는데 다른 지역들은 플럭스타워가 임분 내에 설치되어 있으나 완도(WD)의 경우 바다가 바로 인접해 있어 다른 지역의 생태 환경과 다르다. 삼척을 제외한 다른 지역보다 평균적으로 낮은 총일차생산량을 보이고 여름 기간 동안의 총일차생산량의 변화가 크다. 그러나 이러한 지역적인 특징이 반드시 총일차생산량의 경향과 밀접한 특징이 있는지는 향후 더 연구가 필요해 보인다. 또한 다른 지역과 다른 특징을 가진 자료라 해도 자료의 양이 모델의 예측력을 향상시키는 훈련에 충분히 확보가 된다면 다른 결과를 보일 수 있다고 판단된다. 총일차생산량의 최고점 예측은 단일 지역 모델과 마찬가지로 과소평가 하는 경향을 보였다. 따라서 이 결과만으로는 모델의 개선 여부를 판단하기에는 아직은 데이터의 양이 충분치 않고 특정 이벤트를 모사하기에 정보가 부족한 것으로 판단된다.

Fig. 6. Relationship between EC GPP and SVM GPP utilizing data from unprocessed MODIS input. SVM model trained with pooled data from all eddy covariance tower data. In this case, test years indicated in Table 3 were reduced one year except the Gwangneung site (GDK). Pearson coefficient (R) is given. GDK = Gwangneung, SMC = Seolmacheon, HC = Hongcheon, SC = Samcheok, WD = Wando, and Jeju = Jeju.

Fig. 7. Seasonal variation of Gross Primary Production (GPP) obtained with eddy covariance measurement (closed circle without outline) and with Support Vector Machine (SVM) algorithm (gray circle with outline), which trained utilizing data from the entire study sites, for the years indicated.

Table 3. Statistics results for the different trained models based on the input sets, which are from the single-site or multi-site. R, RMSE, STD and MSE denote correlation coefficient, root mean square error, standard deviation and mean squared error, respectively

그러나 산불과 산림 재생, 기상 조건의 변화나 기후변화 등으로 인한 갑작스러운 이벤트들은 이벤트가 발생하였을 때의 자료와 그 전후의 자료가 모두 훈련 자료에 입력되어야 모델의 예측력이 향상 될 것으로 판단된다. Langhammer and Cesak (2016)의 결과에서도 태풍이나 건조 같이 일반적이지 않은 수문 현상에 대해서는 SVM 알고리즘 모델의 예측력이 떨어지는 것을 보였다. 이것은 SVM 알고리즘 모델이 아닌 다른 기계학습 알고리즘을 사용하더라고 반드시 해결해야 할 과제이다. 또한 기계학습의 특성상 모델을 훈련시킬 때 사용된 기준자료의 불확실성 역시 모델의 예측력에 영향을 준다. 이 연구에서 기준 자료로 사용된 에디공분산 총일차생산성의 추정 과정에서 생태계 호흡량은 야간 온도로 추정하고 있어 야간 시간의 생태계 호흡량이 과소평가 되는 경향이 있다는 것은 여러 연구를 통해 제기되었다(Lavibne et al., 1997; Baldocchi, 2003; Aubinet, 2008; Speckman et al., 2014). 향후 입력자료의 개선과 딥러닝(Deep learning)의 적용 시 보다 개선된 결과가 나올 것으로 기대된다.

4. 결론

이 연구는 기계학습 알고리즘과 인공위성자료를 이용하여 산림의 총일차생산량(GPP)의 추정 가능성과 산림생태계에서의 위성자료 활용 가능성에 대해 논의해보고자 하였다. 이 연구는 농림위성 알고리즘 개발을 위한 기초자료로서의 가치와 위성자료와 기계학습 알고리즘 간의 활용 가능성을 보여주고 산림일차생산성을 추정하는 방법에 대한 보다 접근이 쉬운 방안을 제시하였다. 이 연구를 통해 우리나라의 복잡한 산림 지형에서도 SVM 알고리즘이 적합한 알고리즘 중 하나로 사용될 수 있다는 것을 증명하였다. SVM 알고리즘을 사용해 단일 지역에서 구축된 모델은 다중 지역으로 적용하였을 때 비교적 높은 예측력을 보이나 각 지역에서 발생하는 예상치 못한 이벤트에 대한 예측력은 아직 부족한 것으로 나온다. 다중 지역의 자료를 모두 활용하여 모델을 훈련시키고 그 모델로 총일차생산성을 예측하였을 때 약간의 개선을 보였는데 이는 모델의 훈련자료가 다양해 질수록 예측력이 향상 될 수 있다고 보인다. 장기간 자료가 있는 단일 지역 훈련자료 기반의 모델과 모든 지역의 자료가 훈련자료로 들어간 모델의 예측력은 약간의 차이가 있기는 하지만 매우 향상된 결과를 나타내지는 않는데 이는 각각의 이벤트들도 중요하지만 장기간의 자료가 모델의 예측력에 큰 영향을 미치는 것을 알 수 있다.

위성영상자료가 빅데이터로 분류될 만큼 자료의 양이 많아지고 그 속도가 빨라지고 정확도도 높아지고 있는 가운데 산림 생태계에서의 위성자료의 활용은 인력이나 비용의 한계로 조사되지 못한 전국의 산림에 대한 조사가 가능하게 되었다. 또한 전국 단위의 조사나 북한까지 조사가 가능하게 하여 보다 넓은 관점에서의 산림생태계 이해가 가능하고 이를 산림 관리나 기후변화에 대응한 산림 보전 및 적응 대책을 마련하는데 높은 기여를 할 수 있을 것으로 판단된다. 현재 우리나라에서도 자체 농림위성 알고리즘 개발을 위한 연구들을 집중적으로 추진하고 있으며 이 연구의 결과 또한 우리나라의 농림위성 알고리즘 개발의 자료로 활용될 것으로 기대된다.