孫 盛 劉立露 胡忠文 余 旭

(1.廣東工業(yè)大學(xué)計算機(jī)學(xué)院， 廣州 510006; 2.深圳大學(xué)自然資源部大灣區(qū)地理環(huán)境監(jiān)測重點(diǎn)實(shí)驗室， 深圳 518000; 3.廣東工業(yè)大學(xué)土木與交通工程學(xué)院， 廣州 510006)

0 引言

甘蔗是我國主要的糖料作物，被農(nóng)業(yè)部列為九大主要農(nóng)作物之一，廣東省作為甘蔗的主要生產(chǎn)區(qū)，2019年種植面積達(dá)1 030 065.01 hm2，占據(jù)全國的12.28%[1]。因此，多維度地監(jiān)測甘蔗的生長變化規(guī)律，對廣東甘蔗產(chǎn)業(yè)的可持續(xù)發(fā)展具有較好的推動作用[2]。甘蔗株高作為甘蔗生長情況、產(chǎn)量的重要評估標(biāo)準(zhǔn)之一，選擇株高進(jìn)行反演對于甘蔗的生長監(jiān)測具有較好的代表意義。

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

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

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

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

1 研究區(qū)域和實(shí)驗數(shù)據(jù)

1.1 研究區(qū)域與研究作物概況

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

圖1 研究區(qū)域位置Fig.1 Location of study area

1.2 研究區(qū)域數(shù)據(jù)獲取

通過多次調(diào)研獲取了研究區(qū)域甘蔗生長周期內(nèi)23景時間序列Sentinel-1A數(shù)據(jù)及衛(wèi)星過境時間基本同步的3期實(shí)地數(shù)據(jù)采集；此外，在當(dāng)?shù)貧夂蚺c農(nóng)業(yè)氣象中心的支持下，補(bǔ)充了3個調(diào)研時間節(jié)點(diǎn)之間的甘蔗生長數(shù)據(jù)。數(shù)據(jù)獲取時間為2020年2月底(甘蔗播種期)至2020年12月底(甘蔗成熟期)。

1.2.1SAR數(shù)據(jù)

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

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

1.2.2實(shí)地調(diào)研數(shù)據(jù)

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

表1 實(shí)地調(diào)研時間點(diǎn)及現(xiàn)場圖Tab.1 Field research date and scene pictures

2 株高反演總體方案及雙極化雷達(dá)植被指數(shù)

2.1 總體反演方案

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

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

2.2 時間序列Sentinel-1A 數(shù)據(jù)預(yù)處理

由于SAR衛(wèi)星的成像體制，Sentinel-1A圖像數(shù)據(jù)含有相干斑噪聲[25]。因此首先需對Sentinel-1A原始圖像進(jìn)行預(yù)處理。

利用Python搭建歐洲航天局哨兵系列衛(wèi)星產(chǎn)品的綜合應(yīng)用平臺SNAP(Sentinel application platform)二次開發(fā)環(huán)境；結(jié)合SNAP平臺中的子模塊完成圖像的預(yù)處理。

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

k=[〈SVV〉 〈SVH〉]T

(1)

(2)

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

〈·〉——各向同質(zhì)情況下隨機(jī)散射介質(zhì)的空間平均

上標(biāo)*表示復(fù)共軛。

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

2.3 雙極化雷達(dá)植被指數(shù)DPRVI的計算

2.3.1Cloude-Pottier目標(biāo)分解與特征值歸一化

(3)

其中

Λ=〈|SVV|2〉

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

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

(4)

式中ci——秩為1的獨(dú)立矩陣

λi、ei——C2的實(shí)特征值與特征向量

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

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

(5)

2.3.2極化度(DoP)計算

BARAKAT[24]提出了DoP的表達(dá)式并將其用于全極化數(shù)據(jù)的NN協(xié)方差矩陣。針對雙極化數(shù)據(jù)進(jìn)行改進(jìn)，計算式為

(6)

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

根據(jù)極化度定義可知，它量化了不同散射機(jī)制之間的相對強(qiáng)度，當(dāng)DoP=1時，認(rèn)為目標(biāo)處于完全極化狀態(tài)，由單一散射機(jī)制占主導(dǎo)，當(dāng)DoP=0時，認(rèn)為目標(biāo)處于完全去極化狀態(tài)，由多次散射占主導(dǎo)，多種散射機(jī)制共存。

2.3.3DPRVI參數(shù)計算

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

DPRVI=1-DoP·P1

(7)

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

3 反演方法實(shí)現(xiàn)及實(shí)驗分析

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

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

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

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

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

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

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

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

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

3.2 不同極化參數(shù)對比分析

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

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

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

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

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

圖9 不同甘蔗田的極化參數(shù)(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 結(jié)論

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

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

(3)使用當(dāng)?shù)貧夂蚺c農(nóng)業(yè)氣象中心數(shù)據(jù)進(jìn)行對比分析，決定系數(shù)R2達(dá)0.839，平均絕對偏差為7.4%，驗證了Sentinel-1A時序數(shù)據(jù)對粵港澳大灣區(qū)的甘蔗作物的連續(xù)監(jiān)測是有效的，對于大規(guī)模、低成本的甘蔗作物生長狀態(tài)監(jiān)測，Sentinel-1A數(shù)據(jù)是一種較好的數(shù)據(jù)源。