張海波,汪長城,2，朱建軍，付海強(qiáng)

1. 中南大學(xué)地球科學(xué)與信息物理學(xué)院，湖南 長沙 410083； 2. 中南大學(xué)有色金屬成礦預(yù)測教育部重點(diǎn)實(shí)驗(yàn)室，湖南 長沙 410083

森林覆蓋了地球陸地30%左右的面積，在維持陸地生態(tài)平衡、保護(hù)生態(tài)安全、防止生態(tài)危機(jī)方面起著決定性的作用[1-2]。同時(shí)，森林也是陸地上最經(jīng)濟(jì)的吸碳器和最大的碳儲庫[3],森林生物量表征著森林生態(tài)系統(tǒng)中的碳儲量[4]，因此，精確估測生物量特別是大尺度范圍估測，能使人們更好地理解全球碳循環(huán)和全球變暖。一般來說，生物量包括地上生物量(above-ground biomass，AGB)和地下生物量，由于收集地下生物量的實(shí)測數(shù)據(jù)比較困難，因此，當(dāng)前的研究大多針對地上生物量[5]。

目前森林生物量估算的方法主要有傳統(tǒng)的地面調(diào)查法和遙感監(jiān)測手段。傳統(tǒng)的地面調(diào)查方法可直接獲得單木樹高、胸徑、冠幅等參數(shù)，通過異速生長方程可獲取森林生物量[6]。這種方法需要耗費(fèi)大量的人力、物力和財(cái)力，對森林破壞性大，并且只能局限于小范圍，不適用于區(qū)域、全國乃至全球尺度的生物量估測與動態(tài)變化監(jiān)測[7-8]。與之相比，具有大范圍、快速、準(zhǔn)確的遙感技術(shù)成為了獲取地上生物量的主要途徑，遙感數(shù)據(jù)已廣泛應(yīng)用于區(qū)域尺度植被動態(tài)變化監(jiān)測研究[9-10]。其中，雷達(dá)遙感因其具有全天候、強(qiáng)穿透力、對地表粗糙度敏感等優(yōu)于光學(xué)遙感的特性，更有利于生物量的估測[11]。

近年來，利用雷達(dá)遙感數(shù)據(jù)估測生物量取得了許多進(jìn)展，一些研究表明長波長SAR(L、P波段)數(shù)據(jù)后向散射系數(shù)在評估地上生物量方面有著相當(dāng)大的潛力[12-14]。但是，由于SAR特殊的側(cè)視成像方式，其后向散射容易受復(fù)雜地形的影響，其中明顯的影響是降低了雷達(dá)的頻率[15],而且地形坡度在SAR影像中是很難處理的。目前在降低復(fù)雜地形對SAR影響方面的研究主要是根據(jù)面積積分[16]或利用投影角進(jìn)行地形輻射校正[17-19]。該方法在以單次和二次散射為主的非植被區(qū)域能較好地校正地形引起的后向散射異常，但對于結(jié)構(gòu)復(fù)雜、內(nèi)部存在多次散射現(xiàn)象的森林區(qū)域，難以描述地形變化引起的多次散射的變化，因此對于森林參數(shù)定量分析改善不明顯。在森林AGB估測模型中融入地形因子的方法，不但考慮了地形對后向散射的綜合影響，而且不需要考慮地形對森林散射機(jī)制的具體影響[20]。目前，該方面的研究相對較少。

以往的研究中，缺乏地形因子對估算模型估算森林AGB的能力及穩(wěn)定性影響的分析，地形變化對后向散射強(qiáng)度與森林AGB的響應(yīng)關(guān)系的影響同樣缺乏可靠分析。因此，本文利用機(jī)載P波段全極化SAR數(shù)據(jù)，分析地形變化與不同極化后向散射的關(guān)系，找到相關(guān)性最高的極化方式用于森林AGB估算，建立適用于地形復(fù)雜區(qū)森林AGB估算的模型，并通過不同模型的比較，驗(yàn)證模型的適用性與穩(wěn)定性，同時(shí)利用該模型對研究區(qū)森林AGB進(jìn)行估算。文中的分析方法除了能適用山地、丘陵地區(qū)森林AGB估測之外，對地形復(fù)雜區(qū)的農(nóng)作物估產(chǎn)方面同樣具有一定的借鑒意義。

1 材料與方法

1.1 研究區(qū)概況

研究區(qū)位于瑞典北部西博滕省溫德恩市(Vindeln municipality)的Krycklan流域，處于64°11′E—64°17′E，19°41′N—19°53′N之間，是溫德恩市森林實(shí)驗(yàn)場的一部分，同時(shí)也是瑞典大學(xué)野外森林研究基地，面積約為9390 hm2，屬于亞寒帶針葉林氣候。林區(qū)地形起伏較大，海拔高度范圍為100～400 m。森林類型為自然生長的混交林，其中大多數(shù)為針葉林樹種，包括云杉(Spruce)，松樹(Pine)，落葉闊葉林樹種有樺樹(Birch)，除此之外還有其他少量的落葉喬木如山楊(Aspen)、花楸(Rowan)等。研究區(qū)如圖1所示。

圖1 研究區(qū)及林分分布Fig.1 The location of test site and the distribution of stand

1.2 森林實(shí)測數(shù)據(jù)

森林實(shí)測數(shù)據(jù)來自于研究區(qū)24個(gè)林分采樣區(qū)，樣區(qū)分布見圖1，面積為3.07～24.34 hm2，在每個(gè)林分中根據(jù)林分面積按照一定的間距(50～160 m)布置8到13個(gè)半徑為10 m的樣地，樣地總數(shù)為310個(gè)，樣地調(diào)查時(shí)間為2008年10月13日—17日。在這些樣地中對胸徑大于4 cm樹木，逐一測定單木的樹種、樹高和胸徑，并記錄其他一些參數(shù)，如植被種類、土壤類型等。各樣地中枯死的樹木及落葉不參與生物量計(jì)算，因此不予記錄。實(shí)地測量的單木數(shù)據(jù)根據(jù)不同林分區(qū)域的樹種類型先匯總成樣地尺度的森林參數(shù)，包括每個(gè)樹種的斷面積加權(quán)平均年齡、胸高斷面積加權(quán)平均直徑(cm)、胸高斷面積加權(quán)平均樹高(dm)。不同樹種AGB信息參考H. Perersspm的計(jì)算方式，將每個(gè)樹種的生物量分成3個(gè)組分(樹干、樹枝、樹葉)分別計(jì)算。由于季節(jié)原因，落葉林不計(jì)算樹葉生物量，落葉喬木歸于樺樹一類統(tǒng)計(jì)。

1.3 DEM數(shù)據(jù)

DEM數(shù)據(jù)來源于歐空局BIOSAR項(xiàng)目，不包含植被信息。首先對其進(jìn)行插值和重采樣，以獲取規(guī)則的樣地DEM數(shù)據(jù)，處理之后的數(shù)據(jù)分辨率為1 m×1 m，然后提取出研究區(qū)坡度。由于SAR主要受距離向的坡度影響[21]，因此參照飛行方向，規(guī)定坡面朝向雷達(dá)視線方向的坡度為正，背向雷達(dá)視線方向的坡度為負(fù)。雷達(dá)入射角、局部入射角和地形坡度的幾何特征及關(guān)系如圖2所示，其中θi是雷達(dá)入射角，θl是局部入射角，θs是坡度。

圖2 坡度幾何特征Fig.2 Geometry characteristics of terrain slope

局部入射角可由地形坡度與雷達(dá)入射角獲得，其計(jì)算公式如下

θi=arccosH/SR

(1)

θl=θi-θs

(2)

式中,SR為斜距距離，是斜距方向上近距起點(diǎn)的斜距值加上像素?cái)?shù)乘以斜距向采樣間距，起點(diǎn)斜距為4 294.21 m；H是飛機(jī)的飛行高度。研究區(qū)坡度與局部入射角見圖3。

圖3 研究區(qū)坡度和局部入射角Fig.3 The slope and local incidence angel of the study area

1.4 SAR數(shù)據(jù)

SAR數(shù)據(jù)來源于歐空局BIOSAR項(xiàng)目，是由德國宇航局(DLR)應(yīng)用機(jī)載E-SAR傳感器獲取的P波段全極化數(shù)據(jù)(HH、HV、VH、VV)，數(shù)據(jù)支持平臺高度為4 090.1 m，飛行速度為92 m/s，中心頻率為0.35 GHz，天線俯視角為40°，航向角為-47.2°，脈沖發(fā)射頻率(PRF)為2000 Hz。影像成像時(shí)間為2008年10月14日，數(shù)據(jù)類型為單視復(fù)數(shù)影像(SLC)，影像的重疊度為50%。在經(jīng)過多視處理(2×1)、濾波、地理編碼，重采樣等SAR影像基本處理之后，最終的影像坐標(biāo)系統(tǒng)為UTM WGS-84，像元大小為1 m×1 m，SAR影像的基本處理由DLR完成。根據(jù)雷達(dá)系統(tǒng)的互易性，HV極化方式與VH極化方式近似相等，因此，本文從HH、HV、VV這3種極化方式的影像數(shù)據(jù)中獲取后向散射系數(shù)，其獲取公式如下

(3)

式中，i為P波段不同極化方式的SAR影像(HH、HV和VV)；m為林分中的像元；n為林分；k和c是校正常量，分別為32 769和10-6(文獻(xiàn)[22])；sinθ是幾何投影因子；θ為雷達(dá)視角，式(3)可通過10lgσ0將其轉(zhuǎn)換為dB值。

1.5 估算模型

在對植被生物量的估算中，半經(jīng)驗(yàn)?zāi)Ｐ屯軌蛱峁┣逦慕Y(jié)構(gòu)和思路來理解不同后向散射相互作用的物理機(jī)制，能夠?qū)ρ芯繉ο筇峁┖侠淼慕馕鯷23-24]。本文采用的半經(jīng)驗(yàn)?zāi)Ｐ褪且愿倪M(jìn)的水云模型為前提，公式如下

(4)

(5)

式中，β為需通過經(jīng)驗(yàn)定義的系數(shù)；V為蓄積量。將式(5)代入式(4)可得式(6)

(6)

對式(6)進(jìn)一步參數(shù)化，由于文中估算的參數(shù)為森林AGB，因此，將公式中的V用AGB替換，如下

(7)

本文在上述分析模型的基礎(chǔ)上引入局部入射角，得到最終模型(水云分析模型)如下

σ0cos-1(θi-θs)=β1+β2eβ3AGB

(8)

為驗(yàn)證文中水云分析模型在地形復(fù)雜區(qū)估算AGB的適用性，文中使用顧及地形因子的線性模型、對數(shù)模型、二次模型、指數(shù)模型與水云分析模型對比分析。為方便表述，令γ0=σ0cos-1(θi-θs)，各模型具體形式見表1。

表1 森林AGB估算模型

1.6 模型參數(shù)確定方法

構(gòu)建森林AGB估算模型后，需要通過樣地實(shí)測數(shù)據(jù)與后向散射系數(shù)值確定模型的參數(shù)值。由于使用的估算模型多為非線性模型，難以用傳統(tǒng)方法確定模型的參數(shù)值，而遺傳算法(GA)適用范圍較廣，無論非線性問題能不能轉(zhuǎn)化為線性問題，都可以直接進(jìn)行數(shù)據(jù)擬合，并且對于可以轉(zhuǎn)化為線性問題的非線性數(shù)據(jù)問題，遺傳算法要優(yōu)于其他算法[27]，因此本文使用遺傳算法來確定生物量估算模型的最優(yōu)參數(shù)值。遺傳算法采用了“適者生存”的原則來獲得最佳解決方案，基本操作包括了選擇、交叉和變異，染色體的基因?qū)?yīng)于模型的待求參數(shù)，進(jìn)化過程中的適應(yīng)度函數(shù)為實(shí)測生物量值與模型計(jì)算生物量值之間的誤差，具有最小誤差的參數(shù)組合表示能很好地適應(yīng)環(huán)境。本文方法流程見圖4。

圖4 方法流程Fig.4 Methodology flowchart

2 結(jié)果與分析

2.1 不同極化方式后向散射與森林AGB響應(yīng)關(guān)系

研究區(qū)不同極化方式的后向散射系數(shù)值與對應(yīng)生物量值的相關(guān)性分析結(jié)果見表2。由表2可知，在3種極化方式中HH與HV極化方式后向散射系數(shù)與生物量的相關(guān)性要高于VV極化方式，其中HV極化方式后向散射系數(shù)與生物量呈強(qiáng)相關(guān)性(0.6～0.8)，相關(guān)性達(dá)到0.696，而VV極化方式后向散射系數(shù)與生物量的相關(guān)性最弱，為-0.052。HH和HV極化方式后向散射與生物量呈正相關(guān)，表明在HH和HV極化方式下，后向散射系數(shù)值隨著地面生物量的增加而增加，VV極化方式后向散射與生物量呈負(fù)相關(guān)，表明在VV極化方式下，后向散射系數(shù)隨著生物量的增加呈減少趨勢。這種現(xiàn)象與地面植被覆蓋的情況相關(guān)，在植被覆蓋較高且為大型植被的區(qū)域，交叉極化方式對生物量的敏感程度高于單極化方式，在植被覆蓋較低的區(qū)域，則會出現(xiàn)相反的情況，尤其對于VV極化方式最為明顯，其在高覆蓋的大型植被區(qū)對生物量最不敏感，然而在稀疏低矮的植被區(qū)對生物量的敏感性較好。

表2 后向散射系數(shù)與森林AGB的相關(guān)性

注：* *表示在0.01水平顯著相關(guān)。

從不同極化方式后向散射系數(shù)值的動態(tài)范圍來看，極化方式不同，對應(yīng)的后向散射系動態(tài)范圍不同(圖5)。HV極化方式與HH、VV極化方式相比要集中，動態(tài)范圍為3.94 dB(-16.98～-13.04 dB)，而HH、VV極化方式分別為9.43 dB(-11.41～-1.98 dB)和5.85 dB(-10.79～-4.94 dB)。相同極化方式中不同生物量水平對應(yīng)的后向散射范圍不同，在HV極化方式中，生物量低于120 t/hm2，其后向散射系數(shù)的動態(tài)范圍為3～4 dB，生物量高于120 t/hm2，后向散射系數(shù)動態(tài)范圍僅為1～2 dB。由圖5還可知，HH、HV和VV極化方式的后向散射系數(shù)變化趨勢在森林AGB處于較低值的情況下能與AGB變化趨勢保持一致，但隨著AGB值的增大，這種一致性僅在HV極化方式下繼續(xù)保持，HH和VV極化方式均出現(xiàn)高AGB值對應(yīng)低后向散射系數(shù)值，并且比低AGB值對應(yīng)的后向散射系數(shù)值還要低的情況。

通過以上對比分析可知在不同極化方式中，HV極化方式更適用于復(fù)雜地形區(qū)的森林AGB估算，因此本文后面的分析主要針對HV極化方式展開。

圖5 后向散射系數(shù)與地上生物量散點(diǎn)圖Fig.5 The scatter plot between backscatter coefficient and AGB

2.2 估算模型不同坡度的適用性

為了更深入了解不同地形對估算模型的影響，以及對比不同模型在不同坡度情況下的適用性，本文將坡度分為3個(gè)等級，分別為0～5°、5°～10°、≥10°。每個(gè)等級對應(yīng)的林分樣地?cái)?shù)目分別為24、18和11個(gè)。分級后的樣地具體情況見表3。

表3 不同坡度等級樣地情況

由圖6可知，在坡度小于5°的地區(qū)，5種估算模型決定系數(shù)分別為0.443、0.539、0.552、0.501、0.57，均為3種坡度等級中最高，隨著坡度的增加，決定系數(shù)均有下降，在坡度≥10°的區(qū)域，5種模型的決定系數(shù)達(dá)到最低，分別為0.13、0.197、0.32、0.16、0.424。這表明在復(fù)雜地形區(qū)，地形因子對森林AGB的估算具有極大的影響，后向散射系數(shù)與森林AGB的相關(guān)程度隨地形坡度的增加而減小。這可能是由于在地形起伏較小的區(qū)域，電磁波與植被層的相互作用構(gòu)成的“水云層”在整個(gè)區(qū)域較穩(wěn)定，沒有發(fā)生實(shí)質(zhì)的變化，而隨著地形坡度的增加，降低了雷達(dá)頻率，“水云層”結(jié)構(gòu)發(fā)生改變，電磁波與地面的相互作用增加，這些易使電磁波與植被之間的相互作用出現(xiàn)異常。

在不同坡度等級中，5種模型估算AGB的決定系數(shù)大小排序是相同的，均為水云分析模型>二次模型>對數(shù)模型>指數(shù)模型>線性模型，這表明水云分析模型在復(fù)雜地形區(qū)估測生物量的適應(yīng)能力最強(qiáng)，線性模型適應(yīng)能力最弱。不同坡度等級變化中，各模型估算森林AGB的決定系數(shù)變化值見表4，坡度由0～5°變化至5°～10°中，水云分析模型的決定系數(shù)變化值最小為0.066，線性模型最高為0.153，變化值大小排序?yàn)榫€性模型>指數(shù)模型>對數(shù)模型>二次模型>水云分析模型。坡度由5°～10°至≥10°中，變化值最小的為水云分析模型，為0.08，最大值為對數(shù)模型，為0.206，不同模型變化值大小排序?yàn)閷?shù)模型>二次模型>線性模型>指數(shù)模型>水云分析模型。坡度由0～5°至≥10°中，變化值最小的為水云模型，變化值為0.146，變化值最大的為對數(shù)模型，為0.342，不同模型變化值大小排序?yàn)閷?shù)模型>指數(shù)模型>線性模型>二次模型>水云分析模型。這表明，無論在地形起伏較小的地區(qū)還是地形起伏較大的地區(qū)，水云分析模型估算森林AGB的可靠性和穩(wěn)定性均高于其他4種模型。

表4 不同坡度過渡中決定系數(shù)變化值

注：(a)、(b)、(c)、(d)、(e)分別指代圖6中對應(yīng)的模型

2.3 森林AGB估算結(jié)果

由2.2節(jié)分析可知，不同坡度等級內(nèi)水云分析模型是5種模型中最適用于復(fù)雜地形區(qū)森林AGB估算的。然而在實(shí)際應(yīng)用中，如區(qū)域性森林AGB估算和生物量制圖，一個(gè)林分樣地往往包含了多種坡度等級，其后向散射系數(shù)值是多類坡度下后向散射的結(jié)果。因此為完成研究區(qū)生物量制圖，先利用研究區(qū)24個(gè)完整林分的實(shí)測生物量值與對應(yīng)的HV極化方式后向散射系數(shù)值，通過遺傳算法得到估算模型最優(yōu)系數(shù)值，結(jié)果見表5。

圖6 不同坡度中模型擬合結(jié)果比較Fig.6 Comparison of models fitting results in different slopes

表5 森林AGB估算結(jié)果

由表5可知，在整個(gè)區(qū)域中，水云分析模型估算生物量的決定系數(shù)最高，為0.796，均方根誤差(RMSE)最小，為0.466，說明水云分析模型相對其他4種模型最適用于整個(gè)研究區(qū)生物量反演和制圖，為便于理解及生物量估算結(jié)果檢驗(yàn)，水云分析模型可寫成如下形式

(9)

確定適合復(fù)雜地形生物量估算模型之后，對模型精度檢驗(yàn)是必要的。首先將林分對應(yīng)的HV極化方式后向散射系數(shù)值代入模型中，估算出生物量值，然后計(jì)算出實(shí)測生物量值與模型估算生物量值的均方根誤差，結(jié)果見圖7。實(shí)測生物量與模型估算的生物量值決定系數(shù)為0.597，RMSE為30.876 t/hm2，擬合精度為77.40%。最終的AGB制圖結(jié)果見圖8，研究區(qū)生物量分布范圍大致為：0～260 t/hm2，平均生物量為76.01 t/hm2。

圖7 森林AGB估算結(jié)果精度評價(jià)Fig.7 Accuracy of the estimated forest AGB

3 結(jié) 論

本文分析了復(fù)雜地形區(qū)不同極化方式后向散射系數(shù)與森林AGB的相關(guān)性，首次將地形因子引入改進(jìn)的水云模型，并通過其他模型與之對比分析，采用遺傳算法確定模型的最優(yōu)參數(shù)，并對模型在不同坡度情況下的可靠性與穩(wěn)定性進(jìn)行了對比分析，結(jié)果顯示：

3種極化方式后向散射系數(shù)變化趨勢在森林AGB處于較低值的情況下能與AGB變化趨勢保持一致，但隨著AGB值的增大，這種一致性僅在HV極化方式下繼續(xù)保持，因此相比之下，HV極化方式更適用于復(fù)雜地形區(qū)生物量的估算。

圖8 森林AGB分布圖(t/hm2)Fig.8 Distribution maps of the forest AGB

地形對森林AGB的估算具有極大的影響，后向散射系數(shù)與AGB的相關(guān)性隨著地形坡度的增加而減小。在地形坡度為0～5°、5°～10°和≥10°內(nèi)，5種模型估算生物量的能力大小排序均為：水云分析模型>二次模型>對數(shù)模型>指數(shù)模型>線性模型。

模型在復(fù)雜地形區(qū)估算生物量的穩(wěn)定性隨坡度變化有所改變，當(dāng)坡度由0～5°變化至5°～10°時(shí)，水云分析模型的決定系數(shù)變化值最小,線性模型最高，穩(wěn)定性大小排序?yàn)椋核品治瞿Ｐ?二次模型>對數(shù)模型>指數(shù)模型>線性模型；坡度由5°～10°至≥10°中，變化值最小的為水云分析模型，最大值為對數(shù)模型，穩(wěn)定性排序?yàn)椋核品治瞿Ｐ?指數(shù)模型>線性模型>二次模型>對數(shù)模型；坡度由0～5°至≥10°中，變化值最大的為對數(shù)模型，穩(wěn)定性排序?yàn)椋核品治瞿Ｐ?二次模型>線性模型>指數(shù)模型>對數(shù)模型。雖然模型穩(wěn)定性排序因坡度變化而不同，但對于水云分析模型而言，無論在坡度起伏較小的地區(qū)還是較大的地區(qū)，其穩(wěn)定性都最高。

利用本文所述融入地形因子的水云分析模型估算的森林AGB與實(shí)測AGB值進(jìn)行精度檢驗(yàn)，其決定系數(shù)為0.597，RMSE為30.876 t/hm2，擬合精度為77.40%，比報(bào)告中使用Le Toan方法的精度有所提高(RMSE為35 t/hm2)，估算出研究區(qū)平均生物量為76.01 t/hm2。

本文采用林分實(shí)測數(shù)據(jù)用于生物量擬合及驗(yàn)證生物量估測模型，雖然林分樣地的數(shù)量有限，尤其是將坡度分級后，在坡度較大的區(qū)域，實(shí)測數(shù)據(jù)只有11個(gè)，但是使用實(shí)測數(shù)據(jù)在一定程度上能夠確保模型系數(shù)的正確性及檢驗(yàn)的準(zhǔn)確性。另外，文中試驗(yàn)數(shù)據(jù)為P波段數(shù)據(jù)，其波長長，穿透性強(qiáng)，在植被區(qū)信號能穿透植被直接與地面交互，這種情況下使用本文方法能夠明顯改善地形對后向散射的影響，提高生物量的估算精度。而對于短波波段，由于波段短，穿透能力弱，其信號主要與植被冠層交互，這種情況下，使用在模型中融入地形因子的方法對生物量估算精度的改善可能不明顯。本文研究方法除適用于山地、丘陵地區(qū)森林AGB估測之外，對地形復(fù)雜區(qū)的農(nóng)作物估產(chǎn)方面同樣具有一定的借鑒意義。