基于Sentinel-2影像與PROSAIL模型參數(shù)標定的玉米冠層LAI反演

蘇 偉， 鄔佳昱， 王新盛， 謝茈萱， 張 穎， 陶萬成， 金 添

1. 中國農業(yè)大學土地科學與技術學院， 北京 100083 2. 農業(yè)部農業(yè)災害遙感重點實驗室， 北京 100083

引 言

葉片是玉米光合作用的主要器官， 玉米的干物質積累大多來自葉片， 因此準確描述玉米葉片生長狀態(tài)對于玉米長勢監(jiān)測、 災害脅迫監(jiān)測、 產量預測等具有重要意義[1]。 葉面積指數(shù)是單位土地面積上的植物葉片表面積的一半[2]， 與植被各生理生化過程有密切聯(lián)系， 是表征玉米長勢的一種重要參數(shù)。 遙感數(shù)據(jù)具有空間上的全覆蓋性、 時間上的連續(xù)性和數(shù)據(jù)類型的多樣性， 能夠以不同的時空尺度獲取多種作物冠層信息。 基于遙感影像反演玉米冠層葉面積指數(shù)， 可為玉米長勢監(jiān)測、 病蟲害脅迫監(jiān)測、 產量預測提供有效依據(jù)[3]。

PROSAIL輻射傳輸模型綜合考慮了作物冠層結構、 生長狀況以及遙感觀測環(huán)境的影響， 能夠準確模擬作物冠層反射率， 由于其高敏感性與完善的輻射傳輸機理， PROSAIL模型被廣泛應用于葉面積指數(shù)(leaf area index, LAI)反演中。 近年來， 眾多學者針對PROSAIL模型的LAI反演進行了諸多研究， 反演方法包括查找表法[4]、 經(jīng)驗線性回歸、 機器學習[5]以及模型耦合[6]等， 取得了顯著成果。 其中查找表法是傳統(tǒng)的反演算法， 通過建立大量情況下模型參數(shù)與冠層反射率之間的映射關系， 使用代價函數(shù)進行迭代運算， 將計算所需的時間轉移至反演前， 因此具有廣泛的應用[7]。 機器學習等方法則弱化了模型的機理過程， 訓練模型更為靈活， 近些年發(fā)展迅速。 然而， 由于PROSAIL模型的未知參數(shù)較多且測量費時、 費力， 無論是經(jīng)驗統(tǒng)計方法或機理方法， 均鮮有PROSAIL模型參數(shù)標定的應用， 在反演過程中缺少作物參數(shù)先驗信息， 導致參數(shù)的取值范圍寬泛， 反演的計算量較大， 部分反演結果出現(xiàn)不符合實際的異常值[8]。

所謂參數(shù)標定， 就是調整模型的參數(shù)范圍使其適用特定的研究區(qū)和作物類型[9]。 對于模型參數(shù)標定， 傳統(tǒng)方法是使用試錯法或設置經(jīng)驗值， 設定的參數(shù)取值范圍較為粗略， 難以表達研究區(qū)內作物參數(shù)的時空差異。 此外， 模型中存在著不同參數(shù)可能會帶來同一模擬結果的“異參同效”現(xiàn)象， 在反演計算中將不符合實際情況的參數(shù)作為計算結果， 導致反演結果與實測值的偏差較大。 貝葉斯理論中馬爾可夫鏈蒙特卡洛方法能夠通過構建模型參數(shù)的馬爾科夫鏈， 不斷進化使模型模擬結果達到觀測值的附近， 獲取觀測值不確定性范圍內的模型參數(shù)后驗分布， 其在各類模型的參數(shù)標定中已經(jīng)有了廣泛應用。

為優(yōu)化LAI反演過程提高反演精度， 本研究對如何進行PROSAIL模型參數(shù)標定進行了探索研究。 在利用Sobol法進行PROSAIL模型參數(shù)敏感性分析基礎上， 使用馬爾可夫鏈蒙特卡洛方法(Markov Chain Monte Carlo， MCMC)方法對PROSAIL模型進行參數(shù)標定， 為研究區(qū)內不同地塊的玉米種植區(qū)提供準確的模型輸入?yún)?shù)， 并應用于LAI反演中， 通過參數(shù)標定優(yōu)化前后的LAI反演精度比較分析驗證基于MCMC方法參數(shù)標定的可行性及有效性。

1 研究區(qū)概況與數(shù)據(jù)源

1.1 研究區(qū)概況

研究區(qū)位于河北省廊坊市安次區(qū)， 地理范圍為116°35′—116°51′E， 39°20′—39°33′N， 地理位置如圖1所示。 研究區(qū)地處中緯度地帶， 屬暖溫帶大陸性季風氣候， 主要地形為平原， 平均海拔約13 m， 年平均氣溫8～19 ℃， 年平均降水量550 mm。 主要農作物有冬小麥、 夏玉米、 春玉米、 馬鈴薯等， 研究區(qū)夏玉米的生長時期一般為6月底至9月中下旬。

圖1 研究區(qū)位置及采樣點空間分布

1.2 數(shù)據(jù)源

1.2.1 Sentinel-2A影像

Sentinel-2環(huán)境監(jiān)測衛(wèi)星是歐盟委員會(EC)和歐洲航天局(ESA)共同倡議的全球環(huán)境與安全監(jiān)測系統(tǒng)—“哥白尼計劃”中的第二顆衛(wèi)星， Sentinel-2A、 B星分別于2015年6月23日、 2017年3月7日發(fā)射成功。 Sentinel-2A遙感影像共有十三個波段， 包括可見光、 近紅外、 紅邊、 水汽、 卷云和短波紅外波段， 各波段空間分辨率分別為10， 20和60 m， 影像產品為經(jīng)大氣底部反射率校正的1C級影像， 可以通過ESA的數(shù)據(jù)共享網(wǎng)站下載(https://scihub.copernicus.eu)。 Sentinel-2A衛(wèi)星時空分辨率較高， 具有紅邊波段， 紅邊波段是指植被可見光譜與近紅外光譜之間反射率增大最快的波段， 與植被覆蓋度、 LAI及各類參數(shù)密切相關， 被廣泛用于農作物長勢監(jiān)測中。 本研究用于模型參數(shù)標定及LAI反演的波段為可見光和近紅外波段， 空間分辨率為10 m， 波譜范圍如表1所示。 影像獲取時間2019年8月18日， 影像獲取時太陽天頂角為30.38°， 太陽方位角為144.1°， 平均觀測天頂角為8.15°， 平均觀測方位角為104.08°。

表1 Sentinel-2A影像部分波段信息

1.2.2 地面測量數(shù)據(jù)

與遙感影像對應的地面測量時間為2019年8月16日至22日， 由于不同地區(qū)的種植時間差異， 玉米的生育期為抽穗期至乳熟期。 玉米冠層LAI的測量儀器為美國LI-COR公司生產的植物冠層分析儀LAI-2200C(LI-COR, Lincoln, NE, USA)， 儀器設置為一個天空光、 四個測量目標值， 探頭佩帶45°張角的鏡頭蓋， 四個葉下測量點分別位于壟上、 兩壟之間1/4處、 兩壟中間、 距離壟3/4處。 LAI測量的同時使用中繪i80智能 RTK測量系統(tǒng)進行定位樣點。 在研究區(qū)內不同地塊設置128個采樣點， 每個采樣點采集兩次， 取每一地塊內所有LAI的均值作為該地塊的測量結果。

2 模型與方法

2.1 PROSAIL輻射傳輸模型

PROSAIL模型由PROSPECT葉片光學模型和植被冠層二向反射率SAIL模型耦合而成， 是目前應用最廣泛的植被冠層輻射傳輸模型之一。 PROSPECT模型將葉片結構假設為多層表面粗糙的平板， 假設光線各向同性平行， 通過輸入葉片理化參數(shù)， 模擬光線在平板之間的朗伯散射， 得到400～2 500 nm葉片的反射率與透射率[10]。 SAIL模型假設植被冠層是均一的、 無限延展的混合介質， 并且葉片各向同性， 以輻射傳輸方程為理論依據(jù)， 模擬植被冠層內部輻射傳輸過程， 得到冠層尺度的反射率。 PROSAIL模型將PROSPECT模型輸出的葉片反射率透射率作為SAIL模型的輸入進行模型耦合， 考慮了植被對太陽輻射的吸收、 二向反射、 葉片反射透射率、 葉片結構參數(shù)、 地表狀況等因素的影響， 能較為真實地模擬植被冠層情況， 并且具有很高的穩(wěn)定性， 經(jīng)常用于植被參數(shù)的遙感定量反演中。

PROSAIL模型的輸入?yún)?shù)共包括四類： 葉片理化參數(shù)、 冠層結構參數(shù)、 地面參數(shù)以及角度信息。 用于描述葉片生理生化屬性的參數(shù)有葉片結構系數(shù)N、 葉綠素含量cab(μg·cm-2)、 類胡蘿卜素含量car(μg·cm-2)、 干物質含量cm(g·cm-2)、 等效水厚度cw(cm)； 冠層結構參數(shù)包括LAI和葉傾角分布參數(shù)Lidfa； 地面參數(shù)包括熱點系數(shù)hspot和土壤光譜； 角度信息包括遙感影像拍攝時的太陽天頂角tts、 觀測天頂角tto和相對方位角psi。

2.2 MCMC參數(shù)標定方法

貝葉斯理論可以根據(jù)模型參數(shù)的先驗知識和模型輸出的相應觀測值， 實現(xiàn)模型參數(shù)的后驗估計。 馬爾可夫鏈蒙特卡洛方法在貝葉斯理論框架下， 將馬爾科夫(Markov)過程引入到通過計算機進行模擬的蒙特卡洛方法(Monte Carlo)中， 實現(xiàn)抽樣分布隨模擬的進行而改變的動態(tài)模擬， 構造合適的馬爾科夫鏈進行抽樣而使用蒙特卡洛方法進行積分計算， 即馬爾科夫鏈可以收斂到平穩(wěn)分布。 貝葉斯基本公式如式(1)

(1)

式(1)中，θ為PROSAIL模型的輸入?yún)?shù)，y為模型模擬反射率，p(θ/y)為參數(shù)的后驗概率密度函數(shù)，f(y/θ)為觀測值似然函數(shù)，g(θ)為參數(shù)的先驗分布。

本研究中采用差分進化馬爾科夫鏈算法(differential evolution Markov chain, DE-MC)， 通過多條馬爾科夫鏈的并行更好地搜索參數(shù)空間， 并且為提高采樣效率， Braak等加入snooker更新算法部分代替平行方向上的更新， 克服原DE-MC算法并行鏈數(shù)必需大于空間維數(shù)的限制， 該算法在模型參數(shù)標定中已有應用[9]， 算法的具體原理見參考文獻[11]。

2.3 Sobol模型敏感性分析方法

敏感性分析是指定性或定量地描述模型參數(shù)變化對輸出結果的影響程度， 是模型應用過程中的基礎和重要環(huán)節(jié)。 模型的敏感性分析根據(jù)作用范圍可分為局部敏感性分析和全局敏感性分析， 全局敏感性分析方法除了分析單參數(shù)取值變化對模型結果的影響外， 將各參數(shù)間的相互作用也考慮進來， 分別得到參數(shù)的各階及全局敏感性， 常用的敏感性分析方法有E-FAST， GLUE， Sobol和Morris等。

Sobol法是一種基于方差的典型全局敏感性分析方法， 本研究采用Sobol法進行PROSAIL模型主要參數(shù)的敏感性分析。 Sobol法通過輸出結果的方差分解和在輸入?yún)?shù)范圍內的隨機抽樣取值， 通過輸入?yún)?shù)對輸出方差的貢獻定量測量參數(shù)敏感性[12]。 本研究使用python下SALib敏感性分析環(huán)境提供的Sobol法進行PROSAIL模型的敏感性分析， 目標參數(shù)Y為Sentinel-2A四個波段波長的模擬反射率。

2.4 LAI反演方法

本研究中PROSAIL模型進行LAI反演的方法為查找表法， 查找表法是通過建立離散的參數(shù)與反射率的對應關系表， 使用代價函數(shù)遍歷迭代計算得到目標參數(shù)， 適用于參數(shù)多且復雜時無法求出參數(shù)反函數(shù)的情況。 根據(jù)PROSAIL模型參數(shù)敏感性分析結果， 將敏感參數(shù)作為查找表中的可變參數(shù)， 通過設定可變參數(shù)的取值范圍和變化步長， 其余參數(shù)設為定值， 將模型多次運行的輸入與輸出進行聚合建立LAI反演的查找表。 根據(jù)所建立的模型輸入?yún)?shù)-冠層反射率查找表， 逐像素分析四個波段反射率與查找表中反射率的匹配情況， 尋找像元反射率與模擬反射率最接近的一組參數(shù)。 研究中利用小二乘法構建代價函數(shù)， 計算模擬反射率與影像反射率之間的絕對差值和并最小化， 代價函數(shù)的建立方法如式(2)

(2)

式(2)中，n為參與反演的影像波段數(shù)，Rmodel為模型模擬反射率，RS2A為Sentinel-2A影像反射率， 通過多次迭代計算得到X取最小值時對應的LAI值。

3 結果與討論

3.1 PROSAIL模型的敏感性分析

PROSAIL模型進行敏感性分析， 首先需要明確輸入?yún)?shù)的取值范圍， 模型中的角度參數(shù)根據(jù)影像頭元文件得到， 葉傾角分布函數(shù)通過地基激光雷達數(shù)據(jù)進行構建， 本研究中采用已經(jīng)建立的對應生育期玉米冠層葉傾角函數(shù)[13]， 對余下的六個主要參數(shù)進行敏感性分析， 各參數(shù)范圍根據(jù)經(jīng)驗值及相關文獻獲取如表2所示。

表2 PROSAIL模型主要參數(shù)取值范圍

使用SALib中的Saltelli采樣器， 采樣次數(shù)N為2000， Saltelli采樣器生成的參數(shù)樣本數(shù)由采樣次數(shù)N和參數(shù)個數(shù)D共同決定， 在只輸出一階敏感性和全局敏感性的情況下， 樣本數(shù)n=N×(D+2)。 通過每次在參數(shù)范圍內的隨機采樣， 計算各參數(shù)對模型模擬光譜中波長為490， 560， 665和842 nm的反射率的敏感性大小， 敏感性分析結果如圖2所示。

從圖2(a—b)中可以看出， 各參數(shù)對各波段(波段2， 波段3， 波段4， 波段8)反射率的敏感性各不相同， Sobol法認為敏感性指數(shù)>0.05的參數(shù)對結果表現(xiàn)出敏感性， 敏感性指數(shù)>0.2的參數(shù)為強敏感參數(shù)， 全局敏感性由于考慮了參數(shù)間的相互作用， 一般高于一階敏感性指數(shù)。 從分析結果可以看出： LAI的敏感性最高， 對各波段都表現(xiàn)為強敏感性， 對植被冠層反射率的影響程度最大，N對綠波段和近紅外波段敏感， 且在綠波段表現(xiàn)出強敏感性[圖2(b)]。 cab的敏感性主要體現(xiàn)在綠波段附近， 對其余波段影響程度較小， 而另一色素參數(shù)car對各波段均不敏感。 cm和cw分別反映玉米生長過程中干物質和水分的積累情況， 對可見光波段不敏感， 在近紅外波段[圖2(d)]僅有cm表現(xiàn)出敏感， 且敏感程度較低。 因此， 基于模型敏感性分析結果， 在LAI反演的查找表中， 將LAI，N和cab作為可變參數(shù)， 變化步長分別LAI為1，N為0.1， cab為2， 其余參數(shù)根據(jù)經(jīng)驗值在建立查找表過程中設為定值。

圖2 PROSAIL模型敏感性分析結果(波段2， 3， 4， 8)

3.2 基于MCMC的參數(shù)標定結果與分析

鑒于圖2的PROSAIL模型的敏感性分析， 由于car和cw對研究中所用的波段范圍均不敏感， 在本研究中不進行標定， 僅對LAI，N， cab和cm四個輸入?yún)?shù)進行標定。 除LAI的標定范圍根據(jù)歷年采樣數(shù)據(jù)及生育期經(jīng)驗值設定外， 其余參數(shù)的先驗取值范圍與查找表中相同。 本研究認為冠層反射率服從以影像觀測值為期望的高斯分布， 使用取對數(shù)計算建立似然函數(shù)

logLref=-0.5(x-xobs)TΣ-1(x-xobs)-

0.5Klog(2π)-log(detΣ)

(3)

式(3)中，L為似然函數(shù)的返回值； T為矩陣轉置值； 矢量x和xobs分別表示波段反射率的模型模擬值和觀測值，Σ表示反射率觀測值的協(xié)方差矩陣， 本研究中認為波段反射率觀測值相互獨立， 方差為觀測值的5%；K為空間維數(shù)， 即反射率觀測值的個數(shù)； detΣ表示Σ的行列式值。

算法運行過程中并行鏈數(shù)為4， 每5次采樣進行一次鏈更新， 每10 000次迭代通過比差法進行一次收斂性判斷， 計算診斷指標R， 當連續(xù)超過3次診斷指標R<1.03， 則認為馬爾科夫鏈達到收斂， 收斂后舍棄前期樣本， 保留參數(shù)的后驗樣本。 本研究中的采樣點共分布在25個地塊中， 理論上同一地塊仔的玉米生長狀態(tài)相似， 因此以地塊為單位進行參數(shù)標定， 取地塊內所有玉米冠層反射率的平均值為該地塊參數(shù)標定的輸入， 通過參數(shù)馬爾科夫鏈的進化達到收斂， 得到各參數(shù)的后驗分布與范圍。 圖3中列出了兩次地塊參數(shù)標定后各參數(shù)的先驗后驗分布， LAI和Cab在兩次定標中都得到高斯分布的定標結果。

注： X軸表示參數(shù)取值， Y軸表示概率密度的相對大小， 綠色柱狀圖表示參數(shù)后驗概率分布， 面積和為1

從圖4(a,b)中可以看出， 各參數(shù)先驗分布為取值范圍內的均勻分布(藍色線)， 經(jīng)參數(shù)標定后各參數(shù)的取值及概率分布均發(fā)生了較大變化。 兩次參數(shù)標定后， 各參數(shù)的取值范圍均得到縮小， 的確只有LAI和cab在新的取值范圍內呈現(xiàn)出高斯分布。N和cm兩參數(shù)在實際測量中難以獲取， 給反演中的參數(shù)設置帶來了困難， 參數(shù)標定后兩參數(shù)的取值集中在右邊界， 說明對于玉米及這一生育期而言，N和cm的取值波動較小， 集中在某一定值附近。 參數(shù)標定根據(jù)各地塊的反射率觀測值， 通過馬爾科夫參數(shù)鏈的進化得到了觀測值不確定性范圍內的參數(shù)取值分布， 并提供了每一參數(shù)概率最大的取值， 由此各地塊的玉米生長差異得以表達， 參數(shù)的后驗取值范圍替代原先LAI反演時查找表內參數(shù)較為寬泛的取值， 提高同一地塊內LAI反演的穩(wěn)定性， 避免出現(xiàn)極端值的情況。

3.3 模型參數(shù)標定對LAI反演精度的優(yōu)化分析

根據(jù)2.1中的敏感性分析結果， 首先建立LAI反演的參數(shù)-反射率查找表， 各參數(shù)的取值范圍為默認， 輸入Sentinel-2A影像的觀測反射率， 使用最小化代價函數(shù)進行計算反演得到各采樣點的LAI。 其次， 根據(jù)各地塊參數(shù)標定結果， 使用后驗分布范圍優(yōu)化參數(shù)設置， 進一步縮小各地塊LAI的可能取值范圍建立新的查找表， 并使用每一地塊標定后驗分布中LAI的最高概率密度取值作為該地塊的參考值， 與查找表反演的LAI結果各取50%的可信度進行計算得到新的LAI反演結果。 分別將優(yōu)化前的LAI反演和加入?yún)?shù)標定結果優(yōu)化的LAI反演結果與實測值進行精度評價， 評價指標包括決定系數(shù)R2、 均方根誤差RMSE、 平均偏差Bias和估算精度EA， Bias和EA的計算方法如式(4)和式(5)

(4)

(5)

式(4)和式(5)中，yg和ys分別表示參數(shù)的模型模擬值和實測值， Mean表示實測值的均值， 提取對應實測點位置的LAI反演結果與實測數(shù)據(jù)進行精度驗證， LAI反演的精度驗證結果如圖4(a,b)和表3所示。

圖4 參數(shù)標定前后LAI反演精度比較分析

表3 參數(shù)標定前后LAI反演精度驗證結果

從圖4(a,b)中可以看出， 兩次反演的LAI與實測LAI值擬合效果較好，R2分別達到0.73和0.76， 反演結果沒有明顯的系統(tǒng)偏差， 擬合線與1∶1線較為接近， 證明查找表反演方法的有效性。 傳統(tǒng)的查找表法進行反演時， 查找表中參數(shù)的范圍通常根據(jù)經(jīng)驗值或默認范圍設定， 是一種先驗信息較少的反演方法， 參數(shù)范圍較大導致反演計算量的加大以及異常值的出現(xiàn)。 參數(shù)標定根據(jù)觀測反射率獲取不確定性范圍內的參數(shù)取值， 排除不合先驗的取值， 縮小了各地塊內LAI的可能范圍， 優(yōu)化了原查找表反演中的設定， 使得各地塊玉米生長狀況的差異得以表達， 加快反演速度的同時有效減少了偏差較大的病態(tài)反演出現(xiàn)。 表3列出了加入?yún)?shù)標定結果前后的反演精度差異， RMSE表示反演結果的離散程度， 參數(shù)標定優(yōu)化前后的RMSE分別為0.79和0.32， 有較為明顯的提高， 反演結果均勻地分布在一階擬合線附近。 Bias是反演值與實測值的平均偏差， 描述反演結果偏大或偏小的程度， 在實際應用中這一指標同樣較為重要， 原反演結果的Bias為0.67， 約為實測均值的20%， 參數(shù)標定優(yōu)化后的Bias降低至0.26， 約為實測均值的8%， 同時估算精度提升至90%， 體現(xiàn)了參數(shù)標定在LAI反演中的優(yōu)化效果， 有效地降低了LAI反演值的偏差， 滿足實際應用需求。

4 結 論

以Sentinel-2A影像為數(shù)據(jù)源， 首先使用Sobol方法對PROSAIL模型進行參數(shù)敏感性分析， 量化各參數(shù)對可見光及近紅外波段的敏感性， 為LAI反演提供基礎， 其次基于MCMC方法對PROSAIL模型進行了參數(shù)標定， 目的是提供區(qū)域內玉米各參數(shù)的取值范圍和分布， 獲取更為豐富的參數(shù)信息， 最終對比分析了參數(shù)標定優(yōu)化前后LAI反演精度差異， 探索參數(shù)標定在LAI反演中的應用潛力， 主要研究結論如下：

(1)PROSAIL模型中LAI對可見光和近紅外光均表現(xiàn)出強敏感性， 滿足反演的需求， 同時葉綠素含量Cab對綠波段表現(xiàn)出較強敏感性， 結構系數(shù)N對綠波段和近紅外波段較為敏感， 將這三個參數(shù)作為LAI反演的查找表中的可變參數(shù)， 其余主要參數(shù)對可見光和近紅外波段不敏感， 在反演過程中設為固定參數(shù)。

(2)基于MCMC方法的參數(shù)標定能夠獲取玉米生長區(qū)域內各參數(shù)的后驗概率分布， 表達區(qū)域間的玉米生長狀況差異， 縮小區(qū)域內LAI反演時中各參數(shù)的取值范圍， 為參數(shù)反演提供先驗信息。 參數(shù)標定優(yōu)化后的查找表有效降低了LAI反演偏差， RMSE和Bias分別降低了0.47和0.41， 估算精度由原先的76%提升至90%， 提高了基于PROSAIL模型的LAI反演在實際應用過程中的準確性， 能夠為作物長勢監(jiān)測提供更為準確的數(shù)據(jù)支持。

本研究使用MCMC方法對PROSAIL模型進行過了參數(shù)標定， 對LAI的反演精度有所提高， 但是參數(shù)標定的應用僅限于為查找表及反演提供參考。 在今后的研究中， 將著重挖掘參數(shù)標定在模型中的應用潛力， 與機器學習等方法結合， 進一步優(yōu)化參數(shù)反演過程。