集成多源遙感數(shù)據(jù)與生育期時(shí)序光譜特征的水稻種植面積提取

鄭紫瑞，趙輝杰，位盼盼，方 鵬，王來剛，徐少博

（1.河南世紀(jì)國科空間技術(shù)應(yīng)用有限公司，河南 鄭州 450008；2.河南省農(nóng)業(yè)科學(xué)院農(nóng)業(yè)經(jīng)濟(jì)與信息研究所，河南 鄭州 450002）

水稻作為主要糧食作物之一，為人類社會提供了穩(wěn)定的食物來源。據(jù)統(tǒng)計(jì)，全球有12%的耕地種植水稻，超過50%的人口以水稻為重要食物來源，每人每天消耗的谷物中水稻占19%[1-2]。我國有50%的耕地種植水稻，并且約有65%的人口以水稻為主食[3]。然而，日益增長的全球人口對糧食供應(yīng)提出了更大的需求，城市化影響、環(huán)境問題、氣候變化也對水稻的生長和種植區(qū)域帶來了更大的挑戰(zhàn)[4]。此外，水稻種植也影響著耕地利用強(qiáng)度[5]、農(nóng)業(yè)耗水量[6]和溫室氣體排放[7]，進(jìn)而影響到農(nóng)業(yè)生產(chǎn)管理和生態(tài)環(huán)境。因此，及時(shí)掌握水稻種植面積等農(nóng)情信息對地區(qū)農(nóng)業(yè)調(diào)控和環(huán)境問題決策起到支持作用。

如何更加精確地提取水稻空間分布信息，是諸多研究人員關(guān)注的重點(diǎn)問題。目前從時(shí)間維度來看，水稻面積提取的主要分類算法包括單時(shí)相分類、多時(shí)相分類兩大類。主要有決策樹法[8]、面向?qū)ο蠓╗9]、監(jiān)督和非監(jiān)督機(jī)器學(xué)習(xí)法[10]、深度學(xué)習(xí)算法[11]等。上述方法被廣泛應(yīng)用于水稻提取研究中，但主要側(cè)重對單時(shí)相影像中的農(nóng)作物像素單元的提取，對水稻生育期內(nèi)豐富的物候特征信息利用較少。雖然時(shí)間序列指數(shù)分析法關(guān)注水稻關(guān)鍵生育時(shí)期光譜指數(shù)“波峰”“波谷”變化差異，然而受到水稻生長期內(nèi)多云雨天氣因素的限制，很難獲取高質(zhì)量的光學(xué)影像以實(shí)現(xiàn)水稻種植面積的準(zhǔn)確提取[12]。

為充分利用水稻生育期內(nèi)的物候特征信息，弱化因云雨對關(guān)鍵生育時(shí)期影像質(zhì)量的影響，提出水稻像元時(shí)序頻率約束分類模型（TPFCM），對區(qū)域水稻信息進(jìn)行流程化地快速提取。該流程一方面集成時(shí)間序列光學(xué)、雷達(dá)影像數(shù)據(jù)，綜合考慮水稻生育期內(nèi)多個(gè)關(guān)鍵生育時(shí)期的生長信息，利用機(jī)器學(xué)習(xí)算法進(jìn)行水稻的預(yù)分類；另一方面以實(shí)測統(tǒng)計(jì)數(shù)據(jù)為參考，對分類器輸出的多期結(jié)果進(jìn)行空間分布概率統(tǒng)計(jì)約束，以誤差最小為標(biāo)準(zhǔn)輸出最優(yōu)水稻提取結(jié)果，使水稻分布結(jié)果更符合實(shí)際情況。

1 材料和方法

1.1 研究區(qū)與數(shù)據(jù)源

1.1.1 研究區(qū)概況 研究區(qū)位于江西省南昌市（115°27′～116°11′E，28°09′～29°11′N），地處江西省中部偏北，土地總面積7 194.98 km2，其中耕地2 744 km2（占土地總面積的38.1%），有效灌溉面積1 899 km2（占耕地總面積的69.2%），見圖1。氣候?qū)賮啛釒Ъ撅L(fēng)氣候，2021 年平均氣溫19.7 ℃，最高氣溫37.6 ℃，最低氣溫-4.2 ℃。年降水量1 962.8 mm，年日照1 666.6 h，無霜期312 d。全市主要糧食作物為水稻，2021 年水稻種植面積占糧食作物的98.63%，占全部作物種植面積的70.54%，并且水稻熟制復(fù)雜，多為一年兩熟或者一年三熟。水稻按種植和收割時(shí)間的差異分為早稻、中稻和晚稻。通常，早稻每年4 月份左右開始種植，7 月收割（本研究對象早稻后面均稱水稻）。

圖1 研究區(qū)位置Fig.1 Location of the study area

1.1.2 數(shù)據(jù)來源與預(yù)處理

1.1.2.1 遙感數(shù)據(jù) 采用的遙感數(shù)據(jù)分別有Sentinel-1 地距合成孔徑雷達(dá)數(shù)據(jù)（Ground range detected，GRD）、Sentinel-2 多光譜數(shù)據(jù)（Multi spectral instrument，MSI）Level-2 產(chǎn)品和Landsat-8陸地成像儀（Operational land imager，OLI）Level-2產(chǎn)品數(shù)據(jù)。Sentinel-1 和Sentinel-2 數(shù)據(jù)時(shí)間為2021 年4—8 月，與水稻生育期相一致。本研究在Google Earth Engine 中檢索并調(diào)用Sentinel-2、Sentinel-1、Landsat-8 數(shù)據(jù)，然后進(jìn)行云量篩選、月度合成等預(yù)處理。

1.1.2.2 輔助數(shù)據(jù) 為了更準(zhǔn)確有效地提取南昌市水稻分布信息，剔除提取過程的非水稻信息，同時(shí)使用了以下輔助數(shù)據(jù)：SRTMGL1_003 DEM[13]數(shù)據(jù)、ESA WorldCover 10 m 土地利用數(shù)據(jù)（圖1b）[14]、2021年南昌市統(tǒng)計(jì)數(shù)據(jù)。其中，DEM 數(shù)據(jù)用于生成坡度信息，WorldCover 10 m 數(shù)據(jù)用于獲取耕地信息，此外，研究使用的水稻種植面積驗(yàn)證數(shù)據(jù)來源于南昌市統(tǒng)計(jì)年鑒。

1.1.2.3 樣本數(shù)據(jù) 研究所需樣本來源于2 個(gè)途徑：①Google Earth 高清影像底圖；②2021 年GF1 WFV16 m和Sentinel-2 MSI 10 m多光譜影像。通過目視解譯的方式選取研究所需的水稻樣本、非水稻樣本兩類樣本集。通過對高清影像和多光譜影像的解譯判讀，勾繪了88 個(gè)水稻純像元多邊形、41 個(gè)非水稻像元多邊形。其中水稻像元為41 716 個(gè)，非水稻像元66 227 個(gè)（圖1c）。最后，將獲取的高質(zhì)量樣本數(shù)據(jù)集通過生成偽隨機(jī)數(shù)分成訓(xùn)練數(shù)據(jù)集和測試數(shù)據(jù)集，訓(xùn)練數(shù)據(jù)集用于隨機(jī)森林分類器學(xué)習(xí)，測試數(shù)據(jù)集用于驗(yàn)證分類器輸出的預(yù)測精度。

1.2 研究方法

1.2.1 影像特征月度合成 光譜特征是農(nóng)作物遙感分類的核心要素，具有相對的獨(dú)立性和差異性[15-16]。不同的農(nóng)作物因其植株體內(nèi)具有不同的葉片結(jié)構(gòu)、葉片細(xì)胞含水量、色素濃度、氮含量等，使其反射的磁波譜呈現(xiàn)獨(dú)特特征曲線[17-18]。本研究檢索了水稻生育期內(nèi)大量存檔的Sentinel-2、Landsat-8多光譜影像，生成水稻生育期內(nèi)的月度合成影像，并計(jì)算了歸一化植被指數(shù)（Normalized difference vegetation index，NDVI）[19]、增強(qiáng)植被指數(shù)（Enhanced vegetation index，EVI）[20]、地表水分指數(shù)（Land surface water index，LSWI）[21]3種光譜指數(shù)，擬從水稻全生育期分析其光譜特征的時(shí)序差異性，從而更加有效提取水稻分布信息。此外，檢索了水稻生育期內(nèi)的Sentinel-1 存檔數(shù)據(jù)，并通過中值合成的方式將VV、VH 波段合成為月度數(shù)據(jù)，作為Landsat-8、Sentinel-2的補(bǔ)充數(shù)據(jù)。

1.2.2 隨機(jī)森林算法分類 隨機(jī)森林（Random forest，RF）算法[22]作為一種魯棒性較強(qiáng)的集成學(xué)習(xí)算法，已經(jīng)被廣泛用于遙感地物分類提取領(lǐng)域，包括水稻、小麥、玉米等作物提取。RF 算法本質(zhì)是一個(gè)集成了諸多決策樹的分類器，每棵樹對結(jié)果進(jìn)行投票，選擇得票最多的分類結(jié)果作為最終的分類結(jié)果。此外，由于GINI 指數(shù)的簡單性，其也被選擇作為RF單棵樹的分裂標(biāo)準(zhǔn)。RF集合了許多樹作為分類的基礎(chǔ)，單棵樹越獨(dú)立，其泛化誤差越小。本研究調(diào)用GEE 內(nèi)置的RF 算法，參考前人研究結(jié)果設(shè)置了算法參數(shù)值（tree=500），用于水稻生育期內(nèi)多期預(yù)分類工作。

1.2.3 TPFCM 模型 作物關(guān)鍵生育時(shí)期的影像光譜特征對其種植區(qū)域的提取起到積極作用，而且前人的研究工作也證明了這點(diǎn)[23]。但是在關(guān)鍵生育時(shí)期所獲取的影像質(zhì)量因雨云的污染而難以保證，為弱化云雨因素對水稻提取精度的影響，綜合利用水稻整個(gè)生育期的光譜特性，盡可能監(jiān)測到不同生育時(shí)期的水稻像元，本研究基于多時(shí)相的作物識別和水體遙感領(lǐng)域中的水體頻率監(jiān)測思想[24]，提出了水稻TPFCM 模型，并用于研究區(qū)水稻空間分布區(qū)域的優(yōu)化提取。

式中，TPFCM為水稻像元時(shí)序頻率約束模型；Pixelsrice為水稻預(yù)測像元；Pixelsother為非水稻預(yù)測像元；n為水稻生長季內(nèi)的影像監(jiān)測期數(shù)；m為頻率判斷系數(shù)，取值范圍為0～1；w為權(quán)重系數(shù)，其取值范圍為0～1。

1.2.4 精度評估 使用總體精度OA（Overall accuracy）和Kappa系數(shù)評估水稻每個(gè)生育時(shí)期的預(yù)提取面積精度。使用相對誤差RE（Relative error）評估TPFCM 輸出不同頻率閾值下的水稻提取面積與統(tǒng)計(jì)參考面積誤差大小。另外，使用決定系數(shù)R2評估各區(qū)縣水稻提取面積和統(tǒng)計(jì)參考面積的相關(guān)性。

其中，n為作為對比的行政區(qū)數(shù)量，yi為市級不同頻率閾值下的水稻累計(jì)面積，y為對應(yīng)市級的統(tǒng)計(jì)參考面積均值。xi為不同區(qū)縣提取的水稻面積，為對應(yīng)的面積均值，ki為不同區(qū)縣統(tǒng)計(jì)參考水稻面積為對應(yīng)的統(tǒng)計(jì)參考面積均值。

2 結(jié)果與分析

2.1 水稻影像特征月度變化

與其他作物種植管理不同，水稻在不同生育時(shí)期灌溉需求不同，因此其種植環(huán)境也會產(chǎn)生較大的差異性，這種特性在時(shí)序光譜特征上反映得更加明顯。水稻移栽期通常在4 月份，并且為保證水稻移栽存活，水田始終保持存水的狀態(tài)。此時(shí)NDVI、EVI、LSWI 均呈現(xiàn)開始上升。分蘗-拔節(jié)期也是水稻生長的關(guān)鍵時(shí)期。水稻葉片、莖稈發(fā)育迅速，葉片葉綠素含量、氮含量、含水量等不斷增多，進(jìn)而也提高了NDVI、EVI。植株對水分的需求量加大，也使得LSWI呈現(xiàn)上升趨勢（圖2a）。孕穗-抽穗期，水稻葉片、莖稈與其他器官繼續(xù)生長發(fā)育，該階段的NDVI、EVI、LSWI也呈現(xiàn)上升的趨勢，并且達(dá)到特征曲線的最大光譜值。從灌漿期、乳熟期到成熟期，水稻進(jìn)入生殖生長階段，對水分需求量逐漸降低。同時(shí)，水稻葉片等器官的色素濃度減少、光合作用能力下降，水稻籽粒硬化成熟。對應(yīng)的水稻NDVI、EVI、LSWI 指標(biāo)也表現(xiàn)出顯著的下降趨勢。此外，雷達(dá)數(shù)據(jù)時(shí)序VV、VH 特征曲線在水稻各個(gè)生育時(shí)期的變化與光譜特征曲線一致（圖2b），也證明了集成利用光學(xué)、雷達(dá)極化波段監(jiān)測水稻季節(jié)變化的可行性。

圖2 水稻光譜與VV、VH波段月度特征曲線Fig.2 Monthly characteristic curves of rice spectrums and VV and VH bands

2.2 水稻種植區(qū)域月度分布

圖3a—e 分別是4—8 月份水稻全生育期的月度種植面積提取結(jié)果，每期面積提取結(jié)果與水稻關(guān)鍵生育時(shí)期基本吻合。為評估每期水稻面積提取結(jié)果的精度，利用OA、Kappa 系數(shù)指標(biāo)進(jìn)行分析。水稻移栽期、分蘗-拔節(jié)期、孕穗-抽穗期、灌漿期、成熟期OA 分別為0.98、0.98、0.97、0.97、0.97，Kappa系數(shù)分別為0.97、0.97、0.94、0.93、0.92，移栽期和分蘗-拔節(jié)期提取精度較高，孕穗-抽穗期、灌漿期和成熟期提取精度次之（表1）。對于提取面積的RE絕對值而言，5個(gè)關(guān)鍵生育時(shí)期的排序?yàn)橐圃云冢驹兴?抽穗期＞成熟期＞灌漿期＞分蘗-拔節(jié)期（圖3f）。分蘗-拔節(jié)期的誤差最小，這表明，分蘗-拔節(jié)時(shí)期的水稻生長環(huán)境更加穩(wěn)定，衛(wèi)星傳感器更能夠捕捉水體-植株混合的反射光譜、雷達(dá)后向散射特征，從而提高水稻面積提取精度。

表1 水稻移栽期至成熟期種植面積提取的評價(jià)指標(biāo)Tab.1 Classification indicators of planting area from transplanting stage to mature stage

2.3 基于時(shí)序頻率的水稻提取與分析

圖4a—e 展示了TPFCM 模型輸出的研究區(qū)域水稻空間分布頻率。為精確分析不同頻率間隔的水稻面積提取誤差，將設(shè)置頻率間隔步長為0.2，分析了0.2～1.0 頻率間隔所對應(yīng)的水稻總種植面積。在水稻移栽期至成熟期的時(shí)序像元頻率約束結(jié)果中，頻率為1.0、0.8、0.6、0.4、0.2 的像元占比分別為27.3%、23.7%、20.6%、16.9%、11.4%。頻率值大于等于1.0、0.8、0.6、0.4、0.2 的面積提取RE 分別為67.96%、40.17%、16.01%、-3.83%、-17.19%（圖4f）。根據(jù)上述不同頻率間隔的水稻面積提取RE，分析發(fā)現(xiàn)，頻率大于等于0.4 的水稻面積提取RE 絕對值最?。?.83%），這也表明0.4 為本研究最佳頻率閾值。且與單一生育時(shí)期的最小RE 絕對值（7.77%）相比，該頻率值所對應(yīng)的水稻面積提取RE 絕對值減少了3.94 百分點(diǎn)，有效降低了面積提取產(chǎn)生的誤差。同時(shí)也說明，盡管基于關(guān)鍵生育時(shí)期內(nèi)的單景影像能提取水稻的種植區(qū)域，但綜合水稻全生育期的水稻頻率約束能進(jìn)一步提高水稻面積提取效果。

圖4 水稻像元頻率空間分布Fig.4 Spatial distribution of rice pixel frequencies

圖5a 展示了提取的南昌市9 個(gè)縣級行政區(qū)的水稻面積空間分布與量級大小。其中，水稻種植區(qū)域面積相對較大的行政區(qū)有4 個(gè)，分別為南昌縣（688.47 km2）、新建區(qū)（429.01 km2）、進(jìn)賢縣（346.37 km2）和安義縣（105.88 km2）。這4 個(gè)縣水稻種植面積占總面積的99.50%。水稻種植面積較小的縣區(qū)分別為青山湖區(qū)（3.67 km2）、灣里管理局（3.14 km2）。由于上述地點(diǎn)主要位于城區(qū)，土地利用方式以建筑用地為主，耕地面積較少，因此水稻提取面積偏低。圖5b 展示了各行政區(qū)水稻種植統(tǒng)計(jì)面積與模型提取面積的線性相關(guān)性分析結(jié)果，二者存在顯著的相關(guān)性（R2=0.97），也進(jìn)一步證明了TPFCM 模型提取水稻的準(zhǔn)確性。

圖5 水稻分區(qū)面積統(tǒng)計(jì)與相關(guān)性Fig.5 Rice extraction area statistics and correlation

2.4 與基于RF分類器單一生育時(shí)期水稻提取效果比較

為進(jìn)一步分析TPFCM 模型提取水稻面積的效果，與基于RF 分類器提取的單一生育時(shí)期的水稻空間分布數(shù)據(jù)進(jìn)行對比分析（圖6）。挑選了平原和河谷兩類地形的水稻種植子區(qū)域，以便突出比較RF單一生育時(shí)期水稻面積提取結(jié)果與TPFCM 模型的差異性。結(jié)合GF1 假彩色合成影像可以發(fā)現(xiàn)，在平原區(qū)域（圖6a—b），RF 單一生育時(shí)期的水稻空間分布和TPFCM 模型大致相同（圖6e—h，圖6i—l），但是在細(xì)節(jié)部分存在較大差異?？梢悦黠@看出，單一生育時(shí)期內(nèi)RF 算法的水稻提取結(jié)果在田間小路、撂荒地等區(qū)域造成了誤分（圖6e—f）。雖然單一生育時(shí)期的影像集成了光譜波段和光譜指數(shù)特征，但是在水稻處于分蘗-拔節(jié)期間，RF 分類器難以準(zhǔn)確區(qū)分與水稻光譜相似的田間小路和撂荒地自然生長的雜草等像元，造成同譜異物現(xiàn)象。在河谷地區(qū)（圖6c—d），RF 能相對準(zhǔn)確提取河谷區(qū)域大面積水稻種植區(qū)域，但并未能有效提取河谷區(qū)域小面積的水稻田地塊（圖6g—h）。而TPFCM 模型通過對水稻全生育期的水稻像元進(jìn)行頻率約束控制，確定了水稻種植區(qū)域，更準(zhǔn)確地提取了河谷地區(qū)小區(qū)域水稻田分布空間（圖6k—l）。

圖6 RF單一生育時(shí)期水稻面積提取結(jié)果與TPFCM模型水稻面積提取結(jié)果比較Fig.6 Comparison of RF single-phase rice extraction results and TPFCM model rice extraction results

圖6i—l 為TPFCM 模型提取水稻種植區(qū)域，與單一生育時(shí)期RF 算法提取的水稻區(qū)域相比，TPFCM 模型更能從影像時(shí)序光譜特征上捕獲水稻信號并確認(rèn)為水稻像元，從而有效提取研究區(qū)內(nèi)的水稻空間分布。該模型在平原區(qū)域能準(zhǔn)確區(qū)分水稻和非水稻像元，同時(shí)，也能精細(xì)地監(jiān)測河谷小區(qū)域水稻種植空間分布。

3 結(jié)論與討論

在本研究中，得益于研究區(qū)域內(nèi)大量存檔的多源遙感數(shù)據(jù)的支持，TPFCM 發(fā)揮出較強(qiáng)的分類能力，準(zhǔn)確地提取了南昌市水稻種植區(qū)域。但該模型的提取步驟依然存在著一定的不確定性：（1）水稻全生育期內(nèi)多源影像質(zhì)量。為盡可能多地獲取研究區(qū)內(nèi)全覆蓋無云光學(xué)影像，將多期云掩膜后的Landsat-8、Sentinel-2 進(jìn)行了影像合成，生成了多期全覆蓋的影像數(shù)據(jù)，并計(jì)算了NDVI、EVI和LSWI光譜指數(shù)作為分類特征輸入分類器。由于Landsat-8和Sentinel-2 可見光區(qū)域和短波紅外區(qū)域的波長范圍并非完全一致，2 種星源對應(yīng)波段存在細(xì)微的光譜反射差異，將二者合成影像及其光譜指數(shù)輸入分類器，一定程度上影響了分類效果。此外，對于大尺度的水稻分類區(qū)域，受制于云雨的影響，Landsat-8、Sentinel-2 光學(xué)影像會形成時(shí)間連續(xù)、大面積的“污染”像元，盡管補(bǔ)充了雷達(dá)數(shù)據(jù)作為分類特征，也會影響到水稻的準(zhǔn)確提取。（2）水稻預(yù)提取面積結(jié)果。由于像元時(shí)序頻率約束并不以原始影像為分類對象，而是以水稻關(guān)鍵生育時(shí)期預(yù)提取面積結(jié)果為輸入對象，因此，預(yù)提取結(jié)果的數(shù)量和質(zhì)量直接影響到最佳的水稻分布區(qū)域。像元時(shí)序頻率約束方法通過類似“投票”的方式，評估所有參與提取的影像在同一像元位置是否為水稻像元，并根據(jù)每期水稻提取影像在該像元出現(xiàn)的次數(shù)賦予最終的頻率值。通常水稻生育期內(nèi)預(yù)提取結(jié)果期數(shù)越多，像元時(shí)序頻率約束輸出水稻頻率分布越準(zhǔn)確。

本研究借助Google Earth Engine 云服務(wù)平臺調(diào)用了水稻生育期內(nèi)存檔的Landsat-8、Sentinel-2、Sentinel-1 等多源遙感數(shù)據(jù)，并生成了逐月的特征矩陣，通過TPFCM 準(zhǔn)確提取了研究區(qū)內(nèi)水稻空間分布區(qū)域，RE 絕對值為3.83%，與基于RF 分類器單一生育時(shí)期的水稻提取面積最小RE 絕對值相比減少了3.94百分點(diǎn)。通過研究分析，TPFCM 能在一定程度上通過多期頻率計(jì)算的方式降低云雨和錯(cuò)分的影響，提高水稻提取效果?？傊?，TPFCM 有助于快速、有效提取以南昌市為代表的南方水稻空間分布，為后續(xù)進(jìn)一步的研究提供高質(zhì)量的水稻空間分布數(shù)據(jù)支持。