基于多角度衛星遙感的胡楊林結構參數獲取

楊雪峰 木尼熱·買買提

（新疆師范大學地理科學與旅游學院 新疆干旱區湖泊環境與資源實驗室 烏魯木齊 830054）

森林是一個復雜的統一體，由多種生物和非生物環境組成，其彼此之間相互影響、相互制約，表現出一定的結構特征和相互作用規律。森林結構特征是森林植物與周圍環境條件之間以及森林植物彼此之間相互作用的表現形式，可直觀反映森林生長狀況，既是林學、生態學和地學領域的主要研究對象，也是眾多陸地生態模型的重要輸入參數，獲取大尺度森林結構特征具有重要科學價值（Pedroni, 1997；Schr?teret al.,2019）。

森林結構特征通常采用覆蓋度（fractional vegetation cover，FVC）、林分密度（density）、樹高（height）、冠幅（crown width）、胸徑（diameter at breast height）、葉面傾角分布（leaf angle distribution，LAD）和葉面積指數（leaf area index，LAI）等指標表示（Rogasset al., 2012）。目前，基于遙感手段獲取森林結構參數已取得較大發展，能夠獲取大面積森林結構參數的遙感技術主要有激光雷達（light detection and ranging，LiDAR）、雷達、光學遙感等。機載LiDAR是測量植物冠層結構、形狀和地上生物量的有效探測手段（董立新, 2016；Lefsky, 2010；Omasaet al., 2002；van Leeuwenet al., 2010）；但因LiDAR機載設備觀測范圍較小，使用成本較高，星載LiDAR設備地面采樣點又十分稀疏，限制了LiDAR的應用。干涉合成孔徑雷達（interferometric synthetic aperture radar，InSAR）和極化干涉合成孔徑雷達（polarimetric SAR interferometry，PolInSAR）等技術發展出隨機體散射體/地表二層（random volume over ground，RVoG）、極化相干層析（polarization coherence tomography，PCT）等模型方法，可用來獲取森林結構的精細結構特征（Cloudeet al.,2003；Cloude, 2006；Parket al., 2015；談璐璐等, 2010）；然而，現階段星載極化干涉數據面臨時間失相干、空間失相關等影響，其應用受一定限制，甚至無法取得令人信服的反演結果（李震等, 2014；劉曉萌等, 2007；李廷偉等, 2009；Leeet al., 2012）。隨著高分遙感技術的發展，利用影像高分辨率特點獲取樹木結構參數的研究日益增多；但由于森林陰影的長度、太陽高度角、方位角以及坡度、坡向等多個因子具有不確定性，使用高分影像在測量中容易造成較大的人為誤差，無法推廣到大范圍森林參數測量中（崔少偉等, 2011；韓學鋒, 2008；周志強等, 2012；Keet al., 2011）。與單一方向光學遙感相比，多角度光學遙感不僅利用多光譜觀測提取地物組分的波譜信息, 而且利用衛星前向、后向和星下點不同角度的觀測，既富含地表結構信息，同時較多的觀測角度也為反演求解地表結構參數提供了便利。隨著多角度光學遙感的發展, 多種二向反射分布函數（bidirectional reflectance distribution function，BRDF）模型可用來反演并估計地表植被類型（楊雪峰等，2017），采用物理模型結合多角度觀測數據提取植被結構參數具有很大潛力（Suet al., 2007；王亞楠,2010；金晟業, 2009)。目前的星載多角度傳感器中，多角度光譜輻射計（multi-angle imaging spectro-radiometer，MISR）觀測角度最多，通常采用MISR數據結合幾何光學模型（geometric-optical model，GOM）和基于核驅動模型的AMBRALS （algorithm for model bidirectional reflectance anisotropies of the land surface）算法對目標進行區分并反演推測相應的植被結構參數（Dineret al., 1998）。鑒于MISR獨特的大范圍多角度觀測能力以及免費獲取政策，與 LiDAR和 PolInSAR 相比，采用MISR進行森林結構參數提取，在數據獲取能力、后期數據處理便利性和使用成本方面更具優勢。

在干旱區環境下，當林上冠層郁閉度不高時，直射光通過孔隙直接到達林下灌、草植被，此時林下灌、草以及較亮的土壤背景反射很大程度影響傳感器接收的信號，傳感器接收的信號是可見光與植被（喬、灌、草、土壤）之間非常復雜的相互交互過程；如果不能正確分離樹木、林下植被和土壤的光學信號，則反演森林結構參數的任務很難實現（黃健熙等, 2005；Rahmanet al., 2004）。GOM結構未能表達背景反射的二向異性特性，也未充分利用 MISR 的多角度觀測優勢，因此一些研究對GOM進行修改以使其適用于干旱區環境，如簡單幾何光學模型（simple geometric-optical model，SGM）（Choppinget al., 2005；2008；2009）。SGM在模型中加入Walthall土壤反射模型（Walthallet al.,2000）描述林下土壤背景二向異性反射，Walthall模型參數由定標站點回歸計算得到，SGM理論上更適合森林郁閉度較低的環境下使用，但模型適用性未得到廣泛驗證。

對塔里木河下游而言，胡楊（Populus euphratica）的高度、密度、冠幅和群落覆蓋度既是植被恢復的標志性指標，也是研究胡楊群落生態過程與演化機制、改善優勢種群保護措施的途徑（陳海燕等，2015）。目前，該區域有關胡楊結構參數的研究工作主要包括：利用人工地面調查方式建立胡楊林植被動態變化GIS數據庫（玉米提·哈力克等, 2008）；使用QuickBird、WorldView等高分衛星遙感影像數據，采用人機交互方法提取胡楊冠幅、株數等森林結構數據（牛婷等,2008；萬紅梅等, 2011)；基于Landsat、Hyperion等中等分辨率衛星遙感影像數據，利用像元二分模型分解反演塔河下游植被覆蓋度（黃粵等, 2013）等。總體來看，以往研究評價指標較單一，除覆蓋度外，研究范圍局限于樣點和樣帶等小范圍區域，尚缺乏大范圍、系統的胡楊林結構參數提取研究。鑒于此，本研究以MISR為主要數據源，以二向反射模型為理論基礎，采用SGM模型，在無人機攝影測量技術支持下對塔里木河下游胡楊林主要結構參數（樹高、冠幅、林分密度和覆蓋度）進行反演和評價。

1 研究區概況

塔里木河地處我國西部內陸干旱荒漠地區，既是干旱區生態環境變化的敏感地區，也是我國生物多樣性保護和全球變化研究的關鍵區域之一。在人類對自然水資源時空格局改變為主要形式的擾動影響下，塔里木河下游以胡楊林為主要建群種的自然植被受到嚴重影響，生物多樣性受損，沙漠化過程加劇。自2000年以來，有關部門實施向塔里木河下游應急輸水工程，以恢復下游生態系統，胡楊是應急輸水生態恢復的目標植物之一。

研究區位于若羌縣英蘇附近的一處河段，87°59′36″—88°11′48″ E，40°22′51″—40°27′32″ N，面積約100 km2（圖1），平均海拔840 m，地形較為平坦。屬暖溫帶大陸性干旱氣候，平原地區年均氣溫10 ℃，年均降水量約55 mm，降水匱乏，植被稀少。河岸附近發育有河岸林，主體為胡楊，表層土壤為粉砂土，且鹽堿化程度較高。地下水水位較高區域分布有檉柳（Tamarix chinensis）灌叢、黑枸杞（Lycium ruthenicum）灌叢等植被。

圖1 研究區Fig. 1 Study area

2 數據與方法

2.1 研究數據

2.1.1 MISR數據 MISR提供9個角度的觀測信息，分別為4個前向觀測角AF（26.1°）、BF（45.6°）、CF（60.0°）、DF（70.5°），4個后向觀測角AA（26.1°）、BA（45.6°）、CA（60.0°）、DA（70.5°）以及1個天頂角AN（0°）。每個傳感器均有藍光、綠光、紅光和近紅外4個波段。MISR 9個相機中，AN相機4個波段為275 m空間分辨率數據，其余8個角度相機紅光波段為275 m空間分辨率數據。

MI1B2T、MI1B2GEOP、MIL2ASAE 3種產品是目前最常用的數據。MI1B2T為TOA輻射數據，MIL2ASAE為17.6 km空間分辨率氣溶膠數據，MI1B2GEOP為17.6 km空間分辨率的太陽天頂角、太陽方位角以及觀測天頂角和觀測方位角數據。本研究用數據獲取時間為2014年8月25日，數據軌道號078111。使用MI1B2T 9相機紅光波段和MI1B2GEOP數據構建多角度觀測數據集。

2.1.2 無人機影像 為了對MISR影像數據提取的胡楊對象進行高度定標和數據驗證，在研究區設置3個無人機驗證區。2018、2019、2020年7月，在觀測區使用大疆精靈4 Pro/RTK進行無人機傾斜攝影測量，相機搭載1英寸2 000萬像素CMOS傳感器。基本過程為：規劃航線，設置飛行高度110 m， 航向重疊度80%、旁向重疊度75%，通過5向傾斜攝影方式獲取研究區地面影像數據。

2.1.3 地面實測數據 2018年7月，在研究區進行胡楊結構參數實地調查，位置如圖2。

圖2 地面采樣示意Fig. 2 Schematic diagram of field-measurement

樹高和冠幅調查方法：選取49株不同生態型胡楊，應用激光測距儀的測高功能從3個不同方向測量樹冠最高點，獲取胡楊樹高，取平均值作為實測樹高；使用長卷尺分別獲取相隔60°的3個方向樹木冠幅，取平均值；利用OruxMaps地圖軟件記錄每株樹坐標，并在地面分辨率0.5 m的WorldView-2多光譜影像地圖上進行位置標注，方便后期與無人機影像比對。

林分密度調查方法：利用OruxMaps地圖軟件，以WorldView-2 影像為參照，在地圖中設置 25個100 m×100 m 采樣格網，實地對每個格網中的胡楊進行株樹統計，并在地圖軟件中標記。

樹高、冠幅和林分密度實測值描述統計如表1 所示。

表1 采樣點的描述統計Tab. 1 Descriptive statistics of sampling points

2.2 數據處理與分析

2.2.1 無人機影像處理 采用Pix4d Mapper攝影測量軟件對無人機影像進行自動配準、空三加密等處理后生成點云數據（點云密度為74點·m-3），點云數據經分類得到地表和植被點云，再插值生成分辨率為0.5 m的數字表面模型（digital surface model，DSM）和數字高程模型（digital elevation model，DEM），進而得到冠層高度模型（canopy height model，CHM）（圖3）。CHM數據首先使用eCognition軟件的“Min/Max Filter”最小/最大值濾波算法（主要參數：Mode: Diff. brightest to center；2D kernel Size: 7）生成種子點，根據種子點位置提取樹高，然后使用“Pixel-Based Object Resizing”增長算法（主要參數：mode: Growing；Pixel Layer Constraint: 2）生成植被冠層矢量多邊形對象，繼而通過ArcGIS計算、統計出冠幅和密度。提取與實測49株胡楊樣本位置對應的胡楊對象，二者樹高、冠幅比較發現，樹高R2為0.90，RMSE為0.45 m；冠幅R2為0.84，RMSE為0.68 m；對與25個采樣格網重疊CHM部分分別提取胡楊對象數量，與格網內實測統計株數比較，R2為0.94，RMSE為4.25株·hm-2。這表明，無人機獲取的胡楊林結構數據可代替地面調查數據，作為SGM的定標參數以及最終精度驗證。

圖3 無人機CHM及MISR樣地Fig. 3 CHM obtained by UAV and MISR samples extenta. 無人機采樣區1 1st sample area of UAV b. 無人機采樣區2 2nd sample area of UAV c. 無人機采樣區3 3rd sample area of UAV

2.2.2 MISR影像處理 1） 預處理 基于MI1B2T TOA BRF輻射數據，對9個相機的紅光波段TOA反射率采用SMAC（simplified method for the atmospheric correction）大氣校正方法（Rahmanet al., 1994）轉換為地面反射率。通過提取研究區影像數據并進行投影轉換，與無人機影像配準后提取重疊區域的MISR像元大小（275 m2）樣本47個（圖3），隨機選擇30個作為訓練集，剩余17個作為測試集。

2） AMBRALS核驅動模型反演 AMBRALS（Wanneret al., 1995）是一種半經驗核驅動線性模型，其將場景反射視為各向同性反射、體散射和表面散射的組合，數學形式表示為：

式中：R為反射率；kvol為體散射核；kgeo是幾何光學核；fiso、fvol和fgeo分別為各向同性反射、體散射和表面散射的權重系數；i、v、φ分別為光照天頂角、視天頂角和相對方位角。

fiso、fvol和fgeo是模型的核心參數，可反映出對應場景基本的反射特征，本研究使用這3個系數構建場景土壤背景回歸方程。

根據研究區大部分地區林地灌木稀疏、草本稀薄的特點，反演選擇RossThin核為體散射核，LiSparseModis核為幾何光學核。

構建多角度反射數據集，使用AMBRALS反演得到47個樣本區域對應的fiso、fvol和fgeo數據（圖4）。

圖4 SGM模型反演流程Fig. 4 The inversion flow chart based on SGM model

3） SGM 模型反演 SGM作為一種物理模型，使用不同視角觀測的土壤背景反射和植被反射2部分的組合表示整個場景反射，其數學表達為：

式中：R為反射率；Gwalthall表示采用Walthall模型模擬土壤背景反射；CRoss表示采用RossThin冠層體散射模型模擬植被反射；kG和kc分別表示可見光照地面和可見光照樹冠比例；?i、?v、φ表示光照天頂角、視天頂角和相對方位角。

其余部分數學表示如下：

式中：c1～c4為Walthall模型系數；?i、?v、φ為變換后的光照天頂角、視天頂角和相對方位角；ξ為光照天頂角、視天頂角和相對方位角的三角函數。

式中：λ為林分密度；γ為樹冠水平半徑；ε′為變換后的散射相位角；O為視野和照明陰影部分重疊區域。

SGM模型反演分為3個步驟：

1） 訓練集樣本的多角度反射率觀測數據和森林幾何參數作為已知固定變量輸入SGM模型，使用GRG（Facóet al.,1989）優化算法反演得到Walthall模型系數（圖4A）；

2） 使用訓練集的AMBRALS模型系數和AN相機多光譜反射率數據，建立Walthall模型系數的回歸方程（圖4B）；

3） 使用Walthall模型系數回歸方程獲得測試集土壤背景反射，SGM和GRG算法反演獲取測試集的幾何參數（圖4C）。

3 結果與分析

3.1 背景反射模擬及精度分析

塔里木河河岸林郁閉度低，高亮土壤背景與暗深河水對比反差較大，再加上場景內灌叢和草本植被，SGM模型能否準確模擬土壤背景反射尤為重要。

反演第一階段使用SGM模型將觀測反射強度分解為土壤背景反射和冠層反射2部分，背景反射模擬采用Walthall模型，冠層反射模擬采用RossThin模型。輸入數據包括無人機獲取的訓練集樣地結構參數（樹高、覆蓋度、冠幅、林分密度）、模型固定參數LAI（固定值為2.3）（巴比爾江·迪力夏提等，2012）、HB（固定值為2.0）、BR（固定值為1.0）（圖5）。

圖5 SGM模型樹冠形態參數Fig. 5 The crown shape parameters of SGM model

反演以觀測反射值與模擬反射值的均方根誤差（root mean square error，RMSE）最小化為目標函數（式7）（圖6a），最終輸出Walthall的c1～c4參數值：

圖6 背景反射模擬與評價Fig. 6 Simulation and evaluation of Walthall model parametersa. 背景反射模擬Simulation of background reflection；b. 背景反射評價Evaluation of background reflection.

以上處理過程使用隨機選擇的30塊訓練樣地，模擬反射值與觀測反射值R2最大為0.99，最小為0.72，均值為0.92，大部分樣地的反射值R2為0.90以上（圖6b）。R2較低和RMSE較大的多為河床附近、覆蓋度較大的樣地，原因可能因河水背景影響，以及胡楊呈現出聚生、叢生狀態，在此情形下，SGM模擬誤差較大。

對訓練樣地反演獲取的樹高和覆蓋度與參考值進行比較發現，覆蓋度與參考值的一致性較好，R2達0.99，RMSE為0.73%，樹高的一致性稍差，R2為0.66，RMSE為0.61 m（圖6b）。總體來看，SGM模型能較好模擬樣地各向異性反射特征。在反射率、結構參數等均在適當精度的情況下，獲取的c1～c4參數值描述的Walthall背景反射模型參數也視為能正確反映場景背景反射的各向異性特征。

3.2 Walthall模型參數敏感性分析

Walthall模型的4個參數中，c1描述背景反射的基本強度水平、c2描述前后向反射間的差異程度、c3描述天底角反射強度、c4描述天底角反射與前后向反射的差異程度。為評價Walthall模型參數對最終反演結果的影響，使用無人機獲取的樣地結構參數反演得到的c1～c4值為參考，分別對其中1個參數值進行±20%范圍變化，同時固定其他3個參數，反演得到覆蓋度、樹高、冠幅和林分密度，與無人機測得的結構參數進行平均相對誤差（mean relative error，MRE）計算比較結果見圖7。

可以觀察到，斜率越大的參數對變量影響越大。c1值偏離參照值±10%，覆蓋度誤差增加20%～40%；偏離參照值±20%，覆蓋度誤差增加60%～80%；其次是c4，c2和c3影響最小（圖7a）。由于模型中樹高由冠幅計算得到，c1～c4變化對樹高和冠幅的影響相同。相對覆蓋度而言，樹高、冠幅對c1～c4在正負變化方向上的響應比較復雜（圖7b、c）。總的來說，影響程度最大的還是c1，其次是c4，c2和c3影響最小。c1～c4對林分密度值變化影響相對最小（圖7d）。實際結構參數變化受到c1～c4變化的疊加影響。

3.3 Walthall模型參數回歸及精度分析

通過建立訓練集AN相機的多光譜（Blue、Green、Nir）反射率、AMBRALS模型參數fiso、fvol、fgeo和c1～c4參數的多元線性回歸方程，用于對測試集樣地的背景反射進行預測。

使用多元線性回歸方程預測的訓練集樣地c1～c4參數，與反演值進行對比，c1～c4參數的調整R2分別為0.78、0.98、0.59和0.75（圖8）。

圖8 Walthall參數回歸建模Fig. 8 Regression of Walthall model parametersa. c1；b. c2；c. c3；d.c4.

3.4 測試集森林結構參數反演及精度分析

輸入17個測試集AN相機的多光譜（Blue、Green、Nir）反射率、AMBRALS模型參數fiso、fvol、fgeo，使用Walthall模型參數多元線性回歸方程，計算出測試集樣地的背景反射，作為SGM模型的輸入信息，其他輸入數據還包括模型固定參數LAI、HB、BR值。

以測試集觀測反射值與模擬反射值的RMSE最小化為目標函數（式7），使用GRG優化算法計算輸出17個測試集樣地的樹高、覆蓋度、林分密度和冠幅等參數。

17個測試集樣地反演結果與無人機測量獲取的數值比較如圖9。

圖9 測試集結構參數評價Fig. 9 Evaluation of structural parameters of test dataseta. 覆蓋度FVC；b. 樹高Height；c. 林分密度Density；d. 冠幅Crown width.

在測試樣地（275 m2）尺度上，與無人機獲取的結構參數相比，反演獲取的測試集結構參數使用線性模型擬合時，覆蓋度、樹高、林分密度、冠幅的R2分別為0.54、0.47、0.41、0.24，RMSE分別為3%、0.76 m、112株、0.64 m，MRE分別為24.7%、8.9%、22.5%、10%；使用冪函數模型時，覆蓋度、樹高、林分密度、冠幅的R2分別為0.80、0.53、0.55、0.30，RMSE分別為1%、0.4 m、33 株、0.32 m，MRE分別為10%、5%、6%、6%（圖9）。

4 結論與討論

基于MISR多角度遙感觀測數據和SGM模型，通過模擬塔里木河下游河岸胡楊林的二向反射特性，分解后獲得地面背景反射特征，最終反演得到樣地主要森林結構信息。

研究得到主要結論如下：

1） UAV獲取樹高、冠幅、林分密度與實測值的R2分別為0.90、0.84、0.94，RMSE為0.45 m、0.68 m、4.25株·hm-2，無人機獲取的胡楊林結構數據可用作SGM模型的定標參數以及最終精度驗證數據。

2） SGM對訓練集樣地的模擬反射值R2最大值為0.99、最小值為0.72、均值為0.92，模擬覆蓋度R2為 0.99，RMSE為0.73%，樹高R2為0.66，RMSE為0.61m。SGM模型能較好模擬胡楊林樣地的各向異性反射特征。

3） 對Walthall模型參數進行敏感性分析發現，c1參數對覆蓋度和樹高的影響最大，其次是c4參數，c2和c3參數影響最小，c1～c4參數對林分密度值變化影響很小。

4） 多元線性回歸方程預測的c1～c4參數的調整R2分別為0.78、0.98、0.59和0.75，說明使用AN多光譜反射率和AMBRALS模型參數可以較好描述Walthall模型參數。

5） SGM模型反演獲取的Walthall模型參數精度以及由此建立的參數回歸方程精度對于最終的結構參數反演結果有較大影響。

6） 使用多角度衛星遙感（MISR）反演獲取的胡楊覆蓋度、樹高、林分密度和冠幅的線性模型R2分別為0.54、0.47、0.41和0.24，冪函數模型R2分別為0.80、0.53、0.55和0.30，說明使用多角度衛星遙感（MISR）和反演方法可以實現區域尺度上的森林調查。

研究中存在的問題：

1） SGM反演需要獲取與衛星影像像元尺度對應的森林結構數據作為模型定標數據，這需使用地面實測或無人機LiDAR、無人機攝影測量等技術獲取高精度的森林調查數據。由于MISR像元與森林調查數據存在空間定位誤差，加之分辨率相差較大，造成二者之間空間配準誤差難以確定。干旱區環境下，河岸林植被覆蓋度本身較低且空間分布異質性較大，較小的空間誤差可能造成對應結構參數的大幅度改變，也會間接造成背景反射模型誤差增大，最終影響到反演精度。對于這一問題，目前尚無較好的解決辦法。

2） 通過對Walthall模型參數的敏感性和回歸方程效果分析，發現背景反射模擬精度對反演結果有直接影響。研究區較為復雜的背景反射（高亮土壤、暗色河流和冠層反射與胡楊接近的灌木等）均增大背景反射模擬的難度和不確定性。如何獲得更準確的Walthall模型標定參數以及在此基礎上建立更有效精確的回歸模型，是提高最終反演結構數據精度的關鍵。目前情況下，有限的定標站點不能代表所有實際場景，因此，基于樣地提取的Walthall參數代表性也受到限制，在此基礎上產生的參數回歸方程會進一步增大該誤差。一個可能的改進辦法是采用DART（discrete anisotropic radiative transfer model）等仿真模型，根據場景組分實測光譜，使用不同組分組合方案，重建MISR尺度大小的胡楊林場景。模擬該尺度場景下中等分辨率光學衛星（如Landsat、Sentinel-2）多光譜反射( Mortonet al., 2015; Gastellu-Etchegorryet al., 2015），建立與背景反射模型參數的回歸關系，生成Walthall模型參數的參數矩陣查找表（look up table, LUT），然后使用真實Landsat、Sentinel-2多波段反射率回歸生成背景反射，提高對背景的模擬精度。

3）SGM模型反演結構參數的一個重要部分是優化算法，本研究使用的GRG是處理非線性約束問題的一種有效算法，但是 GRG 無法保證最后的尋優結果是全局最優值。當針對多極值問題時，其最終校準結果是一個臨近初始給定值的局部最優值（李川等，2014）。與常規優化演算法相比，遺傳演算法等近些年來出現的一些模仿自然選擇與進化的隨機搜索演算法只要能夠提供合適的突變概率和運算時間，就能夠確保其很好地遍歷整個參數空間，最后尋優結果也往往能夠落在一個較優的局部參數區域內。因此，使用更好的優化算法是提高反演精度的可行辦法之一。