1. 서론

작물은 기후 등 끊임없이 변동하는 환경 속에서 생산활동을 하며 저온, 고온, 양·수분의 결핍 또는 과다, 병충해, 재해와 같은 스트레스 요인에 항상 노출되어 있다 (Hong, 1998). 농작물의 작황 추정은 생산량 예측을 통한 수급조절, 가격 예측, 농가소득 보전을 위한 정책 수립 등에 중요한 판단 자료로 사용된다(Hong et al., 2012). 우리나라는 작물별 면적조사와 생산량 예측에 많은 시간과 비용을 들여 매년 통계자료를 생산하고 있으며, 그 중 쌀 생산량 통계는 매년 10월경에 발표되고 있다.

쌀은 우리나라 주곡 작물로서 정부는 가격 안정을 위한 수급관리를 위해 공급량 및 수요량을 추정하여 재배면적을 조정·관리하고 있다. 이를 위해서는 벼 재배면적과 생육상황을 모니터링하고, 생산량을 예측하는 것이 필요하다. 따라서 벼 생육의 진척도와 그 현황에 대한 공간분포를 이해할 수 있다면 값으로 표현되는 통계자료보다 더 활용도 높을 것으로 생각된다.

작물의 수량을 결정하는 데는 많은 요인이 관련된다. 우선 일사, 온도, 강수, 바람 등기 후 요인이 광합성 작용의 원동력이다. 또한, 적정한 수분 및 양분의 공급과 관리, 병충해 및 재해에 따른 수량 감소 등에 따라 작물의 실제 수량이 결정된다(Sadras et al., 2015). 이와 같이 작물 품종의 내재적 특성과 수량 결정에 관련된 인자들이 생육 시기별로 미치는 영향이 다르고 지역별로 생산환경이 다르기 때문에 다양한 방법과 기술이 수량 예측에 활용되고 있다.

원격탐사를 이용한 관측은 다양한 파장 영역을 감지하는 센서를 이용하여 넓은 지역을 대상으로 주기적으로 자료를 획득할 수 있다. 즉, 육안으로 감지하기 어려운 물체나 지표면의 상태를 파악할 수 있다는 이점을 가지고 있다 (Campbell, 1996). 지표면 피복체의 분광반사(reflected) 또는 복사의 방출(emitted)을 측정해서 그 신호를 생물 현상으로 해석하여 대상물의 상태를 알아낼 수 있다.

선행연구를 살펴보면, Hong et al. (2012)은 Aqua위성 에서 2002년부터 2010년까지 촬영된 시계열 MODIS(MODerate resolution Imaging Spectroradiometer) NDVI(Normalized Difference Vegetation Index)와 기상자료를 이용한 우리나라 벼 수량 추정 모형을 제안하였다. 위성정보로부터 추출한 논벼의 출수기 NDVI(x1)와기상자 료(x2)를 각각 독립변수로 하는 단순 중회귀 모형을 작성하였다. 모형의 추정 값과 통계 값과의 차이는 5% 이내로 유지되어 왔으나 최근 오차와 연차간 오차 변동폭이 커지는 등의 한계가 나타나고 있다. 같은 개념의 모형이 북한의 벼 수량 추정에도 활용되었다(Hong et al., 2015). Na et al.(2012)도 MODIS NDVI와 다양한 기상요인을 독립변수로 입력하여 단계적 다중회귀 모형을 작성하여 전국과 시군별로 벼 수량을 추정한 결과, 모형의 추정 값과 통계 값과의 오차는 각각 5.08%와 8.46%로 나타났다. 그러나 대상 지역별 독립변수가 일정하지 않아 모형의 오차에 대한 해석이 어려워 개선 방안을 찾기 어려울 것으로 생각된다. 또한 Na et al. (2013)은 MODIS 산출물을 기반으로 작물의 순일차생산성(Net Primary Production, NPP) 추정을 통해 작황을 판단하는 CASA(Carnegie-Ames-Stanford Approach) 모형을 이용하여 전국과 도별로 벼 수량을 추정하였다. 지역별로 모형의 오차는 3%~12%로 나타나 활용성이 있다고 생각되나, 추정에 사용되는 산물출에 대한 불확실성을 정량화하기 힘들고 모수화가 필요하며 방대한 자료가 필요하다는 단점이 있다.

한편, 벼 수량 및 생산량은 기후요인뿐만 아니라 품종 개발, 물·양분 재배관리 등과 학과 기술의 발전에 따라 달라진다. 이러한 특성 때문에 장기적인 시계열 수량 자료는 내재적으로 비선형적이며, 평균과 오차가 고정적이지 않다(Lizumi et al., 2018; Rojas, 2007; Wu et al., 2007). 따라서 작물 수량 예측 모델링에 있어서 기후요인이 수량에 미치는 영향과 기술적인 요인에 의한 변동을 구분해낼 필요가 있다(Lu et al., 2017).

추세 제거(detrending)는 장기적인 시계열 수량 평균의 변동을 통계적으로 제거하는 기법이다. 기존 연구를 살펴보면 상관성에 기반하여 단순 회귀 모형 또는 이차회귀 모형을 이용하는 것이 대부분이다 (Goldblum, 2009; Mishra and Cherkauer, 2010; Wu et al., 2007). Wu et al. (2007)은 ‘추세(trend)’에 대한 정확한 정의나 추출을 위한 논리적 알고리즘은 없고 대개 내재하는 추세를 유도하는 경험적인 방법으로 시간척도에서 데이터의 변동성을 설명한다고 보고하였다.

본 연구에서는 기존에 Hon g et al. (2012)이 개발했던 MODIS NDVI 및 기상자료 기반의 원격탐사 모형의 입력변수 시계열 범위를 확장하고 모형으로 추정한 벼 수량 값의 추세를 분석하여 제거함으로써 모형의 정확도를 높이는 것을 목적으로 하였다.

2. 재료 및 방법

1) MODIS NDVI 자료취득 및 영상 처리

본 연구에서는 Aqua 위성의 h27 v04, h27v05, h28v05 3개 granule 자료를 2002년부터 2019년까지 다운로드한 16일 단위로 작성된 MODIS NDVI(MYD13Q1) 영상을 활용하였다. Fig. 1과 같이 NDVI 레이어와 QC/QA 레이어를 import하고 각 granule을 영상 접합한 후 USGS (United States Geological Survey)의 MRT(MODIS Reprojection Tool)를 이용하여 영상자료 배포 좌표계인  WGS-84 Sinusoidal을 경위도 좌표로 재투영(reproject)하였다. 또한, 보정계수(scale factor)를 적용하여 영상이 데이터 값으로 표현되도록 변환하여 자료 분석에 사용하였다(Hong et al., 2012).

Fig. 1. Flow chart of MODIS data preprocessing.

논벼 지역에 대한 NDVI 추출을 위한 마스크 영상으로는 농림축산식품부에서 2018년에 제공받은 팜맵을 사용하였다. 팜맵은 ‘농경지 전자지도’로써, 최신항공·위성영상과 현장 실제 조사를 통해 논, 밭, 과수원, 시설재배지, 인삼재배지로 구분되어 있다 (http://data. go.kr). 이 중, 논 지역을 추출하여 이 부분에 해당하는 MODIS NDVI 평균값을 추출하여 자료 분석에 활용하였다.

2) 기상자료 및 통계자료 수집

전국 단위 벼 수량 추정을 위한 변수로서 온도, 강수, 일사 등 격자 형태의 기상자료 생산을 위하여 기상청에서 제공하는 89개 지점의 기상관측소 및 599개 지점의 자동 기상관측 장치(Automatic weather Station, AWS)를 통해 관측되는 일별 기상자료(2002~2019)를 수집하였다. Cressman기법을 이용하여 3km 격자 간격의 자료 생산을 위한 객관분석을 수행하였다(Ahn et al., 2010; Cressman, 1959). 경위도로 투영된 3 km 상세 기후자료를 import 하여 IDL(Interactive Data Language)을 이용하여 일별 온도, 강수량, 일사량으로 각각 매핑하여 위성 자료 값과의 비교 및 공간 분석이 가능한 2차 원격 자료로 변환하여 벼 수량과의 상관성 및 회귀분석에 활용하였다(Hong et al., 2012).

전국 벼 생산량에 대한 과거 자료는 통계청 국가통계포털에서 제공하는 농림어업통계자료 중 2002년부터 2019년까지의 벼(정곡) 단위 면적당 수량자료를 수집하여 분석에 사용하였다(Table1).

Table 1. Paddy rice yield statistics of Korea (Statistics Korea, http://kosis.kr) (unit: kg/10a)

3) 추세 분석과 Detrended yield 산출

서론에서 언급한 바와 같이 벼 수량은 기후 요인 외 품종, 관개 상황, 양분, 재배기술 등 다양한 요인에 영향을 받는다. 이를 기술 구성 요인이라 하며, 필지단위 수량 변이를 추정하는 데 중요한 인자로로 알려져 있다. 그러나 기술 구성 요인은 정량적으로 측정이 불가능하고 대상 지역이 국가 단위일 경우 지역 내 공간 변이를 추정하는 것이 매우 어렵다. 따라서 본 연구에서는 벼 수량 이품종 및 재배기술적 요인에 의하여 자연스럽게 증가되는 추세를 제거(detrending)함으로서 기술구성요인이 벼 수량에 미치는 영향을 최소화하였다. 즉 시간에 따라 자연스럽게 증가하는 수량의 추세를 기술의 발전으로 인한 영향이라고 가정하고 연도와 수량간의 회귀분석을 통해 얻어지는 증가추세를 제거함으로서 Eq. 1과 같이 추세가 제거된 수량(Detrended Yield, DY)을 산출하였다.

DY= Y–(a×year) + b       (1)

여기서, Y는 연도별 논벼의 생산량이며, a, b는 연도와 수량의 회귀분석을 통해 산출된 계산 상수이다.

4) 벼 수량 추정 모형 개선

벼 수량과 NDVI 또는 기상요인과의 상관성 분석을 수행하였고 이를 기반으로 Hong et al. (2012)이 제시한 방법에 따라 수량 추정 회귀모형을 작성하였다. 즉, 위성센서로부터 추출한 논 지역의 출수기 NDVI(x1)와 등숙기 기상자료(x2)를 각각 독립변수로 하여 벼 수량 추정을 위한 단순 중회귀 모형(수량=a * 출수기 NDVI + b * 등숙기 기상자료+ c)을 작성하였다. 또한 추세가 제거되기 전 수량(Y)과 추세가 제거된 수량(DY)을 각각의 종속변수로 하여 각각의 모형을 작성하고 통계자료와 비교를 통해 추세 제거가 모형 정확도에 미치는 영향을 평가하였다.

Y= (a×x₁) + (b×x₂) + c       (2)

DY= (d×x₁) + (e×x₂) + f       (3)

여기서, x₁은 출수기 NDVI, x₂는 등숙기 기상자료이며, a~f는 계산 상수이다.

기존 연구에서 나아가 벼 수량 추정을 위한 본 연구의 고도화 방안을 Table 2와 같이 요약하여 살펴보면, 먼저 논 마스크 영상을 기존 지적도에서 최신 농경지 전자지도인 팜맵으로 교체한 점이다. 또한 모형 작성을 위해 사용한 데이터 수집 기간을 18년으로 확장했고, 작황에 대한 품종 육종 및 재배기술의 발전에 따른 수량 변동의 비선형적인 요인을 제거하였다는 점이다.

Table 2. Factors for improving the rice yield estimation model

3.결과 및 고찰

1) 시기별 MODIS NDVI 및 기상자료와 벼 수량과의 상관성

2002년부터 18년간 벼 생육기간 동안 수집한 시기별 MODIS 16-day NDVI 전국 평균값과 벼 수량과 의 상관성을 살펴보았다 (Table3). Hong et al. (2012)이 8년 동안의 자료를 이용한 결과와 비슷한 경향을 보였고, 마찬가지로 수잉기(booting stage) 전후부터의 16일간의 생육상황을 반영하는 8월 21일 NDVI 값과의 상관성이 가장 높게 나타났다. 이는 드론 또는 지상센서를 이용하여 시기별 식생지수와 벼 수량을 비교한 결과 유수 형성기~수잉기 식생지수와 최종 수량과의 상관성이 가장 높았다는 기존의 연구결과와도 일치하였다 (Kimand Hong, 2008; Lee et al., 2020).

Table 3. Correlation coefficient between rice yields and MODIS NDVI products

본 연구 결과와 마찬가지로 벼의 생육기간 중 측정한 NDVI와 최종 수량과의 상관성은 유수 형성기, 개화기, 출수기 사이에 가장 높은 것으로 알려져 왔다(Hong et al., 1997; Kim and Hong, 2008; Lee et al., 2020). 이 시기는 영양생장기의 생육을 통해 벼 이삭의 크기가 결정되고 생식생 장기로 변화하는 기간이다. 이와 같은 결과를 종합하여 8월 21일 MODIS NDVI를 벼 수량 추정 모형의 입력 변수로 사용하였다.

일장(day length)과 온도에 영향을 받는 벼의 생육기간 동안 월별 및 등숙 기온도, 강수, 일사량 자료와 벼 수량과의 상관성 분석 결과 등숙기 일사량과 상관계수가 가장 높게 나타났다(Table4).이는 벼가 익는 등숙기의 기상은 출수기 이후 쌀의 수량을 결정하는 가장 중요한 환경요인이라는 보고(Cock and Yoshida, 1972)를 뒷받침하고, Hong et al. (2012)의 결과와 일치하는 일관성을 보였다. 등숙기는 우리나라에서 벼의 대표 출수기인 8월 20일부터 35일로 하였다. 벼 수량은 8월 일사량, 등숙기 최저온도와도 상관계수가 높게 나타났으며, 강수량 및 평균 최저 온도와는 부의 상관성을 나타내었다. 이를 바탕으로 출수 이후의 기상요인으로서는는등숙기의 일사량·평균·최저 온도, 8월 강수량을 각각 벼 수량 추정을 위한 한 입력 변수로 사용하였다.

Table 4. Correlation coefficient between rice yield and climate factors

2) 벼 수량의 추세분석과 모형의 개선

2002년부터 2019년까지 수집한 자료의 상관성 분석 결과를 토대로 선택한 8월 21일 MODIS NDVI(x1)와 등숙기 기상 자료(x2; 등숙기 평균 온도, 최저온도, 일사량, 8월 강수량)를 이용하여 벼 수량 추정 모형을 각각 작성하였다 (Table5). NDVI와 등숙기 평균온도·최저온도· 강수량·일사량을 각각 이용한 결과 상관계수는 각각 0.63**, 0.68**, 0.65**, 0.61**로 유의하지만, 모형의 결정계수는 0.39*, 0.46**, 0.43*, 0.37*로 Hong et al. (2012)이 보고한 0.70*~0.78*에 비해 매우 낮게 나타났다. 이는 대상 년도의 범위가 확대되면서 벼 수량의 장기변동에 비선형적 요인이 누적되어 추세에 포함되었기 때문으로 생각된다.

Table 5. Rice yield estimation models based on MODIS NDVI and meteorological data

장기간에 걸쳐 수집되는 작물 수량의 추세가 전반적으로 증가하는 까닭은 대개 품종·재배 관리 등 기술의 발전에 기인하고, 해마다 변동폭을 보이는 이유는 기후의 영향이라 설명하고 있다 (Lu and Gao, 2017). 따라서 장기적인 시계열 수량 평균의 변화 추세를 통계적으로 분석하여 제거함으로써 모형의 정확도를 높일 수 있다.

벼 수량 추정 모형의 개선을 위해서 통계로 발표된 벼 수량의 변동 추세를 살펴보았다. Fig. 2는 2002년부터 2019년까지 벼 수량 통계 값과 그 추세(점선), 추세를 제거한 벼 수량 값과 그 추세(실선)를 나타내고 있다. 연도별로 직선적으로 증가하는 단순회귀 모형으로 나타났다. 연차별 시계열 벼 수량 값의 추세를 제거한 후 DY를 산정하였고, 이를 기반으로 NDVI와 등숙기 평균 온도·최저온도·강수량·일사량을 이용하여 각각 벼 수량 추정 모형을 작성하였다.

Fig. 2. Rice yield and detrended yield with their trends from 2002 and 2019.

벼 수량 평균의 변화 추세를 통계적으로 제거하고 보완하여 벼 수량 추정 모형을 작성한 결과 모형의 결정 계수는 각각 0.78**, 0.80**, 0.80**, 0.77*로 원 모형보다 크게 개선되었으며 (Table6), RMSE(Root Mean Square Error)도 20~21 kg/10a에서 12~13 kg/10a로 크게 감소하였다. 개선된 모형을 이용하여 2002년부터 2019년까지 벼 수량을 추정하고 벼 수량(정곡) 통계 값과 비교한 결과, NDVI 및 온도·강수·일사량 조합 모두 원래의 모형(점선)에 비해 개선된 모형(실선)의 추정 값이 벼 수량 통계 값과 차이가 작게 나타났다(Fig. 3).

Table 6. Rice yield estimation models based on MODIS NDVI and meteorological data after detrending of rice yield

Fig. 3. Comparison of yearly rice yield (Y) and detrended yield (DY) estimated, with yield statistics from 2002 to 2019.

모형에 사용된 기상요인 중에서는 결정 계수와 RMSE를 종합적으로 고려할 때, 8월 강수량, 등숙기 최저온도, 평균 온도, 일사량의 순으로 모형의 정확도가 높게 나타났으나 그 차이는 크지 않았다. 따라서 우리나라 벼 수량 추정을 위해서 위성영상과 기상요인을 이용한 통계 모형의 경우에는 하나의 모형을 선택하는 것보다는 온도, 강수량, 일사량 등의 요인으로 다각적으로 추정하여 종합화 해 나가는 것이 정책적 활용에 더 부합할 것으로 생각된다.

3) 개선된 벼 수량 추정 모형을 이용한 우리나라 쌀 수량 지도 작성

MODIS NDVI와 평균·최저기온, 강수량, 일사량 등 각각의 기상요인을 이용하여 우리나라 2019년 벼 수량 지도를 작성하였다 (Fig. 4). 우리나라의 경우 충청남도, 전라북도 및 전라남도와 같은 서쪽 지역이 동쪽 지역과 비교하여 상대적으로 벼 재배지역이 넓게 분포되어 있으며, 태안, 김제, 해남 등 해안가와 인접한 평야지대의 수량이 높게 나타났다. 이와 같이 벼 수량 지도는 각 지역별 작황 및 수량의 변이를 파악하는 데 중요한 자료로 사용될 수 있으며, 연도별 수량 지도를 작성하면 지역에 따른 벼 수량분포의 변화 추이 파악도 가능할 것으로 판단된다.

Fig. 4. Rice yield map of Korea based on MODIS NDVI and each climate factor, (a) average temperature, (b) minimum temperature, (c) rainfall and (d) solar radiation.

4. 결론

본 연구는 Aqua 위성으로부터 얻어진 시계열 MODIS NDVI(x1)와 기상요인(x2)을 독립변수로하여 벼 수량을 추정하는 기존의 단순 다중회귀 모형을 개선하고자 장기적인 시계열 수량 평균 변화의 추세를 분석하고 제거함으로써 모형의 정확도를 높이는 것을 목적으로 하였다.

2002년부터 2019년까지의 MODIS NDVI 영상과 기상청 기상관측소 및 자동기상관측장치를 통해 관측되는 일별 기상자료를 수집하고 MODIS NDVI와 평균 및 최저기온, 강수량, 일사량을 이용하여 벼 수량 추정 모형을 각각 작성하였다. 또한, 벼 수량의 변화 추세를 분석하고 제거한 DY 값을 이용하여 같은 방법으로 수량 추정 모형을 작성하고 결과를 비교하였다. 그 결과, 벼 수량 평균의 변화 추세를 통계적으로 제거하고 보완하여 작성한 벼 수량 추정 모형의 결정 계수는 각각 0.78**, 0.80**, 0.80**, 0.77*로 0.37*~0.46**으로 나타난 원 모형보다 크게 개선된 것으로 나타났다. 또한, 개선된 모형을 이용하여 2002년부터 2019년까지 벼 수량을 추정하고 벼 수량(정곡)통계 값과 비교한 결과, 원 모형의 RMSE가 20~21 kg/10a로 나타났던 것에 비해 12~13 kg/10a로 크게 감소하였다.

이와 같은 결과를 바탕으로 위성영상과 기상요인을 이용한 경험적 통계 모형을 이용한 우리나라 벼 수량 추정을을 위해서는 장기적인 시계열 수량 평균 변화 추세를 제거한 모형을 사용하고, 벼 생육을 결정하는 일장과 온도와 같은 기상요인 중 하나를 선택하는 것보다는 온도, 강수량, 일사량 등의 요인으로 다각적으로 추정하여 종합화 해 나가는 방향으로 관리할 것을 제안하고자 한다. 또한 금후에는 해당 모형을 비롯한 다양한 방법으로 위성영상을 활용한 우리나라 시·군 단위 벼 수량 추정모형의 개발이 요구될 것으로 판단된다.

사사

본 논문은 농촌진흥청 공동연구사업(과제 번호 : PJ 01478701)의 지원을 받았으며, 이에 감사드립니다.