최근 기후 변화는 전 세계의 농업 생산량에 심각한 영향을 미치고 있다(Yuan et al., 2024). 기후 변화 문제에 대한 과학적 평가를 제공하는 국제기구 Intergovernmental Panel on Climate Change (IPCC) 의 6차 보고서에서는 앞으로 몇 십 년 동안 모든 지역에서 기후 변화가 더욱 심화될 것이라고 예측하고 있다(Intergovernmental Panel on Climate Change, 2021). 특히 1990년 이후 온도 상승 속도가 가속화되었으며, 지구 평균 기온은 10년마다 약 0.19°C 상승하고 있고, 이는 기후 변화가 빠르게 진행되고 있음을 보여준다(Samset et al., 2023). 이러한 온도 상승은 옥수수, 밀, 쌀, 대두와 같은 주요 작물의 수확량 감소에 중요한 영향을 미치는 요소로 분석된다(Lobell et al., 2011).

2016년 Food and AgricultureOrganization (FAO) 보고서에 따르면, 현재와 같은 온실가스 배출과 기후 변화 상황이 지속될 경우, 2100년까지 옥수수 수확량은 20~45%, 밀은 5~50%, 쌀은 20~30% 감소할 것으로 전망하였다(Food and Agriculture Organization, 2021). 이는 기후변화로 인한 고온 뿐만이 아닌, 홍수, 가뭄과 같은 극단적인 기상현상들이 작물의 생장에 심각한 영향을 미치기 때문이다(Chen et al., 2018). 그중 쌀 작물은 수잉기, 출수기, 등숙기의 생장 과정을 거치는데, 이 시기의 기후변화는 쌀 생산량에 큰 영향을 미친다. 기온이 1.5°C 상승할 경우 쌀의 출수기와 개화기가 앞당겨지고 생장 기간이 단축된다. 특히 출수기 전까지 생식기의 기간이 단축되어 쌀 생산량에 부정적인 영향을 미친다(Saud et al., 2022).

과도한 강수는 침수를 유발해 벼 뿌리의 산소 공급을 방해하고, 생장과 수확량을 감소시킬 수 있으며(Rahman et al., 2017), 벼의 출수기와 등숙기와 같은 생장 단계에서 침수가 발생할 경우 벼의 생장에 큰 영향을 미친다(Zhen et al., 2024). 중간 생육기나 개화기에 가뭄이 발생하면 벼의 이삭수와 불임률에 영향을 미쳐 수확량 손실이 최대 92%에 이를 수 있다(Hussain et al., 2022). 또한 가뭄으로 인한 벼의 잎 기공 조절을 통해 수분 증발을 줄이려는 방어 기작의 유발이 광합성 속도를 저하시키고 벼의 생장에 부정적인 영향을 미치며, 활성 산소종이 가뭄 스트레스에 의해 증가하면서, 벼 세포에 손상을 줄 수 있다(Sandeep and Godi, 2023). 이러한 기후 변화는 벼 수확량의 변동성을 크게 증가시키며, 세계적인 식량 안보에 심각한 위협이 되고 있다(Rana et al., 2019). 이에 대응하기 위해 정확한 수확량의 예측은 필수적이다.

기존 연구에서는 전통적인 통계 모델을 사용해 왔지만, 이 모델들은 복잡한 기후 요인과 비선형적 관계를 충분히 반영하지 못하는 한계가 있었다(Feng et al., 2023). 하지만 Artificial Intelligence (AI) 기반 예측 모델은 전통적인 통계 모델보다 더 복잡한 상호작용과 비선형적 관계를 처리할 수 있어, 기후 변화로 인한 변동성에 보다 효과적으로 대응할 수 있다(Bali and Singla, 2022). 특히 위성 영상과 AI 모델을 결합하여 넓은 지역의 벼 작황 상태를 실시간으로 모니터링하고, 정밀한 수확량 예측을 제공할 수 있다(Liu et al., 2022). 기상 자료는 벼의 생장에 필요한 조건을 추적하는 데 중요한 역할을 하며, 기후 변화로 인한 수확의 변동성을 효과적으로 예측할 수 있다(Khaki and Wang, 2019).

그러나 이러한 선행 연구는 대부분 특정 지역 또는 단일 변수에 초점을 맞추는 경우가 많아, 다양한 데이터 출처를 통합한 모델링이 미흡했다는 한계가 있었다. 다양한 선행연구에서 AI를 활용하여 위성 데이터 및 기후 데이터를 결합한 벼 수확량 예측의 가능성을 제시했으나(Liu et al., 2022; Gandhi et al., 2016; Noureldin et al., 2013), 미래 기후 시나리오 데이터를 포함하지 않아 장기적인 기후 변동성을 충분히 반영하지 못했다. 또한, 고도와 같은 지형적 요인을 함께 분석한 연구는 드물며, 이는 지역별 특성을 반영하는 중요한 요소로 작용할 수 있다. 따라서 본 연구는 남한에서의 주요 생장단계인 수잉기, 출수기, 등숙기를 고려한 7월부터 9월까지(Cho et al., 2021) 기상청 종관기상관측소(Automated Surface Observing System, ASOS)의 기상 자료와MODIS 및 VIIRS의 식생지수를 활용하였다.

본 연구는 기후 변화에 따른 벼 재배 지역의 수확량 예측 정확도를 높이는 것을 목적으로 한다. 이를 위하여 농림축산부에서 제공하는 팜맵(Farm Map) 자료를 활용하여 논 마스킹(Masking)을 수행하였으며, Variance Inflation Factor (VIF)와 Recursive Feature Elimination withCross-Validation (RFECV) 기법을 적용하여 주요 변수를 선정하였다. 또한, 지역 기후 모델(Regional ClimateModel, RCM)을 활용한 미래 기후 변동성 예측 데이터와 Shuttle Radar Topography Mission Digital ElevationModel (SRTMDEM) v4.1을 통해 각 지역의 고도 데이터를 반영하여 예측 정확도를 더욱 향상시켰다. 이러한 방법론을 기반으로 AutomatedMachine Learning (AutoML)을 활용한 벼 수확량 예측 모델을 제시하였다.Digital ElevationModel (SRTMDEM) v4.1을 통해 각 지역의 고도 데이터를 반영하여 예측 정확도를 더욱 향상시켰다. 이러한 방법론을 기반으로 AutomatedMachine Learning (AutoML)을 활용한 벼 수확량 예측 모델을 제시하였다.

2.1. 연구 지역

본 연구의 연구지역 선정은 농림축산식품부에서 제공하는 팜맵 지도 및 국가통계포털(Korean Statistical Information Service, KOSIS) 에서 제공하는 2000~2023년 벼 조곡량 자료를 사용하였으며, 단위는 ton/ha로 환산하였다. 연구지역의 단위는 시군으로 설정하였으며, 연구지역 선정 과정에서는 다음과 같은 기준을 적용하였다. 본 연구에서는 2000년에서 2023년까지 조곡량 데이터가 없는 지역을 연구 대상에서 제외하였다. 또한, 연도별 편차가 큰 지역은 높은 변동성으로 인해 모델의 정확도를 낮출 수 있기 때문에, 해당 지역의 벼 조곡량 최솟값이 최댓값의 1/3 이하인 경우 연도별 벼 수확량 편차가 큰 지역으로 정의하여 제외하였다. 아울러, 농지 비율이 낮은 지역에서 Relative RootMean Square Error (RRMSE)가 높게 나타나 예측 오차가 증가할 가능성이 있으며, 이는 모델 신뢰도 저하로 이어질 수 있다는 연구 결과가 보고된 바 있다(Sajid et al., 2022). 이에 따라 본 연구에서는 팜맵의 농경지 픽셀만 추출하여, 농지 비율이 1% 미만인 지역도 연구 대상에서 제외하였다. 위 기준으로 인천광역시, 광주광역시, 울산광역시 등 33개 시군이 제외되었으며, 최종적으로 161개 시군 중 128개 시군이 연구 지역으로 선정되었다. Fig. 1은 연구지역으로 선정된 128개 시군을 나타낸다

Figure 1. Study areas in South Korea: 128 selected cities and counties.

2.2. 사용 자료

기상자료는 기상청의 종관기상관측(Automated Synoptic Observing System, ASOS) 자료를 사용하였다. 습도가 낮으면 잎의 성장과 면적 확대가 억제되어 광합성 효율이 떨어져 벼 생장에 부정적인 영향을 미칠 수 있으며(Stuerz and Asch, 2021), 일사량이 적을 경우, 벼의 광합성 활동이 감소하여 생장과 수확량에 부정적이 영향을 미칠 수 있다(Wei et al., 2021). 일교차는 벼의 생육 단계 중 이삭 형성과 같은 중요한 시기에 수확량 증가와 관련이 있으며, 일교차가 일정 수준 이상 유지되었을 때 벼의 수확량이 증가하는 경향을 보여준다(Lobell, 2007). 강수와 기온 또한 앞서 언급한 바와 같이 벼 수확량에 큰 영향을 미치는 기상 조건으로, 따라서 기상변수로 평균 기온(°C), 일교차(°C), 상대습도(%), 일사량(MJ/m²), 강수(mm)를 선정하였다. Table 1은 이러한 기상변수들을 설명한다.

Table 1 . ASOS data used in this study.

Variable Name	Variables Used for Calculation	8–Day Aggregation	Unit
Precipitation		8-Day Total	mm
Relative Humidity		8-Day Average	%
Diurnal Temperature Range	Minimum Temperature, Maximum Temperature	8-Day Average	°C
Solar Radiation		8-Day Total	MJ/m2
Temperature		8-Day Average	°C

위성자료로는 MODIS와 VIIRS 자료를 사용하였다. MODIS는 NASA의 Terra와 Aqua 위성에 탑재된 다분광 센서로, 하루에 1~2회 지구 전역을 주기적으로 관측하여 대기, 해양 및 육지 생태계의 변화를 모니터링한다. 1999년부터 운영 중이며 식생 상태 평가 지수를 제공하여 농업 및 환경 연구에 필수적인 데이터를 제공하고 있다. VIIRS는 National Oceanic and Atmospheric Administration (NOAA)의 Suomi-National Polar-orbiting Partnership (NPP)와NOAA-20 위성에 장착된 최신 센서로, 기후 및 환경 변화를 실시간으로 감지한다. 식생 지수로는 정규화 식생지수(Normalized Difference Vegetation Index, NDVI), 강화 식생지수(Enhanced Vegetation Index, EVI), 광합성 유효 복사 분율(fraction of absorbed PhotosyntheticallyActive Radiation, fPAR), 잎 면적 지수(LeafArea Index, LAI)를 사용했다. NDVI는 식생의 활력도를 측정하는 지수로 근적외선과 적색 빛의 반사율 차이를 사용하여 계산한다. EVI는 식생이 활발한 곳에서NDVI가 높은 값에 대한 변별력이 떨어지는 단점을 개선한 지수이다. fPAR는 흡수된 광합성 유효 복사율로 식물의 광합성 능력과 수확량을 예측하는 데 중요한 지표이며(Sun et al., 2024), LAI는 단위 면적당 식물의 잎 면적의 총합을 나타내는 지수로 잠재적 수확량을 나타내는 중요한 지표로서(Cheng et al., 2022) 상기 4개의 식생지수를 포함시켰다. Table 2는 사용한 위성 데이터의 정보를 나타낸다.

Table 2 . MODIS and VIIRS data used in this study.

Short Name	Long Name	Spatial Resolution (m)	Temporal Resolution	Variable	Operation Start Date
MOD13Q1	MODIS/Terra Vegetation Indices 16-Day L3 Global 250 m SIN Grid V061	250	16-Day	NDVI, EVI	2000.02.18
MYD13Q1	MODIS/Terra Vegetation Indices 16-Day L3 Global 250 m SIN Grid V061	250	16-Day	NDVI, EVI	2002.07.04
MOD15A2H	VIIRS/NPP Vegetation Indices 16-Day L3 Global 500 m SIN Grid V002	500	8-Day	fPAR, LAI	2000.02.18
VNP13A1	VIIRS/NPP Vegetation Indices 16-Day L3 Global 500 m SIN Grid V002	500	16-Day	NDVI, EVI	2012.01.17
VNP15A2H	VIIRS/NPP Leaf Area Index/fPAR 8-Day L4 Global 500m SIN Grid V002	500	8-Day	fPAR, LAI	2012.01.01

또한 본 연구에서는 추가 분석을 위해 EuropeanCentre forMedium-Range Weather Forecasts (ECMWF)에서 제공하는 재분석 기상 자 료인 ERA5와, 기상청의 국지예보모델(Local Data Assimilation and Prediction System, LDAPS) 자료를 활용하여 ASOS 기반 데이터를 일부 변수에서 대체하고 비교 분석을 수행하였다. ERA5 자료는 2000년부터 2023년까지의 기온 변수를 대체하는 데 사용되었으며, 이는 장기간에 걸친 일관된 재분석 자료로 과거 데이터를 안정적으로 제공하기 때문이다. 최근 4개 연도(2020년부터 2023년)에는 LDAPS 자료를 이용해 강수량, 일사량, 상대습도 변수를 대체하였다. LDAPS 자료는 지역적 특성에 맞춘 예보 자료로, 한국 내에서 더욱 높은 지역 예측 정확도를 갖추고 있기 때문이다. 이와 같이 ASOS, ERA5, LDAPS 세 가지 데이터 세트를 구축하고 모델 예측 성능을 비교함으로써, 데이터 출처에 따른 예측 성능의 차이를 검토하였다. Table 3은 ERA5와 LDAPS 각 데이터셋의 주요 특성 및 연구에 사용된 요약을 나타낸다

Table 3 . Summary of key meteorological datasets and variables: ASOS, LDAPS, and ERA5.

Attribute/Variable	LDAPS	ERA5
Data Source	Korea Meteorological Administration (KMA)	European Centre for Medium-Range Weather Forecasts (ECMWF)
Type of Data	Forecast	Reanalysis
Temporal Resolution	3 hour	1 hour
Time Period Covered	2020–2023	2000–2023
Variables Used	Temperature, Precipitation, Solar Radiation, Relative Humidity	Temperature

본 연구에서는 기후 변화에 따른 벼 수확량 예측에 미래 기상 변동성을 반영하기 위해 지역 기후 모델(RegionalClimateModel, RCM) 자료를 추가적으로 활용하였다. RCM 자료는 세계 기후 연구 프로그램(World Climate Research Programme, WCRP)의 Coordinated Regional Climate Downscaling Experiment (CORDEX) 프로젝트에서 제공하는 CORDEX-East Asia 데이터셋을 사용하였다. 이 자료는 독일 GERICS에서 개발한 REMO2015 모델을 기반으로 하며, 영국 MOHC-HadGEM2-ES 모델의 전 지구 기후 시뮬레이션 결과를 다운스케일링 (Downscaling)한 데이터이다. RCM 자료는 2001~2005년 동안의 Historical 시뮬레이션과 2006~2025년 기간 동안의 RCP 8.5 시나리오 기반 데이터를 포함하며, 기후 변화의 극단적 상황에서의 기상 변동성을 반영한다. 본 연구에서는 0.22도, 약 25km의 일평균 기온(°C)과 일사량(MJ/m²) 변수를 활용하여 기후 변화가 벼 수확량 예측에 미치는 영향을 반영하였다. Table 4는 본 연구에서 사용된 RCM자료의 주요 특성을 요약한 것이다.

Table 4 . Summary of RCM data used in this study.

Attribute/Variable	Value
Data Source	CORDEX-East Asia
Model	GERICS-REMO2015
Forcing GCM	MOHC-HadGEM2-ES
Scenario	RCP 8.5
Spatial Resolution	0.22° (~25 km)
Variables Used	Temperature, Solar Radiation
Temporal Resolution	Daily

또한, 지역별 특성을 고려하기 위해 NASA의 SRTM DEM v4.1을 사용하여 각 시군의 고도를 변수로 추가하였다. 고도 변수는 벼 재배환경에서의 미세 기후 차이를 반영할 수 있으며, 이는 벼 생장과 수확량에 영향을 미치는 중요한 지형적 요인으로 판단된다. Table 5는 최종 매치업(Matchup) 파일에 포함된 변수들에 대한 특성을 나타낸다. 변수 뒤의 숫자는 날짜를 나타내며, 위성 데이터의 경우 날짜를 Day of Year (DOY)로 표현하여, 윤년과 평년의 차이를 고려하여 날짜를 일관성 있게 조정하였다.

Table 5 . Composition of matchup variables.

Variable Type	Variable Name	Time Range	Example Variable Name	Source
Meteorological	Temperature (°C)	July 3 – Sept 29	tmp7_3, tmp7_33	ASOS, ERA5, LDAPS
	Relative Humidit (%)	July 3 – Sept 29	rh7_3, rh7_11	ASOS, LDAPS
	Diurnal Temperature Range (°C)	July 3 – Sept 29	dtr7_3, dtr7_11	ASOS
	Solar Radiation (MJ/m2)	July 3 – Sept 29	sr7_3, sr7_11	ASOS, LDAPS
	Precipitation (mm)	July 3 – Sept 29	pre7_3, pre7_11	ASOS, LDAPS
Satellite	NDVI	July 3 – Sept 29	ndvi185, ndvi193	MODIS, VIIRS
	EVI	July 3 – Sept 29	evi185, evi193	MODIS, VIIRS
	fPAR	July 3 – Sept 29	fPAR185, fPAR193	MODIS, VIIRS
	LAI	July 3 – Sept 29	lai185, lai193	MODIS, VIIRS
Regional Climate	Temperature (°C)	July 3 – Sept 29	tas7_3, tas7_11	CORDEX-East Asia
	Solar Radiation (MJ/m2)	July 3 – Sept 29	rsds7_3, rsds7_11	CORDEX-East Asia
Topographic	Elevation (m)	-	Dem	SRTM DEM v4.1
Target Variable	Rice Yield (ton/ha)	Annual	Yield	KOSIS

2.3. 자료 전처리

본 연구에서는 주요 생장 단계인 신장기, 출수기, 등숙기를 고려하여(Cho et al., 2021), 7월부터 9월까지의 기간을 대상으로 8일 단위 시간 해상도를 설정하였다. ASOS 자료는 일별 데이터를 8일 단위로 합성하여 사용하였다. 상대습도, 일교차, 평균기온은 8일간의 평균값으로 계산하였으며, 강수량과 일사량은 8일 동안의 총합으로 산출하였다. 또한, ASOS는 실측 관측치를 기반으로 하므로 결측값이 발생한 경우 가장 가까운 관측 지점의 데이터를 사용하여 결측값을 보완하였다. 추가적으로, ERA5 자료는 단일면 데이터를 데일리로 제공 받아 NC 파일 형태로 다운로드한 후, 8일 간격으로 합성하여 분석에 사용하였다. LDAPS 자료의 경우 3시간 단위의 데이터를 8일로 합성하였다. 이후 마스크 파일의 해상도에 맞춰 리샘플링(Resampling) 작업 및 마스킹을 수행하였고, 각 시군 내의 구역 연산(Zonal statistics)를 적용하여 최종적으로 시군 단위의 통계값을 산출하였다. 또한 ASOS를 대체하기 위한 기상변수들로 ERA5자료는 ASOS자료에 맞춰 켈빈 단위(K)를 섭씨(°C)로 바꿔 계산하였으며, LDAPS의 경우 3시간 간격의 초당 태양복사 총량(W/m²)을 8일 평균으로 계산하여, ASOS 단위에 맞춰 일사량(MJ/m²)으로 변환하였다. 변환 공식은 식(1)과 같다.

위성 자료의 경우MODIS는 Terra와Aqua 위성에서 수집된 자료로, 16일 주기의 데이터를 제공한다. 2002년 이전에는 Terra 위성만이 16일 주기로 데이터를 제공하였기 때문에 전후 시점의 평균값을 계산하여 8일 단위 데이터를 생성하였다. 2002년 이후 Aqua 위성이 운용되기 시작한 후로는 16일 주기의 데이터를 교차하여 시간 해상도 8일을 유지하였다. 2013년 이후에는 후속 위성 센서인 VIIRS가 도입되었고, MODIS와의 상관관계 분석 결과CorrelationCoefficient (CC)를 기준으로NDVI는 0.959, EVI는 0.959, fPAR는 0.924, LAI는 0.92가 나왔다. 이에 따라MODIS와 VIIRS의 상관관계를 분석한 결과를 토대로 VIIRS 데이터를 대체하여 사용하였다. Fig. 2는MODIS 식생지수와 VIIRS 식생지수의 상관관계 분석 결과를 Density scatter plot으로 나타낸 것이다.

Figure 2. Density scatter plot of correlation analysis between MODIS and VIIRS vegetation indices.

위성 자료의 전처리 과정에서는GeoTIFF 포맷으로 데이터를 추출한 후, GDAL 라이브러리OpenCV의 Telea 인페인팅(Inpainting)기법을 사용하여 갭필링(Gap-filling) 작업을 통해 결측값을 보완했다. 갭 필링은 위성 자료에서 발생하는 결측값을 보완함으로써 모델이 완전한 데이터셋을 사용해 학습할 수 있도록 하기 위함이며, Chen et al., (2021)의 연구에서는 Savitzky-Golay 필터를 활용한 갭필링 기법을 통해 Landsat NDVI 시계열 데이터를 복원한 결과, 정확도가 개선되었음을 보여주었다. 이후, 논지 마스킹을 통해 논지에 해당하는 픽셀만 추출하였으며, 구역 연산을 활용하여 시군 폴리곤(Polygon) 내의 래스터(Raster) 픽셀 값의 평균을 계산하여 최종 매치업 DB를 구축하였다. 최종 매치업 DB는 지역, 위경도, 연도, 수확량, 기상 및 위성 설명변수들로 구성되어 있으며, 설명변수의 경우 9개의 위성 및 기상변수와 12개의 시간 단위의 곱으로 총 108개로 구성되어 있다.

또한 연구에서는 모델의 정확도를 향상시키기 위해 국립농업과학원에서 제공한 30 m 해상도의 논지 마스크 파일을 기반으로 리샘플링 및 마스킹 작업을 수행하고, DEM및 RCM변수를 추가하여 매치업 DB를 재구성하였다. 이 과정을 통해 논지의 실제 경계를 반영한 고해상도의 기상 및 위성 데이터를 구축하고, 기후변동성 및 지역의 특성을 반영함으로써, 연구 지역 내 공간적 해상도를 높여 모델의 예측 정확성을 개선하고자 하였으며, 리샘플링 및 변수 추가 전후로 모델의 성능을 비교하였다.

2.4. 변수 선정

변수 선정 과정에서는 RFECV를 사용하여 중요한 변수를 추출하였고, VIF 분석을 통해 다중 공선성을 확인하였다. RFECV는 교차 검증을 통해 모델 성능에 기여도가 높은 변수를 반복적으로 선택하고, 중요하지 않은 변수는 제거하는 방식으로 모델의 복잡도를 줄이면서 성능을 유지하는 방법이다. RFECV 실험에서는 모든 변수를 포함시켰을 때 가장 성능이 좋게 나왔기 때문에, 특정 변수를 제외하지는 않았다. 또한 이러한 변수들이 다중공선성을 가질 경우에는 회귀계수의 분산을 증가시켜 모델의 안정성을 해칠 수도 있기 때문에(O’Brien, 2007), 일반적으로 사용되는 VIF 기준치인 10을 초과하는(Kwon, 2015) 변수들을 조사하였다. 일부 VIF >10인 변수들이 발견되었으나 그 초과하는 정도가 미미하고, RFECV 실험에서 모든 변수로 구성된 모델의 성능이 가장 높았던 것으로 확인되었기 때문에, Table 5의 변수들을 모두 사용하였다.

2.5. AutoML

본 연구에서는 구축된 매치업 DB를 AutoML을 통하여 벼수확량 예측 최적 모델을 구축하였다. AutoML은 머신러닝 모델의 생성, 훈련, 평가, 선택 과정을 자동화하여, RandomForest (RF), GradientBoosting Machine (GBM)등 다양한 AI 모델에 대해 각각 하이퍼파라미터(Hyperparameter) 최적화를 수행한 후, 리더보드(Leaderboard) 상위 N개의 모델을 앙상블하여 최적의 성능을 가진 모델을 선택하는 기술이다(Feurer et al., 2015). 본 연구에서는 AutoML 라이브러리 중 Autokeras, Pycaret, Scikit-learn의 성능을 비교하였다. Fig. 3은AutoML의 파이프라인(Pipeline)에 대한 개요이다.

Figure 3. An overview of AutoML pipeline covering data preparation (Section 2), feature engineering (Section 3), model generation (Section 4) and model evaluation (Section 5) (He et al., 2021).

2.6. 정량적 검증

본 연구에서는 구축한 모델의 성능을 평가하기 위해 5 폴드 교차검증(5-Fold Cross-Validation) 및 Leave-One-Year-Out Cross-Validation (LOYO) 평가를 진행하였다. 5 폴드 교차검증은 모델의 성능을 평가하기 위해 데이터셋을 5개의 하위 집합(fold)으로 나누어 수행하는 검증 기법이다. 전체 데이터셋을 5개의 폴드로 나눈 후, 각 폴드 중 하나를 검증용 데이터로, 나머지 4개를 학습용 데이터로 사용하여 모델을 훈련하고 평가한다(Brownlee, 2023). 이 과정을 5회 반복해 각 폴드가 검증용 데이터로 한 번씩 사용되도록 하며, 마지막에 5회의 검증 결과를 평균하여 모델의 일반화 성능을 추정한다. 5 폴드 교차검증은 과적합을 방지하고 모델의 성능을 안정적으로 평가하는 데 유리하다(Yadav and Shukla, 2016). 본 연구에서는 Pycaret, Autokeras, Scikit-learn의 성능을 비교하였다. 이 중 Pycaret은 다양한 머신러닝 알고리즘을 통합적으로 지원하며, 하이퍼파라미터를 효율적으로 최적화할 수 있는 점에서 우수한 성능을 보였다.

Pycaret을 통해 최종적으로 선택된 모델은 LightGradient Boosting Machine (LightGBM)이다. LightGBM 모델은 Gradient Boosting 알고리즘을 기반으로 하여 높은 예측 정확도와 효율성을 제공하며, 대규모 데이터 처리와 빠른 학습 속도를 지원하는 장점을 갖추고 있다(Ke et al., 2017). 특히, LightGBM은 Leaf-Wise 성장 전략을 사용하여 트리의 각 단계에서 손실 감소가 가장 큰 리프 노드를 우선적으로 분할함으로써 기존의 Level-Wise 방식보다 더 큰 손실 감소를 달성할 수 있다. Leaf-Wise 성장 전략은 학습 과정에서 트리 구조의 깊이를 효과적으로 조정하여 더욱 세밀한 분류를 가능하게 하며, 이를 통해 모델의 예측 성능을 크게 향상시킨다. 또한, Histogram 기반의 분할 방식을 채택하여 계산 효율성을 극대화하며, 연속형 데이터를 효과적으로 처리할 수 있다. 이러한 특성은 대규모 위성 데이터와 복잡한 기후 데이터를 포함한 본 연구의 데이터셋 처리에 적합하다. Fig. 4는 LightGBM의 Leaf-Wise를 나타낸다.

Figure 4. Visualization of the Leaf-Wise growth mechanism in LightGBM.

LOYO 평가는 각 연도를 검증용 데이터로 한 번씩 사용하고, 나머지 모든 연도로 모델을 학습하여 해당 연도에 대한 예측 성능을 평가하는 방식으로, LOYO 평가는 시계열 데이터의 특성을 반영하여 연도별 예측 성능을 검증하는 데 효과적인 방법이다. LOYO 평가를 통해 모델이 특정 연도에 대해 과대 또는 과소 예측하는 경향이 있는지, 다양한 연도에 걸쳐 얼마나 안정적으로 예측을 하는지 평가할 수 있으며, 이는 시계열 데이터에서 모델의 일반화 성능을 평가하는 데 유용하다(Bergmeir and Benítez, 2012). 모델 평가 지표로는 평균 편향오차(Mean Bias Error, MBE), 평균 절대 오차(Mean Absolute Error, MAE), 평균 제곱근 오차(RootMean Square Error, RMSE), 상관계수(CC), 평균 절대 백분율 오차(Mean Absolute Percentage Error)를 사용하였다.

MBE는 예측값과 실제값의 평균 차이를 나타내며, 모델이 과대 또는 과소 예측하는 경향을 파악하는데 유용하다(Willmott and Matsuura, 2005). MBE의 계산식은 식(2)와 같으며, 여기서 P_i는 예측값, A_i는 실제값, n은 총 데이터 수이다.

MAE는 예측값과 실제값의 절대적인 차이의 평균을 측정하여, 예측 오차의 전반적인 크기를 제공하여, 모든 오차를 동일하게 취급하므로 극단적인 오차에 대한 민감도가 낮고, 예측값과 실제값의 간단한 평균 차이를 확인하는 데 유용하다(Hyndman and Koehler, 2006). MAE의 계산식은 식(3)과 같다.

RMSE는 오차의 제곱을 평균하여 제곱근을 취한 값으로, 큰 오차에 대해 더 높은 가중치를 부여하여 전반적인 예측 정확성을 평가하고, 예측값과 실제값 간의 큰 차이에 민감하여, 모델이 큰 오차를 줄이는 데 효과적인지 판단할 수 있다(Willmott and Matsuura, 2005). RMSE의 계산식은 식(4)와 같다.

CC는 예측값과 실제값 간의 상관 관계를 평가하여, 모델이 실제값의 변동 패턴을 잘 따라 가는지 확인할 수 있다. 상관 관계가 높을수록 예측값이 실제값의 경향을 잘 반영하고 있다는 의미이며, 모델의 적합성을 평가하는 데 유용하다(Legates and McCabe Jr, 1999). CC의 계산식은 식(5)와 같으며, 여기서. $\overline{P}$와 $\overline{A}$는 각각 예측값과 실제값의 평균이다.

MAPE는 예측값과 실제값 간의 오차를 백분율로 나타내는 지표로, 예측 정확성을 직관적으로 이해할 수 있게 도와주며, 오차가 실제값에 대해 얼마나 비율적으로 큰지 알 수 있게 한다(Hyndman and Koehler, 2006). MAPE의 수식은 식(6)과 같다

3.1. 다중공선성 변수 제거에 따른 모델 성능 평가

다중공선성 문제를 확인하고 이를 해결하기 위해 수행한 모델 성능 비교 결과(Table 6), sr9_21 변수를 제외한 모델과 sr9_29 변수를 제외한 모델 간에는 유의미한 성능 차이가 나타나지 않았다. 그러나 두 변수를 모두 포함한 모델에서 가장 높은 성능이 확인되었다. 또한, fPAR 7월과 lai 8월 변수를 포함한 모델과 lai 7월과 fPAR 8월 변수를 포함한 모델 간의 성능 차이 역시 크지 않았으나, 7월과 8월의 fPAR 및 lai 변수를 모두 포함한 모델이 전반적으로 가장 우수한 결과를 보였다. 따라서 본 연구에서는 다중공선성에 의해 모델의 안정성이 약화될 가능성을 인식하면서도, 모든 변수를 사용하는 것이 성능 면에서 더 적합하다고 판단하여 모든 변수를 포함하여 진행하였다.

Table 6 . Model performance evaluation with a variable exclusion.

Excluded	Variable	MBE	MAE	RMSE	CC	MAPE
sr9_21	0.005	0.265	0.348	0.754	4.045
sr9_29	0.001	0.258	0.340	0.765	3.939
Lai185, lai193, lai201, lai209, fPAR217, fPAR225, fPAR233, fPAR241	0.000	0.263	0.343	0.760	4.009
FPAR185, fPAR193, fPAR201, fPAR209, lai217, lai225, lai241	-0.002	0.261	0.342	0.762	3.972
Lai185, lai193, lai201, lai217, lai225, lai233, lai241	0.002	0.260	0.343	0.762	3.964
FPAR185, fPAR193, fPAR201, fPAR209, fPAR217, fPAR225, fPAR241	-0.002	0.261	0.342	0.762	3.988

3.2. 5-폴드 교차검증 결과

본 연구에서는 세 가지 라이브러리(Pycaret, Scikit-learn, Autokeras)를 사용하여 5-폴드 교차 검증을 수행하고, 각 라이브러리의 성능을 비교하였다(Table 7).

Table 7 . Performance comparison of AutoKeras, PyCaret, and Scikit-Learn.

Library	MBE	MAE	RMSE	CC
Pycaret	-0.002	0.257	0.337	0.771
Autokeras	-0.053	0.311	0.399	0.666
Scikit-learn	-0.001	0.286	0.377	0.700

5-폴드 교차검증 결과, Pycaret이 Autokeras 및 Scikit-learn보다 우수한 성능을 보였다. Pycaret의MBE는 –0.002로, 예측값과 실제값 간의 편차가 거의 없는 중립적인 예측을 제공하였다. 이는 Pycaret이 과대 또는 과소 예측 경향 없이 안정적인 모델임을 나타낸다. 반면, Autokeras는 MBE가 –0.053으로 실제값에 비해 다소 낮게 예측하는 경향이 있었으며, Scikit-learn은 MBE –0.001로 Pycaret과 유사한 수준을 보였으나 다른 지표에서 성능이 낮았다. MAE와 RMSE를 기준으로 분석한 결과, Pycaret의 MAE와 RMSE는 각각 0.257과 0.337로 가장 낮은 오차율을 보였다. 이는 Pycaret이 극단적인 오차에도 효과적으로 대응하며, 전반적으로 안정적이고 일관된 성능을 유지함을 시사한다. 반면, Autokeras의MAE는 0.311, RMSE는 0.399로 가장 큰 오차를 보였으며, Scikit-learn은MAE 0.286, RMSE 0.377로 Pycaret보다 다소 높은 오차율을 기록하였다.

상관계수(CC) 분석에서는 Pycaret이 0.771로 가장 높은 값을 보였다. 이는 예측값이 실제값의 변동 패턴을 잘 반영하고 있음을 나타내며, 모델의 예측 일관성과 신뢰성이 우수함을 보여준다. 반면, Scikit-learn은CC0.700, Autokeras는CC0.666으로 Pycaret에 비해 낮은 상관성을 보였다. 따라서 Pycaret은MAE, RMSE에서 가장 낮은 오차율과 CC에서 가장 높은 값을 보여, 벼 수확량 예측에서 가장 안정적이고 일관된 성능을 제공하는 것으로 평가되었다. 또한, 리샘플링 및 RCM변수의 추가에 따른 모델 성능 변화를 검증한 결과는 Fig. 5와 같다.

Figure 5. Scatter plots of model Predictions vs. actual Rice Yields: (a) without 30 m resampling, (b) with 30 m resampling using the paddy field mask, and (c) including RCM variables.

리샘플링 및 RCM변수 추가 전후로 모델의 성능이 전반적으로 개선된 것으로 나타났다. 구체적으로 MAE는 0.257에서 0.245로 감소했고, RMSE는 0.337에서 0.326으로 감소하여 오차가 줄어든 것으로 나타났다. 또한, CC는 0.771에서 0.791로 증가하여 예측값과 실제값 간의 상관성이 높아졌음을 확인하였다. 반면, MBE는 –0.002에서 0.001로 변화하여 약간의 과대 예측 경향이 나타났다.

또한, MAPE가 3.918에서 3.734로 감소하여, 전반적인 예측 정확도가 향상되었음을 확인하였다. 리샘플링 작업 및 RCM변수 추가를 통해 모델 성능이 전반적으로 개선된 후, 본 연구에서는 ASOS, LDAPS, ERA5 세 가지 기상 데이터 세트를 활용하여 동일한 5-폴드 교차검증을 수행하고 각 데이터 출처에 따른 모델 성능을 비교하였다. Table 8의 결과에서 ASOS 기반 모델은 MAE 0.245, RMSE 0.323, CC 0.791, MAPE 3.734로, 가장 높은 예측 정확도와 일관성을 보였다.

Table 8 . Comparison of model performance using different meteorological datasets.

Data Source Type	MBE	MAE	RMSE	CC	MAPE
ASOS	-0.001	0.245	0.323	0.791	3.734
LDAPS	0.003	0.249	0.324	0.789	3.794
ERA5	-0.002	0.248	0.324	0.790	3.782

ASOS 데이터는 지상 관측에 기반한 포인트 자료로, 벼 재배지에서 실제 발생한 기상 조건을 잘 반영하는 특성이 있어 높은 예측 성능을 제공한다. 다만, 포인트 관측 특성상 연구 지역 전체를 대표하는 데는 일부 한계가 있을 수 있다. 반면, LDAPS는 격자 형태의 예보 모델로 지역적 변동성을 반영하고 연구 지역을 광범위하게 커버할 수 있는 장점이 있다. LDAPS 기반 모델은MAE 0.249, RMSE 0.324, CC 0.789를 기록하며 ASOS와 유사한 성능을 보였으나, 약간의 과대 예측 경향(MBE 0.003)이 나타났다.

ERA5는 재분석 자료로 일관된 장기 데이터를 제공하여 기후 변화분석에 유리하나, 상대적으로 낮은 공간 해상도로 인해 국내 기상 변동을 세밀히 반영하기 어려운 한계가 있다. ERA5 기반 모델은 MAE 0.248, RMSE 0.324, CC 0.790을 기록했으며, 약간의 과소 예측 경향(MBE -0.002)이 있었다. 따라서 벼 재배지에서의 실제 기상 조건을 가장 잘 반영하며 우수한 예측 성능을 보인 ASOS 데이터를 사용하여 최종 예측 모델을 구축하였다.

DEM 변수를 포함한 모델은MAE와 RMSE에서 각각 0.210과 0.278을 기록하며, 고도 변수를 제외한 모델(0.245, 0.323)에 비해 성능이 개선된 것으로 나타났다. 또한, CC는 0.851로 증가하여 예측값과 실제값 간의 상관성이 높아졌음을 확인하였다(Fig. 7). 이를 통해, 고도 변수를 포함함으로써 모델이 지역적 특성을 더 잘 반영하고, 예측 성능이 향상되었음을 알 수 있다.

Figure 7. Scatter plot of observed vs. predicted rice yield incorporating the DEM variable.

3.2. LOYO 평가 결과

본 연구에서는 실제 2024년 벼 수확량 예측에 앞서, 다양한 연도별 기상 조건과 수확량 변화에 대한 모델의 예측 일관성을 LOYO 평가를 통해 검증하였다.

LOYO평가 결과, 전반적으로 2013년 이후 데이터에서MAPE가 5% 전후로 안정적인 오차율을 보였으며, 이는 VIIRS 식생지수의 도입이 모델의 예측 안정성 향상에 기여했음을 시사한다. Table 9는 MODIS 식생지수를 사용한 연도별 LOYO 평가 결과를, Table 10은 VIIRS 식생지수를 사용한 연도별 LOYO 평가 결과를 나타낸다. MBE 결과를 살펴보면, 대부분의 연도가 -0.2에서 0.5 사이에 위치하여 예측값과 실제값 간의 편차가 크지 않음을 확인할 수 있다. 이는 모델이 과대 혹은 과소 예측 없이 실제 수확량을 비교적 중립적으로 예측하고 있음을 나타낸다. 다만, 2003년과 2007년과 같은 일부 연도에서는 MBE 값이 높아 과대 예측 경향이 관찰되었는데, 2003년은 역대 최대의 태풍으로 기록된 매미가 상륙하여 전국적으로 강풍과 호우를 동반했으며, 2007년에는 태풍 나리로 인해 2000년 이후 가장 높은 태풍 일강수량을 기록하였다. 이러한 기상 이변으로 인해 해당 연도의 수확량이 크게 감소했으며, 그로 인해 모델이 실제 수확량을 과대 예측하는 경향을 보인 것으로 판단된다. LOYO검증 결과는 5-Fold 교차검증 결과에 비해MAPE, MAE, CC 등 모든 지표에서 낮은 성능을 보였다. 5-Fold 교차검증에서 MAPE는 평균 3.216으로 나타난 반면, LOYO검증에서는 이 값을 달성한 연도가 없었다. 이러한 성능 차이는 두 검증 방법의 구조적 차이와 훈련 데이터의 차이로 보이며, 특히 LOYO 검증에서는 테스트 연도가 훈련 데이터에 포함되지 않으므로, 모델이 특정 연도의 새로운 기후 조건, 농업 관행 변화, 또는 예측 변수의 극단적 값에 대응하지 못하는 경우 성능 저하가 발생할 가능성이 크다. 반면, 5-Fold 교차검증에서는 시간적 경계를 넘나들며 데이터가 훈련 및 테스트에 포함되기 때문에 모델의 성능이 다소 과대평가될 가능성이 있다. 따라서 LOYO검증 결과를 통해 모델의 한계를 인지하고, 이를 보완하기 위한 추가적인 데이터 확보 또는 모델 구조 개선이 필요할 것으로 보인다. 또한, 본 연구는 2024년 벼 수확량을 예측하기 전, 2000년에서 2022년까지의 데이터를 통해 정답 자료가 존재하는 2023년 벼 수확량을 예측하여, 실제값과 예측값의 차이를 확인하였다(Fig. 8). 그 결과, 2023년 실제 수확량과 예측 수확량 간 MAPE는 4.091%로, 본 연구의 모델이 높은 예측 정확도를 가지고 있음을 확인하였다. 이를 통해 2024년 벼 수확량 예측의 신뢰성을 높일 수 있을 것으로 기대된다.

Figure 8. Maps of 2023 rice yields: (a) predicted paddy rice yields and (b) observed paddy rice yields.

Table 9 . Performance evaluation using Leave-One-Year-Out Cross- Validation (LOYO) with MODIS data.

Year	No.	MBE	MAE	CC	MAPE
2000	128	-0.159	0.382	0.131	5.656
2001	128	-0.155	0.301	0.498	4.356
2002	128	0.150	0.349	0.561	5.940
2003	128	0.367	0.424	0.498	7.446
2004	128	-0.126	0.296	0.442	4.374
2005	128	0.063	0.238	0.462	3.665
2006	128	0.038	0.283	0.517	4.393
2007	128	0.551	0.556	0.607	9.162
2008	128	-0.098	0.253	0.630	3.726
2009	128	-0.296	0.397	0.570	5.552
2010	128	0.176	0.322	0.692	5.239
2011	128	0.107	0.316	0.681	5.164
2012	128	0.323	0.393	0.599	6.779

Table 10 . Performance evaluation using Leave-One-Year-Out Cross- Validation (LOYO) with VIIRS data.

Year	No.	MBE	MAE	CC	MAPE
2013	128	0.232	0.340	0.647	5.510
2014	128	0.086	0.274	0.629	4.082
2015	128	-0.178	0.299	0.587	4.121
2016	128	-0.354	0.394	0.618	5.392
2017	128	-0.293	0.370	0.640	5.273
2018	128	-0.179	0.292	0.614	4.124
2019	128	0.015	0.314	0.648	4.666
2020	128	0.105	0.303	0.567	4.868
2021	128	-0.208	0.295	0.588	4.114
2022	128	-0.153	0.305	0.553	4.489
2023	128	-0.203	0.305	0.585	4.407

본 연구는 기후 변화로 인한 비선형적 요인과 복잡한 상호작용을 효과적으로 반영함으로써 국내 벼 수확량 예측의 신뢰성을 크게 향상시킬 수 있음을 입증하였다. 특히 NASA의 SRTMDEM과 CORDEX East Asia의 RCM데이터를 통합적으로 활용하여 지역별 특성과 미래 기후 변동성을 모델에 반영하였고, 이를 통해 2023년 예측 실험에서 MAPE 4.091%라는 높은 정확도를 달성하였다. 그러나 일부 해의 기상 현상, 예를 들어 극단적인 폭우와 태풍과 같은 이상기후에서는 상대적으로 높은 예측 오차가 나타났다. 이는 극단적 이상기후가 모델의 예측 성능에 미치는 영향을 보여주며, 향후 이러한 시기의 예측 정확도를 높이는 것이 빈번해질 이상기후에 대비한 수확량 예측의 정밀도를 향상시키는 핵심 과제임을 시사한다. 향후 연구에서는 극단적인 기상 이변에 대한 보정 기법과 다양한 기상 데이터 출처를 비교·분석하여 모델의 예측 안정성과 적응성을 더욱 강화할 필요가 있다. 이러한 발전은 지속 가능한 농업 정책 수립과 기후 변화 대응을 위한 과학적 기반으로 활용되어, 농업 생산성 증대와 식량 안보 강화에 기여할 것으로 기대된다.

Figure 6. Scatter plots of predicted vs. actual rice yields using different data sources: (a) ERA5, (b) LDAPS, and (c) ASOS.

이 연구는 2024년 국립농원과학원 연구용역 “위성영상과 AI를 활용한 벼 생육 시계열 예·관측 체계 구축” 과제의 지원으로 수행되었습니다.

No potential conflict of interest relevant to this article was reported.