基于貝葉斯模型的森林高度極化干涉SAR反演不確定性分析

張庭葦，張王菲，張永鑫，黃國(guó)然

西南林業(yè)大學(xué)林學(xué)院,昆明 650224

1 引言

采用極化干涉合成孔徑雷達(dá)PolInSAR（Polarimetric Interferometry Synthetic Aperture Radar）技術(shù)進(jìn)行森林高度反演是目前森林生物量、陸地碳儲(chǔ)量計(jì)算等相關(guān)研究的主要方法之一，然而采用PolInSAR 技術(shù)進(jìn)行森林高度反演時(shí)仍然存在諸多不確定性因素。森林高度反演結(jié)果的不確定性直接造成陸地碳儲(chǔ)量計(jì)算結(jié)果的不確定性（傅煜等，2014，2015；廖展芒，2019），因此，有必要針對(duì)采用PolInSAR 技術(shù)進(jìn)行森林高度反演的不確定性展開(kāi)研究。

在基于PolInSAR 技術(shù)進(jìn)行森林高度反演的研究中，隨機(jī)體地表散射RVoG（Random Volume over Ground）模型是使用最廣泛的植被散射物理模 型（López-Martínez 和Alonso-González，2014；Ballester-Berman 等，2015）。該模型中，各復(fù)相干值線性的分布在復(fù)平面上，各復(fù)相干值在復(fù)平面上分布的位置由森林高度、消光系數(shù)、地面相位及垂直有效波數(shù)等決定（Treuhaft 等，1996）。RVoG 模型針對(duì)森林在L 波段的散射特征提出，在L 波段中的反演結(jié)果也最優(yōu)（Papathanassiou 等，1998）。采用PolInSAR 技術(shù)和RVoG 模型估計(jì)森林高度是一個(gè)多階段過(guò)程，相干優(yōu)化、地面相位估計(jì)和森林高度反演等階段累積的不確定性會(huì)傳播到每個(gè)后續(xù)階段（吳小丹 等，2014）。各個(gè)階段累積的不確定性通常可以通過(guò)外部驗(yàn)證數(shù)據(jù)來(lái)量化，例如使用激光雷達(dá)得出的高度與RVoG 模型反演的森林高度進(jìn)行比較。該不確定性估計(jì)結(jié)果精度取決于驗(yàn)證數(shù)據(jù)源的不確定性以及不同數(shù)據(jù)集之間的配準(zhǔn)質(zhì)量，不能較好的反應(yīng)出反演結(jié)果的不確定性。Kugler 等（2015）使用蒙特卡羅方法對(duì)PolInSAR 數(shù)據(jù)執(zhí)行重復(fù)且隨機(jī)的RVoG 模型反演，有效的提高了不確定性分析的精度，然而此方法只考慮了觀測(cè)數(shù)據(jù)引起的不確定性，并沒(méi)有充分考慮概率意義上的建模誤差或模型先驗(yàn)知識(shí)引起的不確定性。Simard 和Denbina（2018）在此基礎(chǔ)上加入了模型各輸入?yún)?shù)（如消光系數(shù)與時(shí)間去相關(guān)等）的先驗(yàn)知識(shí)，由此來(lái)考慮先驗(yàn)知識(shí)引起的不確定性。該研究使用外部LiDAR 數(shù)據(jù)來(lái)固定消光參數(shù)或時(shí)間去相關(guān)的值，結(jié)果表明：加入外部輔助數(shù)據(jù)來(lái)估計(jì)模型參數(shù)的方法，可顯著提高森林高度反演結(jié)果的精度，降低反演結(jié)果的不確定性。Riel等（2018）指出采用貝葉斯模型不僅可以確定RVoG 模型輸入?yún)?shù)引起的不確定性，同時(shí)可以將模型輸入?yún)?shù)、理論假設(shè)、觀測(cè)值等引起的不確定性綜合考慮，進(jìn)而可更客觀的評(píng)價(jià)PolInSAR技術(shù)森林高度反演結(jié)果中的不確定性。

SAR 成像時(shí)觀測(cè)到的森林場(chǎng)景中的各參數(shù)變化亦會(huì)影響森林高度的反演結(jié)果，例如土壤含水量、地面粗糙度等地表因子；樹(shù)種、密度、分布等森林結(jié)構(gòu)因子的變化均會(huì)造成森林高度反演結(jié)果的不確定性（黃揚(yáng) 等，1986；Jackson 和Pinter，1981；高元科，2016；陳魯皖 等，2017）。高元科（2016）采用C 波段PolInSAR 數(shù)據(jù)定量反演森林地表土壤水分含量，結(jié)果表明：雷達(dá)圖像的后向散射系數(shù)產(chǎn)品可以反映土壤水分含量，說(shuō)明土壤含水量在一定程度上影響地表散射。羅時(shí)雨等（2017）在假定植被和土壤特征不變的情況下，研究了土壤含水量對(duì)散射矩陣的影響，指出土壤含水量變化會(huì)導(dǎo)致土壤的散射矩陣變化。地表粗糙度也一直是雷達(dá)數(shù)據(jù)建模和反演的關(guān)鍵因素，其對(duì)后向散射系數(shù)的影響往往超過(guò)土壤含水量、土壤質(zhì)地等其他因素（陳思，2019）。除上述土壤各種性質(zhì)的影響外，有研究表明森林密度會(huì)對(duì)森林高度反演結(jié)果精度產(chǎn)生較大影響，而傳統(tǒng)的反演算法中多未考慮森林密度的影響。Wang 等（2016）使用模擬和機(jī)載PolInSAR 數(shù)據(jù)研究了RVoG 模型在不同森林密度的森林高度反演中的適用性，研究結(jié)果表明：森林密度明顯影響森林高度反演精度，RVoG 模型不適用于植被稀疏區(qū)域的森林高度反演，此外，森林密度與地體散射比呈負(fù)相關(guān)的關(guān)系。姜友誼等（2020）使用模擬數(shù)據(jù)分析了森林密度對(duì)傳統(tǒng)森林高度反演算法的影響，提出了一種考慮森林密度影響的相位與幅度聯(lián)合反演算法，其結(jié)果改善了傳統(tǒng)相位與幅度聯(lián)合反演法精度較低的情況，同時(shí)也降低了反演結(jié)果的均方根誤差。

綜上所述，國(guó)內(nèi)外采用模擬數(shù)據(jù)（曹霸等，2016；周筑博，2013）或真實(shí)數(shù)據(jù)（羅環(huán)敏等，2010；范亞雄，2019）開(kāi)展了部分因子引起的森林高度反演結(jié)果的不確定性研究，但較全面的分析模型輸入?yún)?shù)、模型假設(shè)、觀測(cè)數(shù)據(jù)及森林場(chǎng)景因子協(xié)同引起的森林高度反演結(jié)果不確定性的研究則開(kāi)展較少。此外，由于微波全波段真實(shí)SAR 數(shù)據(jù)獲取的困難性，以及真實(shí)森林場(chǎng)景的復(fù)雜性，采用真實(shí)森林場(chǎng)景無(wú)法清晰闡明各影響因子的影響方式及機(jī)理，因此本文選取模擬的森林場(chǎng)景展開(kāi)研究。本文以RVoG 模型為基礎(chǔ)，基于貝葉斯模型，采用L 波段模擬數(shù)據(jù)結(jié)合PolInSAR 技術(shù)，對(duì)RVoG 模型輸入?yún)?shù)、模型假設(shè)、觀測(cè)值、樹(shù)種、森林密度、地面粗糙度和土壤含水量等因子在森林高度反演過(guò)程中協(xié)同引起的不確定性進(jìn)行了探索性分析。

2 森林場(chǎng)景模擬數(shù)據(jù)生成、預(yù)處理

2.1 森林場(chǎng)景模擬數(shù)據(jù)生成

本文利用歐洲航天局ESA（European Space Agency）發(fā)布的PolSARPro 4.2 版本中由Mark L.Williams 博士開(kāi)發(fā)的SAR 相干散射和成像代碼（PolSARpro sim 模塊），獲得了PolInSAR 的模擬數(shù)據(jù)（Pottier 和Ferro-Famil，2012）。模擬數(shù)據(jù)具有與機(jī)載系統(tǒng)相似的性能，且不存在基線、配準(zhǔn)、時(shí)間去相關(guān)及信噪比去相關(guān)等信號(hào)源相關(guān)的不確定性問(wèn)題。

在該模擬器中，本文首先設(shè)置了表1所示的成像參數(shù)；包括平臺(tái)高度、垂直和水平基線長(zhǎng)度、入射角、距離向和方位向分辨率?；诖顺上駞?shù)設(shè)置了4 種樹(shù)種（圖1）、4 種森林密度（圖2）、4 組土壤含水量和4 組地面粗糙度梯度，用來(lái)生成模擬PolInSAR 數(shù)據(jù)。含水量和粗糙度設(shè)置的數(shù)值中，值越大表示含水量或粗糙度程度越高。本研究中共生成256 景PolInSAR 模擬數(shù)據(jù)，模擬器中樹(shù)種、森林密度、土壤含水量及地面粗糙度模擬值的詳細(xì)設(shè)置見(jiàn)表2。

圖1 PolSARprosim模擬的不同樹(shù)種分布圖Fig.1 Forest species and their distribution maps simulated by PolSARprosim

圖2 PolSARprosim模擬的不同密度森林場(chǎng)景圖（以圓形冠層針葉林為例）Fig.2 The simulated Forest scene maps with different densities（coniferous forest（circular）is selected as an example）

表1 PolSARprosim 模擬器中各參數(shù)設(shè)置Table 1 The parameters settings in the PolSARprosim simulator

表2 PolSARprosim 模擬器中森林場(chǎng)景參數(shù)設(shè)置Table 2 The forest scene parameter settings in the PolSARprosim simulator

2.2 模擬PolInSAR數(shù)據(jù)預(yù)處理

本文中模擬的256 對(duì)PolInSAR 數(shù)據(jù)的預(yù)處理包括干涉圖生成（圖3（a））、平地效應(yīng)估計(jì)（圖3（b））、平地效應(yīng)去除（圖3（c））及復(fù)相干估計(jì)（圖3（d））。生成的復(fù)相干數(shù)據(jù)用于結(jié)合RVoG 模型法和分層貝葉斯框架對(duì)森林高度反演的不確定性進(jìn)行分析。

3 研究方法

3.1 RVoG模型

RVoG 模型是目前采用PolInSAR 技術(shù)反演森林高度中使用較為廣泛的一種物理模型，其通過(guò)建立干涉復(fù)相干系數(shù)與森林高度、地體散射比、地相位之間的物理關(guān)系，以實(shí)現(xiàn)采用干涉復(fù)相干系數(shù)對(duì)森林高度的反演。RVoG 模型是一個(gè)雙層模型，其假設(shè)觀測(cè)到的總回波由森林體散射層和地表散射層組成，其中體散射層無(wú)極化依賴性，通常用來(lái)表征森林場(chǎng)景的散射，而地表散射則具有極化依賴性，用來(lái)表征林下地表的散射。RVoG 模型通常由式（1）表示（Cloude 和Papathanassiou，1998；Papathanassiou 和Cloude，2001；張王菲 等，2017；Riel等，2018）：

式中，R(m)p表示對(duì)于給定的極化P的極化干涉復(fù)相干，?0是底層地表地形的相位，μp是極化P下的有效地-體幅度比，γv是體散射復(fù)相干性，由式（2）和（3）表示。

式中，γv是森林高度hv和消光系數(shù)σ的非線性函數(shù)，垂直有效波kz、雷達(dá)平均入射角θ、基線長(zhǎng)度Bn、波長(zhǎng)λ及斜距R為SAR系統(tǒng)參數(shù)。

3.2 貝葉斯統(tǒng)計(jì)原理

3.2.1 貝葉斯模型

貝葉斯定理就是將先驗(yàn)分布p(θ)（在觀測(cè)到數(shù)據(jù)之前對(duì)問(wèn)題的理解）轉(zhuǎn)換為后驗(yàn)分布p(θ|y)（在觀測(cè)到數(shù)據(jù)之后對(duì)問(wèn)題的理解）的過(guò)程，其本質(zhì)上是一種機(jī)器學(xué)習(xí)的過(guò)程（奧斯瓦爾多·馬丁，2018）。其表達(dá)式如式（4）所示。

式中，p(θ)為先驗(yàn)分布，反映的是在觀測(cè)到數(shù)據(jù)之前我們對(duì)待估計(jì)的參數(shù)的了解和認(rèn)識(shí)；p(y|θ)為似然函數(shù)，是對(duì)實(shí)際觀測(cè)數(shù)據(jù)的一種描述；p(θ|y)為后驗(yàn)分布，是通過(guò)貝葉斯定理得到的最終分析結(jié)果，反映的是在給定觀測(cè)數(shù)據(jù)的基礎(chǔ)上，對(duì)于參數(shù)的新的認(rèn)知；p(y)為邊緣概率，在具體應(yīng)用中通常將之作為后驗(yàn)概率p(θ|y)計(jì)算過(guò)程中的標(biāo)準(zhǔn)化常量，并常被省略，即在實(shí)際應(yīng)用中，式（4）通常表示為式（5）：

3.2.2 多參數(shù)貝葉斯模型和分層貝葉斯模型

（1）多參數(shù)貝葉斯模型?；赗VoG 模型進(jìn)行森林高度反演時(shí)，模型輸入?yún)?shù)多于1個(gè)，因此需要構(gòu)建多參數(shù)貝葉斯模型，貝葉斯模型的多參形式見(jiàn)式（6）：

式中，θ1，…，θn表示模型中的不同參數(shù)。

（2）分層貝葉斯模型。分層貝葉斯模型的構(gòu)建就是在先驗(yàn)之上使用一個(gè)共享先驗(yàn)，即超先驗(yàn)（hyper-prior），其參數(shù)稱為超參數(shù)。在加入超先驗(yàn)之后，新的模型相比于原來(lái)多了一層，以反映建模過(guò)程中多種不確定性來(lái)源，例如模型輸入?yún)?shù)、模型理論假設(shè)、觀測(cè)值等引起的不確定性，并且通過(guò)建立分層模型也可以避免模型參數(shù)過(guò)多引起的過(guò)擬合問(wèn)題（Riel等，2018）。

3.2.3 Metropolis-Hastings算法原理

盡管利用貝葉斯框架可獲得參數(shù)的后驗(yàn)分布，但多數(shù)模型都難以獲得封閉的解析后驗(yàn)分布函數(shù)，因此需要采用數(shù)值方法計(jì)算后驗(yàn)。在無(wú)法獲得后驗(yàn)解析式時(shí)，需要使用推理引擎進(jìn)行貝葉斯分析，以得到后驗(yàn)分布。目前，貝葉斯分析主要是通過(guò)馬爾科夫鏈蒙特卡洛MCMC（Markov Chain Monte Carlo）方法隨機(jī)采樣進(jìn)行。MCMC 方法是在貝葉斯理論框架下，將馬爾可夫（Markov）過(guò)程引入到靜態(tài)蒙特卡洛模擬中實(shí)現(xiàn)對(duì)后驗(yàn)分布函數(shù)動(dòng)態(tài)模擬。其算法是通過(guò)構(gòu)造一個(gè)收斂到π 的馬爾可夫鏈來(lái)實(shí)現(xiàn)從目標(biāo)函數(shù)進(jìn)行抽樣。使用MCMC 方法時(shí)，馬爾科夫鏈轉(zhuǎn)移核的構(gòu)造至關(guān)重要，不同的轉(zhuǎn)移核構(gòu)造方法，將產(chǎn)生不同的MCMC 方法。常用的MCMC 方法主要有兩種：Gibbs 抽樣和Metropolis-Hastings 算法（朱新玲，2009），本文中采用Metropolis-Hastings算法進(jìn)行抽樣。

Metropolis-Hasting算法的迭代步驟如下：

（1）給待估參數(shù)賦一個(gè)初始值，通常是隨機(jī)初始化或者某些經(jīng)驗(yàn)值；

（2）根據(jù)先驗(yàn)分布隨機(jī)生成一個(gè)樣本值，先驗(yàn)分布可以是高斯分布或者均勻分布；

（3）根據(jù)接收率準(zhǔn)則（Riel等，2018）計(jì)算接受一個(gè)新的樣本值的概率；

（4）從位于區(qū)間［0，1］內(nèi)的均勻分布中隨機(jī)選一個(gè)值，并與步驟（3）中得到的概率值進(jìn)行比較，若大于該值，則接受新的值，否則接收步驟（3）中的樣本值；

（5）然后返回步驟（2）重新迭代，直到獲取足夠的樣本。

在該過(guò)程中可獲得一連串?dāng)?shù)值，也稱作采樣鏈或跡。由此可知Metropolis-Hasting 算法是把對(duì)后驗(yàn)求積分的過(guò)程轉(zhuǎn)化成了對(duì)采樣鏈所構(gòu)成的向量求和的過(guò)程，進(jìn)而使后驗(yàn)分析變得簡(jiǎn)單。

本文采用Python 平臺(tái)中的PyMC3 模塊實(shí)現(xiàn)文中涉及的Metropolis-Hastings 算法迭代過(guò)程，獲得符合貝葉斯后驗(yàn)分布的采樣點(diǎn)。另外，由于MCMC 方法依賴于模擬的收斂性，本文通過(guò)生成多條馬爾科夫鏈來(lái)判斷采樣過(guò)程是否收斂，若這些馬爾科夫鏈均穩(wěn)定，則說(shuō)明采樣結(jié)果收斂。

3.3 貝葉斯框架應(yīng)用于PolInSAR數(shù)據(jù)

基于3.2 節(jié)所介紹的貝葉斯統(tǒng)計(jì)原理，本文應(yīng)用的貝葉斯框架形式見(jiàn)式（7），其表示為后驗(yàn)概率p(m|γ，R(m))正比于先驗(yàn)概率p(m)與似然函數(shù)p(γ|m，R(m))乘積的形式，式中似然函數(shù)表示基于RVoG模型輸入?yún)?shù)值對(duì)應(yīng)的觀測(cè)數(shù)據(jù)的概率。

在式（7）中增加超先驗(yàn)p(σγ)，使之變?yōu)榉謱迂惾~斯框架（式（8）），超先驗(yàn)p(σγ)是方差σγ的先驗(yàn)分布，以重新定義不確定性（均值已知、方差未知）的分布。通過(guò)這種方式可將觀測(cè)誤差和模型假設(shè)引起的不確定性表示為各個(gè)像元中的隨機(jī)值σγ。式（8）中，后驗(yàn)分布p(m，σγ|γ，R(m))為聯(lián)合分布，通過(guò)該式可同時(shí)得到森林高度hv和不確定性σγ的后驗(yàn)分布采樣鏈?；谠摲椒色@得更全面的森林高度參數(shù)不確定性度量（Wu 等，2010）。

本文使用Python 中PyMC3 模塊將貝葉斯框架應(yīng)用于PolInSAR 數(shù)據(jù)。首先利用貝葉斯框架分析RVoG 模型中輸入?yún)?shù)和模型假設(shè)的不確定性；在此基礎(chǔ)上利用貝葉斯概率框架同時(shí)獲得基于模擬PolInSAR 數(shù)據(jù)反演的森林高度的后驗(yàn)分布采樣點(diǎn)及其不確定性。文中基于分層貝葉斯框架獲得森林高度反演不確定性的主要步驟包括：

（1）首先采用貝葉斯模型確定RVoG 模型中影響最大的輸入?yún)?shù)，然后使用已有觀測(cè)值“固定”該參數(shù)，再使用“固定”該參數(shù)后求得的RVoG 高度hv作為高斯先驗(yàn)分布p(m)的平均值，由于采用RVoG模型估測(cè)的森林高度結(jié)果的標(biāo)準(zhǔn)偏差為10 m，因此該高斯先驗(yàn)分布的p(m)標(biāo)準(zhǔn)偏差設(shè)置為10 m；使用均勻分布作為固定參數(shù)的先驗(yàn)分布，并根據(jù)經(jīng)驗(yàn)確定該參數(shù)的變化范圍（Riel 等，2018），此步驟可結(jié)合關(guān)于模型參數(shù)和任何不確定性的所有先驗(yàn)信息或假設(shè)。

（2）在似然函數(shù)中加入模型各輸入?yún)?shù)方差σγ的先驗(yàn)分布p(σγ)作為超先驗(yàn)，使貝葉斯模型變?yōu)榉謱迂惾~斯模型，似然函數(shù)形式變?yōu)閜(γ|m，R(m)，σγ)。

（3）采用Metropolis-Hastings 從后驗(yàn)分布中提取近似樣本。在采樣過(guò)程中，參數(shù)的極大后驗(yàn)估計(jì)（Maximum a posteriori estimation）點(diǎn)使用PyMC3提供的find_MAP 函數(shù)獲得。文中采樣鏈數(shù)為2，基于每條采樣鏈獲得50000個(gè)采樣樣本，其中剔除最初的10000個(gè)采樣點(diǎn)以增大采樣鏈的穩(wěn)定性；將采樣過(guò)程重復(fù)2 次以減小系統(tǒng)誤差，并選2 條采樣鏈中收斂性較好的結(jié)果進(jìn)行后續(xù)分析。

3.4 不確定性量化

樣本標(biāo)準(zhǔn)差S常用來(lái)反應(yīng)整個(gè)樣本變量的離散程度，但其無(wú)法直接反映樣本平均數(shù)與總體平均數(shù)之間的誤差。因此，本文采用平均數(shù)的標(biāo)準(zhǔn)誤（Standard Error）對(duì)不確定性進(jìn)行量化。下文中平均數(shù)的標(biāo)準(zhǔn)誤簡(jiǎn)稱標(biāo)準(zhǔn)誤，是標(biāo)準(zhǔn)差的（郝拉娣等，2005）。標(biāo)準(zhǔn)誤是用于衡量樣本均值和總體均值的差距，即多個(gè)樣本均值的標(biāo)準(zhǔn)差（Iversen 等，2000；李炳凱，2007），通常使用標(biāo)準(zhǔn)誤來(lái)評(píng)價(jià)樣本平均數(shù)與總體平均數(shù)的誤差，以此評(píng)價(jià)各因子引起的不確定性。

式中，Xi是模擬得到的PolInSAR 數(shù)據(jù)進(jìn)行森林高度反演結(jié)果樣本，n為觀測(cè)樣本中所含元素的個(gè)數(shù)，為各結(jié)果觀測(cè)值的算術(shù)平均值，為標(biāo)準(zhǔn)誤。

4 結(jié)果與討論

4.1 基于貝葉斯框架的RVoG 模型輸入?yún)?shù)不確定性分析

模型的輸入?yún)?shù)是森林高度反演過(guò)程中不確定性來(lái)源之一。為了描述輸入?yún)?shù)引起的不確定性，基于3.1 節(jié)所述RVoG 模型，設(shè)定森林高度為30 m，消光系數(shù)為0.2 dB/m，其他參數(shù)則設(shè)置為沒(méi)有建模誤差的理想情況下的值（即kz=0.1 m、?0=0°、μp=0 和γt=1），以此合成復(fù)相干系數(shù)γv，并通過(guò)Metropolis-Hastings 算法采樣得到森林高度hv（圖4（a））及消光系數(shù)σ（圖4（b））的采樣樣本，并制作二維KDE 圖（圖4（c））。圖4 中紅色的五角星表示樣本的真值，由圖4可知，即使在無(wú)建模誤差影響時(shí)，森林高度反演結(jié)果的不確定性依然存在。

圖4 貝葉斯后驗(yàn)分布采樣Fig.4 Generated samples from the Bayesian posterior distribution

為了進(jìn)一步分析由模型不同輸入?yún)?shù)誤差帶來(lái)的不確定性，文中同時(shí)設(shè)定了圖5所示的兩組合成數(shù)據(jù)，分析由時(shí)間去相關(guān)γt和消光系數(shù)σ分別作為自由參數(shù)時(shí)在RVoG 模型反演中所引起的不確定性。先基于RVoG 模型生成森林高度hv為30 m、消光系數(shù)σ為0.2 dB/m 的復(fù)相干γv（?0=0°、μp=0）。第一組場(chǎng)景模擬中，森林高度hv和消光系數(shù)σ為自由參數(shù)，γt的值設(shè)定為0.75，并獲得貝葉斯框架下hv和σ的采樣點(diǎn)，其結(jié)果如圖5（a）所示，紅星表示消光系數(shù)及森林高度觀測(cè)值，即真實(shí)值。第二組場(chǎng)景模擬中，森林高度和時(shí)間去相關(guān)為自由參數(shù)，σ的值設(shè)定為0.2 dB/m，并對(duì)hv和γt后驗(yàn)分布進(jìn)行采樣，結(jié)果如圖5（b）所示，紅星表示時(shí)間去相關(guān)及森林高度觀測(cè)值，即真實(shí)值。

圖5 兩種RVoG模型貝葉斯后驗(yàn)采樣樣本的KDE熱圖Fig.5 KDE heat maps of samples sampled according to two different RVoG models using the Bayesian posterior distribution

圖5（a）中可見(jiàn)森林高度hv和消光系數(shù)σ的采樣結(jié)果分散，不確定性大，相比圖5（a），圖5（b）中森林高度hv和時(shí)間去相關(guān)系數(shù)γt兩個(gè)參數(shù)的后驗(yàn)分布受到很好的約束，且彼此之間的協(xié)方差較小。因此可得出，在RVoG 模型中增加消光的先驗(yàn)知識(shí)可以大大減少森林高度反演的不確定性，該結(jié)論在Riel 等（2018）、Simard 和Denbina（2018）的研究中也得到了證實(shí)。鑒于本節(jié)研究結(jié)果，本文采用模擬數(shù)據(jù)將成像時(shí)森林高度的值（18 m）作為森林高度hv的先驗(yàn)知識(shí)，基于式（1）和式（2），求出一個(gè)“固定”消光系數(shù)σ，作為RVoG 模型反演中消光系數(shù)的先驗(yàn)概率信息，進(jìn)而研究基于RVoG的PolInSAR森林高度反演中的不確定性。

基于4.1 節(jié)的研究結(jié)論，我們采用L 波段PolInSAR 模擬數(shù)據(jù)，基于RVoG 模型，分析了樹(shù)種、森林密度、地面粗糙度及土壤含水量變化對(duì)反演結(jié)果引起的不確定性。圖6 描述了不同樹(shù)種、不同密度、不同粗糙度及不同土壤含水量組合下森林的原始高度（觀測(cè)值）、原始RVoG 模型反演高度、固定消光系數(shù)后的RVoG 模型反演高度及貝葉斯后驗(yàn)采樣高度2 倍標(biāo)準(zhǔn)差置信區(qū)間的動(dòng)態(tài)變化；圖7進(jìn)一步分析了土壤含水量、地面粗糙度對(duì)固定消光后的RVoG 模型森林高度反演結(jié)果的影響。圖6 和圖7 中，Tree species 1、Tree species 2、Tree species 3、Tree species 4 分別表示樹(shù)種1（針葉林—錐形）、樹(shù)種2（針葉林—圓形）、樹(shù)種3（針葉林—圓形錐形混交）以及樹(shù)種4（闊葉林）4類樹(shù)種。圖6 中，00、03、06…XX 表示地面粗糙度（0、3、6、10）和土壤含水量（0、3、6、10）的不同組合。為了描述方便，在圖7 中，我們用A、B、C、D 分別表示樹(shù)種1、樹(shù)種2、樹(shù)種3 和樹(shù)種4，用Ⅰ、Ⅱ、Ⅲ、Ⅳ分別表示森林密度為150株/hm2、300株/hm2、600株/hm2、及1200株/hm24 種森林分布密度，AⅠ、AⅡ、AⅢ、…、DⅢ、DⅣ表示樹(shù)種和森林密度的不同組合。圖7 中，綠色表示反演結(jié)果位于真實(shí)值的二倍標(biāo)準(zhǔn)差置信區(qū)間，其中綠色越淺表示反演結(jié)果越接近真實(shí)值。由圖6 可知固定消光系數(shù)后，基于RVoG 模型反演的森林高度更接近真實(shí)值。Riel 等（2018）基于RVoG 模型使用LiDAR 數(shù)據(jù)進(jìn)行消光系數(shù)的固定，其研究結(jié)果同樣表明，消光系數(shù)的固定可以大大減小森林高度反演的不確定性。圖7描述了固定消光系數(shù)后L 波段基于RVoG 模型反演的森林高度結(jié)果隨樹(shù)種、密度、土壤含水量及地面粗糙度的變化。

圖6 RVoG模型演結(jié)果與貝葉斯處理結(jié)果比較Fig.6 The comparison between the results from direct RVoG model and combination of RVoG model and Bayesian posterior distribution

圖7 基于貝葉斯后驗(yàn)分布森林高度反演結(jié)果不確定性隨森林場(chǎng)景因子變化的關(guān)系Fig.7 The changes of Forest height inversion uncertainty with the variation of soil moisture content，soil surface roughness，forest species and forest density

4.2 森林高度反演結(jié)果分析

4.2.1 森林高度反演結(jié)果定性分析

圖6橙色線條為分層貝葉斯概率框架后驗(yàn)采樣的結(jié)果，與原始的RVoG模型反演結(jié)果（綠色線條）相比，在各類樹(shù)種、各類密度森林場(chǎng)景中，固定消光的后驗(yàn)采樣結(jié)果均能明顯改善反演結(jié)果，其結(jié)果更接近真值18 m。在4種不同密度的樹(shù)種中，在地面粗糙度小于10時(shí)，反演結(jié)果受地面粗糙度影響不明顯，當(dāng)?shù)孛娲植诙葹?0時(shí)，森林高度出現(xiàn)了明顯的低估現(xiàn)象。由圖6 可知，當(dāng)粗糙度為10時(shí)，4種密度的4種樹(shù)種均有明顯的森林高度低估現(xiàn)象，但在樹(shù)種4中，低估現(xiàn)象最為明顯。這說(shuō)明在PolSARpro sim模擬器中，粗糙度達(dá)到10 時(shí)，地表散射機(jī)制顯著且受到森林密度影響不明顯，這與現(xiàn)實(shí)中森林散射機(jī)制差異較大，其原因需要未來(lái)進(jìn)一步分析。

圖7進(jìn)一步分析了土壤含水量、地面粗糙度對(duì)固定消光后的RVoG 模型森林高度反演結(jié)果的影響。圖7將使用后驗(yàn)采樣的結(jié)果在土壤含水量（0—10）、地面粗糙度（0—10）動(dòng)態(tài)變化范圍內(nèi)進(jìn)行了插值，得到了森林高度隨土壤含水量、地表粗糙度在0—10連續(xù)變化范圍時(shí)的反演結(jié)果，基于此來(lái)進(jìn)一步分析這些環(huán)境因子在森林高度反演中協(xié)同引起的不確定性。其中橫坐標(biāo)表示地面粗糙度，縱坐標(biāo)表示土壤含水量，16 個(gè)子圖表示由4 個(gè)梯度地面粗糙度及4個(gè)梯度土壤含水量?jī)蓛山M合，在二維平面內(nèi)各自插值的森林高度反演結(jié)果，所有結(jié)果均統(tǒng)一量綱為5—25 m，使用綠色表示真值95%置信區(qū)間的范圍，灰色區(qū)域從淺到深表示5—25 m 范圍。觀察圖7中可以發(fā)現(xiàn)，對(duì)于L 波段的模擬數(shù)據(jù)，樹(shù)種、森林密度、地面粗糙度及土壤含水量4個(gè)環(huán)境因子均會(huì)對(duì)森林高度反演結(jié)果造成不確定性：

（1）高度反演結(jié)果不確定性受樹(shù)種影響明顯，圖中表現(xiàn)為置信區(qū)間的動(dòng)態(tài)變化范圍的面積較大。對(duì)比不同樹(shù)種，可以發(fā)現(xiàn)闊葉林反演結(jié)果的不確定性大于針葉林。這可能是由于闊葉樹(shù)葉片較大，散射情況較為復(fù)雜，增加了極化SAR 數(shù)據(jù)復(fù)相干估計(jì)的不確定性，使得闊葉林反演結(jié)果不確定性明顯大于針葉林。針葉林反演結(jié)果整體優(yōu)于闊葉林，其中樹(shù)種2反演結(jié)果最優(yōu)。

（2）森林密度對(duì)高度反演結(jié)果的影響也較明顯，特別是當(dāng)密度為150 株/hm2時(shí)，多數(shù)反演結(jié)果與真值有差異，這可能是由于森林較稀疏時(shí)，有較多的裸露地面（圖2（a）），地面散射影響較大，導(dǎo)致反演結(jié)果整體偏低。值得注意的是，當(dāng)樹(shù)種為闊葉林時(shí)，由于其葉片較大，相比于同樣密度的針葉林，其來(lái)自于地面的散射較少，因此闊葉林即使在森林密度較小時(shí)不確定性也較小。從圖中也可發(fā)現(xiàn)，當(dāng)森林密度較大時(shí)，圖2中綠色區(qū)域覆蓋范圍較大，表明更多的反演結(jié)果接近于真值。但當(dāng)森林密度為600 株/hm2和1200 株/hm2時(shí)，反演結(jié)果差異不大，說(shuō)明當(dāng)森林密度達(dá)到一定密度時(shí)，密度變化引起的不確定性不再明顯。

（3）圖中森林高度反演結(jié)果均在橫軸方向上的變化略大于其隨縱軸的變化，說(shuō)明地面粗糙度的變化對(duì)森林高度反演結(jié)果的影響略大于土壤含水量。在地面粗糙度為6左右時(shí)，整體反演結(jié)果趨于最佳；當(dāng)森林密度為600 株/hm2時(shí)，除了闊葉林外其他3個(gè)樹(shù)種反演結(jié)果誤差均較小。

4.2.2 森林高度反演結(jié)果定量分析

由4.2.1 節(jié)的定性分析可知，4 個(gè)因子均會(huì)造成森林反演結(jié)果的不確定性。為了定量地分析樹(shù)種、森林密度、地面粗糙度及土壤含水量在森林高度反演中協(xié)同引起的不確定性，本文基于3.4 節(jié)所述的標(biāo)準(zhǔn)誤作為衡量指標(biāo)，計(jì)算了不同樹(shù)種森林高度反演中由其他3 個(gè)因子協(xié)同引起的不確定性。其中，表3、表4、表5和表6分別描述了樹(shù)種1、樹(shù)種2、樹(shù)種3和樹(shù)種4的不確定性定量統(tǒng)計(jì)結(jié)果。表中GMC（Ground Moisture Content）表示土壤含水量，SP（Surface Properties）表示地面粗糙度；四種樹(shù)種分別為：樹(shù)種1（針葉林—錐形）、樹(shù)種2（針葉林—圓形）、樹(shù)種3（針葉林—圓形錐形混交）、樹(shù)種4（闊葉林）。

表3 樹(shù)種1標(biāo)準(zhǔn)標(biāo)準(zhǔn)誤S-XTable 3 Standard error of tree species 1/m

表4 樹(shù)種2標(biāo)準(zhǔn)標(biāo)準(zhǔn)誤S-XTable 4 Standard error of tree species 2/m

表5 樹(shù)種3標(biāo)準(zhǔn)標(biāo)準(zhǔn)誤S-XTable 5 Standard error of tree species 3/m

表6 樹(shù)種4標(biāo)準(zhǔn)標(biāo)準(zhǔn)誤S-XTable 6 Standard error of tree species 4/m

由表3 可知，在樹(shù)種1 的森林高度反演結(jié)果中，標(biāo)準(zhǔn)誤的值隨著粗糙度增加而增大，最大與最小標(biāo)準(zhǔn)誤值的差值為0.031 m，其中最大的標(biāo)準(zhǔn)誤的值為0.190 m，最小值為0.159 m。標(biāo)準(zhǔn)誤的值隨著森林密度的增大而減小，兩者的差值為0.174 m，因此可以表明：在樹(shù)種1中，森林密度變化引起的不確定性大于地表粗糙度變化的影響。

表4 描述了樹(shù)種2 的森林高度反演結(jié)果中，森林密度、地表粗糙度和土壤含水量協(xié)同引起的不確定性的定量統(tǒng)計(jì)結(jié)果。地表粗糙度對(duì)樹(shù)種2反演結(jié)果不確定性影響的規(guī)律性略低于樹(shù)種1，但在地表粗糙度較大時(shí)，仍然具有較大的標(biāo)準(zhǔn)誤。與樹(shù)種1類似，森林密度引起的反演結(jié)果的標(biāo)準(zhǔn)誤隨密度增加而減小。由森林密度變化引起的標(biāo)準(zhǔn)誤的差異中，樹(shù)種2的差異大于樹(shù)種1，樹(shù)種2中最大值與最小值的差值為0.246 m，與樹(shù)種1的差異為0.072 m。

表5 為樹(shù)種3 森林高度反演結(jié)果不確定性的定量統(tǒng)計(jì)結(jié)果。樹(shù)種3為模擬的針葉混交林，相比樹(shù)種1 和樹(shù)種2 兩種針葉林純林，森林密度、地表粗糙度對(duì)其森林高度反演結(jié)果的不確定性影響的規(guī)律性略有下降。以地表粗糙度為例，粗糙度變化對(duì)反演結(jié)果的標(biāo)準(zhǔn)誤出現(xiàn)波動(dòng)現(xiàn)象，在粗糙度為3和10 時(shí)，兩者相差不明顯。此外，當(dāng)森林密度由150 株/hm2變?yōu)?00 株/hm2時(shí)，標(biāo)準(zhǔn)誤明顯降低，但隨著密度繼續(xù)增加，標(biāo)準(zhǔn)誤的變化并不明顯。

表6 總結(jié)了樹(shù)種4 的森林高度反演不確定性統(tǒng)計(jì)結(jié)果。模擬的樹(shù)種4為闊葉林，其地表粗糙度對(duì)反演結(jié)果標(biāo)準(zhǔn)誤的影響與樹(shù)種1 和樹(shù)種2 相似，隨著粗糙度的增加，反演結(jié)果的標(biāo)準(zhǔn)誤增加。然而，與前3種針葉林差別較為明顯的是，密度變化對(duì)森林高度反演結(jié)果標(biāo)準(zhǔn)誤的影響并不明顯。

綜合考慮樹(shù)種、森林密度、地表粗糙度及土壤含水量對(duì)L 波段基于RVoG 模型的森林高度反演結(jié)果的影響發(fā)現(xiàn)：樹(shù)種對(duì)森林高度反演中引起的不確定性較大。樹(shù)種1（針葉林—錐形）、樹(shù)種2（針葉林—圓形）、樹(shù)種3（針葉林—圓形錐形混交）3 種針葉林的不確定性大小差別不大，其中針葉林樹(shù)種1 純林不確定性最小，為0.176 m。而當(dāng)樹(shù)種為闊葉林時(shí)，不確定性最大，所有結(jié)果反演標(biāo)準(zhǔn)誤均值為0.296 m，與樹(shù)種1 的差值為0.120 m。森林密度對(duì)針葉林反演結(jié)果的標(biāo)準(zhǔn)誤影響較大，對(duì)闊葉林的影響較小。在針葉林中，當(dāng)森林密度為150 株/hm2時(shí)，樹(shù)種1、樹(shù)種2、樹(shù)種3 的標(biāo)準(zhǔn)誤分別為0.298 m、0.349 m 和0.379 m，當(dāng)森林密度增至300 株/hm2時(shí)，此3 類樹(shù)種標(biāo)準(zhǔn)誤急劇減小，分別降至0.149 m、0.140 m 和0.151 m，該結(jié)論與Wang 等（2016）等研究結(jié)果相似，RVoG 模型不適用于植被稀疏地區(qū)的森林高度反演，森林密度的增加能明顯降低森林反演結(jié)果不確定性。當(dāng)森林密度增加至600 株/hm2時(shí)，標(biāo)準(zhǔn)誤減小不明顯。另外，當(dāng)森林密度增至1200株/hm2，此3類森林標(biāo)準(zhǔn)誤依然有0.124 m、0.103 m 以及0.123 m，說(shuō)明森林高度反演的不確定性不完全由森林密度所決定。其原因可能是當(dāng)針葉林林分的地表植被稀疏時(shí)，裸露的地表會(huì)極大地影響散射，而闊葉林由于其葉片寬大，稀疏林分也不易有較多裸露地表，故散射依然大多來(lái)自植被本身，所以，闊葉林不確定性受森林密度影響較針葉林小。針葉林和闊葉林純林中，地面粗糙度的變化均會(huì)引起森林高度反演結(jié)果較大的不確定性，并且隨著地面粗糙度的增加，標(biāo)準(zhǔn)誤也隨之增加。地面粗糙度從0 增加為最大值10 時(shí)，樹(shù)種1、樹(shù)種2、樹(shù)種3 的標(biāo)準(zhǔn)誤分別從0.159 m 增加為0.190 m、0.174 m增加到0.193 m、0.168 m 增加到0.223 m。較針葉林純林，闊葉林純林受地面粗糙度影響更明顯，當(dāng)粗糙度從0 變?yōu)?0 時(shí)，其標(biāo)準(zhǔn)誤也由0.261 m 增至0.341 m，增幅達(dá)到31%。陳思（2019）對(duì)土壤含水量及地表粗糙度在雷達(dá)數(shù)據(jù)建模和反演中的影響進(jìn)行分析時(shí)，結(jié)果同樣表明地面粗糙度的影響大于土壤含水量。此外，Wang 等（2016）采用模擬數(shù)據(jù)，分析討論了10 m，14 m及18 m等3種樹(shù)高基于RVOG模型反演結(jié)果，其結(jié)果表明，隨著樹(shù)高的增加，其RMSE也相應(yīng)增加。

5 結(jié)論

本文基于L 波段的PolInSAR 模擬數(shù)據(jù)，首先采用貝葉斯模型，明確了RVoG 模型森林高度反演中輸入?yún)?shù)對(duì)反演結(jié)果不確定性的影響，確定了使用固定消光的方法能夠顯著降低RVoG 模型森林高度反演的誤差及不確定性。以此為基礎(chǔ)，本文綜合分析了樹(shù)種、森林密度、地面粗糙度及土壤含水量四個(gè)因子對(duì)森林高度反演結(jié)果的不確定性影響的協(xié)同效應(yīng)。通過(guò)基于256 組模擬PolInSAR數(shù)據(jù)森林高度反演結(jié)果對(duì)比分析得出以下結(jié)論：（1）森林高度反演結(jié)果受樹(shù)種影響較大，闊葉林反演結(jié)果的不確定性高于針葉林；（2）森林密度變化對(duì)基于RVoG 模型的森林高度反演結(jié)果的不確定性影響明顯。特別是在針葉林純林中，密度越大，不確定性越低。森林密度較小時(shí)，RVoG 模型反演的森林高度結(jié)果標(biāo)準(zhǔn)誤大。（3）地面粗糙度的變化與森林高度反演結(jié)果不確定性影響呈正相關(guān)，即粗糙度越大，標(biāo)準(zhǔn)誤越大；當(dāng)?shù)孛娲植诙仍黾訒r(shí)，反演結(jié)果的標(biāo)準(zhǔn)誤也相應(yīng)增加，在不同梯度的森林密度下，粗糙度的增加，導(dǎo)致反演結(jié)果的標(biāo)準(zhǔn)誤急劇增加。（4）與地面粗糙度相比，土壤含水量引起的標(biāo)準(zhǔn)誤變化十分小，幾乎可忽略不計(jì)。由于本文中相關(guān)結(jié)論基于模擬的L 波段PolInSAR 數(shù)據(jù)及模擬器中森林場(chǎng)景的設(shè)置，該數(shù)據(jù)中森林的散射機(jī)制受到模擬器中使用的散射模型、已有森林場(chǎng)景參數(shù)等的限制，其在真實(shí)森林場(chǎng)景中的適用性有待進(jìn)一步研究。此外，PolInSAR 數(shù)據(jù)基線選取、頻段選擇等引起的不確定性在本文中也未考慮，在未來(lái)的應(yīng)用中這些因子引起的不確定性仍然需要進(jìn)一步明確。