基于GF-1與Sentinel-2融合數(shù)據(jù)的地膜識別方法研究

羅 琪，劉曉龍,2*，史正濤，屈 冉，趙文智

(1.云南師范大學(xué)地理學(xué)部，云南 昆明 650500；2.內(nèi)蒙古工業(yè)大學(xué)信息工程學(xué)院，內(nèi)蒙古 呼和浩特 010051；3.生態(tài)環(huán)境部衛(wèi)星環(huán)境應(yīng)用中心，北京 100094；4.北京師范大學(xué)地理科學(xué)學(xué)部，北京 100875)

0 引言

地膜能夠調(diào)節(jié)農(nóng)田土壤光、熱、水、氣等條件，緩解作物遭受低溫凍害、高溫干旱、病蟲害等的影響，進(jìn)而促進(jìn)作物生長發(fā)育、提高作物產(chǎn)量[1,2]，但地膜易破碎且難以回收和降解，大量殘膜造成土壤肥力下降、作物減產(chǎn)及土壤污染等農(nóng)業(yè)生產(chǎn)安全與生態(tài)環(huán)境問題[3-6]。例如，地膜殘膜促使土壤中微塑性物質(zhì)含量增加[7]，其吸收有機(jī)污染物/吸附病原體并通過空氣傳播被生物體吸收，進(jìn)而進(jìn)入食物網(wǎng)對人類健康構(gòu)成威脅[8,9]。隨著我國地膜使用面積的增加，準(zhǔn)確獲取地膜的空間分布信息并進(jìn)行科學(xué)規(guī)劃與管理對提高農(nóng)業(yè)生產(chǎn)與土壤污染防治及保障人類安全具有重要意義[10]。

目前，國內(nèi)外學(xué)者多采用遙感數(shù)據(jù)進(jìn)行大范圍的地膜識別研究。例如：Lu等[11]基于250 m低空間分辨率的MODIS-NDVI時(shí)間序列數(shù)據(jù)，利用光譜特征閾值法實(shí)現(xiàn)了對光譜特征比較穩(wěn)定(田塊尺度較大，作物均一性高)的新疆覆膜棉田的高精度識別，但目前我國大部分地區(qū)種植結(jié)構(gòu)零碎，基于低分辨率數(shù)據(jù)的地膜識別方法存在一定的局限性。中分辨率數(shù)據(jù)則可以實(shí)現(xiàn)種植結(jié)構(gòu)相對均一地區(qū)的地膜識別，如：Lanorte等[12]基于Landsat-8數(shù)據(jù)，利用支持向量機(jī)(SVM)分類法提取了意大利南部阿普利亞地區(qū)的地膜；沙先麗[13]基于Landsat-5 TM數(shù)據(jù)的光譜特征、指數(shù)特征和紋理特征，通過構(gòu)建決策樹對新疆覆膜棉區(qū)的識別精度達(dá)90%以上。但在種植結(jié)構(gòu)破碎的局地地膜精細(xì)識別中，中分辨率遙感影像易受混合像元影響，難以準(zhǔn)確捕捉地膜的細(xì)節(jié)信息，因此，多數(shù)學(xué)者利用高空間分辨率遙感影像提取紋理、形狀等空間結(jié)構(gòu)特征進(jìn)行地膜識別研究。例如：高夢婕等[14]基于GF-2遙感影像探討最佳分割尺度及最佳特征選取，并采用隨機(jī)森林、CART決策樹、支持向量機(jī)分類法提取塑料大棚，結(jié)果顯示隨機(jī)森林法的分類精度最高(89.65%)；趙璐等[15]基于GF-2遙感影像，利用面向?qū)ο笞罱彿椒ㄌ崛∞r(nóng)業(yè)大棚信息，精度達(dá)95.65%；吳錦玉等[16]基于GF-2遙感影像，通過對比分析不同紋理提取算法，發(fā)現(xiàn)利用影像紋理可顯著提高大棚識別精度。

受限于目前的衛(wèi)星成像技術(shù)，高分辨率遙感數(shù)據(jù)大都僅包含可見光波段和近紅外波段，缺少對地物反射敏感的短波紅外波段，而具有短波紅外光譜信息的中分辨率數(shù)據(jù)存在分辨率不足。為此，一些學(xué)者嘗試基于多源數(shù)據(jù)和多特征的組合進(jìn)行地膜識別研究。如哈斯圖亞[17]利用GF-1、Landsat-8和Radarsat-2數(shù)據(jù)的地膜識別精度顯著高于單一GF-1數(shù)據(jù)的地膜識別精度，但Radarsat-2遙感影像重訪周期長且成本高，限制了該方法的實(shí)用性[18]。考慮到GF-1系列衛(wèi)星數(shù)增加及Sentinel-2數(shù)據(jù)的高重訪(小于5 d)特性，本文嘗試將兩種數(shù)據(jù)融合，旨在充分利用GF-1數(shù)據(jù)的空間細(xì)節(jié)信息和Sentinel-2數(shù)據(jù)的光譜信息，探討多源光學(xué)遙感數(shù)據(jù)融合在地膜識別中的可行性。

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

1.1 研究區(qū)概況

本文研究區(qū)為河北省邯鄲市邱縣(圖1)，區(qū)域高程為32.8～40.8 m，屬北暖溫帶亞濕潤大陸性季風(fēng)氣候，年均氣溫13 ℃，年均降雨量524.6 mm，氣候溫和，雨熱同季，無霜期長，以棉花、蔬菜地膜種植為主，其中棉花種植面積達(dá)33萬hm2[19]。

圖1 研究區(qū)位置(GF-1影像RGB顯示)Fig.1 Location of the study area (RGB display of GF-1 image)

1.2 數(shù)據(jù)及其處理

研究區(qū)地膜覆蓋時(shí)間集中于4月中旬，揭(挑)地膜時(shí)間為6月中旬。本文選用的GF-1、Sentinel-2遙感影像獲取時(shí)間均為2019年5月23日穩(wěn)定覆膜期，是地膜光學(xué)遙感識別的理想時(shí)間窗口。

(1)GF-1數(shù)據(jù)。時(shí)間分辨率為5 d，PMS傳感器可獲得空間分辨率為2 m的全色影像和8 m的多光譜(BlueGreenRedNIR波段)影像。首先對影像進(jìn)行輻射定標(biāo)，并利用FLAASH模塊進(jìn)行大氣校正，然后利用GF-1衛(wèi)星數(shù)據(jù)自帶的RPC文件進(jìn)行正射校正。通過Gram-Schmidt Spectral Sharpening算法對影像的多維特征進(jìn)行正交化處理，可消除數(shù)據(jù)冗余且光譜保真度高，能明顯提高影像的空間細(xì)節(jié)信息，融合效果好，尤其是在黑色背景區(qū)域較多的情況下[20-22]。為保證GF-1數(shù)據(jù)的光譜和空間分辨率，本文采用該算法融合GF-1多光譜和全色波段，最終獲得空間分辨率為2 m的多光譜影像，用于對比分析本文方法的地膜識別精度。

(2)Sentinel-2數(shù)據(jù)。Sentinel-2當(dāng)前在軌運(yùn)行A、B兩顆星，覆蓋13個(gè)光譜波段，地面分辨率分別為10 m(BlueGreenRedNIR波段)、20 m(Vegetation Red EdgeSWIR波段)、60 m(Coastal aerosolWater vapourSWIR-Cirrus波段)，單顆衛(wèi)星的重訪周期為10 d，A、B雙星互補(bǔ)后在赤道地區(qū)的重訪周期達(dá)5 d。本文對下載的Sentinel-2-L1C產(chǎn)品用Sen2Cor插件進(jìn)行大氣校正，并根據(jù)研究需求，將數(shù)據(jù)的BlueGreenRedNIRSWIR波段合成并重采樣為空間分辨率為10 m的多光譜影像。

2 研究方法

本文將GF-1全色影像與Sentinel-2多光譜影像進(jìn)行NNDiffuse Pan Sharpening融合后，利用分層分類(HC)法進(jìn)行地膜識別，并與基于GF-1遙感影像/多源融合影像的深度學(xué)習(xí)(DL)法、傳統(tǒng)機(jī)器學(xué)習(xí)分類方法進(jìn)行精度對比，驗(yàn)證方法的可行性。

2.1 融合算法

Gram-Schmidt Spectral Sharpening算法用于GF-1影像融合效果佳[23,24]，但該算法耗時(shí)長，且對不同數(shù)據(jù)源的影像配準(zhǔn)精度要求高，更適用于同一平臺傳感器的影像融合[25],故本文僅用該算法對GF-1數(shù)據(jù)進(jìn)行融合。NNDiffuse Pan Sharpening融合方法將單個(gè)像元光譜值作為融合運(yùn)算的最小元素，通過混合模型產(chǎn)生分辨率增強(qiáng)的融合影像，該方法運(yùn)算速度快，能較好地保留原影像的色彩和紋理，且融合后影像的空間特征和光譜特征不會發(fā)生較大變化[26]，但需注意的是，待融合遙感影像多光譜波段間的波長范圍不能交叉覆蓋，且全色波段能夠基本覆蓋多光譜所有波段的波長范圍[27]，故本文利用該算法對GF-1、Sentinel-2數(shù)據(jù)進(jìn)行融合。首先，以預(yù)處理后的GF-1全色影像為參考影像，Sentinel-2多光譜影像為待校正影像，采用6個(gè)地面控制點(diǎn)和一階多項(xiàng)式變換模型進(jìn)行幾何精校正，減少非系統(tǒng)因素產(chǎn)生的誤差，校正后影像的幾何配準(zhǔn)偏差小于0.5個(gè)像元；再以GF-1全色影像為高分辨率影像，Sentinel-2多光譜影像為低分辨率影像，進(jìn)行NNDiffuse Pan Sharpening融合，得到空間分辨率為2 m的多光譜(BlueGreenRedNIRSWIR波段)影像。

2.2 地膜識別方法

對地膜識別的多特征信息進(jìn)行可分離性分析有利于區(qū)分地膜和其他地物。本文通過目視解譯獲取190個(gè)樣本，計(jì)算研究區(qū)不同土地類型(地膜、植被、水體、裸土和休耕地、建筑物)在融合數(shù)據(jù)RedGreenBlueNIRSWIR波段的反射率均值，并基于NDVI[28]、改進(jìn)的歸一化差異水體指數(shù)(MNDWI)[29]、歸一化建筑指數(shù)(NDBI)[30]、建筑用地指數(shù)(IBI)[31]計(jì)算不同土地覆蓋類型的反射率(圖2)，進(jìn)行地膜與其他地物的可分離性分析。由圖2可知，地膜在Blue、Green、Red、NIR、SWIR波段的反射率均高于裸土和休耕地，且地膜在Blue、Green、Red、SWIR波段的反射率明顯高于水體、植被，但其與建筑物的光譜反射率較為接近。由于利用歸一化差異水體指數(shù)(NDWI)[32]提取水體時(shí)存在許多非水體信息，且建筑物等陰影在Green和NIR波段的波譜特征與水體相似[29]，故利用SWIR波段數(shù)據(jù)，采用MNDWI提取水體信息，再利用NDVI提取植被，實(shí)現(xiàn)地膜與水體、植被的分離；同時(shí)，利用NDBI和IBI提取建筑物，但這兩種指數(shù)難以實(shí)現(xiàn)地膜與建筑物的分離。像元形狀指數(shù)(PSI)[33]是利用中心像元及其相鄰像元的關(guān)系，計(jì)算具有光譜相似性相鄰像元組的維數(shù)，并拉大不同方向鄰域像元間的特征值，以增強(qiáng)影像的空間細(xì)節(jié)信息。利用PSI能提取高分辨率影像的形狀和結(jié)構(gòu)特征，可彌補(bǔ)光譜特征的不足，在規(guī)則紋理人工地類提取中具有較強(qiáng)優(yōu)勢[16]，故本文利用PSI提取研究區(qū)建筑物信息。PSI的計(jì)算步驟為：1)定義圍繞中心像元的方向線，不同方向線的長度通過相鄰像元間的光譜同質(zhì)性測度值(式(1))和閾值確定，每條方向線從中心像元向兩邊擴(kuò)展，擴(kuò)展條件包括光譜約束條件(PHi(x,y)

圖2 不同土地覆蓋類型的反射率與遙感指數(shù)值(均值)Fig.2 Reflectance and remote sensing index values of different land cover types (mean value)

(1)

(2)

(3)

分層分類(HC)法是根據(jù)各類地物的光譜特征，按照一定的準(zhǔn)則對每一層級進(jìn)行分解，以達(dá)到高精度分類或?qū)ｉT提取某一地類的目的。針對待分類地物的光譜特征選擇不同的特征參數(shù)或分類方法，逐層提取地物信息并進(jìn)行信息掩膜處理，從而消除對其他地物提取產(chǎn)生的影響，然后將逐層分類結(jié)果疊加得到最終分類結(jié)果[34]。首先，基于融合影像利用MNDWI、NDVI提取水體和植被，并對提取結(jié)果進(jìn)行掩膜處理，然后利用PSI提取建筑物，并再次進(jìn)行掩膜處理。本文是在覆膜穩(wěn)定時(shí)期進(jìn)行的地膜識別研究，適逢農(nóng)作物生長發(fā)育期，裸土和休耕地分布較少。通過分析地膜與裸土和休耕地在Blue、Green、Red、NIR、SWIR波段的光譜反射率，可知地膜的反射率高于裸土和休耕地的反射率(圖2)，利用反射率差異，基于最大類間方差法(OSTU)[35]實(shí)現(xiàn)地膜與裸土和休耕地的區(qū)分(圖3)。OSTU計(jì)算簡單，自適應(yīng)性強(qiáng)，對類間方差為單峰的影像分割效果好[36-39]，故本文通過該方法確定水體、植被、建筑物及地膜的提取閾值。

圖3 分層分類框架Fig.3 Framework of hierarchical classification

2.3 對比驗(yàn)證方法

為驗(yàn)證本文方法在地膜識別中的可行性，將其與基于深度學(xué)習(xí)、傳統(tǒng)機(jī)器學(xué)習(xí)的分類方法進(jìn)行比較。其中，深度學(xué)習(xí)法通過分層處理機(jī)制，構(gòu)建含有多個(gè)隱藏層的學(xué)習(xí)模型，進(jìn)行影像特征的自動提取和選擇，具有很強(qiáng)的特征表達(dá)、容錯(cuò)能力和魯棒性[40]，但其需要大量標(biāo)記的訓(xùn)練樣本，在小尺度遙感影像數(shù)據(jù)集情況下會降低分類性能，且需要對分類結(jié)果進(jìn)行優(yōu)化處理，以提高地物邊界的細(xì)化程度，在遙感影像分類應(yīng)用領(lǐng)域存在局限[41]；另外，該方法對計(jì)算機(jī)性能要求較高，需要使用GPU或TPU加快學(xué)習(xí)模型的訓(xùn)練速度，成本較高[42,43]。目前，基于深度學(xué)習(xí)的遙感影像分類以監(jiān)督分類為主，標(biāo)記樣本數(shù)據(jù)(影像塊、語義分割)的數(shù)量和質(zhì)量也將影響模型應(yīng)用的精度[44]。支持向量機(jī)能較好地解決小樣本統(tǒng)計(jì)學(xué)習(xí)問題，對高維數(shù)據(jù)不敏感，具有較高適應(yīng)性和泛化能力，且能避免過學(xué)習(xí)問題[45]。隨機(jī)森林法的計(jì)算效率高，適合高維數(shù)據(jù)處理，對數(shù)據(jù)噪聲和小樣本數(shù)據(jù)不敏感，具有較好的魯棒性和穩(wěn)定性，并且在運(yùn)算過程中人工干擾少，不易出現(xiàn)過擬合問題[46]。

3 結(jié)果與分析

3.1 融合結(jié)果及精度評價(jià)

本文采用定性和定量相結(jié)合的方法對GF-1與Sentinel-2的融合影像進(jìn)行評價(jià)，以驗(yàn)證本文融合方法是否有益于地膜識別。1)為便于定性分析融合影像的顏色、清晰度和空間紋理特征等細(xì)節(jié)信息，在研究區(qū)裁剪出地膜分布集中且與其他地物邊界明顯的局部區(qū)域(圖4)。由圖4可知，相對于Sentinel-2多光譜影像，融合影像的空間分辨率明顯提高，黑色圓圈中地物的紋理特征、邊緣細(xì)節(jié)更清晰，但其相比GF-1影像清晰度有所降低。2)利用目前常用的平均梯度、信息熵和相關(guān)系數(shù)[47,48]指標(biāo)對融合影像和GF-1、Sentinel-2影像的各波段數(shù)據(jù)進(jìn)行評價(jià)(表1)。由表1可知，融合影像各波段的信息熵均接近7.0，保留了較多原始影像的信息量，且平均梯度與GF-1接近，但低于Sentinel-2。由于受多源數(shù)據(jù)獲取方式以及衛(wèi)星姿態(tài)、擾動等影響，多源遙感影像的幾何精校正存在偏差，使融合方法存在不確定性[49]；另外GF-1數(shù)據(jù)與Sentinel-2數(shù)據(jù)的空間分辨率存在差異，且受地表復(fù)雜性等因素影響，難以獲取理想?yún)⒖加跋?，加之GF-1融合后的數(shù)據(jù)與Sentinel-2數(shù)據(jù)存在不可避免的偏差，導(dǎo)致利用GF-1、Sentinel-2數(shù)據(jù)和融合數(shù)據(jù)進(jìn)行的指標(biāo)評價(jià)存在較大的不確定性[50]。融合影像各波段相關(guān)系數(shù)高于0.8，且與Sentinel-2數(shù)據(jù)的相關(guān)系數(shù)達(dá)0.98以上，說明融合數(shù)據(jù)的保真度較高。

圖4 原始影像與融合影像對比Fig.4 Comparison of original image and fusion image

表1 融合質(zhì)量評價(jià)結(jié)果Table 1 Fusion quality evaluation results

3.2 地膜識別結(jié)果比較

本文設(shè)計(jì)不同分類方案(表2)，將分層分類法的識別結(jié)果與基于深度學(xué)習(xí)、支持向量機(jī)和隨機(jī)森林方法的識別結(jié)果進(jìn)行比較。為更好地體現(xiàn)不同分類方法的可比性，將研究區(qū)地類分為地膜和其他(建筑物、裸土和休耕地、植被、水體)兩個(gè)類別。

表2 不同分類方案Table 2 Different classification schemes

本文利用混淆矩陣法，選用總體精度(OA)、Kappa系數(shù)、生產(chǎn)者精度(PA)、用戶精度(UA)4個(gè)指標(biāo)[51]對地膜識別結(jié)果進(jìn)行定量分析(表3)，精度驗(yàn)證的參考數(shù)據(jù)為人工目視選取的173個(gè)樣本。

表3 不同方案識別精度Table 3 Accuracy of different identification schemes

由表3可知，F(xiàn)us+HC方案的分類精度最高，總體精度達(dá)94.87%，Kappa系數(shù)達(dá)0.89；其次是GF-1+DL和Fus+RF方案，總體精度分別為93.14%、92.49%，Kappa系數(shù)分別為0.85、0.83。與GF-1+DL方案相比，F(xiàn)us+DL方案總體精度降低1.16%，Kappa系數(shù)降低0.03。多次實(shí)驗(yàn)表明，基于深度學(xué)習(xí)的識別方法易受影像分辨率的影響，用于混合像元較多的中、高分辨率影像時(shí)，其識別效果有限[41]。Fus+RF方案比GF-1+RF方案總體精度提高5.71%，Kappa系數(shù)提高0.11；Fus+SVM方案比GF-1+SVM方案總體精度提高1.59%，Kappa系數(shù)提高0.03。整體上，除基于深度學(xué)習(xí)的地膜識別方法外，基于多源融合數(shù)據(jù)的地膜識別精度均高于基于GF-1遙感數(shù)據(jù)的識別精度。

從識別結(jié)果的空間分布看(圖5)：GF-1+SVM和GF-1+RF方案分類結(jié)果中地膜與建筑物存在較多混淆，尤其是在地物復(fù)雜區(qū)域；Fus+RF和Fus+SVM方案將建筑物區(qū)域內(nèi)的植被誤分為地膜，導(dǎo)致分類結(jié)果出現(xiàn)“椒鹽”現(xiàn)象；GF-1+DL和Fus+DL方案對“椒鹽”現(xiàn)象進(jìn)行了改善，地膜識別精度均高于90%，但其識別結(jié)果中地物邊界過于平滑，需引入影像后處理方法(如形態(tài)學(xué)濾波[52]、全連接CRFs[53]等)優(yōu)化分類結(jié)果，使分類后的地物邊界更接近真實(shí)地物邊界。Fus+HC方案識別出的地物邊界輪廓與真實(shí)地物邊界基本一致，空間細(xì)節(jié)特征表現(xiàn)明顯，分類效果更好，故本文利用Fus+HC方案提取研究區(qū)地膜的空間分布(圖6)。

圖5 不同方案分類結(jié)果Fig.5 Classification results of different schemes

圖6 基于分層分類法的地膜識別結(jié)果Fig.6 Identification results of plastic mulch based on hierarchical classification

4 結(jié)論與討論

針對單一遙感數(shù)據(jù)源在地膜識別中存在混合像元、光譜特征不足等問題，本文提出一種基于GF-1與Sentinel-2融合數(shù)據(jù)的地膜識別方法，并與機(jī)器學(xué)習(xí)、深度學(xué)習(xí)方法提取結(jié)果進(jìn)行精度對比，結(jié)論如下：1)基于大樣本空間的深度學(xué)習(xí)、傳統(tǒng)機(jī)器學(xué)習(xí)分類方法雖然能夠獲得高地物識別精度，但其對于有限樣本空間的地物識別不如閾值法，本文基于多源融合數(shù)據(jù)的分層分類法通過同時(shí)相GF-1數(shù)據(jù)和Sentinel-2數(shù)據(jù)融合，在地膜識別應(yīng)用中獲得了較為可靠的結(jié)果(總體精度為94.87%，Kappa系數(shù)為0.89)，明顯高于Fus+DL、Fus+RF和Fus+SVM方案地膜識別的總體精度(分別為91.98%、92.49%和87.50%)和Kappa系數(shù)(分別為0.82、0.83、0.74)。2)地膜具有較強(qiáng)的方向反射特性，多角度成像是實(shí)現(xiàn)地膜識別的關(guān)鍵。本文未考慮不同衛(wèi)星數(shù)據(jù)源多角度特征對識別精度的影響，后續(xù)將進(jìn)一步分析識別精度與成像幾何問題，并探討不同空間分辨率(8 m、10 m、20 m)多光譜影像組合對識別精度的影響，以期獲得更佳的識別效果。隨著Sentinel-2數(shù)據(jù)和我國GF-1系列數(shù)據(jù)的發(fā)展，這兩種數(shù)據(jù)的結(jié)合可實(shí)現(xiàn)多源遙感數(shù)據(jù)的優(yōu)勢互補(bǔ)，為解決地膜識別以及積雪和冰川識別等問題提供了新思路。