基于PROSAIL 模型的玉米冠層參數(shù)敏感性分析

李 青

（天津市測(cè)繪院有限公司，天津 300381）

植被葉面積指數(shù)等特征參數(shù)是表征植被生長(zhǎng)發(fā)育狀態(tài)的重要量化指標(biāo)，被廣泛應(yīng)用于農(nóng)業(yè)管理、環(huán)境檢測(cè)、生態(tài)保護(hù)和氣候變化等領(lǐng)域。植被冠層輻射傳輸物理模型建立了植被特征參數(shù)與植被光譜之間的數(shù)學(xué)關(guān)系，能夠通過遙感手段觀測(cè)光譜，進(jìn)而反演植被特征參數(shù)。PROSAIL 模型是一種基于輻射傳輸理論建立的植被冠層反射率模擬模型，其輸入?yún)?shù)包括葉片結(jié)構(gòu)和生理特征、冠層結(jié)構(gòu)、土壤反射率、天空光散射比和觀測(cè)幾何等。為分析不同參數(shù)對(duì)應(yīng)的光譜特征，并探究不同輸入?yún)?shù)在不同波段對(duì)模型輸出結(jié)果的影響，本文以玉米冠層為研究對(duì)象，利用OAT 和EFAST 敏感性分析方法，分析各輸入?yún)?shù)對(duì)冠層反射率變化的貢獻(xiàn)。

1 模型與方法

1.1 PROSAIL 模型

PROSAIL 模型是建立在PROSPECT 葉片模型和SAILH 冠層模型之上的。在葉片尺度上，PROSPECT提供了400～2 500 nm 光譜上葉片的方向半球反射率和透過率，并將其作為SAILH 的輸入?yún)?shù)。PROSPECT 模型以Allen 的平板模型（plate model）為基礎(chǔ)，將葉片描述為一個(gè)或多個(gè)疊加起來的表面粗糙的各向同性平板[1]。這個(gè)模型假設(shè)葉片是由N 層同性層疊加起來的，每?jī)蓪又g被空氣間隔，葉片結(jié)構(gòu)參數(shù)N 就是用來描述葉片中的平均空氣/細(xì)胞間隔的。在葉子內(nèi)部，光線被認(rèn)為是各向同性的[2]。PROSPECT 模型的輸入?yún)?shù)包括：葉片結(jié)構(gòu)參數(shù)，葉綠素含量、類胡蘿卜素含量、褐色素含量、葉片含水量和葉片干物質(zhì)含量等。

SAILH 模型是Kuusk 在SAIL（Scattering by Arbitrarily Inclined Leaves）模型的基礎(chǔ)上引入熱效應(yīng)得到的改進(jìn)模型。SAIL 模型的數(shù)學(xué)形式是一個(gè)基于輻射傳輸方程的4 方程9 參數(shù)的線性微分方程組，是Verhoef[3]和Bunnik 根據(jù)Suits 模型擴(kuò)充得到的。由于輻射傳輸方程微分-積分的形式無法求得其嚴(yán)格的解析解，應(yīng)用比較廣泛的解法是KM 理論，即用一組微分方程近似輻射傳輸方程，只考慮下行漫散射和上行漫散射。Suits模型在此基礎(chǔ)上加入了觀測(cè)方向的輻亮度[4]。SAIL 模型則以Suits 為基礎(chǔ)，加入葉片傾角和分布函數(shù)，可以計(jì)算任意天空散射光比例、太陽高度、觀測(cè)方向下的冠層反射率[3]。由于SAIL 模型不能很好地對(duì)熱點(diǎn)效應(yīng)進(jìn)行模擬，因此Kuusk[5]在SAIL 模型的基礎(chǔ)上加入熱點(diǎn)效應(yīng)從而發(fā)展出了SAILH 模型，引入熱點(diǎn)效應(yīng)因子，熱點(diǎn)效應(yīng)因子為葉片長(zhǎng)度和冠層高度的函數(shù)。SAILH模型的輸入?yún)?shù)包括：葉面積指數(shù)、平均葉傾角（葉傾角分布函數(shù)）、熱點(diǎn)參數(shù)、土壤反射率、天空光散射比、太陽天頂角、觀測(cè)天頂角和相對(duì)方位角等。

1.2 模型參數(shù)設(shè)置

本文以玉米冠層為研究對(duì)象，設(shè)定PROSAIL 模型輸入?yún)?shù)見表1。根據(jù)玉米的特征，玉米葉片大多平展且狹長(zhǎng)，在生長(zhǎng)的過程中葉片長(zhǎng)度與植株高度的比值從接近1 逐漸減小到0.25 左右。因此將平均葉傾角（ALA）的默認(rèn)值設(shè)為50°，取值范圍為30～60°；熱點(diǎn)參數(shù)（hspot）默認(rèn)值設(shè)為0.25，取值范圍為0.25～1。除此之外，由于褐色素含量與葉綠素含量相沖突，且其含量不易測(cè)量，因此將褐色素含量設(shè)為0，即衰老的植株不予考慮。

表1 PROSAIL 模型輸入?yún)?shù)及取值范圍

1.3 敏感性分析方法

敏感性分析可以確定輸入?yún)?shù)對(duì)模型模擬結(jié)果影響的大小，分析方法主要分為局部敏感性分析和全局敏感性分析。

局部敏感性分析是研究單個(gè)輸入?yún)?shù)在局部范圍內(nèi)變化時(shí)對(duì)模型輸出響應(yīng)的影響，局部敏感性分析方法簡(jiǎn)單、計(jì)算量小，但其缺點(diǎn)在于無法充分描述模型參數(shù)的空間分布形態(tài)，并且忽略了參數(shù)之間的相互作用[6]。本文采用一次一個(gè)變量法（OAT）計(jì)算每個(gè)參數(shù)在其默認(rèn)值附近進(jìn)行變化導(dǎo)致的模型輸出的變化率，其變化率的絕對(duì)值代表了該參數(shù)的敏感性。由于之后的全局敏感性分析會(huì)得出個(gè)輸入?yún)?shù)在各波段的一階敏感性，因此這一部分不做定量分析，僅定性展示各輸入?yún)?shù)的局部敏感性。

全局敏感性分析則能夠分析整個(gè)參數(shù)空間對(duì)模型結(jié)果的影響，考慮參數(shù)之間的相互作用。常用的全局敏感性分析方法包括回歸分析法、Morris 篩選法、Sobol’方法和EFAST 方法等[7]。PROSAIL 模型數(shù)據(jù)間呈現(xiàn)非線性關(guān)系且各輸入?yún)?shù)間相互作用，不宜采用回歸分析法。Morris 篩選法可以簡(jiǎn)單快速地篩選出敏感性高的參數(shù)，但無法量化影響，Sobol 方法可以定量分析各參數(shù)敏感性大小，但對(duì)于復(fù)雜模型計(jì)算成本高。因此本文選用能定量分析且計(jì)算量較小的EFAST 方法進(jìn)行敏感性分析。

EFAST 方法是由Saltelli 等結(jié)合Sobol’法和傅里葉振幅靈敏度檢驗(yàn)法（Fourier amplitude sensitivity test，F(xiàn)AST）的優(yōu)點(diǎn)提出的基于方差分解的全局敏感性分析方法[7]，可以計(jì)算得到每個(gè)參數(shù)及參數(shù)相互作用對(duì)模型結(jié)果的方差的貢獻(xiàn)率。該方法的核心思想是通過搜索函數(shù)對(duì)模型的每個(gè)參數(shù)引入一個(gè)共同獨(dú)立參數(shù)的函數(shù)，并給每個(gè)輸入?yún)?shù)定義一個(gè)整數(shù)頻率，最后使模型成為獨(dú)立參數(shù)的周期函數(shù)，對(duì)模型的輸出結(jié)果進(jìn)行傅里葉變換，獲得傅里葉級(jí)數(shù)的頻譜曲線，通過該頻譜曲線分別得到參數(shù)的敏感性。

式中：Vi是參數(shù)xi引起的方差，Vij、Vijm、…、V1，2，…，k是各參數(shù)相互作用貢獻(xiàn)的方差。那么，通過歸一化處理，參數(shù)xi的一階敏感性指數(shù)可定義如下

一階敏感性指數(shù)反應(yīng)的是參數(shù)對(duì)模型輸出總方差的直接貢獻(xiàn)率。同理，參數(shù)xi的二階及三階敏感性指數(shù)可定義為

參數(shù)xi對(duì)模型輸出總方差的直接貢獻(xiàn)和通過參數(shù)間相互作用對(duì)模型輸出方差的貢獻(xiàn)可以用總敏感性指數(shù)STi來表示，即各階敏感性指數(shù)之和

該方法采樣的總個(gè)數(shù)為

式中：Nr是搜索曲線的重復(fù)次數(shù)，可取2；M 是干涉因子，一般取4 或更高；ωmax為最大采樣頻率；n 為參數(shù)個(gè)數(shù)[7]。

EFAST 敏感性分析借助專業(yè)敏感性分析軟件Simlab（Version 2.2.1）[8]完成。方法實(shí)現(xiàn)過程包括樣本生成、模型選擇、Monete Carlo 模擬及模擬結(jié)果的不確定性和敏感性分析。具體步驟如下所示。

1）對(duì)PROSAIL 模型的所有輸入變量在取值范圍內(nèi)按均勻分布采樣；

為確保在行洪安全的前提下，充分合理地開發(fā)利用沂沭河汛期洪水資源，實(shí)現(xiàn)洪水資源化，臨沂市通過新建和改建河道蓄水工程，對(duì)河道來水進(jìn)行節(jié)節(jié)攔蓄、梯級(jí)開發(fā)，汛期調(diào)蓄洪水，非汛期攔蓄來水，從而有效緩解流域內(nèi)工農(nóng)業(yè)用水矛盾，改善區(qū)域水環(huán)境，實(shí)現(xiàn)區(qū)域水循環(huán)。目前，已在祊河、沂河沿城區(qū)段建成7座橡膠壩和8座攔河閘，在臨沂城形成了88.5 km長(zhǎng)、面積達(dá)48.5 km2的水面，一次性蓄水可達(dá)1.62億m3。

2）使用MATLAB 語言編寫的程序讀入樣本，運(yùn)行PROSAIL 模型進(jìn)行模擬；

3）將PROSAIL 模型輸出的反射率按Simlab 軟件樣本輸出格式將寫入文件；

4）在Simlab 中選擇結(jié)果文件，執(zhí)行Monete Carlo模擬和EFAST 敏感性分析。

2 結(jié)果與分析

2.1 局部敏感性分許

局部敏感性分析采用的方法是OTA，對(duì)PROSAIL模型的葉綠素含量（Cab）、類胡蘿卜素含量（Car）、干物質(zhì)含量（Cw）、葉片水分含量（Cm）、葉片結(jié)構(gòu)參數(shù)（N）、葉面積指數(shù)（LAI）、熱點(diǎn)參數(shù)（hspot）和平均葉傾角（ALA）8個(gè)冠層參數(shù)進(jìn)行敏感性分析，每次只改變一個(gè)輸入?yún)?shù)的值，其他參數(shù)取表1 中的默認(rèn)值，觀察輸入?yún)?shù)的改變對(duì)冠層反射光譜在400～2 500 nm 波段的影響，結(jié)果如圖1 所示。

圖1 冠層參數(shù)局部敏感性O(shè)TA 分析結(jié)果

由圖1（a）（b）可知，對(duì)于葉片中的色素來說，Cab和Car對(duì)冠層反射光譜的影響主要集中在可見光波段，其中對(duì)Cab最敏感的波段為550～670 nm，對(duì)Car最敏感的波長(zhǎng)為525 nm。Cab對(duì)反射率的影響程度遠(yuǎn)遠(yuǎn)大于Car。隨著色素含量的增加，冠層的反射率減小，且減小的幅度越來越小。這可能是由于色素對(duì)可見光特定波段的強(qiáng)吸收造成的。由圖1（c）可知，水分子在1 200、1 400、1 950 nm 處存在吸收谷，因此Cw在1 200～2 500 nm 的波段范圍內(nèi)對(duì)植被反射率均有較大影響，反射率隨著葉片水分的增加而減小。對(duì)于冠層葉面積指數(shù)來說，當(dāng)LAI＞2 時(shí)，冠層反射率隨著LAI 的增大而增大；當(dāng)LAI≤2 時(shí)，冠層反射率在400～700 nm 和1 250～2 500 nm 處有明顯的抬高，這是因?yàn)榈乇碇脖幌∈柙斐赏寥辣尘奥懵?，模擬得到的反射光譜中包含了土壤的光譜特性。ALA 對(duì)反射率的影響主要集中在700～1 400 nm 波段，葉傾角越大（葉片越平）反射率越小。除此之外，干物質(zhì)含量（Cm）、葉片結(jié)構(gòu)參數(shù)（N）和熱點(diǎn)因子（hspot）沒有表現(xiàn)出在特定波段的敏感性，對(duì)反射率的影像不大。葉片結(jié)構(gòu)參數(shù)和熱點(diǎn)因子與冠層反射率成正相關(guān)，干物質(zhì)含量則相反。

針對(duì)觀測(cè)幾何參數(shù)，進(jìn)行如下分析：①保持觀測(cè)天頂角和相對(duì)方位角均為0°，改變太陽高度角，發(fā)現(xiàn)隨著太陽高度角的減小，植被冠層的反射率增加。分析原因主要是由于觀測(cè)角度和太陽角度逐漸靠近，并且植被陰影逐漸減少。②保持太陽高度角為30°，改變觀測(cè)天頂角和相對(duì)方位角，計(jì)算觀測(cè)半球上的冠層反射率，得到結(jié)果如圖2 所示。從圖中可以發(fā)現(xiàn)，整體上隨著觀測(cè)天頂角的減小和相對(duì)方位角的增大，冠層反射率減小，但相對(duì)方位角為0°時(shí)除外，在這個(gè)角度可以觀察到明顯的熱點(diǎn)效應(yīng)。

圖2 觀測(cè)半球的反射率（黑圈為太陽入射位置）

2.2 全部敏感性分析

對(duì)于表1 中的所有冠層生理結(jié)構(gòu)參數(shù)進(jìn)行EFAST全局敏感性分析，對(duì)400～2 000 nm 的波段范圍結(jié)果進(jìn)行輸入輸出均勻采樣，得33 組的敏感性分析結(jié)果，同時(shí)計(jì)算得到了各輸入?yún)?shù)的一階敏感性（局部敏感性）和全局敏感性。實(shí)驗(yàn)結(jié)果表明，各輸入?yún)?shù)在各波段的一階敏感性和全局敏感性趨勢(shì)基本一致。全局敏感性較強(qiáng)的輸入?yún)?shù)包括：LAI、Cab、Cw和ALA。在400～500 nm、750～950 nm 和1 900～2 000 nm 波段，LAI 對(duì)反射率方差的貢獻(xiàn)最大，全局敏感性指數(shù)均超過0.5。Cab在500～750 nm 波段占據(jù)絕對(duì)優(yōu)勢(shì)，對(duì)反射率方差貢獻(xiàn)高達(dá)0.7。而Cw則是在1 150～1 850 nm 成為敏感性最高的輸入?yún)?shù)，全局敏感性大多超過0.6。

針對(duì)能夠反應(yīng)植被光譜特征的藍(lán)光、綠光、紅光、近紅外波段，對(duì)輸入?yún)?shù)的敏感性進(jìn)行了排序和具體分析。如圖3 所示，在藍(lán)光波段和近紅外波段LAI 敏感性最高，分別超過0.8 和0.4；在綠光和紅光波段，Cab的敏感性最高，能夠解釋超過一半的反射率方差變化，這是由于葉綠素吸收紅光反射綠光的特性造成的。

圖3 植被光譜特征波段EFAST 敏感性分析結(jié)果

3 結(jié)論

本文對(duì)PROSPECT 模型的輸入?yún)?shù)進(jìn)行了定性的局部敏感性分析和定量的EFAST 全局敏感性分析，得到以下結(jié)論：

全局敏感性較強(qiáng)的輸入?yún)?shù)包括：冠層葉面積指數(shù)LAI、葉片葉綠素含量Cab、葉片水分含量Cw和平均葉傾角ALA。

1）對(duì)于葉片葉綠素含量Cab，其對(duì)冠層反射光譜的影響主要集中在綠光和紅光波段，能夠解釋超過一半的反射率方差變化。隨著葉綠素含量的增加，冠層的反射率減小。

2）水分子在1 200、1 400、1 950 nm 處存在吸收谷，葉片水分含量Cw在1 150～1 850 nm 波段為敏感性最高的輸入?yún)?shù)，反射率隨著葉片水分的增加而減小。

3）對(duì)于冠層葉面積指數(shù)來說，在400～500 nm、750～950 nm 和1 900～2 000 nm 波段，LAI 對(duì)反射率方差的貢獻(xiàn)最大，全局敏感性指數(shù)均超過0.5。當(dāng)LAI＞2時(shí)，冠層反射率隨著LAI 的增大而增大；當(dāng)LAI≤2時(shí)，地表植被稀疏造成土壤背景裸露，模擬得到的反射光譜中包含了土壤的光譜特性，導(dǎo)致冠層反射率在400～700 nm 和1 250～2 500 nm 處有明顯的抬高。

4）隨著太陽高度角的減小，植被冠層的反射率增加；保持太陽高度角為30°，觀測(cè)半球上的冠層反射率整體上隨著觀測(cè)天頂角的減小和相對(duì)方位角的增大，冠層反射率減小。相對(duì)方位角為0°時(shí)可以觀察到明顯的熱點(diǎn)效應(yīng)。