GEE環境下融合主被動遙感數據的冬小麥識別技術

馬戰林 劉昌華 薛華柱 李靜茹 房 旭 周俊利

(1.河南理工大學測繪與國土信息工程學院， 焦作 454003； 2.河南省遙感測繪院， 鄭州 450003)

0 引言

隨著我國糧食生產的市場化，冬小麥的種植面積呈現年際變化大的特點，及時獲取冬小麥精確的種植面積信息及空間分布，對促進我國農業生產發展和糧食安全具有重要意義[1]。

遙感觀測技術在大范圍作物種植面積和結構的監測上已被廣泛應用[2]。但由于技術條件限制和工作原理的不同，單一的傳感器難以全面反映復雜的地表覆被特征。光學遙感數據受光譜分辨率、空間分辨率、光譜波段的因素及限制，導致同物異譜、異物同譜的現象存在，影響了地物識別的精度。已有學者通過不同方法提高地物分類精度[3-9]。上述方法需要足夠的光學數據支持[10]。但在實際應用中，受云雨、光照等因素的影響，光學遙感數據源的質量和數量無法保障，一定程度上限制了地面作物信息的準確提取[11]。

合成孔徑雷達(Synthetic aperture radar,SAR)不受云雨天氣影響，具有全天時、全天候監測等優點。然而，SAR數據的信號易受相干斑噪聲的干擾，影響對目標地物的提取，單一時相的SAR數據很難對地物進行精確提取。有研究表明，使用多時相或多極化SAR數據比單時相或單極化SAR數據能獲得更好的分類結果[12]。若能將光學和微波遙感數據結合，利用兩者數據的優點進行地物分類，則分類精度會有較大提升。文獻[13-16]利用最大似然分類(Maximum likelihood classification,MLC)、人工神經網絡(Artificial neural network, ANN)、支持向量機(Support vector machine,SVM)和監督分類等方法，對比光學遙感、微波遙感、光學遙感與微波遙感數據結合3種方式的作物分類精度，得到光學和微波遙感數據相結合的分類精度最高。

相比于單期光學數據和SAR數據的融合，時序數據進行融合應用時，作物分類精度進一步提高。文獻[17]結合Landsat-8和Sentinel-1A數據構建時間序列遙感數據以提取作物物候特征，利用卷積神經網絡算法進行不同作物的分類，并取得了較高的分類精度。文獻[18]以多時相SAR(Sentinel-1A)和Landsat-8光學影像對南京市高淳區的冬小麥種植面積進行提取，結果表明添加紋理特征的時序SAR數據和光學數據結合時分類精度最高。

除數據源外，分類算法的選擇也會直接影響分類結果[19]。近年來隨著機器學習分類算法的發展，隨機森林算法(RF)計算精度高，模型訓練時間少，能夠確定變量在模型中的相對重要性，同時對訓練樣本數量和質量的敏感度低，而被廣泛應用于作物識別中[20-21]。

以上研究中均需將所需遙感數據下載至本地進行預處理，所需時間較長，而谷歌地球引擎(Google Earth Engine，GEE)是一個提供全球尺度地球觀測數據儲存和用戶友好界面的數據訪問平臺，具有強大的數據處理能力，僅需內置的少數代碼即可完成數據獲取、預處理及算法的應用。因此，本文在GEE平臺上，充分利用SAR極化數據所包含的地物結構信息和光學數據的光譜信息，結合RF算法，探究時間序列SAR數據、融合時間序列SAR和光學數據的不同特征值組合對冬小麥識別精度的影響。

1 研究區域與數據來源

1.1 研究區概況

河南省襄城縣、鄢陵縣和臨潁縣位于河南省中部地區(圖1)，是中國北方典型的冬小麥種植區和河南省高標準冬小麥種植示范區。三縣地跨113°2′～114°2′E，33°41′～34°15′N，總面積為2 612.60 km2，屬北溫帶季風氣候區，熱量資源豐富，雨量豐沛，光照充足；年平均氣溫14.3～14.7℃，年平均降水量為671.1～736.0 mm，無霜期為217.5 d；地形地貌以平原為主，占75.81%，土壤肥沃，具有較高的自然肥力；農作物以冬小麥、夏玉米為主，其中，小麥播種面積高達夏糧總種植面積的90%以上。該區域的冬小麥10月中旬種植，到次年6月上旬收割。鄢陵縣種植結構復雜，地塊細碎，與襄城縣和臨潁縣的連片種植冬小麥形成鮮明對比，在分類精度的評價上具有較好的區域對比性，可驗證融合Sentinel主被動遙感數據在提取冬小麥方面的應用潛力。

1.2 GEE及應用數據

GEE云平臺存儲了Sentinel數據、MODIS數據集、降水數據、The National Agriculture Imagery Program(NAIP)數據、海洋表面溫度數據、Landsat數據、CHIRPS(Climate hazards group infraRed precipitation with station data)氣候數據和海拔數據等海量數據，容量超過5PB，且每天增加約5 000幅影像，可以解決大面積土地覆蓋制圖方面最重要的數據存儲下載問題[22]。用戶可以使用基于網絡的集成開發環境代碼編輯器分析所有可用的遙感圖像，而無需將這些數據下載到本地機器上。這樣，用戶可以輕松訪問、選擇和處理大研究區域的大量數據。除了快速處理之外，它提供了幾個包含大量算法的軟件包，可以簡化專家和非專家對遙感工具的訪問。GEE也允許用戶上傳自己的柵格和矢量數據(例如GeoTIFF或Shape文件)進行分析，完全控制訪問[23]。

Sentinel-1主動微波遙感衛星由兩顆極軌衛星A星和B星組成，搭載C波段的合成孔徑雷達(SAR)傳感器，重訪周期小于10 d，本文采用分辨率為10 m、極化方式為VV和VH的后向散射系數數據。Sentinel-2由Sentinel-2A和Sentinel-2B兩顆高分辨率衛星組成，單顆衛星的重訪周期為10 d，兩顆互補，重訪周期為5 d。其搭載的MSI傳感器有可見光、近紅外和短波紅外波段，空間分辨率在不同的波段有10、20、60 m。本研究主要選取該衛星分辨率10 m的藍光、綠光、紅光和近紅外波段。選用的Sentinel-1微波遙感數據和Sentinel-2 MSI多光譜遙感數據均在GEE云平臺上進行調用、處理。在GEE云平臺上使用函數直接獲取不同極化方式的后向散射系數，減輕處理SAR數據的難度。根據圖像的大氣校正狀態，Sentinel-2數據可分為兩個不同的級別:1C級和2A級。1C級為大氣頂部圖像，需要大氣校正。2A級數據為大氣底部反射數據，已進行大氣校正。選用2A級數據。同時考慮云霧對數據質量的影響，在GEE中設置云量低于10%。遙感影像參數如表1所示。

表1 遙感影像參數Tab.1 Remote sensing image parameters

本文使用高分六號(GF-6)數據，該衛星配備2 m全色/8 m型多光譜高分辨率相機以及16 m型多光譜中分辨率寬幅相機。應用2020年4月6日和2020年5月17日兩景GF-6融合后的2 m多光譜數據，分析評價分辨率10 m Sentinel主被動數據在碎部地區的分類效果。

1.3 特征變量與數據集

結合研究區植被的生態環境和物候特點，選取Sentinel-2影像的光譜反射率及相關植被指數和熵值紋理(Entropy,ENT)共11個特征變量。Sentinel-1 SAR數據的特征變量包含極化特征變量和紋理特征變量，為避免多紋理特征帶來的特征信息冗余[24]，在保留最大信息量、便于計算和較高分類精度的前提下，選取灰度共生矩陣(Gray-level co-occurrence matrix,GLCM)生成的角二矩陣(Angular second moment,ASM)、對比度(Contrast,CONTRAST)、相關性(Correlation,CORR)和熵值4個紋理特征變量，共6個特征變量[25]。研究中共17個特征變量，如表2所示。

表2 特征變量Tab.2 Characteristic variable

歸一化差異植被指數(Normalized difference vegetation index,NDVI)研究較高吸收和葉綠素反射的波段來識別光合作用植被。歸一化水體指數(Normalized difference water index,NDWI)用來識別水體，NDVI和NDWI的結合可用于提取一些沿河畔、坑塘生長的植物的信息[26]。增強型植被指數(Enhanced vegetation index,EVI)包括近紅外、紅色和藍色波段[27]，可在高生物量區域獲得更好的靈敏度，并通過解耦冠層背景對不同物種的冠層結構變化做出更好的響應[28]。綠色歸一化差異植被指數(Green normalized difference vegetation index,GNDVI)已被開發用于估計葉片葉綠素含量[29]。綠葉指數(Green leaf index,GLI)為顏色植被指數，反映植被的生長至衰老顏色變化信息。優化土壤調節植被指數(Optimization soil-adjust vegetation index,OSAVI)能夠反映冬小麥的生長狀況信息[30]。

受制于研究區作物種植結構和田塊零碎的影響，研究中未使用Sentinel-2數據中分辨率為20 m的紅邊波段。紋理特征的提取能夠較好地描述紋理的細節和隨機性，適于地物分布復雜的SAR影像[31]。

為獲取襄城縣、臨潁縣、鄢陵縣主要地物類型分布，結合Google Earth Engine，在2020年3月15日到2020年5月12日冬小麥生長期間進行實地采樣。依據研究區的農田、道路、房屋、水體占土地覆被總面積90%左右，將覆被類型分為建筑用地、冬小麥、水體、林地、其他植被、其他用地6大類別。建筑用地包括房屋道路等建設用地，其他植被包含蔬菜、花卉、瓜果等，其他用地包括裸土、裸巖以及采石場等。應用中海達V30 型GPS儀器進行實地采樣，共獲得1 051個樣本點(表3)，其中田地采樣點均在大面積(邊長100 m以上)田塊中心采集。驗證數據集應用GPS在研究區內隨機均勻地對不同類別的數據進行區域采集，生成矢量數據，應用矢量范圍數據對遙感圖像進行裁剪得到驗證數據集像元數量。

表3 樣本數據集數量Tab.3 Sample data sets

2 研究方法

2.1 隨機森林算法

隨機森林算法的預測結果是在對組成森林的多棵決策樹的決策結果進行眾數求解得出的。每次采用Bootstrp采樣技術從原始的訓練數據集中有放回地抽取大約2/3的樣本用于對當前決策樹模型進行訓練，故原始訓練集中的部分樣本可能被同時用來對決策樹進行訓練，也有可能從未被任意一棵決策樹使用。抽取剩余的約1/3的樣本構成對隨機森林模型進行性能評估的袋外數據，計算該模型的預測錯誤率，該錯誤率被稱為袋外誤差。隨機森林分類的基本思想：使用Bootstrap抽樣從原始訓練集中提取k個樣本，每個樣本的容量與原始訓練集相同；為k個樣本建立k個決策樹模型，以獲得k個分類結果；對每個記錄進行表決以確定其最終分類。雖然隨機森林分類器計算速度比單個決策樹慢，需要自身算法推斷出超出范圍的獨立變量或非獨立變量，但其能處理高維特征，不易產生過擬合，模型訓練速度比較快，可以產生高準確度的分類器，分類效果比較好。但隨機森林決策樹的數量對隨機森林算法的效率有重要影響。當決策樹數量偏少時，分類精度降低；決策樹偏多時，分類精度會趨于穩定，但運算速度緩慢。有研究表明，在1～100范圍內，當決策樹大于50時，再繼續增大隨機森林決策樹數量已無意義，因此本文選取決策樹數量為50[32]。

2.2 精度評價

為有效評估不同數據對不同覆被類型的提取精度，本文從分類總體精度、混淆矩陣、Kappa系數、制圖精度和用戶精度5方面進行比較分析[33]。

3 實驗結果與分析

3.1 融合月平均時序Sentinel-1 SAR數據分類

在氣候濕潤、多霧多云的地方，光學數據難以獲取。Sentinel-1衛星搭載的C波段波長遠高于云粒子的典型尺寸，且自身發出能量脈沖，故可穿透大部分云霧并不受光照條件的影響，被廣泛應用于農作物提取研究。單利用某一時相的SAR極化特征數據有極大的局限性。有研究表明，將月尺度上觀測數據進行均值合成后，能夠降低降水對分類精度的影響，較大程度提高禾谷類作物的分類精度[12]。因此本研究選取冬小麥生長期內2019年11月到2020年5月Sentinel-1 SAR數據，在月尺度上進行均值合成，用于冬小麥識別研究。

融合7個月均值Sentinel-1 SAR數據的VV、VH極化雷達后向散射系數和其通過GLCM計算的紋理特征，應用RF算法進行分類。在GLCM計算紋理特征共生矩陣大小的選擇上，選取4、8、16鄰域數值進行計算，分類精度最高為4鄰域。此時，融合多時相Sentinel-1極化特征和紋理特征的分類總體精度為91.00%，Kappa系數為0.84。冬小麥的制圖精度為95.66%，用戶精度為96.58%，結果如表4和圖2所示。從碎部地區分類結果與GF-6的融合多光譜分辨率2 m數據(圖3)對比發現，加入紋理特征導致覆被結構復雜的區域呈現斑塊現象，碎部覆被細節體現不明顯，這是共生矩陣計算紋理特征時應用周圍像素參與計算和破碎地塊像素變化不規律導致的結果。

表4 多時相Sentinel-1 SAR數據的極化特征和紋理特征分類結果Tab.4 Classification results of multi-temporal Sentinel-1 SAR data combining texture features and polarization features

在遙感影像分類中，特征變量權重能夠增強對特征變量之間的理解。在RF算法中，基于袋外數據，可以對輸入數據的權重進行驗證，獲得每個特征的權重值，用平均精度減少值來表示，其原理是將某特征參數數值變為隨機數，計算其對模型準確率的影響，根據得到的精度減少值來計算此參數的權重，該值越大表明該變量的權重越大[24]。GEE中通過explain()函數可輸出各特征的權重，特征變量權重的總和與變量的數目相關，并非固定值，其權重在自身變量組中具有相對意義。圖4表示在分類過程中各變量的特征值權重，各特征變量在不同時期的重要性差異較大。對比VH、VV與紋理特征的特征值，在每個月平均值中，除去2019年11月和2020年2月，極化數據特征對分類的貢獻率相對紋理特征較大。有研究表明，分類精度和特征重要性與特征值之間的差異性呈正比例關系[13]。從圖5中發現，2019年11月和2020年2月重要性降低的原因是冬小麥與其他用地、其他植被、林地的極化值差異較小。

去除紋理特征，僅使用月平均的VH和VV極化數據進行冬小麥提取，結果如表5和圖6所示。此時，融合月平均VH和VV極化數據的分類總體精度達到85.93%,Kappa系數為0.75，相對于添加紋理特征，總體精度降低5.07個百分點，Kappa系數降低10.71%。由此表明，綜合極化數據和紋理特征使分類精度較高，說明多源特征量的綜合有利于總體分類結果的改善。受冬小麥地塊連片種植的影響，冬小麥的制圖精度和用戶精度變化僅為0.01、0.4個百分點，識別精度受紋理特征影響較小。對于其他覆被類型，雖單使用融合極化月平均數據識別精度較低，但破碎地塊細節體現明顯，其原因有待進一步研究。

表5 多時相Sentinel-1 SAR的極化特征分類結果Tab.5 Classification results of multi-temporal Sentinel-1 SAR data polarization features

3.2 融合Sentinel-1 SAR和Sentinel-2光學數據分類

SAR數據的地物后向散射特性異于光學遙感影像。光學數據反映目標體的光譜特性，SAR數據的穿透性不僅能夠獲取植被的表面信息，對植被的葉、莖、枝干等信息也有一定的反映，獲取的是不同于光學數據的地物信息[34]。圖7為不同覆被在冬小麥生長期間的Savizky-Golay濾波NDVI變化圖。越冬期過后，2—3月溫度上升，冬小麥進入返青期和拔節期，葉綠素含量逐漸升高，而林草與其他植被在此時生長較為緩慢，導致NDVI特征值差異較大。3—4月為冬小麥拔節期/孕穗期，整體植被指數較高，NDVI差異值減小。在5月中旬以后，冬小麥在灌漿期、收獲期葉片含水率、葉綠素含量降低，葉黃素含量升高，植被指數呈陡然下降趨勢，但其他植被的植被指數開始升高，此時冬小麥NDVI會與其他覆被的NDVI存在交叉。基于特征值差異值越大，RF算法分類精度越高的原則，選取返青期作為光學數據的輸入，數據日期為2020年2月22日。在GEE中，通過get(‘CLOUDY_PIXEL_PERCENTAGE’)函數，可獲取該地區2020年2月22日Sentinel-2 光學影像的含云量，輸出結果為0.45%。

首先，對返青期的光學數據進行分類研究，結果如表6和圖8所示。返青期Sentinel-2光學數據的分類總體精度為94.01%,Kappa系數為0.89，冬小麥制圖精度和用戶精度分別為96.34%、97.67%。相對于加紋理特征的融合Sentinel-1 SAR數據分類精度，返青期Sentinel-2的總體精度提高3.01個百分點，Kappa系數提高5.95%。但冬小麥的制圖精度和用戶精度差距為1個百分點左右。單期光學影像相對于時序SAR數據，在識別冬小麥過程中優勢明顯。

表6 冬小麥返青期Sentinel-2光學數據分類結果Tab.6 Sentinel-2 optical data classification results of winter wheat at greenback period

圖9為冬小麥返青期Sentinel-2光學數據特征值的權重。光學數據熵值紋理特征此時為0，植被指數特征值權重相對均衡，主要是在選取各類指數過程中，不同指數反映了作物生長階段的一些生長特征。但遙感信息反映的只是地表或作物群體瞬間的物理狀況，而作物生長是一個隨時間變化的連續過程。在進一步研究中，應考慮光學遙感信息與物候信息的結合，突出作物生長-衰老交替引起的季節性植被變化，降低同物異譜和異物同譜現象，特別在混合像元及破碎地塊對冬小麥的提取上具有很強的優勢。但本文主要研究時序SAR數據和光學數據的互補性，開發一種應用少量光學數據或無光學數據的分類方法。

融合時序Sentinel-1 SAR、SAR紋理特征和返青期Sentinel-2光學數據，進行冬小麥提取研究。表7和圖10結果表明，當時序Sentinel-1 SAR加入紋理特征時，與返青期Sentinel-2光學數據融合的分類總體精度為95.31%，Kappa系數為0.92，冬小麥的制圖精度及用戶精度為97.71%和97.23%。融合Sentinel主被動遙感數據的總體精度和Kappa系數，相對于加入紋理特征的融合月平均Sentinel-1 SAR數據提高4.31個百分點和9.52%，冬小麥制圖精度和用戶精度提高2個百分點左右，相對于返青期的Sentinel-2光學數據提高1.3個百分點和3.37%。融合Sentinel主被動遙感數據能夠提高分類總體精度和冬小麥的識別精度。

表7 多時相Sentinel-1 SAR數據的極化特征、紋理特征融合單期光學數據分類結果Tab.7 Classification results of multi-temporal Sentinel-1 SAR data polarization features and texture features integrating a single optical data

在融合Sentinel主被動遙感數據特征值權重方面，如圖11所示，光學數據的權重高于SAR數據的權重，這也是應用RF算法時，單基于光學數據高于融合月平均時序SAR數據的原因。紋理特征值的權重相對較弱。時序Sentinel-1 SAR不加入紋理特征時，與返青期Sentinel-2光學數據融合的分類總體精度為95.78%，Kappa系數為0.92，相對于未添加紋理特征的融合月平均Sentinel-1 SAR數據提高9.85個百分點和22.67%，冬小麥的制圖精度及用戶精度為98.56%和97.92%，如表8和圖12所示。相對于添加紋理特征的聯合數據總體精度增加0.47個百分點，兩者精度差別不明顯。但在冬小麥及其他地物提取方面，添加紋理特征的融合數據冬小麥識別精度降低0.8個百分點左右，特別是其他用地的制圖精度降低18.86個百分點。綜合對比表4和表5，表7和表8，紋理特征對冬小麥的識別精度影響程度較小。可見，融合Sentinel-1 SAR和Sentinel-2光學數據分類與融合Sentinel-1 SAR數據同樣表現為，添加紋理特征的融合數據在覆被類型復雜地塊的分類精度上，并不隨特征值的增加而占優勢。

表8 多時相Sentinel-1 SAR極化特征融合單期光學數據分類結果Tab.8 Classification results of multi-temporal Sentinel-1 SAR data polarization features integrating a single optical data

4 結論

(1)在缺乏光學數據的情況下，融合冬小麥生長期時序月平均SAR數據對冬小麥的提取精度能夠達到95%，受冬小麥地塊連片區域的影響，紋理特征對冬小麥識別精度影響程度較小。

(2)時序月平均SAR數據與單期光學數據聯合時，添加紋理特征的分類總體精度與未添加紋理特征的總體精度差異微小，都達到95%左右。融合Sentinel主被動遙感數據，在破碎地塊，添加紋理特征會導致分類精度降低，但對冬小麥的提取精度同樣達到98%左右。