張新樂,官海翔,劉煥軍,2,孟祥添,楊昊軒,葉 強,于 微,張漢松※
(1.東北農業大學公共管理與法學院,哈爾濱150030:2.中國科學院東北地理與農業生態研究所,長春130012)
中國是一個農業大國,也是世界上受風災影響較為嚴重的國家[1]。風災倒伏面積的準確提取,能夠為災后的農業生產、政府決策以及農業保險理賠,提供數據和技術支持[2-3]。傳統的作物風災倒伏信息是通過田間實測的方式獲取,工作量比較大;衛星遙感技術雖然能夠進行全方面、大范圍地提取倒伏面積信息,但時效性相對較差,且影像易受云層干擾[4-5],而且一般衛星影像的空間分辨率難以滿足地塊尺度下的倒伏面積信息精準提取。近些年來,無人機遙感憑借精度高、受地形約束小、成本低、操作簡便等優勢,迅速占領了遙感市場,并逐漸成為農業災害信息獲取的重要方式[6-7]。與衛星遙感相比,無人機遙感不僅具有更高的時間和空間分辨率,而且能夠獲取田塊尺度作物的空間變異信息,能反映地塊內部不同位置和生長期作物的長勢差異[8-11]。因此,無人機遙感技術非常適用于田塊尺度的作物倒伏面積提取。
目前國內外研究主要使用無人機圖像或衛星遙感影像,選取對倒伏作物識別度較高的特征對作物倒伏面積信息進行提取。Liu 等利用無人機獲取的彩色圖像和熱紅外數據,構建了倒伏水稻的識別模型,并對模型的識別精度進行了驗證[12]。李廣等對無人機彩色圖譜進行了二次低通濾波分析,并將兩種單特征的線性擬合作為特征集,由此提取了小麥倒伏的空間分布[13]。李宗南等利用Worldview-2 多光譜影像,分析了灌漿期倒伏玉米地塊的紋理和光譜特征,基于優選的3 波段紋理和光譜特征估算了倒伏玉米面積[14]。Yang 等將無人機圖像、紋理特征和數字表面模型進行了不同組合,并計算了不同組合特征中各類作物形態的分類精度,選取了彩色圖像和數字表面模型組合特征作為倒伏水稻識別的指標[15]。李宗南等將研究區內的作物形態劃分為倒伏玉米與未倒伏玉米兩類,通過對該兩類作物形態的紋理特征對比,選取了均值紋理特征作為倒伏玉米識別的指標[16]。Han 等使用無人機獲取了倒伏玉米的多光譜和可見光圖像,并提取紋理、冠層結構、植被指數等參數,以其為輸入量構建了2 種倒伏提取的Logistic 模型,并驗證了模型的精度[17]。Chu 等利用超空間分辨率的無人機影像構建了玉米的冠層高度模型,通過設立高度閾值,實現了倒伏與未倒伏玉米的分類[18]。毛智慧等通過無人機傾斜測量的方式獲取了倒伏作物的數字表面模型(DSM),并結合無人機圖像的色彩特征,進行了玉米倒伏信息提取[19]。劉良云等基于倒伏前后的兩景TM 影像,利用NDVI 值監測了小麥抽穗初期的倒伏發生程度[20]。王立志等利用多時相的多光譜衛星影像,分析了抽雄期玉米倒伏前后的多種植被指數變化,采用RVI差值監測了倒伏受災范圍[21]。但以上研究主要集中于作物生長旺盛期的倒伏信息提取,而未對完熟期的作物倒伏提取進行研究,并且所使用的倒伏提取特征主要以單類紋理、光譜與植被指數特征或不同類型特征間的組合為主,弱化了紋理特征對于倒伏作物的空間域形態信息的反映,也未對紋理特征提取倒伏作物的最優維數進行探討,這在一定程度上造成了倒伏識別特征的信息冗余和分類結果的過擬合效應。
正常生長的玉米在完熟期葉綠素含量較低,其葉片呈現淡黃色,與倒伏玉米葉片的顏色很相似,兩者的反射光譜特征差異較小[22]。而處于生長發育旺盛期的玉米,其倒伏前后的光譜特征差異較大,因此基于生長旺盛期的玉米倒伏識別特征并不適用于完熟期玉米倒伏面積識別。本文利用5 波段的無人機多光譜數據,分析研究區內不同作物形態的光譜、植被指數及紋理特征,優選分類誤差最低的特征組合,提取完熟期玉米地塊內的倒伏面積,定量分析和評價倒伏面積提取精度,以期尋求最佳的完熟期玉米倒伏面積識別的特征組合和方法,為農業災損的精準評估提供技術參考。
友誼農場位于黑龍江省雙鴨山市友誼縣西部,隸屬于黑龍江農墾總局紅興隆管理局,研究區位于該農場西北部,常年種植玉米,占地面積為172 975 m2,地理坐標為131.6°E,46.7°N,具體位置如圖1所示。此地塊玉米正處于完熟期,并存在4 類典型作物形態,主要包括:長勢較差、葉片呈現綠色的未倒伏A 類玉米;長勢較好、葉片呈淡黃色的未倒伏B 類玉米;葉片呈現淡黃色的倒伏玉米;輻射亮度較低、分布零散的陰影區域。該地區位于三江平原地區,氣候屬于溫帶大陸性季風氣候,雨熱同期,冬季干燥,夏季溫熱多雨。平均年降雨量在540 mm 左右,年日照時數約為2 500 h,有效積溫達2 300 ℃,海拔高度在64~70 m 之間,地勢平坦,土壤類型為暗棕壤,肥力較高。主要種植玉米和水稻,一年一熟。
本實驗采用大疆600 M Pro 六旋翼式飛行器作為遙感平臺,該平臺質量為10 kg,最大可承受風速8 m/s,能夠在-10~40 ℃的溫度環境下正常工作,最大承重上限5.5 kg,最高續航時間是38 min。根據飛行器的荷載重量,搭載MicaSense RedEdge?3 多光譜相機,該相機有紅、綠、藍、紅邊、近紅外5 個波段,配備光強傳感器和一塊白板。在飛行過程中,光強傳感器可以校正太陽光線的變化對影像所產生的影響,白板則具有固定的反射率信息,利用白板對影像進行校正即可生成反射率圖像。數據采集時間為2018 年9 月12 日,天氣晴朗無云,風速小于1 級,適合飛行。為保證影像的完整和準確,通過地面站設置飛行的航向重疊為80%,旁向重疊為79%,主航線角度99°,邊距-1.8 m,航線9 條,飛行高度為110 m,實際高度109.5 m,圖像的空間分辨率為0.074 m,航拍完成后采用自動返航的方式降落。

圖1 研究區位置及驗證點分布Fig.1Locationandverificationpointdistributionofresearcharea
首先對無人機采集的照片進行篩選,去除姿態角度異常、成像存在問題的圖像。將篩選后的圖像和POS 數據輸入Pix4D mapper 軟件中,根據相機的配置調整處理參數。運行軟件后,應用會自動生成連接點并與POS 數據參與空三運算,由此得到每一張航拍影像的準確外方位元素和加密點的坐標。經點云加密后,軟件會自動生成數字表面模型,以此作為數字微分糾正的依據對影像進行正射校正[23-24]。
1.4.1 研究技術路線
研究首先基于無人機多光譜特征圖譜分別選取4類作物形態的統計區域,提取并分析每一類統計區域的光譜、植被指數與紋理特征。然后經過特征篩選后分別構建適于倒伏提取的光譜、植被指數、單類紋理與多類紋理特征,利用以上4種特征組合分別結合最大似然法進行分類得到4種倒伏玉米面積的提取結果。最后驗證每種結果的Kappa系數和面積誤差,分析倒伏玉米面積提取精度,優選倒伏玉米面積提取的最佳方法。技術路線如圖2所示。
1.4.2 特征統計樣區選擇
根據目視解譯,依照玉米葉片的不同顏色將未倒伏玉米分為未倒伏A類和未倒伏B類;垂直特征不明顯,壟向性較差的玉米歸為倒伏玉米;黑色區域歸為陰影。從多光譜影像上均勻選取未倒伏A類2 113.76 m2、未倒伏B類2 189 m2、倒伏區域2 164.07 m2,3類作物形態分別選5處,由于圖上陰影區域分布零散且面積較小,因此陰影區域選擇50處,共計253 m2。4類統計樣區的選擇應均質且不能包含其他類別。其中,未倒伏A類主要為地塊東部受大風影響嚴重,葉表呈翠綠色,長勢較差的未倒伏區域;未倒伏B 類為地塊西部受大風影響較小,葉片表現為淡黃色特征,長勢較好的未倒伏區域;倒伏區域為地塊東部邊緣,受大風影響而發生倒伏的區域;陰影區多為未倒伏玉米在倒伏玉米區域上的投影以及南側防護林的陰影。

圖2 倒伏玉米面積提取技術路線圖Fig.2 Technology roadmap of maize lodging area
1.4.3 光譜特征選擇
分別計算4 類特征統計樣區中各個波段的反射率均值,對應各波段波長的中心位置,繪制每個類別的光譜反射率曲線,統計4 類區域在各波段上反射率的平均差異。雖然本研究使用的多光譜數據具有定量的光譜信息,但不同類別作物形態之間的反射率差異過小仍會產生較多的錯分像元[25-27],進而影響倒伏面積提取的精度,因此為增強各類區域之間的反射率差異,減少錯分像元數量,本研究根據計算結果,選取平均差異相對較大的藍、綠、紅3個波段作為倒伏面積提取的光譜特征指標。
1.4.4 植被指數法
歸一化植被指數(normalized difference vegetation index,NDVI)通過非線性拉伸的方式增強了紅波段與近紅外波段反射率之間的差異,是反映農作物生長狀況和植被覆蓋密度的重要指標[28-29];而紅邊歸一化植被指數(red edge normalized difference vegetation index, NDVIR-edge)則更能夠反映葉冠層的細小變化、植被片斷和衰老等信息,這些信息通過紅波段與紅邊波段反射率的非線性拉伸而得以體現[30-31]。NDVI與NDVIR-edge的計算公式:

式中BRed為無人機多光譜影像中的紅波段;BNir為近紅外波段;B717為紅邊波段,其中心波長為717 nm。
基于無人機獲取的多光譜影像分別提取NDVI 和NDVIR-edge,提取每類特征統計樣區中的植被指數,分別統計2種植被指數中每類作物形態的均值。
1.4.5 紋理特征選擇
對影像進行灰度共生矩陣的紋理濾波,得到各波段的均值、方差、協同性、對比度、相異性、信息熵、二階矩、相關性等8類,共40項紋理特性,并統計各類區域中每項紋理特征的均值[32-33]。以統計結果作為分析的原始數據,分別進行基于單類和多類紋理特征的倒伏面積提取特征選擇。
1)單類紋理特征選擇:在只選取單類紋理特征作為倒伏面積提取指標的條件下,為使提取的倒伏面積更加準確,實驗將基于8 類紋理特征,分別建立4 類作物形態的訓練樣區,利用最大似然法進行分類,根據分類結果選取適于區分4 類作物形態,倒伏面積提取精確的紋理特征。通過反復試驗,最終選取均值紋理特征。
2)多類紋理特征選擇:為了盡可能地降低紋理特征的維數,去除維數冗余所產生噪聲信息,本文首先選取均值、對比度2類紋理描述,作為倒伏玉米面積提取的基礎特征池,利用該特征池結合最大似然法進行分類,均勻選取倒伏和未倒伏玉米驗證樣本各50 個,計算提取結果的Kappa系數。然后,以該特征池為基礎,以一類紋理向量為步長,對該特征集進行維數疊加,驗證每次紋理向量疊加后分類結果的Kappa 系數,從而確定最優的分類維度,驗證過程如圖3 所示。根據圖3 的結果,本研究最終確定多類紋理特征組合的最優維數為15 維,共計3 類紋理特征。

圖3 不同維數下的Kappa系數變化Fig.3 Variation of Kappa coefficients at different dimensions
在確定用于分類的紋理特征之前,本研究將先對劃分的未倒伏A類、未倒伏B類、倒伏以及陰影區的各項紋理特征進行對比分析,并統計4 類作物形態之間各類紋理特征的相對差異度。為了使所選紋理特征保有最大類間差異和最低類內離散度,進而提高分類特征的魯棒性,本文選取了每2 類區域對比之下差異度較大的4 類紋理特征,分別為均值、對比度、二階矩和協同性紋理特征。任選3 類特征進行組合,驗證每種特征組合下倒伏玉米提取結果的Kappa 系數,選出Kappa 系數最高的紋理特征組合。經反復試驗,最后確定用于分類的紋理特征分別為均值、對比度和協同性。
1.4.6 倒伏面積提取和精度驗證
分別對NDVI、NDVIR-edge、單類紋理特征、多類紋理特征、光譜特征,共5 類特征集進行分類,由此獲得每一類特征集下的倒伏玉米提取結果,共計5種。分別統計5種分類結果中倒伏玉米的柵格數量,然后將柵格數與單個像元的面積作乘求得每種結果中的倒伏玉米面積。
為驗證不同特征組合的分類精度,根據目視解譯,選取4類典型倒伏作物形態的驗證樣區各4處,其中未倒伏A類3 599.26 m2、未倒伏B類3 563.09 m2、倒伏3 512.17 m2、陰影區域3 843.3 m2。選取的驗證樣區與特征統計樣區不能有重疊,驗證樣區的選取結果如圖4b,樣區選取后目視勾畫并統計各類別驗證樣區中的目標類別面積,以此作為各類別驗證的實測面積。
在影像中,倒伏玉米的線性種植結構和垂直生長特征不明顯,并且其葉片呈現淡黃色特征[34],結合目視解譯勾畫地塊內部所有的倒伏玉米邊界,解譯結果如圖4c 所示。經統計邊界內的倒伏面積為21 611.91 m2,以此作為實測值,驗證本文5種特征組合的倒伏面積提取精度。

圖4 樣區分布及實測倒伏區域Fig.4 Sample area distribution and measured lodging area
空間一致性是將指定位置的分類結果和驗證樣本所對應的類別進行比較,多采用混淆矩陣來度量[35]。而Kappa 系數能夠揭示景觀的空間變化信息,說明數量和位置信息的變換,全面反映分類結果的精度[36]。因此該研究采用Kappa 系數評價倒伏提取結果的空間一致性。計算公式如下:

式中N 是驗證樣本的總數;n 為混淆矩陣中的總列數;xii為混淆矩陣中第i行、第i列上的樣點數;xi+和x+i分別是第i行和第i列的總樣點數;K是Kappa系數。
計算每類特征統計樣區在各波段的反射率均值,繪制4 類作物形態的光譜反射率曲線,如圖5 所示,根據圖5,與未倒伏A 類相比,倒伏玉米由藍波段至近紅外波段的反射率平均上升0.04;與未倒伏B 類相比,倒伏玉米各波段的反射率平均上升0.01。前者的反射率上升幅度遠大于后者,這主要是未倒伏B 類玉米的各波段反射率大于未倒伏A 類玉米的原因。未倒伏A 類與B類之間在綠波段和紅邊波段上反射率數值差異要高于其他波段。陰影區域在藍、綠、紅波段上與其他3 類區域的反射率差異較大,在紅邊和近紅外波段上反射率差異較小。

圖5 4類典型作物形態的光譜反射率曲線Fig.5 Spectral reflectance curves of 4 typical crop morphology
由藍波段至綠波段,陰影的光譜反射率變化趨于平緩,而其他3 類作物形態的反射率變化趨勢相近,且變化速率較快。未倒伏A類與陰影在綠波段至紅波段上變化趨勢較緩,與未倒伏B類和倒伏形態的反射率變化趨勢差異較大。在紅波段至近紅外波段上,未倒伏A類、未倒伏B類和倒伏形態的光譜反射率變化趨勢大體相同,均呈現較快的變化速率。但陰影在紅波段至紅邊波段的反射率變化要大于其他3類作物形態,而在紅邊波段至近紅外波段上與其他3類作物形態的光譜曲線的變化基本相似。
計算特征統計樣區內各類作物形態的植被指數如表1,在兩類植被指數中陰影與其他3 類作物形態之間的差異最大。在NDVI 中未倒伏B 類與倒伏的數值差異最小,而在NDVIR-edge中,差異最小的是未倒伏A 類與未倒伏B 類。對于未倒伏的玉米而言,相機探測的反射率以葉片為主,而玉米在倒伏后,相機所探測的反射率以莖稈和葉片為主。由于莖稈中的葉綠素含量低于葉片,因此玉米發生倒伏后植被指數也會下降,其中NDVI 平均下降了0.08,NDVIR-edge下降0.07。

表1 4類典型作物形態的植被指數Table 1 Vegetation index of 4 typical crop morphology
計算特征統計樣區內每類作物形態的紋理特征均值,分析不同形態之間的相對差異度。根據表2,倒伏與2類未倒伏玉米的紋理特征之間有不同的相對差異度。在倒伏與未倒伏A 類玉米的紋理特征對比中,差異最大的是均值紋理,在該特征空間內,2 類作物形態的的相對差異度達1.11,比其他紋理特征平均高0.59;而在倒伏與未倒伏B類玉米的紋理差異分析中,相對差異最大的是對比度,在該特征中,2類作物形態的差異度為2.76,比其他紋理向量平均高0.8;與其他紋理特征相比,二階矩紋理特征在陰影與倒伏玉米、陰影與未倒伏A類玉米、陰影與未倒伏B類玉米的差異分析中,具有最高的區分度,其平均相對差異為31.87,比其他紋理特征平均高29.63。

表2 紋理特征的相對差異統計Table 2 Relative difference statistics of texture features
本研究在使用植被指數、多類、單類紋理特征以及光譜反射率特征集進行倒伏提取時,均采用最大似然法。圖6a 為利用多類紋理特征集合(均值、對比度以及協同性紋理)提取倒伏玉米的結果;圖6b 為使用單類紋理特征集(均值紋理)所得到的倒伏玉米結果;圖6c 為利用藍、綠、紅波段的光譜特征識別倒伏玉米的結果;圖6d 和6e分別為NDVI和NDVIR-edge特征的倒伏提取結果。

圖6 不同特征的玉米倒伏面積提取結果Fig.6 Extraction results of maize lodging area under different features
利用每類作物形態的驗證樣區,驗證不同特征組合的分類誤差,統計結果如表3 所示。根據表3 可知,使用多類紋理特征識別陰影和未倒伏A類形態的誤差稍高于植被指數和單類紋理特征識別方法,但在對其他作物形態進行分類時,多類紋理特征的分類誤差顯著低于光譜、植被指數和單類紋理特征識別方法。基于5 種特征組合所提取的倒伏結果如圖6 所示,根據此結果計算倒伏面積提取誤差,利用100 個倒伏和未倒伏驗證樣本,計算提取結果的Kappa系數,計算結果如表4所示。
由表4 可知,紋理類的特征組合提取倒伏結果的Kappa 系數最高,平均為0.61,比光譜和植被指數特征結果的Kappa系數分別高0.28、0.46。分析原因是玉米發生倒伏后其主要的變化是形態結構的變化,而光譜和植被指數特征變化相對較小,特別是完熟期未倒伏B 類玉米和倒伏玉米間的光譜和植被指數特征的物理差異很小,因此利用這兩類特征進行倒伏提取會產生較多的錯分,而紋理特征能夠很好地反映未倒伏和倒伏玉米間的空間形態差異,因此利用該類特征進行倒伏提取產生的錯分較少,空間一致度較高。
在5 種特征組合中,基于多類紋理特征進行倒伏提取的誤差最低,4 類典型作物形態的識別平均誤差為9.82%,Kappa系數最高,與單類紋理特征提取結果相比,其Kappa系數提高了0.47,倒伏提取誤差降低了27.57 個百分點。分析原因是單類紋理特征雖然能夠在一定程度上反映作物的空間域形態信息,但研究區內存在不同倒伏程度和不同生長狀況的玉米,它們均具有不同的形態特點,單類紋理的空間映射并不能很好地反映倒伏和未倒伏玉米間的形態差別,而多類紋理特征中所含有的多方向的空間映射信息能夠較好地體現倒伏和未倒伏玉米間的形態差異,因此具有較好的空間一致性。
利用光譜特征和NDVI 所提取的倒伏面積過大,其原因主要是大量的未倒伏B 類形態被劃分為倒伏,從而使倒伏像元數增多;而利用NDVIR-edge所得的倒伏面積偏大的主要原因是,未倒伏B 類形態和受邊緣效應影響的未倒伏玉米被錯分為倒伏;基于單類紋理特征進行分類,部分未倒伏B 類形態也被錯分為倒伏,但提取的倒伏面積反而過小,分析其原因是更多的倒伏玉米被劃分為未倒伏玉米。

表3 不同特征下4類典型作物形態的分類誤差Table 3 Classification error of 4 classes of typical crop morphology under different features

表4 不同特征下玉米倒伏面積提取誤差Table 4 Extraction error of maize lodging area under different features
本文利用無人機多光譜影像的植被指數、紋理和光譜特征結合最大似然法分類,對完熟期倒伏玉米面積的提取進行了定量分析,得到了以下結果:
1)地塊內土壤和作物的時空異質性會使完熟期的玉米處于不同的生長狀態,其中倒伏玉米與長勢較好的未倒伏玉米的葉片特征和光譜特征很相似,兩類玉米的葉片均呈現淡黃色,各波段的反射率平均差異僅為0.01,對這兩類作物形態進行分類會產生較多的錯分像元,降低倒伏面積提取的精度。單一反射光譜特征或植被指數特征不能準確地區別完熟期倒伏玉米區域。
2)利用高維多類紋理特征能夠顯著區別地塊內不同作物形態之間的差異,增強類內一致性,從而提高倒伏玉米的識別精度。針對完熟期的倒伏玉米識別,該特征組合具有較高的實用價值。
3)陰影區域的判斷是完熟期倒伏玉米精準識別的一個難點。未倒伏玉米地塊內部的陰影主要包括葉片間隙與防護林陰影兩部分,而倒伏發生后地塊內部不僅存在以上兩部分陰影,還包括未倒伏玉米投影到倒伏玉米上的陰影,這部分陰影集中分布在小面積倒伏區域內部以及片狀倒伏邊緣,應同樣算作倒伏面積。研究發現,除NDVI 和NDVIR-edge外,在選取合適的紋理濾波窗口的條件下,單類、多類紋理特征也能夠較好地剔除葉片間隙的陰影,但后者的提取精度更高。對于地塊邊緣的防護林陰影,本文的5 種特征組合均不能對其進行剔除,考慮可通過實地調查,確定防護林陰影為未倒伏區域之后,根據分類結果建立矢量文件,對陰影區進行手動去除。
本研究利用完熟期玉米的無人機多光譜圖像,分析了圖像中各類區域的植被、光譜以及紋理特性,針對光譜和紋理特征提取方法分別選取了最優的分類特征,然后利用優選的光譜、植被指數、單類、多類紋理特征向量集結合最大似然分類方法提取了倒伏玉米面積,通過分析分類和提取結果得到以下結論:倒伏玉米與未倒伏A 類玉米的光譜特征差異較大,倒伏玉米與未倒伏B 類玉米的光譜特征差異較小,光譜特征組合無法準確區分倒伏玉米和未倒伏B 類玉米,該特征組合所提取的倒伏玉米面積偏差較大,因此不適于完熟期倒伏玉米面積提取;植被指數特征會使受邊緣效應影響的未倒伏玉米錯分為倒伏玉米,在該類特征下,倒伏玉米提取結果的空間一致性較差。因此也不適于完熟期倒伏玉米面積提取;對比不同特征分類和倒伏面積提取結果發現,利用均值、對比度及協同性的多類紋理特征組合提取的倒伏玉米面積精度最高,在該類特征下分類的平均誤差為9.82%,倒伏面積提取誤差為3.40%,Kappa 系數為0.84,比其他特征組合的Kappa 系數平均高0.59,因此多類紋理特征組合可以實現完熟期倒伏玉米面積的準確提取;與光譜特征相比,利用植被指數和紋理特征進行分類能夠更好地剔除玉米葉片間隙的陰影,但仍無法準確判斷防護林陰影下的玉米是否為倒伏玉米;本研究初步證明了多類紋理特征能夠增強倒伏與未倒伏玉米間的空間形態差異,結合不同維數下倒伏提取結果的空間一致性分析能夠降低特征信息冗余,防止分類結果的過度擬合,進而提高倒伏提取結果的精度。但受當前無人機續航能力與通訊距離的制約,本研究還無法實現大范圍的倒伏面積提取。因此在今后的研究中應嘗試利用無人機圖像與衛星影像結合的方式來擴大倒伏信息的獲取范圍。