梁 好，徐維新*，段旭輝，張 娟，代 娜，肖強智，王淇玉

1. 成都信息工程大學(xué)資源環(huán)境學(xué)院，四川 成都 610225 2. 青海省氣象科學(xué)研究所，青海 西寧 810001

引 言

青藏高原獨特的生態(tài)環(huán)境與自然氣候條件是推動全球變化研究與認(rèn)識深化的重點區(qū)域。 其中，高寒草地是青藏高原最主要的植被覆蓋類型以及我國大規(guī)模生態(tài)環(huán)境保護與建設(shè)的對象。 冬季是草原過度采食與生態(tài)破壞的主要時期，冬季牧草即枯草的存量則是保護生態(tài)平衡與畜牧業(yè)生產(chǎn)與災(zāi)害防御的關(guān)鍵參量。 已有工作表明，枯草是生物量、水分、葉綠素等參數(shù)持續(xù)協(xié)同變化的結(jié)果，對其關(guān)鍵參數(shù)的認(rèn)識具有重要的科學(xué)意義與直接的應(yīng)用需求。 然而，當(dāng)前國內(nèi)外對冬季枯草的研究報道極少，而青藏高原地區(qū)冬季枯草監(jiān)測空白的現(xiàn)狀，則與枯草性狀參數(shù)認(rèn)識的不足有直接關(guān)系。

傳統(tǒng)草地監(jiān)測主要依靠人工采樣完成，不僅費時費力，更難以完成大面積的監(jiān)測。 近年來，輻射傳輸模型PROSAIL由于可以定量描述色素含量、水分、葉面積指數(shù)等性狀參數(shù)對反射光譜的影響，而被廣泛用于植被反演監(jiān)測和機理研究中[1]。 多項研究指出PROSAIL模型在現(xiàn)實環(huán)境、模型復(fù)雜性、準(zhǔn)確性和計算時間之間達(dá)成了良好的統(tǒng)一[2-3]，使其適用于多視角高光譜數(shù)據(jù)的分析和反演。 Berger等[4]對1992年到2017年間281篇與PROSAIL相關(guān)的論文分析發(fā)現(xiàn)，利用PROSAIL模型針對多種農(nóng)作物(以小麥、玉米、甜菜為主)、森林及天然植被(草地，灌木)展開的研究，均取得了成功應(yīng)用的同時，獲得并提升了對各類植被反射光譜特征與植被性狀參數(shù)作用的認(rèn)識。 梁順林等[5]指出準(zhǔn)確的參數(shù)運行區(qū)間是基于PROSAIL模型開展相關(guān)研究與反演活動的基礎(chǔ)。 然而，目前針對衰落期植被參數(shù)及其光譜特征的研究鮮有報道，個別研究僅作為生長期對照組而涉及衰落期[6-7]，缺乏對衰落過程中的各參數(shù)變化與機理狀態(tài)的系統(tǒng)研究。

研究將利用PROSAIL模型，結(jié)合實測光譜與性狀參數(shù)，通過模擬枯草冠層反射率曲線，探究PROSAIL模型對枯草的適用性，完成模型敏感性與不確定性分析。 在此基礎(chǔ)上率定高寒枯草各項性狀參數(shù)及其取值范圍，最終得到LAI，LAD和Cm等關(guān)鍵參數(shù)取值區(qū)間參考表。 為利用模型對枯草進行研究和反演提供有利的基礎(chǔ)數(shù)據(jù)支撐，并期望推動對枯草性狀參數(shù)的認(rèn)識與理解，為冬季枯草的遙感監(jiān)測提供基礎(chǔ)與方法支持。

1 實驗部分

1.1 樣本

樣本采集于青海省海北藏族自治州西海鎮(zhèn)中國氣象局海北牧業(yè)氣象試驗站。 該試驗場位于100°51′33″E，36°57′33″N，海拔3 140 m，地處青藏高原北部青海湖畔，屬典型高寒天然草原。

在該試驗場沿南北方向平行設(shè)置5行23個50 cm×50 cm枯草觀測固定樣方。 在牧草自然衰落并呈完全枯干狀態(tài)的2020年11月27日進行觀測采樣。 采樣時晴朗無風(fēng)，枯草反射光譜及性狀參數(shù)采樣時間限定于9:30—15:30。 最終獲得反射光譜230條，參數(shù)數(shù)據(jù)207條。

1.2 數(shù)據(jù)采集

1.2.1 光譜數(shù)據(jù)

使用ASD公司的FieldSpec4便攜式地物光譜儀進行枯草反射光譜采集。 波長范圍為350～2 500 nm，在可見光(VIS)和近紅外(NIR)波段采樣間隔為3 nm，短波紅外(SWIR)采樣間隔為8 nm，單次采集時間為0.2 s，探頭視場角θ為25°。 按操作規(guī)范預(yù)熱儀器并進行暗電流測試和參考板校正，為保證視場范圍半徑處于樣方區(qū)內(nèi)，觀測時將探頭高度保持在距地面約80 cm[8]進行觀測。 每個樣方每次同步采集10條冠層反射光譜，以保證觀測數(shù)據(jù)可靠性與穩(wěn)定性。

1.2.2 枯草性狀參數(shù)

使用KonicaMinolta公司的SPAD502葉綠素儀進行枯草葉綠素含量同步觀測。 該儀器測量窗口2 mm×3 mm，精度±1.0 SPAD, 采樣間隔2 s，重現(xiàn)性±0.5 SPAD內(nèi)。 測量時每個樣方隨機選擇3株枯草，在枯草的上中下三部分各進行一次測量，取平均值作為該株枯草的葉綠素相對含量SPAD值。

使用直尺測量冠層高度L、平均葉片長度D，并使用量角器量取葉傾角。 隨后用剪刀齊根剪下樣方中所有牧草并裝袋送回實驗室。 同步記錄觀測時間、天氣等狀況。 根據(jù)式(1)計算得到熱點參數(shù)Hspot。

(1)

使用天平(可讀性0.01 g)測量枯草鮮重。 隨后將同一樣方的葉片均勻、分散平鋪于多張A3紙上，經(jīng)拍照及后期圖像識別，獲取樣方葉面積指數(shù)(LAI)。 隨后放入烘箱以70 ℃恒溫烘干，烘至重量不再改變時作為該樣方干重。 鮮重與干重相減得到枯草水分含量。

1.3 光譜模擬

1.3.1 PROSAIL模型

PROSPECT葉片光學(xué)傳輸模型與SAIL冠層雙向反射率模型耦合組成PROSAIL模型。 Féret等[9]指出，最新推出的PROSPECT-D優(yōu)于之前的版本，不僅降低了模型預(yù)測的不確定性、更好地恢復(fù)光合色素，同時能在可見光范圍內(nèi)以最小的誤差模擬真實的葉片光學(xué)特性。 4SAIL模型引入了熱點效應(yīng)和土壤的二向反射特征, 提高了模型的數(shù)值穩(wěn)定性和魯棒性。 因此, 本研究采用PROSPECT-D與4SAIL耦合的PROSAIL模型進行枯草光譜模擬。

根據(jù)實測枯草性狀數(shù)據(jù)，同時參考LOPEX93數(shù)據(jù)庫與已有關(guān)于人工落葉林和冬小麥衰落期的研究成果[6-7]預(yù)設(shè)了模型運行初值(表1)。

1.3.2 模型敏感性分析

敏感性分析是定量或定性分析模型參數(shù)對模型結(jié)果產(chǎn)生的影響程度，對模型輸入?yún)?shù)進行敏感性分析是進行參數(shù)率定的前提。 本研究采用全局敏感性分析方法EFAST法來完成PROSAIL輸入?yún)?shù)敏感性分析。 EFAST法可以通過方差分解得到各參數(shù)的方差貢獻(xiàn)比重，并以敏感性指數(shù)值來衡量各參數(shù)的敏感性。 其優(yōu)點在于所得敏感性指數(shù)充分考慮了參數(shù)間耦合作用[10]。 其計算式為

(2)

式(2)中，V(Y)為全部參數(shù)的總方差，k為參數(shù)個數(shù)，Vi為第i個參數(shù)自身的方差，Vij為第i和j參數(shù)耦合作用的方差，V1…k為所有參數(shù)耦合作用的方差。

1.3.3 模型不確定性分析

不確定性是模型運行中參數(shù)間存在耦合作用無法完成全局最優(yōu)判斷導(dǎo)致的輸出結(jié)果離散[11]。 不確定性分析可以識別模型不確定性的來源并量化參數(shù)對輸出不確定性的貢獻(xiàn)。 進行不確定性分析有助于提高枯草各參數(shù)組分間相互作用的認(rèn)識和枯草衰落機理的理解。

采用全局敏感性指數(shù)與局部敏感性指數(shù)的差值來量化不確定性[12]。 如果某個參數(shù)與其他參數(shù)間存在耦合作用，造成了模型不確定性，那么它的全局敏感性指數(shù)值會大于局部敏感性指數(shù)值。 由式(2)知，兩者的差值可寫作Vij+V1…k，即為耦合關(guān)系指數(shù)，作為量化參數(shù)不確定性的指標(biāo)。

1.3.4 枯草高敏感性參數(shù)率定

枯草反射光譜受不同草品種、倒伏狀態(tài)、枯黃程度等影響，在光譜波形相似的情況下反射率絕對值波動較大，因此選取注重維度差異而不是數(shù)值差異的率定函數(shù)是必要的。

余弦距離[13]又稱余弦相似度，衡量維度間取值方向的一致性，同時修正了變量間可能存在的度量標(biāo)準(zhǔn)不統(tǒng)一問題。 其不同于以距離為測度的歐式距離，也不同于以單位變化相似程度為依據(jù)的相關(guān)系數(shù)，在具有時序特性并對趨勢(方向)敏感的曲線相似度評價中更有優(yōu)勢。 本研究采用余弦距離為參數(shù)率定評價函數(shù)。 通過對比不同輸入?yún)?shù)取值下的模型模擬反射光譜與實測反射光譜間的余弦距離完成關(guān)鍵參數(shù)的率定。

2 結(jié)果與討論

2.1 枯草光譜模擬結(jié)果

2.1.1 光譜數(shù)據(jù)預(yù)處理

對野外實測的固定樣方光譜數(shù)據(jù)，逐一檢查每個樣方同步觀測得到的10條光譜，剔除因儀器、吹風(fēng)等因素引起的明顯異常波動數(shù)據(jù)后，得到枯草原始光譜數(shù)據(jù)序列，每個樣方隨機挑選1條原始光譜展示于圖1中。 由于PROSAIL模型模擬的有效范圍為400～2 500 nm，本研究處理與分析均基于此區(qū)間完成。

從圖1可以看出，枯草光譜在400 nm附近差異較小，隨著波長的增大，枯草光譜間的反射率差異逐漸變大，并在1 300 nm附近達(dá)到反射率峰值(0.25～0.46)。 從1 300 nm開始，光譜反射率差異始終保持在0.1～0.25間，各樣方曲線分布區(qū)間分散，區(qū)分度明顯。

圖1 2020年11月27日海北試驗場實測枯草原始反射光譜Fig.1 Measured withered grass spectrum (MW) on November 27, 2020 in Haibei

為了更好凸顯枯草光譜特征，消除不同光柵間的系統(tǒng)性“跳躍”偏差及個別光纖波動導(dǎo)致的“毛刺”性異常值。 利用ASD光譜儀附配的ViewSpecPro軟件進行光譜預(yù)處理。 完成拋物線修正及同次觀測樣本的均值合成輸出，實現(xiàn)數(shù)據(jù)校正的同時降低測量隨機誤差。

2.1.2 模型運行

根據(jù)1.3.1設(shè)定的模型模擬參數(shù)初始值，設(shè)定了較寬的模型運行參數(shù)閾值范圍(表1)，以保證模擬結(jié)果處于可能的波動范圍區(qū)間。

其中LAI，LAD，Cw，Cm，Cab和Hspot根據(jù)實測數(shù)據(jù)確定。 例如：Cm實測值在0.003～0.020之間，平均值為0.010，因此設(shè)定Cm參數(shù)運行范圍為0.002～0.050。Cbr，Car，Cant，N初值參考文獻(xiàn)[6-7]設(shè)定，OZA，SZA和RAA由試驗點坐標(biāo)和時間確定。

利用MATLAB對PROSAIL模型編程運行。 根據(jù)表1參數(shù)范圍，以蒙特卡洛采樣[14]方式獲得在參數(shù)范圍內(nèi)隨機分布的15 000組參數(shù)組合對，模型運行輸出得到對應(yīng)的15 000組模擬反射光譜。

表1 模型運行初值及范圍設(shè)定Table 1 Initial value of PROSAIL model

2.1.3 基于光譜特征的枯草初級篩選

圖2為試驗場夏季(8月30日)綠草和冬季(11月27日)枯草樣方經(jīng)預(yù)處理后的平均光譜。 從圖中可以看出，高寒冬季枯草與綠草光譜在可見光與近紅外波段存在著顯著差別。 綠草光譜表現(xiàn)為典型的綠色植被光譜特征，紅光波段(650～680 nm)的吸收谷與近紅外波段(760～1 300 nm)的高反射特征清晰。 然而，枯草的光譜特征在短波紅外波段前則明顯不同于綠草，其反射光譜曲線自400～1 300 nm基本呈明顯的線性增加趨勢，綠色植被的藍(lán)光與紅光波段的吸收谷特征完全消失，近紅外波段較一致的高反射率特征被逐漸增加的線性分布特征代替。 反映了牧草干枯、葉綠素流失、葉片內(nèi)部結(jié)構(gòu)破壞后的獨特光譜特征。

圖2 海北試驗場8月30日和11月27日樣方平均光譜Fig.2 The average spectrum for all samples on August 30and December 30 respectively in Haibei

由于模型設(shè)定的參數(shù)閾值范圍較寬且輸入?yún)?shù)隨機組合，PROSAIL模型生成的15 000組模擬光譜序列中，存在著反映綠草特征以及由綠至枯過渡階段特征甚至“錯誤”的光譜。 因此，需要依據(jù)枯草的光譜特征，進行非枯草特征光譜數(shù)據(jù)的篩選排除。

觀察圖2可以發(fā)現(xiàn)，枯草在可見光波段與綠草有顯著光譜差異，尤其是紅光波段吸收谷的存在與否，可作為綠/枯草區(qū)分與判別劃界依據(jù)。 由于綠草紅光波段吸收谷與綠光波段小反射峰的存在，綠草綠光波段反射值必定大于紅光波段，而枯草則是紅光波段大于綠光。 因此，模擬結(jié)果中枯草光譜的篩選可通過式(3)判別

Rred-Rgreen>0

(3)

式(3)中：Rred為紅光波段，本研究選取660 nm反射率；Rgreen為綠光波段，選取560 nm反射率。

當(dāng)該式計算值大于0時，視該模擬光譜值反映了枯草光譜特征。 以此完成15 000組模擬數(shù)據(jù)的初級篩選, 并得到枯草相關(guān)光譜序列。

2.1.4 基于線性擬合的枯草光譜二級篩選

根據(jù)圖2枯草光譜特征，枯草與綠草光譜在1 300～2 500 nm波段雖然數(shù)值差異較大，但波譜特征相似，難以區(qū)分，而在400～1 300 nm波段二者間數(shù)值與分布特征均有顯著差異。

此外，圖1中1 300～2 500 nm波段枯草光譜分布區(qū)間分散，最大值與最小值的極值分布區(qū)間寬，與非枯草地物光譜區(qū)間重疊范圍較大。 而400～1 300 nm波段，枯草的極值分布區(qū)間窄而收斂性強，波形特征與其他地物具有顯著差異性。 同時，該波段處于大氣窗口，也不屬于水汽的強吸收帶，受水汽的影響較小。 因此從枯草400～1 300 nm波段波形特征出發(fā)，進一步區(qū)分枯草光譜。

從經(jīng)過數(shù)據(jù)預(yù)處理后400～1 300 nm波段實測枯草反射光譜可以發(fā)現(xiàn)(圖3)，這一區(qū)間枯草光譜呈準(zhǔn)線性分布，波形線性特征顯著。 對這一區(qū)間內(nèi)實測光譜數(shù)據(jù)序列的線性擬合結(jié)果表明，絕大多數(shù)觀測樣本線性擬合方程的決定系數(shù)R2達(dá)到了99%置信區(qū)間，所有23個觀測樣本的線性擬合R2均通過95%置信區(qū)間。 因此，以400～1 300 nm區(qū)間模擬光譜數(shù)據(jù)序列線性擬合方程系數(shù)R2>0.95作為枯草光譜二級篩選的方法與閾值。

圖3 經(jīng)預(yù)處理的400～1 300 nm波段實測枯草光譜Fig.3 The measured spectrum of withered grass from400 to 1 300 nm after data preprocessing

2.1.5 枯草光譜分布區(qū)間模擬

根據(jù)式(3)及2.1.4中提出的線性擬合方法，對PROSAIL模擬輸出的15 000組數(shù)據(jù)分步進行枯草光譜篩選，最終得到1 799條完成兩級篩選符合枯草特征的光譜數(shù)據(jù)序列。 挑出這些序列中的最大值與最小值，構(gòu)建完成400～2 500 nm枯草模擬光譜最大值、最小值與實測平均值光譜數(shù)據(jù)序列(圖4)。

該最大與最小值指示了指定波段枯草反射光譜響應(yīng)的閾值區(qū)間。 在400 nm處光譜反射率在0.01～0.02波動，差距微??；1 300 nm處其極值變動區(qū)間在0.24～0.42，閾值變幅0.18，而在1 450 nm水分強吸收谷波動在0.12～0.27。

從400～2 500 nm的全光譜序列變化特征看，模擬光譜與平均實測光譜顯著相似。 對篩選得到的模擬光譜與平均實測光譜進行擬合評價，其R2值在0.904～0.994之間，通過了0.01顯著性水平檢驗。 說明PROSAIL模型模擬冬季高寒牧草的效果準(zhǔn)確可信，篩選出的模擬光譜及對應(yīng)參數(shù)可用于枯草特征分析。

圖4 二級篩選得到的枯草模擬光譜序列極值

2.2 模型敏感性分析

采用EFAST方法進行全局敏感性分析，將2.1.2中的參數(shù)輸入組合對和產(chǎn)生的模擬光譜結(jié)果輸入EFAST分析得到各輸入?yún)?shù)的敏感性指數(shù)值。

圖5為10個枯草性狀參數(shù)在400～2 500 nm的全局敏感性指數(shù)百分比面積堆疊曲線。 從圖5可以看出： 在VIS波段(400～700 nm)，LAI，LAD和Cab對枯草光譜影響較大。 在NIR波段(700～1 350 nm)，Cm和LAD平均葉傾角影響顯著，影響幅度均接近50%，而在這一波段N也有持續(xù)少量影響，值得注意的是，Cbrown僅在700～800 nm之間有近40%的影響。 在SWIR波段(1 350～2 500 nm)，主要影響因素依次為Cm，LAI，Cw及LAD。 其中Cm除了在兩個水分吸收帶指數(shù)值較低外，其余波段均維持在50%左右，Cw則僅對兩個典型水分吸收谷(1 450，1 940 nm)有突出影響，LAD的影響程度則隨波長逐漸變小。

需要注意的是，雖然枯草內(nèi)Cab含量非常低，是枯草低敏感參數(shù)，但由于其在牧草衰落過程中具有直接指示意義，在此也將其劃入高敏感性參數(shù)。

根據(jù)以上分析結(jié)果，將輸入?yún)?shù)分為高敏感性參數(shù)(LAI，LAD，Cm，Cab，Cw)和低敏感性參數(shù)(Cbr，Car，Cant，N，Hspot)。 通過優(yōu)化高敏感性參數(shù)的取值范圍，獲得參數(shù)合理閾值區(qū)間。 由于低敏感性參數(shù)對模型輸出的影響較小，幾乎可以忽略不計，將低敏感性參數(shù)的取值固定，取值為篩選出的模擬枯草光譜序列對應(yīng)的參數(shù)平均值。

參照圖4枯草光譜極值分布范圍，基于敏感性分析結(jié)果，分別統(tǒng)計5個枯草敏感參數(shù)的光譜值區(qū)間內(nèi)數(shù)值分布，準(zhǔn)確界定其數(shù)值分布范圍，逼近理想閾值區(qū)間的同時，實現(xiàn)了模型初始運行時相對寬泛閾值區(qū)間的縮小與優(yōu)化。

枯草參數(shù)敏感性分析優(yōu)化的結(jié)果如表2所示。 從表2與表1的對比可以看出，LAI由0.2～2.5變動區(qū)間縮小至0.2～1.5，說明枯草葉面積指數(shù)最高閾值界限低于綠色植被的同時，其閾值范圍明顯偏窄。 同樣，LAD的下限閾值由0提高到11，Cab的閾值范圍由0～35縮窄聚集至0.03～25。 枯草敏感的5個參數(shù)閾值區(qū)間得到了較精準(zhǔn)化的界定，為最終實現(xiàn)枯草參數(shù)率定進一步縮小了范圍。

2.3 模型不確定性分析

由1.3.3可知，耦合關(guān)系指數(shù)可指示模型各參數(shù)間耦合作用的強弱，籍此判斷各輸出參數(shù)隨其他參數(shù)變化而波動的不確定性。

表2指出LAI，LAD與Cw耦合關(guān)系指數(shù)值均超過了100，表明這3個參數(shù)的耦合關(guān)系強，其取值易受其他參數(shù)的影響。 因此，在反演時不能簡單地將三者看作獨立變量進行分析，需考慮該參數(shù)和其他參數(shù)的耦合關(guān)系帶來的不確定性影響。 其他7個參數(shù)的耦合關(guān)系指數(shù)值較低，表明模型獲得的參數(shù)的不確定性較小、可靠性高。

2.4 枯草關(guān)鍵參數(shù)閾值率定

基于表2枯草參數(shù)閾值范圍，對5個高敏感性參數(shù)再次進行OFTA[15](one factorata time)局部敏感性分析。 即： 通過固定模型內(nèi)所有其他參數(shù)取值，以均勻步長對表2中該參數(shù)閾值區(qū)間一一取值，逐一判別單個參數(shù)變化對模擬結(jié)果的貢獻(xiàn)。 并以余弦距離作為評價函數(shù)，以99%置信度為標(biāo)準(zhǔn)，對比不同取值下的模擬反射光譜與實測反射光譜間的余弦距離，識別出單個參數(shù)與枯草取值方向的一致性程度及數(shù)值敏感區(qū)間，實現(xiàn)高敏感性參數(shù)的進一步率定，結(jié)果如表3所示。

表3 枯草關(guān)鍵參數(shù)取值區(qū)間參考表Table 3 Reference table for 10 main parameters threshed of withered grass

3 結(jié) 論

利用PROSAIL輻射傳輸模型，結(jié)合實測光譜與枯草性狀參數(shù)，進行枯草反射光譜模擬，以余弦距離作為評價函數(shù)，率定了10個枯草性狀參數(shù)閾值。 模擬結(jié)果顯示，400～2 500 nm的全光譜模擬序列與實測光譜R2通過了0.01顯著性水平檢驗，證實PROSAIL模型對高寒冬季枯草具有很好的模擬能力。 并依據(jù)紅光與綠光波段反射率差值和顧及水汽吸收影響下枯草在400～1 300 nm的準(zhǔn)線性波形特征提出了一種有效篩選枯草光譜的方法。

通過對模擬結(jié)果的兩級敏感性分析與不確定分析，篩選出葉面積指數(shù)、葉傾角、干物質(zhì)、等效水厚度、葉綠素含量等5個枯草變化高敏感參數(shù)，指出了枯草參數(shù)不確定性響應(yīng)程度。 進一步，給出了不敏感性參數(shù)的推薦數(shù)值，界定了枯草敏感的LAI的閾值介于0.2～0.89、LAD為11°～90°、Cab為0～1.29 μg·cm-2、Cw為0.000 1～0.005 cm，Cm為0.008～0.05 g·cm-2，完成了高寒枯草參數(shù)率定。 為提高對高寒冬季枯草的基礎(chǔ)認(rèn)識及遙感反演等研究提供基礎(chǔ)數(shù)據(jù)與依據(jù)。