Introduction
Materials and Methods
시험지역 및 실험 설계
드론 영상 취득
바로미2 가루쌀 수량 조사
토양 시료 채취 및 분석
흙토람 경지 현황 조사
기상 자료 취득
Sentinel-2 NDVI 시계열 NDVI 및 OffsetV 산출 절차
데이터 전처리 및 통계 분석
Results and Discussion
지역별 토양 특성의 차이
지역별 바로미2 벼 수량 변동성
단순/다중회귀 분석의 한계 및 Random Forest 분석의 필요성
기상 및 토양 요인과 수량의 관계
Sentinel-2 기반 NDVI와 수량의 관계
연구의 의의와 한계
Conclusions
Introduction
벼 (Oryza sativa L.)는 우리나라에서 가장 중요한 식량 작물로, 안정적인 생산은 기후변화와 토양 열화가 심화되는 현대 농업에서 가장 중요한 과제이며, 국가 식량 안보와 직결된다 (Statistics Korea, 2024; Kim et al., 2025). 최근에는 기능성 및 가공용 쌀 품종에 대한 수요 증가에 따라 가루쌀 품종의 재배가 확대되고 있으며, 그중 조생종의 바로미2 품종은 제분 적합성 갖추어 기존의 품종에 비해 가공에 용이한 것 장점으로 주목받고 있다 (Kwak, 2025). 그러나 2019년도에 등록된 바로미2 벼는 다른 품종에 비해 수량에 영향을 미치는 인자에 관한 연구가 전무한 실정이며 (Ha et al., 2022), 이는 안정적인 원료 확보와 가공 산업 활성화에 제약 요인으로 작용할 수 있다.
벼 수량은 기상, 토양 화학성, 지형, 관리 등 복합적 요인에 의해 결정되며, 이 중에서도 토양 양분 균형 및 지형적 조건은 생육과 수량의 공간적 변동성을 유발하는 핵심 요인으로 지적되어 왔다 (Deng et al., 2015; Franch et al., 2021; Quille-Mamani et al., 2023). 그러나 이러한 토양 요인들은 기상 및 기타 환경요인들과 복합적으로 작용하기 때문에, 단순 선형 회귀 분석으로는 벼 수량을 충분히 설명하기 어렵다. 최근 농업 환경 데이터의 비선형성과 상호작용 효과를 효율적으로 분석할 수 있는 머신러닝 기법의 활용이 확대되고 있다 (Breiman, 2001; Becker-Reshef et al., 2010; Shook et al., 2021; Li et al., 2023). 특히 Random Forest (RF) 모델은 다중공선성에 대한 내성이 높고, 변수 간 비선형 관계와 상호작용을 포착할 수 있다는 장점이 있다 (Cao et al., 2021; De Clercq and Mahdi 2025). 이에 따라 RF 기반 특징 선택 (recursive feature elimination with cross-validation, RFECV) 과 순열중요도 (permutation importance) 분석은 농업 환경요인의 상대적 기여도를 정량화하고 모델 해석력을 높이는 효과적인 도구로 활용되고 있다 (Paudel et al., 2023).
집약적인 농경지에서는 장기간의 비료 및 농약 사용으로 인해 양분 불균형이 자주 나타난다 (Hwang et al., 2025; Lee and Kim, 2025). 이러한 불균형은 근권에서 유효한 양이온의 흡수를 억제하고 세포막 투과성을 저하시켜 광합성 효율을 감소시킴으로써 벼 생육과 수량에 부정적인 영향을 미친다 (White, 2011). 특히 간척지 토양은 일반 논에 비해 Mg 과잉, 낮은 유기물 함량, 배수 불량, 염분 축적 등의 특성을 보여 벼의 생육과 양분 균형을 저해한다 (Lee et al., 2013; Jeong et al., 2020; Seo et al., 2021). 또한 지형적으로 저지대에 위치해 과잉 수분과 환원 상태가 지속되며, 이는 생육 부진과 수량 불안정으로 이어진다 (Ock et al., 2022; Park et al., 2025). 이러한 특성 때문에 간척지에서는 동일 품종이라도 지역 간 수량 변동이 크게 나타나 생산성 변동의 핵심 지표로 작용한다.
한편, 최근에는 원격탐사 기술과 토양 및 기상인자를 통합해 활용하여 작물의 생육 상태와 수량을 예측하려는 연구가 활발히 진행되고 있다. 위성이나 드론 영상으로부터 산출되는 정규화식생지수 (normalized difference vegetation index, NDVI)는 작물의 엽록소 함량, 군락 밀도, 생체량을 반영하여 생육 수준을 간접적으로 평가하는 지표로 널리 사용되고 있다 (Rouse et al., 1974; Tucker, 1979; Ji et al., 2021; Soriano-González et al., 2022). 특히 ESA의 Sentinel-2 위성 기반 NDVI는 10 m 수준의 높은 공간 해상도를 제공하여, 필지 단위의 미세한 생육 변화를 포착하는 데 효과적이다 (Deng et al., 2018).
본 연구의 목적은 국내 주요 가루쌀 재배단지를 대상으로 (1) 바로미2 벼 수확량 변동에 영향을 미치는 주요 환경 요인을 규명하고, (2) 선형 및 비선형 모델의 예측 정확도를 비교하며, (3) 핵심적인 토양·기상·지형 요인의 수량 반응 특성을 규명하고, (4) Sentinel-2 위성를 활용하여 환경 요인에 의해 나타난 수량 차이가 생육 반응으로 어떻게 반영되는지를 보완적으로 평가하는 것이다. 이를 통해 바로미2 벼 생산성 안정화를 위한 토양 관리 및 정밀농업 적용 기반을 제시하고자 하였다.
Materials and Methods
시험지역 및 실험 설계
본 연구는 2023년 500 ha 규모의 부여(Buyeo, BY) 등 10개 지역의 바로미2 가루쌀 벼 (Oryza sativa L. cv. Baromi2) 재배단지를 내에서 선정되었다 (Table 1). 지역별 대표성을 확보하기 위해 3개의 독립된 필지를 선정하였으며, 필지는 동일한 이앙 및 수확 일정으로 관리되는 직사각형 논으로 정의하였다. 각 필지 내에서 무작위로 2 m × 2 m 크기의 소구획 3개를 선정하여 수확 조사구로 설정하였고, 총 90개 조사구 (10지역 × 3필지 × 3반복)를 확보하였고, 각 필지별 수확량 (Yield, kg 10a-1)을 목표 변수로 설정하였다. 모든 지역의 논은 2023년 6월 중순에서 7월 중순 사이에 이앙이 이루어졌고, 수확은 10월 12일 또는 13일에 수행되었다. 각 조사구의 중심점 GPS 좌표는 Leica Captivate GNSS로 기록되었으며, 모든 공간 데이터는 EPSG:4326 좌표로 정규화하여 관리하였다.
Table 1.
Geographic locations of the ten study sites and cultivated areas of the Baromi2 rice fields.
드론 영상 취득
선정된 필지의 공간적 분포 및 조사구 위치를 제시하기 DJI Mavic 3 Multispectral 드론 (지상 샘플링 거리 ≤ 3 - 5 cm)을 이용해 수확기 전후 RGB 영상을 획득하였다 (Fig. 1). 본 연구에서 드론 영상은 선정된 대표 필지의 시각화 용도로만 사용되었고, NDVI값 추출을 통한 정량 분석은 Sentinel-2 위성을 활용해 수행되었다.
바로미2 가루쌀 수량 조사
수량 조사는 각 필지 내에서 무작위로 2 m × 2 m 크기의 소구획 3개를 조사구로 선정하여 실시되었다. 각 조사구에서 수확은 2023년 10월 12-13일에 걸쳐 이루어졌다. 이후 탈곡·정조 과정을 거친 후 15% 수분함량으로 건조하여 천립중을 측정하였으며, 재식밀도 측정의 절차를 거쳐 10a 환산 수량 (kg 10a-1)으로 표준화하였다 (IRRI, 2013; Ma et al., 2019). 각 필지의 수량은 조사구 (n = 3) 평균으로 산정하고, 지역 단위 비교 시 필지의 평균값을 사용하였다 (총 조사구 = 90).
토양 시료 채취 및 분석
토양 시료의 채취, 전처리 및 분석은 농업과학기술 연구조사분석기준 (NAAS, 2010)에 따라 수행되었으며, 모든 토양분석은 한국농업기술진흥원 (KOAT)에 위탁하여 수행되었다. 토양 시료는 벼 수확일 (2023.10.12 - 13)에 각 필지별 3개 조사구 (2 m × 2 m)의 0 - 20 cm 깊이에서 채취하였으며, 시료는 풍건 후 2 mm 체거름을 하였다. pH (pH1:5)와 전기전도도 (EC1:5, dS m-1)는 1:5 토양 현탁액을 통한 초자전극법, 유기물 (OM, g kg-1)은 Tyurin법을 통한 적정법, 유효인산 (Avail-P2O5, mg kg-1)은 Lancaster법을 통한 비색법, 교환성 양이온 (Ca, Mg, K, cmolc kg-1)은 초산암모늄 침출법으로 추출 후 ICP-OES 정량, 유효규산 (SiO2, mg kg-1)은 NaOAc 침출법을 통한 비색측정법으로 분석을 실시하였다. 필지의 고도 (Altitude)는 Google Earth (earth.google.com/web)를 활용하여 조사구 중심점 좌표를 기준으로 기록되었고, 재식밀도는 수확한 조사구의 실제 식재 간격 (이식 폭 × 주간 거리, m) 조사로 산정하였다.
흙토람 경지 현황 조사
필지의 주소를 기반으로 농촌진흥청 흙토람-비료사용처방 (https://soil.rda.go.kr/sibi/sibiExam.do) 데이터베이스를 참조하여 필지의 토양통, 지형, 심토토성, 배수등급, 토양특성을 확인하였다 (Supplementary Table 1).
기상 자료 취득
성장기간 (이앙-수확) 동안의 일별 기온 (°C), 강수량 (mm), 강우일수 (day), 일조시간 (hrs day-1), 일사량 (MJ m2) 자료는 각 필지가 위치하는 행정구역 내의 기상관측소에서 취득하였다. Table 3은 기상인자의 구분 및 정의에 대한 설명이다.
Sentinel-2 NDVI 시계열 NDVI 및 OffsetV 산출 절차
본 연구에서는 토양 요인에 의해 발생하는 벼 수량의 공간적 변동성을 가시화하기 위해 Sentinel-2 위성 기반 NDVI 값을 보조지표로 산출하였다.
Google Colab 환경에서 Google Earth Engine (GEE)을 이용하여 Sentinel-2 Level-2A (S2_SR) 지상반사 영상으로부터 조사구 (2 m × 2 m) 단위의 NDVI 시계열을 구축하였다. 조사구 중심을 포함하는 10 m × 10 m 공간해상도의 Sentinel-2 픽셀을 기본 분석 단위로 설정하였으며 (Fig. 1), 연구 필지의 경계는 GeoPackage (.gpkg) 형식으로 업로드하였다. 시계열 범위는 2023년 Day of Year 155 (6월 4일)부터 327 (11월 23일) 까지로 제한하여 이앙기부터 수확기를 포괄하였다.
Sentinel-2 위성은 약 5일 주기로 동일 지역을 재관측하지만, 실제 활용 가능한 영상은 구름, 연무, 또는 그림자 등으로 인해 불규칙하게 확보된다. 특히 본 연구 기간 중 7월에는 강수와 구름이 빈번하여 해당 시기의 NDVI 자료 대부분이 마스킹되어 누락되었다. 이러한 관측 간격의 불균일성과 누락값 (missing data)은 시계열 분석 시 심각한 왜곡을 초래할 수 있으므로, 노이즈 제거 및 시계열 복원을 위해 평활화 과정을 수행하였다.
구름 및 그림자 마스킹은 Sentinel-2의 QA60 품질밴드와 s2cloudless 확률층 (Cloud Probability ≤ 60%)을 병행 적용하였으며, Scene Classification Layer(SCL)의 3, 8, 9, 10, 11 클래스 (구름·그림자·적설 등)는 모두 제외하였다. 이러한 마스킹 이후에도 남는 센서 또는 마스킹 오류 의심치는 사분위 범위(IQR) 기반 이상치 탐지 규칙으로 제거하였다.
최종적으로 클리핑된 영상에 대해 10 m 공간해상도를 가지는 Red 밴드 (ρRed)와 NIR 밴드 (ρNIR)의 반사도를 일별로 추출하고, zonal statistics 를 통해 날짜별 NDVI를 계산하였다 (Rouse et al., 1974; Tucker, 1979). NDVI는 -1 에서 1 사이의 값을 갖으며, 값이 높을수록 식생의 활력이 높음을 의미한다.
NDVI 시계열은 Savitzky-Golay (SG) 필터 (창 = 5, 차수 = 3) 를 이용해 평활화하였다. 이 과정은 노이즈를 줄이면서도 생육 곡선의 전환점 (기울기 및 곡률 변화)을 보존하기 위함이다 (Savitzky and Golay, 1964; Chen et al., 2021). 평활화된 NDVI 곡선으로부터 ‘생육 전개 (녹음기) - 정점 - 노화 (등숙/고엽)’의 세 구간을 분할하고, 각 구간의 전환시점(T)과 해당 시점의 NDVI 값(V)을 도출하였다.
본 연구에서 OffsetT는 생육 종결 시점 (10월 12일 또는 13일)으로 정의하였으며, OffsetV는 단순히 수확 직전 NDVI가 아니라, SG 평활화된 NDVI 생육곡선에서 데이터 기반으로 결정된 생육 종결 시점의 NDVI 값을 의미한다. OffsetV는 수확기 수관의 지속성 (stay-green) 및 광합성 유지력 (canopy persistence) 을 정량화하는 현상학적 지표로서, 벼 수량 변동의 생리학적 신호를 포착하는 후기 생육 특성이다 (Ghaley 2012; Moldenhauer et al., 2021).
Table 2.
Description of meteorological variables and their definitions. All variables were aggregated by phenological stage (vegetative, reproductive, and ripening).
| Variables | Definition |
|
Precip_Mean_VE, Precip_Mean_RE, Precip_Mean_RI | Mean precipitation (mm) during VE1, RE2, RI3 |
|
Precip_Max_VE, Precip_Max_RE, Precip_Max_RI | Maximum daily precipitation (mm) during VE, RE, RI |
|
Precip_Ac_VE, Precip_Ac_RE, Precip_Ac_RI | Total precipitation (mm) during VE, RE, RI |
|
Precip_frieq40_VE, Precip_frieq40_RE, Precip_frieq40_RI | Days with precipitation 40 - 80 mm during VE, RE, RI |
|
Precip_frieq80_VE, Precip_frieq80_RE, Precip_frieq80_RI | Days with precipitation over 80 mm during VE, RE, RI |
|
Precip_frieq_VE, Precip_frieq_RE, Precip_frieq_RI | Days with precipitation during VE, RE, RI |
|
Temp_Mean_VE, Temp_Mean_RE, Temp_Mean_RI | Mean air temperature (°C) during VE, RE, RI |
|
Temp_Min_VE, Temp_Min_RE, Temp_Min_RI | Minimum air temperature (°C) during VE, RE, RI |
|
Temp_Max_VE, Temp_Max_RE, Temp_Max_RI | Maximum air temperature (°C) during VE, RE, RI |
|
Temp_Ac_VE, Temp_Ac_RE, Temp_Ac_RI | Total accumulated air temperature (°C) during VE, RE, RI |
|
Insol_Mean_VE, Insol_Mean_RE, Insol_Mean_RI | Mean solar irradiance (MJ/m2) during VE, RE, RI |
|
Insol_Ac_VE, Insol_Ac_RE, Insol_Ac_RI | Total accumulated solar irradiance (MJ/m2) during VE, RE, RI |
|
Sun_Mean_VE, Sun_Mean_RE, Sun_Mean_RI | Mean of daily sunshine duration (hrs/day) during VE, RE, RI |
|
Sun_Ac_VE, Sun_Ac_RE, Sun_Ac_RI | Total accumulated sunshine duration (hrs/day) during VE, RE, RI |
데이터 전처리 및 통계 분석
벼수량·토양화학성·기상자료, Sentinel-2 위성기반 NDVI (OffsetV)는 단일 데이터베이스로 통합하였다. 모든 데이터는 필지별 고유 ID와 위치좌표를 기준으로 결합되었으며, 동일 필지에서의 반복 측정값은 평균값으로 대표하였다. NDVI 시계열의 이상치 (센서/마스킹 오류 의심)는 사분위 범위 (interquartile range, IQR) 기반 이상치 탐지법을 적용하고 제거하였다. 모든 수치형 변수는 표준화 (Z-score scaling)하여 척도의 차이에 따른 편향을 방지한 후 모델 입력자료로 사용하되었다. 이상의 전처리 절차를 통해 생성된 통합 데이터는 Random Forest (RF) 및 Multiple Linear Regression (MLR) 모델의 입력자료로 사용하였다.
벼 수량과 개별 예측변수 간의 일차적 관계를 파악하기 위해 Pearson-Spearman 상관분석을 수행하였으며, 상관성이 유의한 변수를 선별하여 단순선형회귀 (simple linear regression) 를 적용하였으며, 다중 회귀 관계 검증을 위해 MLR 모델을 구축하였다. 모든 모델은 분산팽창계수 (variance inflation factor, VIF) ≤ 5를 기준으로 다중공선성이 낮은 변수를 최종 입력값으로 채택하였다. 이는 상호 상관된 변수의 중복 효과를 제거하여 회귀계수의 신뢰성을 확보하기 위함이다 (O’Brien, 2007).
벼 수확량 예측의 비선형 관계를 반영하기 위해 Random Forest (RF) 모델을 구축하였다. 모델의 하이퍼파라미터는 n_estimators = 200, max_features = sqrt, random_state = 42로 설정하였다. RF 모델의 검증은 지역 간 공간 종속성 (spatial dependence)을 방지하기 위해 5-fold Group 교차검증으로 수행하였으며, 전체 데이터는 학습:검증 = 80:20 비율로 분할하였다. 동일 지역 내 조사구를 하나의 fold로 묶어 학습·검증 세트를 분리하였으며, 예측 잔차의 공간적 독립성은 Moran’s I 공간자기상관 통계 (PySAL, 거리 500 m) 로 검증하였다. 잔차의 공간 자기상관이 0에 근접하여 공간 누출 (spatial leakage)이 발생하지 않았음을 확인하였다. 각 fold에 대해 학습-검증을 반복 수행한 후, 5회 반복 결과의 평균 R2, RMSE, MAE, 및 Bias를 최종 성능 지표로 제시하였다. 변수 중요도는 순열중요도 (permutation importance, PI) 로 산출하였고, 모델별 fold 결과를 평균하여 변수별 기여도를 정량화하였다. 이는 각 변수를 무작위로 섞은 뒤 모델 성능 저하 폭으로 변수의 상대적 기여도를 측정하는 방법이다 (Breiman, 2001). Recursive feature elimination with cross-validation (RFECV) 을 추가로 적용하여 불필요한 변수를 제거하고 최적 R2 및 RMSE 변화율 < 1% 를 기준으로 최적 변수 조합 (n_feature)을 결정하였다 (Guyon et al., 2002). 모델의 반응함수를 시각적으로 검증하기 위해 부분종속도곡선 (partial dependence plot, PDP) 을 생성하였다. 이는 특정 변수의 값을 변화시키면서 예측값의 변화를 관찰하는 방법으로, 변수 수준별 수확량 반응을 해석하는 데 활용하였다. Random Forest 모델의 예측결과를 300회 부트스트랩 재표집을 통해 95% 신뢰구간 (confidence interval, CI)을 계산하고, 이를 PDP에 음영 리본 형태로 추가하였여 불확실성을 확인하였다.
보조 지표로 활용된 OffsetV 는 Sentinel-2 NDVI 시계열에서 Savitzky-Golay 필터 (창 = 5, 차수 = 3) 를 적용해 평활화한 후, 생육곡선의 종결 시점 (offsetT)에서의 NDVI 값으로 정의하였다. 이 지표는 등숙기 수관의 광합성 지속성과 토양 제약 조건을 반영하는 보완적 생육 특성으로 활용하였다.
모든 통계분석은 Google Colab 환경 (Python 3.10; scikit-learn 1.3, pandas, NumPy, statsmodels)에서 수행하였다.
Results and Discussion
지역별 토양 특성의 차이
본 연구에서 조사된 10개 지역의 토양 화학성은 간척지와 비간척지 간 뚜렷한 차이를 보였다. 간척지 토양은 비간척지에 비해 전기전도도 (EC)와 토양 Mg 농도가 높고, OM 함량과 필지의 고도 (altitude)가 낮은 경향을 나타냈다 (Table 3). 특히 GJ를 제외한 간척지 지역 (DH, HN, MA)에서는 교환성 토양 Mg이 과도하게 축적되어, 교환성 Ca/Mg 비율이 3 이하로 떨어졌다. 이러한 값은 양분 불균형이 발생할 수 있는 수준으로, 작물 생산성을 제한할 수 있음을 시사한다 (Jeong et al., 2020; Seo et al., 2021; Hailemariam et al., 2023).
Table 3.
Soil chemical properties, altitude, and planting density across the ten study regions (mean ± standard deviation, n = 3).
| Region | BY | DJ | GJ | GS | HN | IS | JE | MA | NW | SC |
Range (paddy)1 |
|
Soil properties | |||||||||||
| pH1:5 | 4.9 ± 0.1 | 5.6 ± 0.4 | 4.9 ± 0.1 | 5.1 ± 0.1 | 5.4 ± 0.1 | 4.7 ± 0.0 | 4.7 ± 0.1 | 5.3 ± 0.1 | 4.9 ± 0.1 | 4.8 ± 0.1 | 5.5 - 6.5 |
|
EC1:5 (dS m-1) | 1.0 ± 0.0 | 3.7 ± 1.8 | 0.5 ± 0.0 | 0.2 ± 0.1 | 4.5 ± 0.5 | 0.3 ± 0.0 | 0.5 ± 0.2 | 4.1 ± 0.3 | 0.2 ± 0.1 | 0.2 ± 0.0 | <2.0 |
|
Ca (cmolc kg-1) | 5.8 ± 0.2 | 5.0 ± 0.7 | 3.4 ± 0.3 | 4.5 ± 0.2 | 4.4 ± 0.2 | 4.0 ± 0.1 | 4.4 ± 0.3 | 3.9 ± 0.5 | 2.1 ± 0.1 | 1.8 ± 0.1 | 5.0 - 6.0 |
|
Mg (cmolc kg-1) | 1.3 ± 0.1 | 3.1 ± 0.9 | 1.7 ± 0.1 | 0.8 ± 0.1 | 7.2 ± 0.1 | 1.0 ± 0.0 | 1.2 ± 0.2 | 4.8 ± 0.2 | 0.4 ± 0.0 | 0.4 ± 0.0 | 1.5 - 2.0 |
|
K (cmolc kg-1) | 0.4 ± 0.0 | 0.6 ± 0.2 | 0.3 ± 0.0 | 0.2 ± 0.0 | 1.6 ± 0.1 | 0.1 ± 0.0 | 0.2 ± 0.0 | 1.0 ± 0.0 | 0.2 ± 0.0 | 0.1 ± 0.0 | 0.2 - 0.3 |
|
OM (g kg-1) | 0.04±0.01 | 0.02±0.00 | 0.03±0.00 | 0.05±0.00 | 0.03±0.00 | 0.04±0.00 | 0.04±0.00 | 0.03±0.00 | 0.03±0.00 | 0.02±0.00 | 0.02 - 0.03 |
|
SiO2 (mg kg-1) | 144 ± 23 | 193 ± 117 | 73 ± 10 | 339 ± 90 | 295 ± 5 | 91 ± 6 | 400 ± 89 | 136 ± 5 | 89 ± 16 | 67 ± 6 | >157 |
|
P2O5 (mg kg-1) | 67 ± 3 | 27 ± 16 | 165 ± 17 | 83 ± 8 | 21 ± 1 | 42 ± 3 | 34 ± 4 | 41 ± 8 | 204 ± 69 | 134 ± 6 | 80 - 120 |
|
Plant density (plant 3.3 m-2) | 73.3 ± 2.5 | 78.6 ± 3.2 | 72.2 ± 7.0 | 78.6 ± 5.9 | 61.1 ± 5.1 | 84.6 ± 5.0 | 68.1 ± 1.2 | 78.6 ± 2.9 | 78.6 ± 7.8 | 68.8 ± 9.0 | - |
|
Altitude (m) | 6.6 ± 0.7 | 3.5 ± 0.8 | 4.5 ± 0.5 | 95 ± 0.2 | 0.1 ± 0.0 | 6.3 ± 0.3 | 8.2 ± 1.0 | 1.9 ± 0.4 | 100 ± 0.6 | 51.7 ± 0.2 | - |
| Ca/Mg ratio | 4.3 ± 0.2 | 1.8 ± 0.7 | 2.0 ± 0.0 | 5.7 ± 0.3 | 0.6 ± 0.0 | 4.1 ± 0.0 | 3.6 ± 0.2 | 0.8 ± 0.1 | 5.2 ± 0.3 | 4.6 ± 0.2 | - |
지역별 바로미2 벼 수량 변동성
전체 10개 지역에서 바로미2 벼의 평균 수량은 약 371 kg 10a-1였으나, 지역 간 변동 폭은 매우 컸다 (Fig. 2). 간척지 지역은 평균 수량 (kg 10a-1)이 낮을 뿐만 아니라 변동 폭이 크고 불안정한 수량 분포를 보였다. 반면 비간척지 지역 (kg 10a-1)에서는 비교적 안정적인 수량이 확보되었다. 이는 간척지 토양의 특징이 벼 수량 변동성에 중요한 원인으로 작용했음을 시사한다 (Lee et al., 2013; Ma et al., 2020; Seo et al., 2021; Ock et al., 2022).
단순/다중회귀 분석의 한계 및 Random Forest 분석의 필요성
상관 분석 및 단변량 회귀 분석 결과, 여러 토양 및 기상 요인이 벼 수확량과 유의미한 상관관계를 나타냈다. 생식생장기 동안의 최저 기온과 일조 시간, 그리고 유기물 (OM), 마그네슘 (Mg), 전기전도도 (EC)와 같은 토양 특성은 중간 정도의 선형 관계 (R2 최대 0.46)를 보였다. 그러나 어떤 단일 예측 변수도 수확량 변동의 절반 이상을 설명하지 못했다 (Supplementary Table 2). 예측 변수 간의 다중공선성을 해결하기 위해 분산팽창계수 (VIF) 분석을 적용해 모든 변수가 VIF ≤ 5 기준을 충족할 때까지 순차적으로 제거했고, 그 결과 얻어진 예측 변수 하위 집합을 사용하여 다중 선형 회귀 분석을 수행했다. 그러나, 다중 선형 회귀 모델 역시 낮은 결정계수 (R2 = 0.53)와 높은 예측오차 (RMSE = 33.7 kg 10a-1)를 보여 모델의 설명력이 제한적이었다 (Supplementary Table 2). 이러한 결과는 벼 생산성이 다차원적 요인에 의해 영향을 받으며, 선형적 접근만으로는 변수 간 비선형적 상호작용을 충분히 반영하기 어려움을 시사한다.
반면, 랜덤 포레스트 (random forest, RF) 모델은 R2 = 0.62와 RMSE = 27.8 kg 10a-1로 예측 성능이 현저히 향상되었으며, 이는 비선형적 관계와 변수 간 상호작용을 효과적으로 포착했기 때문으로 해석된다. RF 모델은 또한 변수 중요도 평가 (순열중요도) 및 교차검증 기반 특징선택 (RFECV)이 가능하여, 모델 해석력 측면에서도 선형 회귀보다 우수하였다 (Supplementary Table 2).
기상 및 토양 요인과 수량의 관계
순열중요도 (permutation importance) 분석 결과, 고도 (altitude)가 수확량 변동에 상대적으로 높은 중요도를 보였으며 (평균 중요도 = 0.154 ± 0.28), 그 다음으로 생식생장기의 평균 강수량 (Precip_Mean_RE), 평균 기온 (Temp_Mean_RE), 영양생장기의 극한 강우빈도 (Precip_frieq_80_G), 그리고 토양 Mg 순으로 중요도가 높았다 (Table 4). 이러한 결과는 지형적 요인과 기상 스트레스가 수량 변동의 1차적 결정요인으로 작용하고, Mg 등 토양 인자는 이러한 반응을 매개하는 조절요인으로 기능함을 보여준다.
Table 4.
Top five variables ranked by cross-validated permutation importance with corresponding mean and standard deviation values.
| Feature | Mean | Standard deviation |
| Altitude | 0.154 | 0.280 |
| Precip_Mean_RE | 0.149 | 0.186 |
| Temp_Mean_RE | 0.107 | 0.368 |
| Precip_frieq_80_VE | 0.099 | 0.214 |
| Mg | 0.072 | 0.152 |
교차검증 기반 재귀적 특징제거 (RFECV) 결과, 토양 Mg과 Altitude 두 변수만 포함된 모델이 평균 R2 = 0.62, RMSE = 27.8 kg 10a-1로 최적의 성능을 보였으며, 추가 변수의 포함은 예측 정확도를 유의하게 개선하지 않았다 (Table 5). 이는 Mg 축적 (화학적 제약)과 저지대 지형 (물리적 제약)이 수량 변동의 핵심 제약 요인임을 입증한다 (Table 5).
Table 5.
Results of the recursive feature elimination with cross-validation (RFECV) analysis showing the top five selected predictors for Baromi2 rice yield. Variables are ranked by higher mean R2 and lower RMSE values.
부분 의존도 플롯(PDP) 분석 결과, 토양 Mg 농도가 약 1.5 cmolc kg-1을 초과하는 지점부터 예측 수량이 뚜렷하게 감소하는 경향을 보였다. 고농도의 Mg 구간에서는 예측 수량이 일관되게 낮게 나타났으나, 해당 구간의 표본 수가 적어 CI 폭이 넓어지는 등 예측 불확실성이 상대적으로 높게 나타났다 (Fig. 3a). 이러한 경향은 과도한 토양 Mg으로 인해 Ca 및 K 흡수 저해와 삼투 스트레스 증가가 이루어져 수량이 감소 (Ahn et al., 2016; Ma et al., 2020; Park et al., 2021; Hailemariam et al., 2023)할 가능성을 시사하지만, 통계적 불확실성을 함께 고려할 필요가 있다. 고도(altitude)에 대한 PDP에서도 수량이 낮은 구간(< 10 m)에서 예측치가 상대적으로 낮게 나타났으며 (Fig. 3b), 이는 토양 배수성 저하 및 염류 축적의 영향을 반영하는 것으로 해석된다 (Jeong et al., 2020; Ock et al., 2022; Seo et al., 2021; Park et al., 2025). 순열중요도 (permutation importance)와 교차검증 기반 재귀적 특징 제거 (RFECV) 분석에서 모두 핵심 변수로 식별된 토양 Mg과 고도 (altitude)는, 토양 영양 불균형과 지형적 고저가 수확량 변동성을 결정하는 데 중요한 역할을 함을 시사한다.
Fig. 3
Partial dependence plots (PDP) showing the marginal effects of (a) soil Mg concentration (cmolc kg-1) and (b) Altitude (m) on predicted rice yield (kg 10a-1), derived from the Random Forest model. Shaded ribbons indicate 95 % bootstrap confidence intervals (n = 300), representing the predictive uncertainty around the mean marginal effects.
더불어, 토양 Mg 단독보다는 Ca/Mg 비율 (Ca/Mg ratio)이 수량 변동성을 더 잘 설명하는 것으로 나타났다. 단순 선형회귀 분석 결과, 토양 Mg과 수량 간 관계 (R2 = 0.37, p < 0.05, n = 90)보다 Ca/Mg 비율과 수량 간 관계에서 더 유의한 양의 상관이 확인되었다 (R2 = 0.63, p < 0.001, n = 90) (Fig. 4). 이는 Ca/Mg 불균형이 벼 수량 변동의 주요 토양 화학 요인 중 하나임을 시사한다 (Liu et al., 2023).
또한, 농촌진흥청 흙토람 데이터베이스를 활용해 대상 필지의 토양통을 확인한 결과, 낮은 수량을 기록한 필지들은 대체로 만경통과 광활통 토양으로 분류되었다 (Table 3). 이들 토양은 해안 저지대에 위치하고 배수성과 염분집적 가능성이 높은 특징을 보였다. 이러한 토양 물리·화학적 제한이 기상 변동보다 수량에 더 큰 영향을 미쳤을 가능성이 있다 (Ahn et al., 2016; Park et al., 2021; Hailemariam et al., 2023).
Sentinel-2 기반 NDVI와 수량의 관계
Sentinel-2 L2A 영상으로부터 산출된 NDVI 지표인 OffsetV는 벼 수량과 유의한 양의 상관관계를 보였다 (R2 = 0.52) (Fig. 5). 특히 토양 Mg 농도가 높은 지역에서는 OffsetV값이 낮게 나타나, 토양 양분 불균형이 수관활력과 광합성 지속성을 저하시키고, 결국 수량 감소로 이어졌음을 유추할 수 있다. NDVI 지표인 OffsetV는 토양 요인의 직접적 설명 변수라기보다, 토양 제약이 작물 생육 반응으로 투영된 결과를 나타내는 보완 지표로 해석된다. OffsetV와 수량 간의 유의한 양의 상관관계 (R2 = 0.52)는 등숙기의 수관활력 지속성이 최종 생산성과 밀접히 연계됨을 의미하며 (Wang et al., 2014; Onojeghuo et al., 2017), 이는 Sentinel-2 기반 NDVI가 토양-작물 상호작용을 모니터링하는 비용 효율적인 지표로 활용될 수 있음을 시사한다.
연구의 의의와 한계
본 연구는 간척지에서 바로미2 벼 생산성을 제약하는 주요 원인이 토양 Mg 과잉 축적, Ca/토양 Mg 비율 불균형, 그리고 저지대 조건임을 명확히 규명하였다. 더불어, Sentinel-2 기반 NDVI가 이러한 토양 요인의 결과를 가시화하는 보조 지표로 활용될 수 있음을 제시하였다. 하지만, 조사 필지별 시비 내역이나 비료 사용 이력에 대한 정보가 부재한 한계점이 있어, 양분 관리가 수량 변동성에 미친 영향을 정량적으로 평가하지 못했다. 이는 비료 시비량 및 토양관리 이력의 미반영으로 인한 오차 요인이 모델 설명력 (R2 = 0.62)을 제한했을 가능성을 시사한다.
Conclusions
본 연구는 국내 주요 가루쌀 품종인 바로미2의 수량 변동성을 규명하기 위해 토양·지형·기상 자료와 Sentinel-2 기반 NDVI를 통합 분석한 결과, 토양 Mg 축적과 저지대 지형 (Altitude < 10 m)이 생산성을 결정하는 핵심적인 물리·화학적 제약요인임을 확인하였다. Random Forest 분석에서 두 변수만으로도 예측력이 높았다는 점은, 토양 Mg 과잉과 배수 취약성이 지역 간 바로미2 벼 수량 격차를 유발하는 구조적 원인임을 의미한다. 또한 NDVI 기반 OffsetV 지표는 이러한 토양 제약이 벼의 수관활력으로 어떻게 반영되는지를 효과적으로 포착하여, 현장 토양 정보가 제한된 상황에서도 수량 위험을 조기 진단할 수 있는 실용적 보완지표임을 보여주었다. 따라서 본 연구는 (1) Mg 축적을 완화하는 Ca/Mg 비율 관리, (2) 저지대 배수개선, (3) 위성 기반 생육 모니터링과의 통합 활용이 바로미2 생산성 안정화를 위한 핵심 관리 전략임을 제시한다. 이러한 결과는 간척지와 유사한 환경·토양 조건을 가진 지역에서도 적용할 수 있으며, 향후 정밀농업 기반의 토양-생육-수량 통합진단체계 구축에 기여할 수 있을 것이다.
