基于Sentinel-1A雙極化時序數據的甘蔗株高反演方法

孫 盛 劉立露 胡忠文 余 旭

(1.廣東工業大學計算機學院， 廣州 510006; 2.深圳大學自然資源部大灣區地理環境監測重點實驗室， 深圳 518000; 3.廣東工業大學土木與交通工程學院， 廣州 510006)

0 引言

甘蔗是我國主要的糖料作物，被農業部列為九大主要農作物之一，廣東省作為甘蔗的主要生產區，2019年種植面積達1 030 065.01 hm2，占據全國的12.28%[1]。因此，多維度地監測甘蔗的生長變化規律，對廣東甘蔗產業的可持續發展具有較好的推動作用[2]。甘蔗株高作為甘蔗生長情況、產量的重要評估標準之一，選擇株高進行反演對于甘蔗的生長監測具有較好的代表意義。

傳統反演作物參數主要通過實地調研，需要耗費較多人力物力與時間，遙感技術以其在時空尺度上的優勢，逐漸在幾種大宗糧食作物如水稻[3-4]、玉米[5-6]、小麥[7-8]、油菜[9-10]等得到了廣泛應用，但應用于甘蔗的研究仍然較少。

目前，學術界較多采用光學遙感技術完成甘蔗種植信息的識別[11-12]、生物物理參數反演與估產[13]。但是，光學遙感圖像受天氣和氣候影響較大，我國甘蔗種植多在中國南部廣東、廣西地區，這些地區空氣濕度大，常年多云多雨，應用光學遙感技術進行甘蔗生長監測的效果較差。合成孔徑雷達(Synthetic aperture radar, SAR)遙感[14]具有全天時、全天候的特點，可以較好地適用于廣東地區農作物生長監測。近年來，基于SAR后向散射系數對甘蔗進行作物類型鑒定[15-16]和長勢參數株高[17]等生長監測的研究有所增加。但是，后向散射系數過早飽和的特點限制了對甘蔗生長后期的監測，因此，需進一步挖掘SAR的極化特征在甘蔗株高反演中的潛力。

目前，已有不少利用極化特征開展其他作物的研究，如CHANG等[18]通過對全極化PALSAR(L波段)衛星數據進行極化特征分解，推得極化雷達植被指數(Polarimetric radar vegetation index,PRVI)定量反演灌木叢生物量，該研究認為灌木叢生物量與PRVI的相關性(R2=0.75)優于雷達植被指數(Radar vegetation index,RVI)[19]的相關性(R2=0.5)。WANG等[20]基于全極化無人機合成孔徑雷達(Uninhabited aerial vehicle synthetic aperture radar，UAVSAR)數據的極化分解方法，為不同的散射機制定義SAR監測作物生長的植被生長指數。結果表明，散射特性可以隨作物類型和物候發展階段而變化，該植被生長指數與作物株高和生物量的地面測量值具有很好的相關性。盡管這些雷達植被指數可以很好地表征植被情況，但僅限于全極化或簡縮極化SAR衛星數據；另一方面，這些研究都是基于經濟成本較高的商業SAR衛星，不太適用于大范圍區域作物監測的應用。

本文采用歐洲航天局(European Space Agency, ESA)提供的免費科研模式的雙極化SAR Sentinel-1A時序數據，分析雙極化雷達植被指數(DPRVI)[21]對于株高反演的可行性；分析討論DPRVI對甘蔗不同生長期的響應，采用4種經典的經驗回歸模型定量反演甘蔗株高并對最佳反演模型進行驗證。將DPRVI與3種經典的極化參數的時間模型進行性能評價，以驗證本文提出的株高反演方法的效果和性能。

1 研究區域和實驗數據

1.1 研究區域與研究作物概況

研究區域位于廣東省廣州市南沙區，中心坐標東經113°56′，北緯22°78′，年降雨量為1 623.6～1 899.8 mm，2017年、2018年與2019年甘蔗種植面積分別達到6 600.00、5 973.34、4 913.34 hm2，是該區域種植面積最大的農作物[22]。在廣州市南沙區大崗鎮廟貝村進行研究區域的數據采集，如圖1所示。甘蔗生長周期一般分為發芽期、幼苗期、分蘗期、伸長期和成熟期5個時期。在2、3月播種，至當年11月或次年1、2月收獲結束。其中分蘗期5月至6月上中旬，伸長期6月下旬至10月，成熟期11月至次年1、2月。因甘蔗伸長期較長，將6月下旬至8月下旬分為伸長前期，8月下旬至10月分為伸長后期。

圖1 研究區域位置Fig.1 Location of study area

1.2 研究區域數據獲取

通過多次調研獲取了研究區域甘蔗生長周期內23景時間序列Sentinel-1A數據及衛星過境時間基本同步的3期實地數據采集；此外，在當地氣候與農業氣象中心的支持下，補充了3個調研時間節點之間的甘蔗生長數據。數據獲取時間為2020年2月底(甘蔗播種期)至2020年12月底(甘蔗成熟期)。

1.2.1SAR數據

為避免不同傳感器觀測配置對觀測結果的影響，選取Sentinel-1A衛星的IW模式一級產品(Level-1)中的單復視(Single looking complex, SLC)圖像，其重返周期為12 d、空間分辨率為5 m(距離向)×20 m(方位向)，運行軌道為升軌右視；單次成像時，完整條帶長約250 km，距離向上由3個子條帶構成，每個子條帶在方位向上有9個bursts并按方位向順序排列(圖2)。

圖2 Sentinel-1A IW 圖像Fig.2 Sentinel-1A IW image

1.2.2實地調研數據

Sentinel-1A時序數據覆蓋了研究區甘蔗的整個周期，地面實地調研時間分別為與衛星過境采集基本同步的3個時間節點：2019年11月24日、2020年9月29日與2020年12月22日(表1)。為保證足夠的純凈像元，減少其他作物的影響，即降低混雜像元，實驗區域要求種植面積不小于20 hm2，并通過北斗高精度分米級移動平臺Qmini A5(亞米級差分信號)采集經緯度等數據。每次調研獲取實驗區域12塊代表性地塊的甘蔗株高、伸長期、種植密度與現場天氣情況等；每個代表性地塊大小為2 m×2 m，以米尺測量地面至甘蔗自然狀態下最高點的距離作為甘蔗株高，將測量得到的每個地塊的株高平均值作為該地塊的甘蔗株高。

表1 實地調研時間點及現場圖Tab.1 Field research date and scene pictures

2 株高反演總體方案及雙極化雷達植被指數

2.1 總體反演方案

首先對時間序列Sentinel-1A數據進行預處理，將SAR散射矩陣轉換為協方差矩陣；再通過Cloude-Pottier目標分解方法[23]求得特征值與極化度(Degree of polarization, DoP)[24]，算出雷達雙極化植被指數DPRVI；對甘蔗株高h進行反演，得出甘蔗生長最佳反演模型。對于模型采用決定系數R2與均方根誤差RMSE進行性能評測，并將DPRVI與3種經典極化參數進行對比。總體設計方案見圖3。

圖3 株高反演總體方案Fig.3 Scheme of plant height inversion

2.2 時間序列Sentinel-1A 數據預處理

由于SAR衛星的成像體制，Sentinel-1A圖像數據含有相干斑噪聲[25]。因此首先需對Sentinel-1A原始圖像進行預處理。

利用Python搭建歐洲航天局哨兵系列衛星產品的綜合應用平臺SNAP(Sentinel application platform)二次開發環境；結合SNAP平臺中的子模塊完成圖像的預處理。

將圖像從散射矩陣轉換為2×2的協方差矩陣(C2)，該矩陣是一個二階散射信息量，由雙極化散射矩陣目標向量k[26]的空間平均生成，計算式為

k=[〈SVV〉 〈SVH〉]T

(1)

(2)

式中SVH、SVV——雙極化數據VH、VV極化通道的散射矩陣

〈·〉——各向同質情況下隨機散射介質的空間平均

上標*表示復共軛。

在距離向和方位向上以4∶1進行多視處理[27]，減弱相干斑噪聲影響和減少后續數據處理量；采用refined Lee filter 方法進行極化濾波處理；最終，將這些像素地理編碼到UTM投影坐標系中。

2.3 雙極化雷達植被指數DPRVI的計算

2.3.1Cloude-Pottier目標分解與特征值歸一化

(3)

其中

Λ=〈|SVV|2〉

α=〈|SVV|2〉/〈|SVH|2〉

由式(3)可知，協方差矩陣C2為半正定的Hermite矩陣，故Cloude-Pottier分解可表達為

(4)

式中ci——秩為1的獨立矩陣

λi、ei——C2的實特征值與特征向量

從式(3)、(4)得到協方差矩陣的兩個特征值為λ1=Λ,λ2=Λα。

由定義可知，特征值量化了主要的散射機制。λ1與VV極化通道相關，代表奇偶次散射強度；λ2與VH極化通道相關，代表多次散射(體散射)強度。當目標由奇偶次散射占主導時，有λ1?λ2。將λ1歸一化處理得到歸一化參數P1，計算式為

(5)

2.3.2極化度(DoP)計算

BARAKAT[24]提出了DoP的表達式并將其用于全極化數據的NN協方差矩陣。針對雙極化數據進行改進，計算式為

(6)

式中 Tr()——矩陣的跡操作符

根據極化度定義可知，它量化了不同散射機制之間的相對強度，當DoP=1時，認為目標處于完全極化狀態，由單一散射機制占主導，當DoP=0時，認為目標處于完全去極化狀態，由多次散射占主導，多種散射機制共存。

2.3.3DPRVI參數計算

隨著作物冠層的發育，散射機制的類型也變得復雜。作物的發育初期，主要為土壤表面的漫散射(奇次散射類型之一)；發育晚期則由來自于冠層和土壤的多次散射占主導,散射類型的隨機性增加。因此，從作物發育初期到晚期，DoP與P1均呈下降趨勢。MANDAL等[21]發現這兩個參數對作物生長狀態的表征存在一定差異性，基于這兩個參數得到一個新的參數，即雙極化雷達植被指數(DPRVI)[8，28]，并成功用于反演油菜花與小麥的葉面積指數(Leaf area index,LAI)和植物含水率。另一方面，文獻[28-32]均證明了葉面積指數與株高有較好的相關性，對于甘蔗，通過對比何亞娟等[13]的甘蔗不同生育期LAI分析與本文的實驗結果，可發現甘蔗伸長前期與株高具有很好的正相關性，即株高越高，葉面積越大。但隨著后期葉片衰老、葉面積減少，株高并不降低甚至在增高，二者的正相關性不明顯。綜上，本文將DPRVI用于甘蔗伸長前期狀態的監測并反演株高；同時對該參數在甘蔗伸長后期的株高反演進行初步分析和討論。DPRVI的計算式為

DPRVI=1-DoP·P1

(7)

從定義可知，DPRVI∈[0,1]，DoP和P1的乘積對應于奇偶次散射的比例，通過單位1減去該比例可以推得目標生長過程的多次散射變化，即散射隨機程度。理論上，當目標完全隨機散射時，DoP=0(完全去極化)，λ1=λ2，即P1=0.5，此時DPRVI=1，等價于作物冠層完全發育狀態下的散射特征；在光滑的裸露表面，散射現象主要為布拉格散射時，有DoP=1(完全極化)，λ1=λ1+λ2，λ2=0。此時DPRVI=0，等價于裸土的散射特征；數值在0～1之間的變化反映了作物生長的過程。

3 反演方法實現及實驗分析

3.1 雙極化雷達植被指數(DPRVI)對甘蔗株高(h)的反演

通過實地數據采集和咨詢當地氣候與農業氣象中心及蔗農，得到不同生長期甘蔗平均株高，見圖4。分析發現，株高隨生長時間變化明顯，通過比較，以Sigmoid函數進行擬合效果最佳，決定系數為0.980。從擬合曲線可看出，從發芽開始到伸長期前期甘蔗呈加速生長趨勢，株高迅速增加，這是因為此期間光合作用逐漸增強，在此階段，甘蔗生長特點為葉面積迅速增大，節間伸長增粗；進入伸長后期，甘蔗拔節變緩，株高增加變緩，這是由于此期間甘蔗吸收的營養轉向蔗莖的增粗與糖分的積累。

圖4 甘蔗不同生長期株高的Sigmoid函數擬合模型Fig.4 Sigmoid function fitting model for plant height of sugarcane at different growth stages

基于2—12月的時間序列Sentinel-1A數據，將株高擬合曲線與對應研究區域3個甘蔗田(甘蔗田1、甘蔗田2、甘蔗田3)的植被指數DPRVI繪制于圖5。其中，甘蔗田1、甘蔗田2與甘蔗田3的種植面積分別約為3.33、4.00、3.67 hm2。通過實地測量，3塊甘蔗田的平均行距約為100 cm，平均株距約為30 cm。

圖5 甘蔗h與DPRVI隨播期的動態變化曲線Fig.5 Dynamic changes of h and DPRVI with sowing date of sugarcane

由圖5可知，在2、3月DPRVI的值較低，這是由于南沙甘蔗播種采取覆膜技術，田地經整土覆膜后以有序的奇偶次散射λ1為主導；隨后，甘蔗進入幼苗期，體散射分量λ2增加，DoP降低，但仍以奇偶次散射為主導(λ1?λ2)，DPRVI從3月到4月下旬的緩慢上升也驗證了這一變化；而后DPRVI加速上升，于6月中旬達到峰值，這時甘蔗到了分蘗期，分蘗加快，冠層的快速發育使得體散射分量λ2進一步增加，此時田地表現為多種散射類型并存，散射隨機程度最高；之后DPRVI下降，是由于甘蔗進入伸長前期，甘蔗分蘗變慢，節間加速分節與增粗，葉面積迅速增大，使得奇偶次散射分量λ1增加；到8月下旬后，DPRVI緩慢上升則是由于甘蔗進入伸長后期，分節與增粗變緩，葉片數積累增多，甘蔗冠層逐漸完全發育，體散射分量λ2增加，奇偶次散射分量λ1減少，DoP變小；進入成熟期后，DPRVI或因甘蔗葉枯萎掉落減低體散射分量λ2而有所降低，數值不再有較大波動。綜上，圖5中DPRVI的變化規律較好地反演了甘蔗的物候周期。在甘蔗生長期的分蘗期之前與伸長后期之后兩者呈正相關關系，伸長前期兩者呈負相關關系。鑒于以上分析，采用分段函數與4種經典經驗回歸模型(線性、二次多項式、指數、對數)針對不同生育期進行株高回歸模型反演，以決定系數R2和均方根誤差(RMSE)評價反演性能，結果如表2所示。

由表2可以看出，DPRVI與株高有較顯著相關性，綜合考量計算量、R2與RMSE性能，以二次函數模型反演效果最好，決定系數分別為0.882、0.625與0.716；將不同生育期的回歸曲線繪制如圖6所示，可以看出，幼苗期、分蘗期的相關性最佳，其次為伸長后期與伸長前期。

為進一步驗證模型的可靠性，使用甘蔗田4(面積大約為20 hm2，平均行距約為100 cm，平均株距約為30 cm，位置見圖1)作為驗證數據，選擇預測效果最佳的幼苗期、分蘗期二次函數模型進行比較(圖7)。可以看出，相對于甘蔗田1、甘蔗田2和甘蔗田3，甘蔗田4的數據離散程度略大，這是由于甘蔗田4較少的耕種面積，受到了更多混雜像元的影響。接著進行預測值與實測值的比較(圖8)，求得決定系數R2為0.839，平均絕對偏差為7.4%，說明利用該模型依然可以較好地反演甘蔗株高與DPRVI的關系，基本滿足DPRVI對伸長前期(分蘗期之前)的甘蔗株高監測的動態需求，保證利用雙極化雷達植被指數預測甘蔗伸長前期株高的準確性與可靠性。

表2 甘蔗各生長期h與DPRVI的關系Tab.2 Relationship between h and DPRVI in sugarcane growth period

3.2 不同極化參數對比分析

繪制2—12月的時間序列Sentinel-1A數據，求得3個甘蔗田研究區域的DPRVI、交叉極化與同極化通道的后向散射系數強度(σVH°、σVV°)、不同極化比值(σVH°/σVV°)、雷達植被指數(RVI)的時間模型，如圖9所示。

由圖9可知，3種經典的極化參數σVH°、σVV°、σVH°/σVV°與RVI的時間模型在3個甘蔗田的伸長前期均表現出一定的生長響應。其中σVH°、σVV°與h呈正相關性，σVH°/σVV°、RVI與h呈負相關性。但是，從甘蔗進入伸長期開始，該3種參數均已經趨向穩定，難以反演伸長期之后株高的生長變化。這是由于甘蔗生長早期散射機制相對簡單，而到了后期，由于隨著甘蔗冠層的發育分蘗，散射機制復雜性增加，后向散射系數達到飽和，基于后向散射系數的σVH°、σVV°、σVH°/σVV°與RVI也無法較好地表征株高，而基于極化特征求得的DPRVI在甘蔗不同生育期仍然保持了與株高較好的相關性。因此，在反演甘蔗作物時，DPRVI比Sentinel-1A后向散射系數及其推導得出的參數性能更加優異，對于指導甘蔗生產具有一定應用價值。

圖6 甘蔗不同生育期株高最佳反演模型Fig.6 Optimal inversion model of plant height in different growth stages of sugarcane

圖7 伸長前期株高反演模型與實測甘蔗株高的關系Fig.7 Relationship between inversed plant height at early growth stage and measured plant height of sugarcane

圖8 伸長前期甘蔗株高實測值與預測值比較Fig.8 Comparison of measured and predicted height of sugarcane plant in early growth period

圖9 不同甘蔗田的極化參數(DPRVI、σVH°、σVV°、σVH°/σVV°、RVI)時間模型Fig.9 Comparison of polarization parameter (DPRVI,σVH°,σVV°,σVH°/σVV° and RVI) time models in different sugarcane fields

4 結論

(1)將時序數據Sentinel-1A求得的DPRVI與3種經典極化參數進行對比。結果顯示DPRVI能較好地反演甘蔗生長前期的株高變化。

(2)采用4種經典的經驗回歸模型對甘蔗株高反演建模。結果顯示分蘗期之前兩者相關性最高，二次多項式模型擬合效果最好，決定系數R2達0.882，均方根誤差達0.118 cm。

(3)使用當地氣候與農業氣象中心數據進行對比分析，決定系數R2達0.839，平均絕對偏差為7.4%，驗證了Sentinel-1A時序數據對粵港澳大灣區的甘蔗作物的連續監測是有效的，對于大規模、低成本的甘蔗作物生長狀態監測，Sentinel-1A數據是一種較好的數據源。