基于Sentinel-2影像與MaxEnt模型識別云南華坪縣芒果種植區

趙管樂 ， 劉勤 *， 彭培好

（1.成都理工大學地球科學學院，成都 610059；2.中國科學院、水利部成都山地災害與環境研究所，成都 610041；3.成都理工大學生態資源與景觀研究所，成都 610059）

云南省麗江市華坪縣自20世紀60年代開始種植芒果，是我國高緯度芒果產地之一。華坪縣生產的晚熟芒果憑借其優異的品質為當地農業、旅游業帶來了的巨大經濟效益，是助農脫貧的重要產品，華坪芒果也因此成為國家地理標志產品。華坪縣自20世紀90年代開始規?；N植后芒果種植面積逐年擴大，至2020年全縣種植面積達到2.5萬hm2。隨著芒果種植面積的不斷擴大，果園人工管理成本也在不斷增加。因此，利用現代化的信息技術進行果園管理成為芒果種植業發展的必然趨勢。遙感技術作為一種新興的信息獲取手段，因其具有高時效性、大面積同步觀測等特點被廣泛應用于農業生產中。利用遙感技術準確地識別芒果種植區域為大范圍地監測果樹長勢、病蟲害防治、估產等提供了可能。

國內外許多專家學者利用遙感數據對農林作物進行了分類。任傳帥等[1]利用單期SPOT影像通過面向對象分類方法，結合植被覆蓋度與坡度因子對海南省三亞市某地芒果林地進行了提取，且分類效果十分理想。李章程等[2]利用多時相高分一號衛星通過決策樹分類法對四川德陽地區水稻種植區進行了識別，制圖精度達到96.5%。楊其坡[3]利用多時相Sentinel-2影像構建多種植被指數結合多種地形因子通過MaxEnt模型對重慶市江津區花椒種植區域進行提取，得到了較好的分類效果。Kobayashi等[4]利用單期Sentinel-2A影像構建了若干光譜指數，并利用隨機森林算法對日本北海道某地的主要農作物進行了分類，最終基于4個波段反射率和8個光譜指數的分類精度達到了93.1%。以上研究區域多集中在平原地區，對山地丘陵地區的農作物遙感分類研究還較少，其原因在于山區地貌復雜、地物類型破碎、山體具有陰影等因素極大制約了遙感影像對作物識別的精度。此外，部分研究使用的高分辨率遙感影像雖然能極為準確地對農作物進行分類，但并不能免費獲取要進行大范圍、多時相的農作物，識別與監測需要高昂的成本。

Sentinel-2系列衛星作為歐空局“哥白尼計劃”的重要組成部分，其影像具有3個紅邊波段，為遙感技術在植被與農作物等方面的應用帶來了巨大的應用潛力。本研究基于多時相Sentinel-2遙感數據對華坪縣已投產的成熟芒果種植區域進行識別，規避了使用單一時相遙感數據對地物進行識別分類時由“同譜異物”等現象所帶來的弊端，充分考慮了植被的光譜特征隨季節變化而變化的特點，此外，在利用傳統植被指數的前提下還引入了紅邊波段數據對芒果種植區進行識別分類。利用主成分分析將Sentinel-2數據3個紅邊波段的信息進行綜合，使其簡便高效地參與模型預測。

MaxEnt模型也稱為最大熵模型，是一種基于最大熵理論的物種空間分布生態位模型[5]。該模型通過已知的物種分布區域或明確的無分布區域，在給定的環境變量中尋找限制或促進物種分布的變量[6]，從而建立目標物種空間分布情況與環境變量之間的關系，以此關系作為依據對未知區域目標物種的潛在分布情況做出預測。本研究通過MaxEnt模型預測華坪縣已掛果投產的芒果種植區分布，并按不同閾值進行二值化分割，確定該研究條件下合理的分割閾值，以期為研究方法的應用和推廣提供科學參考。

1 材料與方法

1.1 研究區概況

華坪縣隸屬于云南省麗江市，位于云南省西北部的金沙江中段北岸，全縣總面積約2 200 km2，是云南省最大的芒果生產縣。華坪縣為典型的南亞熱帶低熱河谷氣候，無霜期長達300 d以上，光熱資源豐富，降水集中，年平均溫度19.8℃，年平均降水量870 mm。全縣地勢高低呈西北-東南走向，平均海拔1 160 m，最高海拔3 198 m，最低海拔1 015 m。

1.2 數據來源與處理

1.2.1 地形數據 本研究使用ALOS衛星(https://earthdata.nasa.gov/)所采集的數字高程模型（digital elevation model，DEM）數據（圖1），分辨率12.5 m。對獲取的多幅DEM數據進行鑲嵌后再重采樣至10 m分辨率，最后用華坪縣行政區劃矢量文件進行裁剪得到全縣DEM數據。

圖1 華坪縣地理位置與高程信息Fig.1 Geographical location and elevation of Huapin county

1.2.2 遙感影像數據 本研究選取的Sentinel-2影像數據來源于歐空局(https://scihub.copernicus.eu/)。根據芒果生長周期的典型節點以及影像云量分布等條件，分別選取了2019年12月16日、2020年3月15日、2020年6月8日、2020年8月17日共4個時相的影像。

所獲取的影像為Level-1C級別的產品，即只經過輻射定標與幾何校正的產品，需要使用歐空局提供的Sen2cor模塊進行大氣校正。通過大氣校正后得2A級別的產品，然后對數據統一重采樣到10 m分辨率，再利用華坪縣行政區劃邊界矢量文件對影像進行裁剪得到全縣的各期影像。

1.2.3 芒果種植區調查點數據 本研究于2020年3月對華坪縣已投產芒果種植區域進行野外調查。借助天地圖高分影像和手持GPS接收機選擇果樹葉片密集、冠幅1 m以上、已進入掛果投產階段的芒果種植園，并在園地內獲取經緯度坐標。此外還獲取了其他非芒果種植區地類的坐標信息，用于對后續分類精度的評價?？偣搏@取芒果種植區域樣點482個，非芒果種植區域樣點1 667個（圖2），其中非芒果種植區域的地類主要包括農田、園地、林地、裸地、建設用地等。

圖2 樣本點分布Fig.2 Sample distribution

1.3 MaxEnt模型

本研究利用Phillips等[7]基于最大熵模型通過Java語言所開發的MaxEnt軟件，將物種分布數據與環境變量數據輸入該軟件后即可得到可視化的物種分布預測結果、不同環境變量對物種分布的貢獻率以及相應的模型性能評價等信息。輸入的環境變量數據為*.asc格式，物種分布點位數據為*.csv格式。隨機選擇芒果種植區樣本點中的75%（362個）作為訓練樣本，25%（120個）作為測試樣本，其他參數選擇默認設置。

最大熵模型采用受試者特征曲線（receiver operating characteristic curve，ROC）對預測結果精度進行評價。ROC曲線以假陽性率（1-特異性）為橫坐標，真陽性率（敏感度）為縱坐標。ROC曲線下與坐標軸所圍成的面積（area under curve，AUC）的大小即作為評價模型精度的指標。AUC值處于0.5～0.6之間表示精度極差；處于0.6～0.7表示精度較差；處于0.7～0.8表示精度一般；處于0.8～0.9表示精度較好；處于 0.9～1.0表示精度極好。

Maxent軟件運行后的結果包含1個與輸入環境變量分辨率相同的圖像，該圖像每個像元值介于0～1之間，代表目標物種在此處存在的概率。該概率表示通過環境變量與已知物種分布點進行預測得到的未知區域該物種存在的可能性。為了更好地表示模型預測結果，將概率發生圖通過合理的閾值進行分割轉換為二值圖像，即1個像元只表示物種存在或不存在2種情況。閾值的選擇對于物種概率發生預測圖的二值化結果影響極大，但目前沒有一種全自動的方法能確定合理的閾值[8]。MaxEnt軟件運行結果中給出了一些常用閾值規則下的對應閾值，本研究基于部分給出的閾值對概率發生圖進行二值化，并利用芒果種植區分布點與其他非芒果種植區地物分布點對不同閾值分割下得到的二值化分類圖像進行分類精度評價，將分類精度最高的二值圖像所對應的閾值作為最佳閾值。

1.4 植被指數與地形因子的提取

1.4.1 植被指數的提取 植被指數是指將遙感影像的不同波段按照線性或非線性組合來表現植被光譜特征的方法，目前其被廣泛應用于植被生長情況的定性或定量監測，是植被與生態遙感領域的重要研究手段[9]。大多數植被指數是利用遙感影像的紅光與近紅外波段組合得到的，也有一些利用綠光波段，但隨著遙感技術與傳感器的不斷發展，發現紅邊波段對于植被遙感具有廣闊的應用前景。本研究所使用的Sentinel-2影像具有3個紅邊波段，這為植被遙感研究帶來巨大的應用價值。

結合華坪縣芒果的生長情況，本研究選取6種植被指數以及Sentinel-2影像的3個紅邊波段進行主成分分析得到第1主成分分量，并分別計算了上述4期影像的7個指數。

①歸一化植被指數（normalized difference vegetation index，NDVI）。 NDVI是使用最為廣泛的植被指數，其能在一定程度消除多種因素所導致輻射誤差所帶來的影響，較好地反映植被覆蓋度與生長情況等[10]，計算公式如下。

式中，NDVI為取值-1.0～1.0之間的無量綱；ρNIR為近紅外波段反射率，對應原始影像數據第8波段，取值范圍在0～100%之間；ρR為紅光波段反射率，對應原始影像數據第4波段，取值范圍在0～100%之間。

②差值植被指數（difference vegetation index，DVI）。該指數對土壤背景的變化十分敏感，但當植被覆蓋度較高時其對植被的敏感度就會下降[11]，因此，該指數常用于中低覆度區域的植被監測。

③比值植被指數（ratio vegetation index，RVI）。與差值植被指數相反，RVI在植被覆蓋度較高時敏感度更高，而當植被覆蓋度降低時敏感度也相應地降低[12]。該指數受大氣條件的影響較為明顯，因此必需進行大氣校正后才能進行計算，取值大于0。

④土壤調節植被指數（soil adjusted vegetation index，SAVI）。SAVI引入了土壤調節系數來減弱土壤背景帶來的干擾，能夠從一定程度上描述土壤與植被之間的關系[13]。

式中，L為土壤調節系數，取值范圍在0～1.0之間，本研究中L取值為0.5。

⑤權重差值植被指數（weighted difference vegetation index，WDVI）。WDVI是在差值植被指數的基礎上改進得到的[14]，其為了減弱土壤背景帶來的影響引入了土壤調節參數，計算公式如下。

⑥色素特異性簡單比（pigment specific simple ratio，PSSRa）。 PSSRa是針對植物葉片色素含量的植被指數，能體現出植被冠層光合色素含量的差異[15]。

式中，ρRedEdge3為第3個紅邊波段反射率，對應原始影像數據的第7波段，取值范圍在0～100%之間。

⑦紅邊波段。其介于紅光波段和近紅外波段之間，植被在此波段范圍內的光譜反射率會出現顯著的突變，植被的種類、生長情況、環境脅迫等在此波段范圍內具有豐富的信息[16]。因此本研究選擇3個紅邊波段的信息參與芒果種植區的識別。但為了減少數據的冗余，最大程度地綜合3個紅邊波段的信息，利用主成分分析PCA（principal component analysis）[17]對 數 據 進 行 降維，第1主成分（PCA1）貢獻率大于85%，說明PCA1分量包含了3個紅邊波段的絕大部分信息，可以代表3個紅邊波段參與芒果種植區的識別。

1.4.2 地形因子的提取 除了DEM數據外，還利用重采樣后10 m分辨率DEM數據生成坡度（slope）與坡向（aspect）2個常用地形因子。由于農業生產中地塊的坡度與坡向影響該區域耕種作物的類型，因此，選擇這2個地形因素有利于輔助遙感數據對地物進行分類識別，提高分類精度。

1.4.3 模型評價 利用選取的4期Sentinel-2影像分別計算出7種植被指數，共計28個；并通過DEM數據計算出坡度與坡向，共計3個地形因子。將上述環境變量輸入MaxEnt軟件，得到基于多時相Sentinel-2影像植被指數與地形因子的芒果種植區ROC曲線與AUC值。

2 結果與分析

2.1 模型性能分析

從圖3可以看出，訓練數據的AUC值為0.962，測試數據的AUC值為0.955，兩項AUC值都接近于1，表明模型預測的精度非常好，預測結果可信度高。

圖3 基于環境變量的芒果種植區預測結果Fig.3 Prediction result of mango planting regions based on environmental variables

2.2 主導因子確定

31個環境變量對芒果種植區分布預測模型的相對貢獻率如表1所示。

表1 環境變量對芒果種植區預測模型的貢獻率Table1 Contribution rate of environmental variables to the prediction model of mango planting regions

本研究選取累積貢獻率達到95%的環境變量作為影響芒果種植區預測準確程度的主導因子，共選擇11個環境變量作為主導因子，累積貢獻率為95.2%，按照相對貢獻率的大小分別為：DEM（27.8%）、20191216DVI（15%）、20191216WDVI（10.6%）、20200817PCA1（9.8%）、20200608WDVI（9%）、20200608SAVI（6.2%）、20200315DVI（5.6%）、aspect（4%）、20200817WDVI（3.3%）、20200608PCA1（2.6%）、20200315PSSRa（1.3%）。

2.3 不同閾值分割下分類精度對比

利用不同閾值分別對結果中的概率發生圖進行二值化，并通過混淆矩陣得到總體分類精度與Kappa系數，如表2所示?？梢钥闯?，所選取的7種閾值規則對芒果種植區進行二值化分類的分類精度都達到了很高的水平，其中，只有EETOD規則的分類精度低于90%為89.57%，其他6種規則的分類精度均在90%以上，最高的PTP規則分類精度最高達到了93.72%，而7種閾值規則下分類結果Kappa系數均在0.70以上。ETESS、ETRSS以及PTP 3種閾值規則下分類結果的Kappa系數高達0.81及以上，PTP規則下的Kappa系數最高為0.82，處于幾乎完全一致性水平。因此，結合總體分類精度與Kappa系數最終選擇PTP為本研究的最佳閾值規則，對應的閾值為0.257。

表2 不同閾值規則下芒果種植區分類精度與Kappa系數Table 2 Classification accuracy and Kappa coefficient of mango planting regions under different threshold rules

2.4 華坪縣芒果種植區預測

根據確定的最佳閾值（0.257）對華坪縣已投產芒果種植區預測分布圖進行二值化分類（圖4）。圖中芒果種植區所占像元數為1 067 121個，每個像元所代表的實際面積為100 m2，因此已投產芒果種植區預測面積折合約1.07萬hm2。據有關報道，2019—2020年華坪縣全縣已掛果投產的芒果種植面積約1.13萬～1.20萬hm2。

圖4 華坪縣已投產芒果種植區分布預測Fig.4 Predicted distribution of mango planting regions putting into production in Huaping country

將最終二值化分類的結果調整到一定透明度后疊加在天地圖高分影像上，并選取了3處具有代表性的區域進行展示，如圖5所示。圖中紅線為預測出已投產芒果種植區矢量邊界，邊界內部淡紅色的區域為疊加在高分影像上預測結果中的已投產芒果種植區，高分影像上芒果種植區呈現明顯的均勻點狀分布。

圖5 華坪縣已投產芒果種植區分布預測結果局部展示Fig.5 Partial display of predicted distribution results of mango planting regions in Huaping county

預測的芒果種植區范圍與高分影像中實際的芒果種植區范圍基本吻合，且能夠較好地與散落在芒果種植區內的耕地區分開。但由于本研究所使用的數據空間分辨率僅為10 m，因此預測結果的邊界不能很好地與面積較小的地物邊界擬合，呈現明顯的鋸齒狀。預測的芒果種植區也能與周邊的其他林地區分開。零星分布的房屋、蓄水池等地物能夠與大面積的芒果種植區區分開，連續的芒果種植區能夠被很好地預測識別出，但與其他類型地物相鄰很近的芒果種植區存在被錯分的現象。由于遙感數據的空間分辨率有限，在不同地物類型交界處單個像元的地物光譜信息混合嚴重，因而在這類區域利用光譜信息識別物體時存在明顯不足。

綜上可以看出，本研究所預測出的華坪縣已投產芒果種植區種植面積和空間分布吻合度都到達了較好的精度水平。因此，本研究所采取的方法與最終得到的結果對芒果種植區的空間規劃與管理具有一定的參考價值。

3 討論

本研究選擇樹齡成熟的芒果種植園為研究對象是由于華坪縣芒果種植均為定植嫁接苗，苗木出圃時通常很矮小、葉片稀少、葉面積指數小，通常需要至少3年以上才能發育為成熟植株掛果投產。而剛定植的芒果幼苗由于葉片稀少，在光學遙感影像上幾乎體現不出植被所具有的光譜特征，很難被識別出來，且新植的芒果樹幼苗的園間管理成本較低、病蟲害發生較少。

本研究在利用遙感數據對已投產芒果種植區進行識別時引入了地形因子作為輔助。這是因為調研發現，華坪縣是緯度較高的芒果產地之一，在選擇芒果種植園時海拔是重要的指標，全縣海拔1 015～1 500 m處均可商品化種植生產芒果[18]，因此引入了高程、坡向等地形因子，這與楊其坡[3]在利用Sentinel-2數據和MaxEnt模型對重慶市江津區花椒種植區的識別所采用的方法類似。本研究結果（表1）表明，DEM對模型預測的貢獻率居于首位，達到27.8%，這一結果與客觀實際相符，即海拔影響華坪縣芒果種植區的分布。

此外，付發群[19]和王克曉等[20]在利用Sentinel-2影像對浙江省金華市城市綠地與重慶市合川區的山地油菜分布進行提取時，都引入了Sentinel-1雷達數據，其研究結果均表明，使用雷達數據參與分類的方案均比不使用雷達數據的分類方案精度高。因此，在數據選擇方面還可以進一步引入雷達遙感數據。

本研究采用MaxEnt模型對已投產芒果種植區進行分布預測，方法操作簡便、效果較好，但分類精度還存在提高的空間，且二值化分類閾值的選擇決定著最終結果的準確性，同時也因閾值選擇的不確定性而帶來誤差。趙娟[21]利用國外高分影像數據與深度卷積神經網絡模型基于FCN語義分割框架對四川省攀枝花市仁和區芒果種植區進行了精確的識別，且該分類識別方法的最終結果不依賴于外部閾值的選擇，同時該方法還能夠識別出處于各生長階段的芒果種植區，這優于本研究只能識別出已投產的芒果種植區。因此，在后續的研究中可以繼續探索將不同分類方法相結合，優勢互補，以期獲得更理想的分類效果。