A Study on the Calculation of Evapotranspiration Crop Coefficient in the Cheongmi-cheon Paddy Field

청미천 논지에서의 증발산량 작물계수 산정에 관한 연구

Abstract

In this study, crop coefficients were calculated in two different methods and the results were evaluated. In the first method, appropriateness of GLDAS-based evapotranspiration was evaluated by comparing it with observed data of Cheongmi-cheon (CMC) Flux tower. Then, crop coefficient was calculated by dividing actual evapotranspiration with potential evapotranspiration that derived from GLDAS. In the second method, crop coefficient was determined by using MLR (Multiple Linear Regression) analysis with vegetation index (NDVI, EVI, LAI and SAVI) derived from MODIS and in-situ soil moisture data observed in CMC, In comparison of two crop coefficients over the entire period, for each crop coefficient GLDAS K_c and SM&VI K_c, shows the mean value of 0.412 and 0.378, the bias of 0.031 and -0.004, the RMSE of 0.092 and 0.069, and the Index of Agree (IOA) of 0.944 and 0.958. Overall, both methods showed similar patterns with observed evapotranspiration, but the SM&VI-based method showed better results. One step further, the statistical evaluation of GLDAS K_c and SM&VI K_c in specific period was performed according to the growth phase of the crop. The result shows that GLDAS K_c was better in the early and mid-phase of the crop growth, and SM&VI K_c was better in the latter phase. This result seems to be because of reduced accuracy of MODIS sensors due to yellow dust in spring and rain clouds in summer. If the observational accuracy of the MODIS sensor is improved in subsequent study, the accuracy of the SM&VI-based method will also be improved and this method will be applicable in determining the crop coefficient of unmeasured basin or predicting the crop coefficient of a certain area.

본 연구에서는 두 가지 방법으로 작물계수를 산정하고, 그 결과를 평가하였다. 첫 번째 방법에서는 GLDAS 자료를 청미천 플럭스타워의 증발산량 실측값과 비교하여 적정성을 평가한 뒤 GLDAS 기반 실제증발산량을 잠재증발산량으로 나눠 작물계수(GLDAS K_c)를 산정하였으며, 두 번째 방법에서는 MODIS기반 식생지수(NDVI, EVI, LAI, SAVI)와 플럭스타워에서의 토양수분 실측치를 이용해 다중선형회귀분석으로 작물 계수(SM&VI K_c)를 산정하였다. 전체기간에 대한 두 가지 작물계수(GLDAS K_c, SM&VI K_c)를 통계(mean, bias, RMSE, IOA)를 통해 비교해 본 결과 평균값은 각각 0.412와 0.378, bias는 0.031과 -0.004, RMSE는 0.092와 0.069, 적합도 지표(IOA)는 0.944와 0.958로 두 방식 모두 전반적으로 실측값과 유사한 패턴을 보여주었다. 그라나 SM&VI 회귀모형 방식이 더 우수한 것으로 나타났다. 또한, 벼의 생장 단계별로 GLDAS K_c와 SM&VI K_c에 대한 통계적 평가를 수행해본 결과 초기와 중기에는 GLDAS 기반의 K_c가 더 우수했으며, 후기에는 SM&VI 기반의 K_c가 더 우수한 것으로 나타났다. 이는 봄철에는 황사, 여름철에는 비구름으로 MODIS 센서의 정확성이 감소했기 때문인 것으로 판단된다. 향후 연구를 통해 MODIS 센서의 관측 정확성이 향상된다면, SM&VI 기반 작물계수 산정방식의 정확성 역시 향상될 것으로 판단되며, 미계측 유역의 작물계수 산정이나 작물계수의 예측에 사용될 수 있을 것으로 판단된다.

1. 서론

증발산은 증발과 증산을 통틀어 이르는 말로서 수문 순환을 구성하는 주요한 현상 중의 하나로서, 육지에서 일어난 강수의 약 64%를 대기로 환원시키는 역할을 한다(Baumgartner and Reichel, 1975). 증발산은 지표와 대기 사이의 에너지 및 수자원의 교환과 직접적인 관련이 있기 때문에 증발산량은 지표상의 수자원의 손실량을 파악하는데 필수적으로 고려되어야 하는 인자라고 할 수 있다. 정량화된 증발산량 자료는 수자원 관리, 가뭄 예측, 농업 계획 등을 위한 주요한 지표로 사용될 수 있다.

현재 증발산량은 다양한 방식으로 계측되고 있는데, 그 중 대표적으로 에디공분산(eddy covariance) 측정 방법, 라이시미터(lysimeter) 측정 방법 등의 직접측정방법과 Food and Agriculture Organization에서 제시한 PenmanMonteith 증발산 산정식(FAO-PM식), Priestley-Taylor 공식을 통해 기상자료를 바탕으로 증발산량을 산정하는 간접측정방법이 대표적이다. 직접 관측하는 방법 중 가장 대표적으로 쓰이는 에디공분산 방법은 플럭스를 관측을 핵심적으로 측정하는데, 이는 플럭스 타워에서 관측한 대기 중의 난류 플럭스의 수증기 농도와 연직 풍속 사이의 공분산을 계산함으로써 실제 증발산량을 관측하는 방식이다. 최근에는 이러한 플럭스타워를 이용해 증발산량을 비롯하여 이산화탄소, 수증기, 에너지 교환 등을 관측하는 FLUXNET(AsiaFlux,Ameri-Flux,EuroFlux, Oz-Flux, Canada-Flux)과 같은 증발산량 관측 네트워크가 전 세계적으로 구축되고 있으며(Baldocchi et al., 2010), 국내에서도 물, 탄소, 지표와 대기 간의 에너지 교환 및 순환에 대한 정확한 해석의 필요성이 대두됨에 따라 2001년부터 KoFlux가 운영되고 있다(Kwon and Kim, 2010). 에디공분산 방법의 가장 큰 장점은 직접적인 관측이 가능하다는 것이지만 고가의 장비가 필요하며, 민감한 측정을 하는 센서에 물 맺힘 등 오염이 발생하면 오차가 크게 발생한다는 단점이 있다.

위와 같은 이유로 인해 증발산량을 간접 측정하는 방법을 이용하는데 FAO-PM 방법, Priestley-Taylor 방법과 같은 증발산량을 간접 측정 방법이 이용하기도 한다. 이러한 방법들은 기상요소 및 주변 환경에 대한 관측 값을 이용해 물리식 기반으로 표준 조건 증발산량을 계산하는 방법이다. FAO-PM 방법은 전 세계적으로 에너지 수지와 공기동역학적을 고려한 항을 포함하는 식을 기반으로 표준 조건 증발산량을 산정하며(Allen et al., 1994), Priestley-Taylor 방법은 Penman-Monteith 방법에서 공기동역학적 항을 간략화 시킨 방법으로써 표준 조건 증발산량을 산정한다. 이러한 방식들은 통합기상장비를 활용하여 비교적 저렴하게 측정이 가능하나, 표준 조건의 증발산량을 산정하기 때문에 주변 환경을 반영하는 작물계수와 별도의 계수들을 곱해야 한다. 세분하자면, 증발산량은 표준 조건 증발산량(ET₀)과 실제증발 산량(ET_a)으로 구분할 수 있는데, 표준 조건 증발산량은 기준 작물에 수분 제한 없이 공급되는 조건에서의 증발산량을 의미하고, 실제증발산은 물과 실질적인 수분공급, 주변 환경 등을 모두 고려한 증발산량을 의미 한다. 이는 실제증발산량은 농업에서 관개를 정확하게 계획하고 실행하는데 에는 실제 농작물의 증발산량이 이용되며, 작물계수를 이용해 이를 산정하는 방식이 널리 이용되고 있다(Doorenbos and Pruitt, 1977; Jensen and Allen, 2000).

작물계수를 산정하는 방법으로는 주로 실제증발산에서 표준 조건 증발산량으로 나누어 직접적으로 산출하는 방법이 있다. 두 증발산량 값에 오차가 포함되면 작물계수의 값이 부정확해지는 치명적인 단점이 있다. 이를 극복하기 위해 여러 인자들과 경험식들을 이용하여 작물계수를 산정하는 방법이 연구가 되어왔다(Allen et al., 2005; Ko et al., 2009). 최근 인공위성 등의 원격탐사 기반의 활용도가 높아짐에 따라 이용할 수 있는 여러 인자들이 늘어나 작물계수를 산정하는 방법에도 다양하게 활용되고 있다. 특히 Park et al. (2017)은 인공위성 기반으로 측정한 토양수분과 식생지수를 통해 동아시아에서 작물계수를 평가하였으며, Kamble et al. (2013)에서 정규식생지수(Normalized Difference Vegetation Index, NDVI)를 이용하여 작물계수를 산정하고자 하였다.

이와 같이 많은 작물계수 산정 방법이 존재하고 있으나 이를 정량화 하거나 산정된 식에 대한 검증에 대한 연구들이 부족한 실정이다. 본 연구에서는 최근 자료동화(assimilation)를 활용한 지형모델 자료의 정확성이 높아짐에 따라 관측의 타당성과 검증 연구에 높은 연구 자료로써 활용되고 있는 자료와 인공위성을 활용한 식생 지수 및 토양수분량을 활용하여 논지에서의 작물계수 산정방법에 대해 분석을 수행하였다.

2. 연구 지역

연구지역으로는 한국수자원조사기술원(Korea Institute of Hydrological Survey, KIHS)에서 관리하는 청미천 플럭스 타워 플럭스 타워를 선정하였다(Fig. 1). 청미천 플럭스 타워(CMC)는 경기도 여주시 점동면 농업기술원 종자관리소의 인근에 위치하고 있으며, 그 주변의 약 41%는 임야지, 약 49%는 농경지로 구성된 농경지 및 산지형 하천의 모습을 보인다. 임야는 침엽수와 활엽수가 주를 이루는 혼합림으로 구성되어 있으며, 토양은 양토 및 사양토가 주를 이루고 있다(Byun et al., 2013). 연평균 기온은 11.6°C, 연평균 강수량은 1,220 mm로서, 한반도의 일반적인 기후 및 기상과 유사한 성향을 보인다(Kim et al., 2016). CMC에는 에디공분산 기반 증발산량의 관측을 위해 순복사계(CNR2, Campbell Scientific), 3차원 초음파 풍향풍속계(CSAT3, Campbell Scientific), 폐회로 기체분석기(EC-155, Campbell Scientific) 등이 설치되어 있다(Park et al., 2015).

Fig. 1. Location of Study Area in CMC.

3. 연구방법 및 자료

1) 자료

(1) 증발산량 자료

① 에디공분산 자료

농경지(논)에서 증발산량을 관측하기 위해 한국수자원조사기술원에서는 청미천 유역에 에디공분산 장비를 설치·운영하고 있다. 에디공분산 장비는 공기중의 난류 플럭스에서 수증기 농도와 연직 풍속 사이의 공분산 값으로 증발산량을 산출한다. 에디공분산 장비를 통에 증발산량을 산정하기 위해서는 공기중의 난류 플럭스에서 수증기 농도와 연직 풍속사이의 공분산을 이용하여 산출을 한다. 에디공분산을 활용한 증발산량 산정은 주요 가정이 만족되어야 하는데, 이를 위해선 플럭스 관측 타워의 주변에 경사가 거의 없고, 식생이 균질하며, 난류장의 형성을 방해하는 방해물이 관측소 주변에 없어야 한다. 하지만, 이러한 가정이 모두 만족되는 이상적인 환경은 흔치 않기 때문에 정확한 증발산량을 산정하기 위해서는 관측결과에 대한 품질관리가 반드시 이루어져야 한다. 그러한 이유로 한국수자원조사기술원에서는 KoFlux 표준화 프로그램을 통해 좌표변환, 밀도보정, spike 제거 등의 품질관리 프로세스를 수행하고 있으며, 좌표변환에는 평면 맞추기 회전(Planar Fit Rotation), 밀도보정에는 Webb-Pearman-Leuning(WPL) 방법을 사용하여 현장의 관측 자료가 에디공분산의 기본가정을 만족시킬 수 있도록 하고 있다(Hong et al., 2010). 해당 프로그램을 통해 불량 자료로 판별된 관측 구간은 제거되어 결측 구간(missing data)으로 남게 되는데, 한국수자원조사기술원에서는 FAO-PM(Food and Agriculture Organization Penman-Monteith) 방법, 평균 일 변동(Mean Diurnal Variation, MDV) 방법, 칼만 필터 (Kalman filter) 등을 이용해 이러한 결측 구간을 보완하고 있다. 결측 구간을 메우는 방법은 Kim et al. (2016)에 자세히 기술되어있다.

Fig. 2. The methodology of gapfilling flowchart (Kim et al., 2016).

② Priestley-Taylor

잠재증발산량(ETp)을 추정하는 Penman(1948) 공식을 기반으로 하여 공기동력학적 항이 에너지 항의 약 30%의 비율을 가진다는 경험 하에 입력 자료를 간소화하여 다음과 같은 식 (1)으로 표현하였다.

\(\lambda E=\alpha_{\varepsilon} \frac{\Delta}{(\Delta+\gamma)}\left(R_{n}-G\right)\)       (1)

여기서, λE는 잠열 플럭스, R_n은 순복사량(W/m² ), G은 토양열플럭스(W/m²), Δ는 포화증기압곡선(kPa °C^-1), γ는 건습계 상수(kPa °C^-1)이다. 또한 α_e는 Priestley-Taylor 상수로서 주로 습윤한 지역에서 통상적으로 1.26, 건조한 지역에서는 1.75 정도의 값을 주로 이용하고 있으며, 본 연구지역은 농경지로서 항시 습윤한 상태를 유지하고 있기 때문에 상수 1.26을 이용하여 계산하였다.

위에서 산정한 Priestley-Taylor 공식을 통한 잠재증발 산량은 식 (2)과 같이 토양의 토양수분 응력(soil water stress, K_s)와 작물계수(crop coefficient, K_c)를 곱하여 실제증발산량을 산정할 수 있다.

\(E T_{a}=K_{s} \times K_{c} \times K_{p}\)       (2)

즉, 본 연구지역에서는 청미천 논지임으로 토양수분응력은 1과 가깝기 때문에 식 (3)와 같이 표현 할수있으며 본 연구 취지에 맞게 FAO-PM에서 산출된 기준 작물 증발산량을 활용하여 작물계수를 산정할 수 있다.

\(E T_{a}=K_{c} \times E T_{p}\)       (3)

③ Global Land Data Assimilation System(GLDAS)

전 지구 지표동화시스템(Global Land DataAssimilation System, GLDAS)은 Mosaic, Noah, Community Land Model (CLM)의 3가지의 지표모델(Land surface models)을 기반으로 한 자료동화 데이터 시스템이며, 다양한 기상 자료(GDAS, GEOS, ECMWF, TRMM, PERSIANN 등)와 지표면 정보를 입력 자료로 이용한다. GLDAS는 최근 위성 데이터와 모델 데이터를 검증하는 연구에서 활발히 사용되고 있는 신뢰도 높은 자료이다. 특히 결측이 빈번히 발생하는 곳과 미계측 지역 등에 주로 유용하게 쓰이는 자료이다. 본 연구에서는 GLDAS 2.1 NOAH 모형을 사용하여 제공되는 실제 증발산을 추출하였으며, 또한 GLDAS에서 제공되는 다양한 기상인자들을 활용하여 Priestly-Taylor 방식을 활용하여 잠재증발산을 산출하였다. 약 25 km의 공간해상도와 3시간의 시간해상도를 가진 자료를 활용하였으나 MODIS 데이터와 비교를 위해 시간해상도를 16일 단위로 자료를 취합하였다.

2) 식생 지수

(1) NDVI(Normalized Difference Vegetation Index

NDVI는 식생의 유무와 상태에 관한 정보를 나타내는 지표로서, 원격탐사 기법에서 식물의 상대적 분포, 활동성, 엽면적지수(LAI) 등과 관련하여 널리 쓰이고 있다. NDVI는 위성영상을 통해 얻은 근적외선 파장대의 반사도(ρ_NIR)와 적색 파장대의 반사도(ρ_red)를 이용해 계산할 수 있는데, 이는 일반적인 건강한 식물이 지닌 클로로필(chlorophyll)이 가시광선의 적색영역 파장을 많이 흡수하고, 근적외선 파장은 대부분을 반사하는 특성을 이용한 것이다. NDVI 식생지수는 아래의 식 (4)과 같이 계산할 수 있다.

\(N D V I=\frac{\rho_{N I R}-\rho_{red}}{\rho_{N I R}+\rho_{red}}\)       (4)

이렇게 계산된 NDVI는 일반적으로는 -1에서 +1사이의 값을 가지며, 그 값이 음수인 경우는 보통 식생이 없거나 적은 수역, 도시지역으로 판단할 수 있으며, NDVI 값이 양수인 경우는 농지, 초지, 산림 등 식생피복 지역으로 판단할 수 있다. 또한, NDVI 값은 식물의 클로로필 농도와 연관이 있기 때문에 병충해 식생이나 노화식생이 많은 지역은 그렇지 않은 지역에 비해 낮은 NDVI값이 산출되며, NDVI값이 +1에 가까울수록 식물의 양이 많고 건강한 지역임을 의미한다. 본 연구에 서는 TERRA위성의 MODIS 센서에서 산출되는 16일 주기의 NDVI 값을 사용하였다.

(2) EVI(Enhanced Vegetation Index)

기존의 식생지수들은 식생의 밀도가 높은 지역에서는 비선형적인 관계가 나타나 식생지수의 산출이 불안정한 측면이 있었으며, 대기 산란과 감쇠효과에 민감한 특성이 있었다. EVI는 캐노피로 인한 노이즈와 대기 효과를 모두 고려하여 그러한 문제점을 보완한 식생지수로서, MODIS 센서를 사용하는 Terra 위성과 Aqua 위성 의식생에 관한 표준 산출물로서 개발되었다(Jiang et al., 2008). EVI는 아래의 식 (5)를 통해 산출된다.

\(E V I=G \frac{\rho_{\mathrm{NIR}}-\rho_{\mathrm{red}}}{\rho_{\mathrm{NIR}}+C_{1} \times \rho_{\mathrm{red}}-C_{2} \times \rho_{\mathrm{blue}}+L}\)     (5)

여기서 G는 gain 계수, C₁과 C₂는 에어로졸 저항계수, ρ_blue는 청색 파장대의 반사도이며, L은 토양보정계수와 에어로졸 저항에 관한 함수로 표현되는 캐노피 보정계 수이다. MODIS를 통한 EVI의 산출에는 일반적으로는 L=1, C₁=6, C₂=7.5, G=2.5의 값을 사용한다(Jiang et al., 2008). 본 연구에서는 TERRA위성의 MODIS 센서에서 산출되는 16일 주기의 EVI값을 사용하였다.

(3) LAI(Leaf Area Index)

LAI는 식생들의 잎 면적(한쪽)의 총합을 단위면적으로 나눈 값으로서, 캐노피의 양이나 식물의 생산성들을 특징짓는 지수이다. LAI의 추정에는 다양한 방법이 존재하지만, 합리적인 비용으로 신속하고 객관적인 결과물을 얻을 수 있다는 측면에서 인공위성을 이용해 LAI를 산출하는 방식이 가장 널리 이용되고 있다. 인공위성 기반 LAI 산출은 Terra/Aqua, LANDSAT, GLASS, ERS-1, Sentinel-2 등 다양한 위성을 통해 이뤄지고 있는데, 그 중 Terra/Aqua를 이용하는 MODIS-LAI가 가장 널리 쓰이며, 본 연구에서도 이를 활용하였다.

(4) SAVI(Soil-Adjusted Vegetation Index)

식생 밀도가 높은 지역에서는 토양의 색상, 수분량, 포화 정도에 따라 달라지는 토양밝기에 따라 NDVI 값이 부정확하게 산출되는 문제가 있는데, 이를 보완하기 위해 고안된 식생지수가 SAVI이다(Huete, 1988). SAVI 는 기존의 NDVI 산출식에 토양보정계수 L을 추가해 아래의 식 (6)와 같이 계산함으로써 토양 밝기에 의한 영향을 최소화한 식생지수이다.

\(S A V I=\frac{\rho_{NIR}-\rho_{red}}{\rho_{NIR}+\rho_{red}+L}(1+L)\)       (6)

일반적으로 L=0.5 일 때 토양 밝기의 영향이 최소가 되기 때문에 캐노피 보정계수를 0.5로 고정하여 SAVI 를 산출하는 방식이 널리 쓰이지만, 이는 시간에 따라 변화하는 식생을 잘 반영하지 못하기 때문에 회기분석을 통해 구한 최적의 L값이나 함수로 표현되는 L이 SAVI의 산출에 쓰이기도 한다(Qi et al., 1994). 본 연구에서는 TERRA 위성의 MODIS 센서에서 산출되는 16일 주기의 밴드(band) 값을 사용하였다.

3) 토양 측정

청미천 관측소에는 Flux tower 아래 토양수분(CS616), 온도(TCAV) 및 토양열 플럭스(HFT)를 측정하는 센서가 설치되어있다. 토양 측정센서들은 모내기철 및 추수기간에 일시적으로 관측을 중지하여 센서의 손상을 최소화 하고 있다. 토양은 논지 특성상 점토질로 이루어져있으며, 여름철 관개용수가 지속적으로 공급되어 포화 토양수분상태를 유지하고 있다. 포화 토양수분 량은 보통 0.65 m³ m^-3 내외이고 강우사상에 대해서도 유효한 범위 내에서 반응하였다.

4) Crop coefficient(Kc)

식량농업기구(Food and Agriculture Organization, FAO)에서 제안하는 Penman-monteith 식은 표준조건에서의 증발산량(ET₀)을 산정 하고 있으며, 실제 증발량 (ET_a)을 산정하기 위해서는 작물계수(K_c)를 이용해야 한다(식 7, 8). 이렇게 산정된 식을 역산해 보면 ET₀ 대비 ET_a을 통해 K_c를 산정할 수 있음을 내포하고 있다.

\(E T_{0}=\frac{0.408\left(R_{n}-G\right)+\gamma \frac{37}{T+273} u_{2}\left(e_{s}-e_{a}\right)}{\Delta+\gamma\left(1+0.34 \mathrm{u}_{-} 2\right)}\)       (7)

\(E T_{c}=K_{c} \times E T_{0}\)       (8)

여기서, ET₀는 기준 증발산량(mm h^-1), ETa는 실제 증 발산(mm h^-1), K_c)는 작물계수를 나타낸다. T는 평균온도 (°C), u₂는 지상으로부터 2 m에서 측정된 평균풍속(m/s), e_s는 포화수증기압(kPa), e_a는 실제 수증기압(kPa)으로 (e_s – e_a)는 포차(kPa)를 나타낸다. 작물계수는 작물의 증산(K_cb)과 토양의 증발(K_e) 영향을 통합한 작물계수로서 식물과 토양의 상태에 대한 특성이 중요한 요소로 작용한다(식 9).

\(K_{c}=K_{c b} \times K_{e}\)       (9)

본 연구에서는 논지(rice paddy)에서 이루어지므로 FAO Irrigation and Drainage paper No. 56에서 제안한 작물계수범위(0∼1.2)를 고려하여 산정하였다. 또한 Penman 식에서 파생된 Priestly-Taylor 방식으로 얻은 잠재증발산을 3.1.1.2 절에서 설명한 바와 같이 표준조건 증발산량으로 설정하고 각각 에디공분산 방법과 GLDAS에서 산정된 실제증발산과의 비율을 통해 작물 계수를 산정하였다. 두 번째 토양 수분과 식생지수를 이용하여 다중선형회귀모형을 통해 K_c를 산정하였으며, 이는 3.5절에 자세히 기술되어있다.

5) Multiple Linear Regression(MLR)

다중선형회귀모형은 2개 이상의 독립변수(independent variable)를 가지고 종속변수(dependent variable)를 예측하기 위한 선형회귀모델이며 식 (10)과 같이 표현할 수 있다.

\(\left[\begin{array}{c} y_{1} \\ y_{2} \\ \vdots \\ y_{n} \end{array}\right]=\left[\begin{array}{cccc} 1 & x_{11} & \cdots & x_{p 1} \\ 1 & x_{12} & \cdots & x_{p 2} \\ \vdots & \vdots & \ddots & \vdots \\ 1 & x_{1 n} & \cdots & x_{m} \end{array}\right]\left[\begin{array}{c} \beta_{1} \\ \beta_{2} \\ \vdots \\ \beta_{n} \end{array}\right]+\left[\begin{array}{c} \epsilon_{1} \\ \epsilon_{2} \\ \vdots \\ \epsilon_{n} \end{array}\right]\)         (10)

여기서, y₁는 종속변수를 나타내며, p개(p≥2)의 독립 변수 x를 통해 모집단의 회귀계수(regression coefficients, β)를 추정하며, є은 추정할 때 발생하는 오차이다. 추정량을 구하기 위해서는 잔차의 제곱 합을 최소로 하는 최소자승법(least squares method)을 이용하였다. 본 연구 에서는 Kc값을 산정하기 위해 종속변수인 y에 설정하였으며, K_c값에 큰 영향을 미치는 토양수분 값과 NDVI, SAVI, EVI, LAI를 독립변수 x로 설정하여 다중선형회 귀모형을 수행하였다.

4. 결과

1) 작물계수 산정

본 절에서는 GLDAS에서 제공하는 자료와 토양수분 및 식생지수 자료(SM&VI)를 활용하여 작물계수를 산정하였다. 우선 GLDAS 자료의 적정성을 비교하기 위해 실제증발산과 잠재증발산을 산정한 값은 Flux tower 에서 실측한 값과 비교평가를 하였다(Fig. 3).

Fig. 3. (a) Actual & potential evapotranspiration between Flux tower and GLDAS, (b) Time series of actual evapotranspiration in 2016, 2017.

Fig. 3(a)를 확인해보면 실제증발산량과 잠재증발산량은 GLDAS 값이 전체적으로 Flux tower 값보다 비교적 과대산정(overestimation)이 되었으며, 특히 잠재증발 산량 부분에서는 높은 일 증발산량이 발생할수록 과대산정이 심화되었다. 이러한 GLDAS의 과대산정은 선행 연구에서도 확인된 바 있다(Baik et al., 2019). 또한 일 실제증발산량 최댓값은 Flux tower에서 5.88 mm, GLDAS는 6.26 mm으로 나타났으며, 잠재증발산량은 각각 7.96 mm, 10.86 mm으로 나타났다. 이는 여름철 국지성 호우가 빈번히 발생하여 나타나는 것으로 보인다. 또한 Flux tower에서는 10 Hz 측정을 하고 있으나 GLDAS는 인공위성 등 3시간 이상의 시간해상도를 갖는 입력 자료를 이용하기 때문에 자료의 평활화(smoothing) 현상으로 볼 수 있다. 또한 GLDAS의 특성상 공간해상도가 다르기 때문에 발생하는 오차로도 볼 수 있다. Fig. 3(b)는 GLDAS와 Flux tower의 시계열 경향을 나타낸 그림으로서, GLDAS에서 산출한 실제증발산 시계열 경향은 Flux tower의 실측값과 유사한 범위와 경향을 나타내고 있다. 이는 작물의 생장의 생장곡선과 유사하며 두 작물계수를 산정하기에 적절할 것으로 보인다. 이렇게 평가된 결과를 바탕으로 GLDAS를 통해 산정한 실제증발산과 잠재증발산을 나누어 작물계수로 나타내었다.

SM&VI 회귀식 방법에서는 토양수분과 식생지수를 입력 자료로 하여 다중선형회귀모형을 이용하여 작물 계수를 산정하였다. 이는 식 (9)과 같이 주어지는 작물과 토양의 수분 조건에 따라 정해지는 단순 증가조건으로서 다중선형회귀모형을 사용하여 산출함에 적합한 것으로 판단하였다.

Fig. 4에서는 식생지수의 시계열 그래프로서 여름철 식생의 증가와 가을철 추수로 인한 식생의 감소가 발생하는 증감패턴을 확인할 수 있다. 본 연구에서 사용되는 식생지수는 Aqua/Terra에 탑재되어있는 분광반사계 (spectroradiometer)에서 관측된 밴드(band) 값을 이용하고 있으나 대기의 상태와 주변 환경에 따라 영향을 받기 때문에 잦은 호우사상이 발생하는 기간에서 식생지수 값에 영향을 미친다. 따라서 NDVI 뿐만 아니라 대기 및 토양 보정된 값을 모두 이용하여 신뢰도를 높였으며, 다양한 식생지수와 토양수분을 활용하여 다중선형회 귀모형을 통해 산출한 결과 식 (11)과 같이 도출되었다. 상관계수와 결정계수는 각각 0.91과 0.82로 높은 근사치를 나타내었으며 토양수분과 식생지수를 다중선형회 귀모형을 활용하여 K_c에 적용 가능성을 확인하였다.

\(K_{c}=\begin{array}{l} 0.063 \times N D V I+0.001 \times S M-0.31 \times E V I \\ +0.003 \times L A I+0.386 \times S A V I+0.18 \end{array}\)       (11)

Fig. 4. Time series of soil moisture and various vegetation index.

2) 작물계수 평가

위의 절에서 산정된 자료는 Flux tower에서 실측한 자료를 기준으로 평가하였다. Fig. 5는 Flux tower, GLDAS 및 SM&VI 회귀식을 활용하여 산정한 K_c를 산정하여 시계열로 나타낸 그림이다. 자료간의 비교를 위하여 식생 지수의 관측시간 범위인 16일에 맞추어 Flux tower와 GLDAS로 계산한 K_c는 16일 단위로 산술평균 계산하여 비교하였다. 전반적으로 시계열 패턴은 Flux tower와 유사함을 보였으며, 모내기가 시작되는 4월부터 추수기간인 9월까지 증감 패턴이 잘 나타났다. 이러한 결과는 FAO Irrigation and Drainage Paper No.56(Allen et al., 1998)에서 제시한 작물계수의 증감 패턴과도 유사한 패턴을 보여준다.

Fig. 5. Time series of crop coefficient

Flux tower는 전체기간에서 평균값이 0.382로 나타났으며, GLDAS는 0.412, SM&VI 회귀식은 약간 낮은 0.378의 값을 가졌으며, Bias는 GLDAS와 SM&VI 회귀 식에서 각각 0.031, -0.004를 보였다. 평균 제곱근 오차 (Root Mean Square Error, RMSE)는 GLDAS와 SM&VI 회귀식에서 각각 0.092, 0.069로 SM&VI 회귀식이 더 낮은 값을 보였다. 적합도 지표(Index Of Agree, IOA)에서도 GLDAS가 0.944, SM&VI 회귀식이 0.958로 두 데이터 모두 높은 IOA를 가졌으며 그 중 SM&VI 회귀식으로 산정한 K_c가 더 높은 값을 나타냈다(Table 1). 하지만 이러한 결과는 GLDAS로 산정하는 경우와 SM&VI 회귀식으로 산정하는 경우 결론적으로 많은 차이를 나타내지 않아 미계측 지역에서의 활용성을 나타낼 수 있을 것으로 보인다.

Table 1. Statistical comparison of Crop Coefficient (Kc) using GLDAS and SM&VI Regression

Fig. 6는 GLDAS와 SM&VI 회귀식을 이용한 K_c를 각각 벼의 생장을 단계별로 나누어 통계치를 Taylor diagram으로 도식하여 평가하였다. 1월부터 모내기가 시작되는 4월 중순까지를 초기(initial)단계, 4월 중순부터 벼가 성장하고 9월 중순에 추수하기 이전까지를 중간(mid)단계, 9월 중순부터 12월까지 마지막(late)단계로 설정하여 분석을 수행하였다. 기간별 통계치 살펴본 결과 GLDAS를 이용한 Kc는 초기(상관계수: 0.730, 정규화된 표준편차: 0.994), 중간(상관계수: 0.722, 정규화된 표준편차: 1.089), 마지막(상관계수: 0.651, 정규화된 표준편차: 0.910)순으로 높게 나왔으며, SM&VI 회귀식을 이용한 K_c는 마지막(상관계수: 0.775, 정규화된 표준편차: 1.076), 중간(상관계수: 0.683, 정규화된 표준편차: 0.965), 초기(상관계수: 0.583, 정규화된 표준편차: 1.000)순으로 높은 상관계수를 보였다. 이러한 결과는 GLDAS를 활용한 K_c가 K_c의 패턴을 나타내는데 더 적합하며, SM&VI 회기식은 정량적인 값을 잘 표현할 수 있는 것으로 해석할 수 있다. 특히, 초기와 중간 단계에서 전반적으로 SM&VI 회귀식 보다 GLDAS가 논지의 관개와 식물의 생장에 대한 반영을 하는 것으로 나타냈는데, 이는 SM&VI 회귀식이 MODIS 센서의 한계로 인해 봄철의 황사와 미세먼지, 여름철 장마와 잦은 비구름으로 식생 지수의 정확도가 떨어져 SM&VI 회기식의 정확도가 떨어진 것에 의한 것으로 추측된다. 관측에 정확도가 떨어지는 현상의 결과라고 추측할 수 있다. 따라서 식생 지수를 관측할 수 있는 센서들의 정확성이 확보(결측 자료나 악성자료들이 발생하는 한계를 극복)하는 연구들이 수행된다면 작물계수 정확성의 상승을 기대할 수 있을 것으로 보인다(Jung et al., 2013).

Fig. 6. The statistics of the crop coefficient (Kc) both Flux tower and SM&VI Regression, GLDAS using Taylor diagram.

5. 결론

본 연구의 목적은 미계측 유역에서의 작물계수를 산정하기 위하여 자료동화 지형모델(GLDAS)에서 생산된 자료와 인공위성을 통한 식생지수와 지점 토양수분을 활용하여 다중회귀모형을 적용하여 작물계수를 제시하였다. 또한 제시한 작물계수를 전체 기간과 작물의 생장 기준으로 기간을 나누어 청미천 유역에서 평가하였다. 그 결과를 정리하면 다음과 같다.

1) GLDAS의 작물계수를 산정하기 위해 우선적으로 지점의 Flux tower 자료와 실제증발산과 잠재증발산을 비교하였다. 이에 대한 결과는 높은 경향성을 나타내었으며 정량적으로도 유효한 범위에 위치하였다. 하지만 GLDAS의 값이 전반적으로 약간 과대산정 되었으며 특히 여름철 일 증발산량이 증가할수록 과대산정 폭이 증가하였다.

2) MODIS 위성을 활용한 식생지수와 지점의 토양수분에 다중선형회귀모형을 이용하여 작물계수를 산정한 결과, 일부분 기상조건에 의해 경향에서 왜곡된 데이터들이 존재하였다. 그러나 대부분의 작물계수가 계절별 경향을 잘 나타낸 것으로 보아 향후 작물계수를 산정함에 있어 다중선형회귀모형을 활용하는 것이 타당할 것으로 판단된다.

3) 두 방식을 통해 산정된 작물계수는 전반적으로 실측값과 유사한 패턴을 나타내었다. GLDAS, SM&VI 회귀모형의 Bias는 각각 0.031과 -0.004를 보였으며, RMSE 는 0.092, 0.069, IOA는 0.944, 0.958의 값을 나타내었다. 도출된 통계량적인 측면에서 보면 SM&VI 회귀모형 방식이 GLDAS 방식보다 더 우수한 결과를 보였다.

4) 계절적 특성을 확인하였을 때 GLDAS와 SM&VI 회귀모형을 이용한 작물계수의 통계값은 각각 초기(상관계수: 0.730, 0.583, 정규화된 표준편차: 0.994, 1.000), 중간(상관계수: 0.722, 0.583, 정규화된 표준편차: 1.089, 0.965), 마지막(상관계수: 0.651, 0.775, 정규화된 표준편차: 0.910, 1.076)의 특성을 보였다. 초기와 중간 기간에 서는 GLDAS가 우수한 결과를 나타내었으나 마지막 기간에서는 SM&VI 회귀모형이 우수한 결과를 보였다.

5) 본 연구에서 활용한 SM&VI 회귀모형 작물계수 산정방법은 인공위성 센서의 한계점이 보완된다면 본 연구에서 수행되었던 방식을 통해 미계측 유역 또는 결측이 발생한 관측소에 실제 현장에 활용될 수 있을 것으로 판단된다.