王新闖 吳金汝 陸鳳連 焦海明 張合兵

(河南理工大學(xué)測繪與國土信息工程學(xué)院 焦作 454000)

森林冠頂高是評價森林生態(tài)系統(tǒng)功能的核心參數(shù)之一，全球碳循環(huán)、森林生產(chǎn)力、氣候變化等研究均需要量化森林冠頂高。傳統(tǒng)的森林冠頂高參數(shù)測定方法需要耗費大量的人力、物力和時間，不利于大范圍地測定森林參數(shù)，而且在實際應(yīng)用中還存在測量儀器誤差、樣地空間分布代表性不足等問題(婁雪婷等， 2011)。近年來，被動光學(xué)遙感技術(shù)在森林參數(shù)估測領(lǐng)域應(yīng)用廣泛，但多用于獲取森林水平結(jié)構(gòu)信息，很少能準確反映森林垂直結(jié)構(gòu)信息，限制了遙感在估測森林冠頂高等方面的精度(邢艷秋等， 2008； 李立存等， 2011； 孫國清等， 2006)。激光雷達技術(shù)具有很強的穿透能力，可以高精度地獲取森林垂直結(jié)構(gòu)信息，搭載于冰、云和陸地高程衛(wèi)星(the ice, cloud and land elevation satellite，ICESat)上的地學(xué)激光測高系統(tǒng)(the geoscience laser altimetry system, GLAS)星載激光雷達ICESat-GLAS數(shù)據(jù)在森林垂直結(jié)構(gòu)參數(shù)中得到了廣泛應(yīng)用。

大光斑激光雷達信號是入射激光脈沖與森林、地面相互作用的結(jié)果，激光回波信號是由發(fā)射脈沖和激光光斑內(nèi)的森林、地表參數(shù)共同決定的，前者包括工作的波長、脈沖寬度、脈沖能量、光斑尺寸和記錄回波脈沖的時間間隔，后者包括激光光斑內(nèi)每株樹的位置、高度、樹種、樹冠大小和形狀、冠層的反射率以及地表的反射率、坡度、坡向等，這些因素都會對基于GLAS數(shù)據(jù)的森林冠頂高反演精度產(chǎn)生影響(龐勇等， 2006)。此外，以森林分布為主要特征的激光腳點處的陡坡、枯落木等地物信息造成波形數(shù)據(jù)失真、疊加，也增加了森林冠頂高的反演難度(劉經(jīng)南等， 2005)。因此，在基于GLAS數(shù)據(jù)提取森林冠頂高的研究中，首先要對原始波形進行去噪處理。近年來，國內(nèi)外學(xué)者在GLAS回波波形去噪方面做了很多努力。Iqbal等(2013)基于傅里葉變換對GLAS波形進行去噪并建立森林冠頂高估算模型，該方法能解釋研究區(qū)79%的森林冠頂高變化。邱賽等(2012)在基于小波變換的基礎(chǔ)上，選擇常用的小波基函數(shù)對GLAS波形進行去噪，結(jié)果發(fā)現(xiàn)sym7小波基比db1小波基去噪效果更好，均方根誤差為0.912 325，證明sym基函數(shù)更適合GLAS波形去噪。王愛娟等(2016)以吉林汪清區(qū)為例，提出了基于窗函數(shù)的林區(qū)GLAS數(shù)據(jù)去噪方法，結(jié)果表明布萊克曼窗函數(shù)對林區(qū)波形數(shù)據(jù)去噪效果最好，預(yù)測冠頂高度與實測冠頂高度的回歸精度由0.725 增至0.820 。以上研究表明，這些去噪算法在一定程度上降低了背景噪聲對GLAS回波波形的影響，但在地形復(fù)雜的地區(qū)適應(yīng)性不強，因此有必要對不同去噪算法進行研究，選取去噪效果更好的算法。目前關(guān)于高斯窗函數(shù)對波形去噪的研究較少，由于GLAS激光發(fā)射器發(fā)射的信號為高斯脈沖信號，因此本研究在前人研究的基礎(chǔ)上，對高斯窗函數(shù)的波形去噪效果進行研究，以進一步提高森林冠頂高反演精度。

由于GLAS按照較大尺度離散采樣，地面的坡度和粗糙度混淆了來自植被冠層的返回信號，且地面的坡度和粗糙度會使脈沖寬度變寬，從而影響GLAS估算植被冠頂高度的精度(Pangetal.， 2006)。邢艷秋等(2009)研究顯示，在0°～15°地形坡度范圍內(nèi)，GLAS 完整波形數(shù)據(jù)能夠合理地反演森林冠頂高度，解釋能力在51% ～90%之間，當坡度超過15°時，地形對GLAS 波形的影響很大，冠頂高估算模型對變量的解釋能力顯著下降。為了減弱地形的影響，Lefsky等(2005)開發(fā)了“地形指數(shù)法”，以波形長度和地形指數(shù)為自變量，提高了坡度較大地區(qū)森林最大冠頂高的提取精度，最終所建估測模型能解釋59% ～68%的森林最大冠頂高變化，并進一步利用森林冠頂高估測了森林生物量，所建模型R2達到0.73。Xing 等(2010)在長白山地區(qū)應(yīng)用全波形GLAS 數(shù)據(jù)，以脈沖寬度的對數(shù)和地形指數(shù)為獨立變量，對Lefsky等(2005)所建模型進行改進，改善了傾斜地面最大冠頂高估算精度，此模型能解釋在0°～30°坡面上56% ～92%的最大冠頂高變化。Hayashi等(2013)分析了不同因素對冠頂高估算的影響，認為GLAS 數(shù)據(jù)的獲取時間和信噪比以及坡度均會對森林冠頂高的估算結(jié)果產(chǎn)生顯著影響。在剔除因獲取時間和信噪比低導(dǎo)致估算誤差較大的數(shù)據(jù)后，采用地形指數(shù)法分坡度區(qū)間(小于10°和大于等于10°)進行建模，提高了最大冠頂高的估算精度。董立新(2008)改進了“地形指數(shù)法”，提出“質(zhì)心-地形指數(shù)法”，對坡度大于20°的地區(qū)應(yīng)用“質(zhì)心-地形指數(shù)法”，提高了坡度較大地區(qū)的森林最大冠頂高提取精度。湯旭光(2013)在地形指數(shù)法的基礎(chǔ)上，通過引入波形半能量高這一參數(shù)，提高了最大冠頂高的估算精度。在上述研究中，基于GLAS 估算森林冠頂高均需要區(qū)域的DEM 支持，為了消除對DEM 的依賴，緩解地形的影響，從而得到更加準確的冠頂高，Lefsky等(2007)在7個研究區(qū)基礎(chǔ)上利用3個波形指數(shù)(波形寬度、下端邊緣糾正系數(shù)和上端邊緣糾正系數(shù))修正了光斑內(nèi)森林冠頂高的變化和地形的影響，定量估算得到了平均冠頂高，并且估算精度較以往研究有所提高。前人研究表明，森林冠頂高模型參數(shù)的選取對森林冠頂高反演精度具有重要影響。盡管基于地形指數(shù)和僅利用波形參數(shù)建立冠頂高反演模型都在一定程度上減弱了坡度對冠頂高估算精度的影響，但目前結(jié)合地形指數(shù)和多種波形參數(shù)建立冠頂高模型估算冠頂高的研究較少。鑒于此，本研究嘗試基于地形數(shù)據(jù)和GLAS全波形數(shù)據(jù)提取的多種波形參數(shù)，對各林型分別建立波形參數(shù)模型、地形因子模型和全模型，檢驗不同模型對森林冠頂高的解釋能力，確立反演精度較高的冠頂高反演模型，提高冠頂高估算精度，為森林生物量估測等研究奠定數(shù)據(jù)基礎(chǔ)。

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

1.1 研究區(qū)概況

研究區(qū)設(shè)置在露水河林區(qū)及其附近區(qū)域，該林區(qū)(127°29′—128°02′E, 42°24′—42°49′N)位于吉林省撫松縣境內(nèi)，長白山西北麓。全局經(jīng)營區(qū)東西寬約40 km，南北長約50 km，總經(jīng)營面積121 295 hm2。東南部地勢起伏不大，比較平坦，西北部地勢起伏較大，南部位于起伏的熔巖高臺地，坡度為0°～53°。林區(qū)屬中溫帶大陸性氣候，氣溫較低，降水充沛，年平均氣溫0.9～1.5 ℃。研究區(qū)森林屬長白山頂級植物群落區(qū)系，主要為闊葉林和紅松(Pinuskoraiensis)闊葉混交林，大部分為天然林和天然次生林，其中成熟林居多，森林高大茂密，郁閉度較高，也分布有部分人工林。該區(qū)森林類型有闊葉林、針闊混交林和針葉林，三者面積比例大致為3∶1∶1。本研究樣地在不同坡度區(qū)間(0° ～5°、5° ～10°、10° ～15°、15° ～20°、20° ～25°、25° ～30°、>30°)均有分布，樣地林型及不同坡度區(qū)間的分布情況如表1所示。

表1 各林型樣地不同坡度區(qū)間分布Tab.1 Distribution of different slope intervals for each forest type

由表1可知，各林型樣地在不同坡度區(qū)間均勻分布，闊葉林在研究區(qū)分布范圍最廣，地形最為復(fù)雜?？傮w來說，所有實測數(shù)據(jù)在平緩地區(qū)和坡度較大地區(qū)均勻分布，在10°～25°坡度范圍內(nèi)分布最多，樣地地形復(fù)雜。

1.2 ICESat-GLAS數(shù)據(jù)

星載激光雷達ICESat-GLAS數(shù)據(jù)來自搭載于美國2003年發(fā)射的科學(xué)實驗衛(wèi)星ICESat上的第1臺星載激光雷達傳感器GLAS。星載激光雷達ICESat-GLAS數(shù)據(jù)免費獲取，可用來測量冰蓋高及其隨時間的變化、云層和氣溶膠的高度、植被的高度以及海冰的厚度等，是第1個用于連續(xù)全球觀測的星載激光測高系統(tǒng)(邢艷秋等， 2009)。GLAS搭載3個激光儀器，以91天重復(fù)軌道方案來執(zhí)行觀測任務(wù)，相鄰軌道間距在赤道附近為15 km，在緯度80°處間距為2.5 km，激光光斑直徑為60～70 m，光斑間隔為170 m。GLAS通過測量脈沖從衛(wèi)星到地面的往返時間來測量地面高度，標準的脈沖寬度為4 ns，相當于地面高度0.6 m，解壓后的1 000幀數(shù)據(jù)采樣間隔為1 ns，對應(yīng)地面高度為0.15 m，可以通過數(shù)據(jù)幀數(shù)來計算地面高度。本文使用研究區(qū)及周邊附近地區(qū)的2006年的L3F、L3G和2007年的L3I 的GLA01數(shù)據(jù)進行研究，數(shù)據(jù)從美國國家冰雪數(shù)據(jù)中心(http:∥nsidc.org/data/ice-sat/)下載。

1.3 樣地實測數(shù)據(jù)及其處理

1.3.1 樣地設(shè)置與調(diào)查 利用已經(jīng)搜集到的林區(qū)GLAS數(shù)據(jù)、數(shù)字林相圖、DEM 和森林資源二類調(diào)查資料進行樣地設(shè)置。露水河林區(qū)的有效GLAS光斑數(shù)據(jù)主要分布在研究區(qū)西北部，所有林型和坡度區(qū)間均有光斑分布?？紤]林型、坡度對估算冠頂高的影響和樣地的代表性，所設(shè)樣地涵蓋了3種不同林型(闊葉林、針葉林和針闊混交林)和不同坡度區(qū)間(0°～5°、5°～10°、10°～15°、15°～20°、20°～25°、25°～30°、>30°)的森林。由于GLAS光斑為60 m×70 m的橢圓形狀，因此在實際樣地測量中，為與GLAS光斑相對應(yīng)，以光斑中心為圓心，利用全站儀設(shè)置半徑為30 m的圓形樣地。在樣地內(nèi)，以樣地最大冠頂高樹木為中心設(shè)置半徑為7.5 m的圓形樣方，并在樣地不同冠頂高區(qū)間(以10 m為間距)設(shè)置半徑為7.5 m的圓形樣方。為避免樣方光譜信息互相干擾，樣方之間距離盡可能大。本文所選用的47個樣地數(shù)據(jù)于2015年8—9月實測得到，樣地空間分布如圖1所示。

樣地調(diào)查內(nèi)容主要包括GPS經(jīng)緯度、高程、森林類型、坡度、坡向、樣地情況、樹種、胸徑、樹高、郁閉度、葉面積指數(shù)和林齡等。部分森林參數(shù)詳細測量方法如下： 用高精度GPS 定位樣地和樣方中心坐標及其高程。對樣地和樣地內(nèi)的樣方中胸徑大于2 cm 的喬木進行每木檢尺，利用激光測高儀獲得樹高，通過胸徑尺測量距地面 1.3 m 處的立木直徑。測量胸徑時，胸徑尺與樹干垂直且緊貼樹干周圍。本文利用樣地內(nèi)與最高樹木相差2 m的所有樹高的平均值作為該樣地樹木實測最大冠頂高。

圖1 野外樣地分布Fig.1 Distribution of sample plots

1.3.2 樣地數(shù)據(jù)處理 由于樣地調(diào)查時間和獲取的GLAS數(shù)據(jù)時間不一致，為了減弱該期間因林木自然生長、森林采伐或自然災(zāi)害等因素導(dǎo)致的森林變化，首先與林業(yè)局管理人員進行溝通，查閱其2006—2014年有GLAS 數(shù)據(jù)分布區(qū)域的采伐記錄，對比與GLAS 數(shù)據(jù)時相一致的SPOT5 影像并實地調(diào)查期間的SPOT5 影像，剔除那些位于采伐區(qū)或受自然災(zāi)害影響區(qū)域的GLAS 數(shù)據(jù)； 對于因林木自然生長引起的冠頂高變化，根據(jù)當?shù)亓铸g與樹高和胸徑的關(guān)系予以校正，獲取與GLAS 數(shù)據(jù)時間一致的森林胸徑、樹高等基本數(shù)據(jù)。然后根據(jù)冠頂高數(shù)據(jù)，獲取各樣地和樣方內(nèi)的最大冠頂高。

2 數(shù)據(jù)處理與研究方法

2.1 技術(shù)路線

整個研究的開展主要包括遙感數(shù)據(jù)獲取與處理、樣地數(shù)據(jù)獲取與處理以及森林冠頂高模型建立3大部分。首先對GLAS波形進行預(yù)處理，將GLA01產(chǎn)品中記錄的激光雷達回波原始數(shù)據(jù)提取出來，由于激光脈沖受大氣環(huán)境和激光器本身性能的影響，其接收能量會發(fā)生變化，且地形也會改變波形能量，因此需要將波形數(shù)據(jù)進行標準化，即用每刻的接收能量除以回波總能量，標準化波形對應(yīng)的面積都為1，便于對比分析不同脈沖的波形； 其次采用布萊克曼窗和高斯窗2種窗函數(shù)對GLAS波形進行去噪處理； 然后利用SNR和RMSE 2種精度衡量指標，對2種波形去噪方法進行定量比較，優(yōu)選出去噪效果最優(yōu)算法； 最后對去噪后的波形進行波形參數(shù)提取，并根據(jù)DEM和提取后的波形參數(shù)建立冠頂高反演模型。具體方案及路線如圖2所示。

圖2 技術(shù)路線Fig.2 Technical roadmap

2.2 GLAS波形去噪算法

2.2.1 高斯窗函數(shù) 高斯窗函數(shù)是一種指數(shù)窗，主瓣較寬，因此其頻率分辨率低，無負的旁瓣。第一旁瓣衰減達-55 dB，常用來截短一些非周期信號，如指數(shù)衰減信號。

假設(shè)f(t)∈L2(R)，以g(t)作為窗函數(shù)的窗口傅里葉變換定義為：

(1)

式中：b為時間位移；ω為頻率位移；f(t)為在離散時間t的信號； WFR(ω,b)為t=b附近的頻譜特性；j為虛數(shù)單位。

令：

gω,b(t)=g(t-b)e-jωt，

(2)

則：

(3)

若g(t)的有效窗口寬度為Dt，則WFR(ω,b)給出的是f(t)在[b-Dt/2,b+Dt/2]時間區(qū)間內(nèi)的頻譜信息。Dt越小，對信號的定位能力越強。

假設(shè)f(t)的傅里葉變換為F(η),gω,b(t)的傅里葉變換為Gω,b(η),根據(jù)Parseval定理：

WFR(ω,b)=〈F(η),Gω,b(η)〉/2π。

(4)

若G(η)的有效窗口寬度為Dω，則WFR(ω,b)給出的是f(η)在[ω-Dω/2,ω+Dω/2]時間區(qū)間內(nèi)的頻譜信息。Dω越小，對信號的定位能力越強。

從公式可以看出，根據(jù)不確定原理，傅里葉變換的時域窗口和頻域窗口寬度不可能同時提高。在所有窗函數(shù)中，只有高斯函數(shù)Dt和Dω的乘積最小，因此為了盡可能減小信號通過濾波器失真的同時保證濾波器有效地進行濾波，選擇高斯窗函數(shù)對信號進行去噪是不錯的選擇。

2.2.2 布萊克曼窗函數(shù) 布萊克曼窗的時域形式可以表示為：

(5)

頻域特性為：

(6)

式中：t為離散時間；ω為頻率；WR(ω)為矩形窗函數(shù)的幅度頻率特性函數(shù)；N為窗函數(shù)長度。

布萊克曼窗函數(shù)的最大旁瓣值比主瓣值低57 dB，但主瓣寬度是矩形窗函數(shù)主瓣寬度的3倍，為12 π/N，雖其頻率識別精度低，但幅值識別精度高，有更好的選擇性。

2.3 窗函數(shù)設(shè)計FIR濾波器去噪

本文基于Matlab軟件，用2種窗函數(shù)設(shè)計FIR濾波器對波形進行去噪處理，具體過程如下。

1) 根據(jù)過渡帶寬和阻帶衰減要求，選擇窗函數(shù)的類型并估計窗口長度N。本文選擇高斯窗函數(shù)和布萊克曼窗函數(shù)。經(jīng)多次對比試驗，窗函數(shù)長度N取值為21。然后，調(diào)用Matlab中的窗函數(shù)求出中心值歸一化為1的窗函數(shù)序列w(n),n為離散時間。

2) 根據(jù)待求濾波器的理想頻率響應(yīng)求出理想單位脈沖響應(yīng)hd(n)。根據(jù)線性疊加原理，對于其理想單位脈沖響應(yīng)hd(n)，可以根據(jù)給定的窗函數(shù)長度N由Matlab語句得到。

3) 計算濾波器的單位脈沖響應(yīng)h(n)，即濾波器的系數(shù)向量。h(n)是理想單位脈沖響應(yīng)和窗函數(shù)的乘積,在Matlab中用點乘命令表示為h(n)=hd(n)×w(n)。

4) 將原始波形數(shù)據(jù)分別與各窗函數(shù)相對應(yīng)的濾波器的系數(shù)向量相互作用，進行波形去噪處理。

2.4 GLAS波形參數(shù)提取

本研究采用4組回波波形參數(shù)嘗試建立和改進森林冠頂高反演模型： 波形長度L、波形半能量高(HOME，height of median energy)H0、波形前緣長度l1和波形后緣長度l2。在噪聲估計的基礎(chǔ)上，利用噪聲均值與其4倍的標準差之和作為波形信號的始末閾值。對于波峰位置的確定，首先計算信號的一階導(dǎo)數(shù)，一階導(dǎo)數(shù)由正變負的零點位置即為波峰，以此檢測波形中的所有波峰。由于信號穿過冠層時，穿透能力減小，能量產(chǎn)生衰減，因此以信號最大強度的1/5作為檢測波峰的閾值，從信號起始位置向后搜索，第1個超過波峰閾值的波峰視為波形的第1個波峰位置，波形前緣長度l1即為第1個波峰位置與信號起始位置的差值； 從信號結(jié)束位置向前搜索，第1個超過波峰閾值的波峰即為地面波峰位置，波形長度L為地面波峰位置與信號起始位置的差值，波形后緣長度l2為地面波峰位置與信號結(jié)束位置的差值。波形半能量高是回波中植被部分能量一半的高度，波形的質(zhì)心位置即波形面積一半處的位置，該位置與地面回波位置的差值即為波形半能量高H0。各波形參數(shù)如圖3所示。

圖3 GLAS波形參數(shù)Fig.3 GLAS waveform parameters

2.5 GLAS森林冠頂高反演模型建立

在基于GLAS數(shù)據(jù)估算森林冠頂高的研究中，直接法和統(tǒng)計法較為常用(Iqbaletal.， 2013； Lefskyetal.， 2005； 2007； Rosetteetal.， 2008； Popescuetal.， 2011； Chen， 2010)。直接法利用GLAS回波信號開始位置和地面回波位置的垂直高度差來估算冠頂高，統(tǒng)計法則基于統(tǒng)計分析來估算冠頂高。在地形較為平坦的區(qū)域，基于GLAS完整波形，利用直接法提取波形長度即可得到較為準確的森林冠頂高，如圖4a中，冠頂高即為信號起始位置和地面回波位置的差值； 但當坡度增大時，回波信號受到影響，植被冠層部分會與地面回波信號混雜，且地面回波位置及最大冠頂位置會發(fā)生展寬，精確提取冠頂高較為困難。如圖4b，f為GLAS光斑大小，θ為地面坡度，ftanθ為地面回波信號向上延伸的高度，即坡度對信號影響的垂直距離。地面回波位置隨著坡度的增大其延伸也越大，僅僅依靠波形長度對冠頂高進行估算受地形影響較大，因此，采取一定方法減弱地形影響，對冠頂高的精確提取具有重大意義。

Chen(2010)研究表明，與直接法相比，統(tǒng)計法可以提高森林冠頂高估算精度，通過引入地形指數(shù)和GLAS波形參數(shù)可在一定程度上減弱地形的影響(Xingetal.， 2010； Lefskyetal.， 2007； Nieetal.， 2015)。前人研究中普遍采用地形起伏度(高程最大和最小值之差)(Iqbaletal.， 2013； Lefskyetal.， 2005； Hayashietal.， 2013； Chen，2010)和坡度對回波信號的垂直影響距離(湯旭光， 2013)作為地形指數(shù)。由表1可知，本研究區(qū)并非平坦區(qū)域，且各林型在不同坡度區(qū)間(0°～5°、5°～10°、10°～15°、15°～20°、20°～25°、25°～30°、>30°)均有分布，因此本研究引入地形指數(shù)，基于統(tǒng)計法建立森林冠頂高反演模型。

圖4 坡度對GLAS回波信號的影響Fig.4 Effects of slope on GLAS echo signal

本研究地形指數(shù)的具體選取流程為：首先提取GLAS激光光斑所處地理位置多種采樣窗口(3×3，5×5，7×7，9×9)下的地形起伏度并實測樣地坡度對回波信號的垂直影響距離，然后計算波形長度與地面實測最大冠頂高的差異，通過建立地形起伏度和坡度對回波信號影響的垂直距離與高度差異之間的相關(guān)性，從而評估哪一個地形指數(shù)最適合估算因地形坡度引起的差異，最后選擇其作為地形指數(shù)。本研究中地形起伏度從研究區(qū)10 m分辨率的DEM影像中提取，坡度對回波信號的垂直影響距離根據(jù)樣地實測坡度計算得到。2種地形指數(shù)與高度差異之間的相關(guān)性分析如表2所示。

表2高度差異與各地形指數(shù)的Pearson相關(guān)系數(shù)
Tab.2Pearsoncorrelationcoefficientsbetweenforestcrownheightdifferenceandterrainreliefofdifferentmovingwindows

變量Variable各采樣窗口下地形起伏度Terrainreliefofdifferentmovingwindows3×35×57×79×9坡度對回波信號的垂直影響距離Influencedverticaldistanceoftheslopeontheechosignal高度差異Heightdifference-0217-046-0192-0164-0011

由表2可知，5×5采樣窗口下地形起伏度與高度差異之間的相關(guān)性最大，因此本文以5×5采樣窗口下地形起伏度作為地形指數(shù)g參與森林冠頂高模型建立。

森林冠頂高模型參數(shù)選取對森林冠頂高反演精度有重要影響。相比于前人僅利用波形長度結(jié)合地形因子或僅利用波形參數(shù)建立的模型，本研究嘗試建立3類模型： 采用線性回歸方法，首先建立為以波形長度為參數(shù)的波形參數(shù)模型及以波形長度和地形指數(shù)為參數(shù)的地形因子模型，然后在地形因子模型基礎(chǔ)上，逐步引入波形半能量高、波形前緣長度、波形后緣長度等參數(shù)分林型建立全模型。引入波形半能量高，是利用其對冠層的垂直分布敏感而對坡度不敏感的特性作為冠頂高反演的重要參數(shù)(董立新， 2008)； 波形前緣長度可以反映植被冠層和地形復(fù)雜起伏對回波信號的綜合影響，引入該參數(shù)可有效抵消由于植被自身特征對回波造成的影響(Lefskyetal.， 2007)； 波形后緣長度可以反映地形坡度和粗糙度對回波信號的影響(Lefskyetal.， 2007)，引入該參數(shù)可在一定程度上修正由于DEM對地形表達不準確對冠頂高提取的影響。所建冠頂高反演模型及參數(shù)如表3所示。

表3 森林冠頂高反演模型及參數(shù)①Tab. 3 Forest crown height inversion models and parameters

①H為森林冠頂高估測值；L為波形長度；g為地形指數(shù)；H0為波形半能量高；l1為波形前緣長度；l2為波形后緣長度；其他均為系數(shù)。H: Estimated value of forest crown height;L: Waveform length;g: Terrain index;H0: Height of median energy;l1: Leading edge extent;l2: Trailing edge extent.

通過各林型實測數(shù)據(jù)的回歸分析確定不同模型的參數(shù)。

3 結(jié)果與分析

3.1 波形去噪結(jié)果對比與分析

本文選擇研究區(qū)內(nèi)47個脈沖點分別基于高斯窗函數(shù)和布萊克曼窗函數(shù)對GLAS波形進行去噪處理，其中1號點經(jīng)過去噪處理后的波形如圖5所示。

圖5 去噪前后波形對比Fig.5 Comparison of waveform before and after denoising

由圖5可知，1號點經(jīng)高斯低通濾波和小波去噪后，波形均得到了一定程度的平滑，保留了原始波形的大部分細節(jié)信息。為定量檢驗2種方法的去噪效果，本文選用均方根誤差(root mean square error,RMSE)和信噪比(signal to noise ratio,SNR)2個指標對去噪效果進行評價。計算公式如下：

；

(7)

(8)

式中：s(i)表示原始信號強度(V)；f(i)表示去噪后信號強度(V)；i表示信號位置；L表示信號長度。

根據(jù)式(7)和(8)，對所選47個脈沖點的布萊克曼窗去噪和高斯窗去噪效果分別進行定量評價，結(jié)果表明： 由高斯窗去噪得到的信號RMSE大部分都小于布萊克曼窗去噪后的RMSE，SNR大部分都大于布萊克曼窗去噪后的SNR。經(jīng)布萊克曼窗去噪處理后波形的信噪比平均值為20.958 81，均方根誤差0.000 285； 經(jīng)高斯窗去噪處理后波形的信噪比平均值為23.360 704，均方根誤差0.000 233。信噪比越高，均方根誤差越低，說明信噪分離效果越好。由此可得，高斯窗去噪方法信噪分離效果優(yōu)于布萊克曼窗去噪方法，能夠較好地保留原始信號的有用信息。

3.2 冠頂高模型精度評價與分析

本研究共利用47組實測冠頂高數(shù)據(jù)及地形指數(shù)和去噪前后的不同波形參數(shù)進行回歸分析，回歸分析在SPSS軟件中進行。各林型在不同模型下的去噪前后模型精度如表4所示。

采用R2、調(diào)整后的R2和RMSE 3個指標對波形去噪前后所建模型對各林型冠頂高擬合效果進行評價。由表4可知，去噪前所建模型的R2范圍為0.533 ～0.848，RMSE范圍為2.475 0～3.836 4 m；而去噪后各林型所建模型的R2顯著提高，范圍為0.686 ～0.972，RMSE普遍降低，范圍為1.004 1～2.756 2 m。這說明，對波形進行去噪，可以提取較為精確的波形參數(shù)，進而可以提高冠頂高估算精度。

由表4可知，去噪后5種模型均能較好地對復(fù)雜地形條件下的森林冠頂高進行估算。與波形參數(shù)模型相比，引入地形指數(shù)所建地形因子模型中，RMSE有所降低。針葉林和針闊混交林模型模擬效果最佳，R2高達0.853，闊葉林相對較差，R2為0.697 。

在引入波形半能量高、波形前緣長度和波形后緣長度等波形參數(shù)后，所建模型改善了各林型的地形因子模型。董立新(2008)研究發(fā)現(xiàn)，在坡度較大地區(qū)，引入波形半能量高對坡度的敏感性較低，引入該參數(shù)對森林冠頂高的提取結(jié)果比較穩(wěn)定。由于本研究區(qū)各林型在各坡度區(qū)間均有分布，因此本文引入波形半能量高建立模型3，減弱坡度較大地區(qū)對冠頂高估算的影響。模型3在各林型中總體較為穩(wěn)定，該模型在地形指數(shù)校正模型的基礎(chǔ)上引入波形半能量高，R2范圍為0.721～0.915，RMSE范圍為2.031 0～2.554 9 m； 模型4在引入波形前緣長度后，冠頂高估算精度在不同林型中總體最高，但調(diào)整后的R2普遍降低。Lefsky等(2007)研究指出，波形后緣長度與地形最為相關(guān)，波形前緣長度與冠頂高變化有關(guān)，適用于估計平均樹高而不是最大樹高，這可以解釋本文在估算最大冠頂高時引入波形前緣長度模型調(diào)整后的R2普遍降低的問題。進一步引入波形后緣長度建立模型5，可以發(fā)現(xiàn)，模型的R2有了進一步提高，范圍為0.728 ～0.972，RMSE范圍為1.001 4～2.216 4 m； 效果最好的仍為針葉林和針闊混交林，其中針闊混交林的R2提高到0.972，RMSE為1.001 4 m，調(diào)整后的R2為0.902 。對于各森林類型，全模型的RMSE范圍在1.001 4～2.554 9 m之間，地形因子模型的RMSE范圍在2.337 3～2.689 8 m之間。其中地形因子模型中針葉林效果最好：R2=0.853， RMSE=2.519 7 m； 全模型中針闊混交林效果最好：R2=0.972，RMSE=1.001 4 m。當不分森林類型時，全模型R2范圍為0.761～0.790，RMSE范圍為2.559 7～2.734 9 m，總體效果較好，這說明當缺乏森林類型數(shù)據(jù)時，利用統(tǒng)一的森林冠頂高，基于GLAS數(shù)據(jù)建立冠頂高反演模型也能準確地估算森林冠頂高。

表4 冠頂高反演模型精度比較Tab. 4 Accuracy comparison of forest crown height inversion models

由以上分析可知，基于去噪后的GLAS波形數(shù)據(jù)估算森林冠頂高時，全模型總體優(yōu)于波形參數(shù)模型和地形因子模型，針葉林優(yōu)于針闊混交林，闊葉林最差。

4 討論

基于GLAS波形數(shù)據(jù)提取森林冠頂高時，對GLAS波形的去噪效果在一定程度上影響波形參數(shù)的準確提取，而建立冠頂高反演模型所選取的模型參數(shù)可直接影響到所建模型對冠頂高的模擬效果，因此選擇合適的波形去噪算法和模型參數(shù)至關(guān)重要。

王愛娟等(2016)探討了漢寧窗、布萊克曼窗、矩形窗、三角窗和凱撒窗5種窗函數(shù)對林區(qū)GLAS波形數(shù)據(jù)去噪的效果，得出布萊克曼窗去噪效果最好。本文考慮GLAS激光發(fā)射器發(fā)射的信號為高斯脈沖信號，可假設(shè)激光回波由高斯分量構(gòu)成(邢艷秋等， 2009)，因此采用高斯窗函數(shù)對波形進行去噪處理，結(jié)果發(fā)現(xiàn)高斯窗去噪效果優(yōu)于布萊克曼窗，且從基于高斯窗去噪后的波形對各林型建立冠頂高反演模型的模擬結(jié)果來看，所建模型R2最大為0.972，高于王愛娟等(2016)所建模型R2最大為0.91的研究結(jié)果。這說明相比其他窗函數(shù)，高斯窗函數(shù)在由高斯分量構(gòu)成的GLAS激光回波去噪處理中更有優(yōu)勢。在基于去噪后的GLAS波形直接估算森林冠頂高時，Iqbal 等(2013)通過對GLAS波形進行去噪處理，利用直接法，可以解釋79%的冠頂高變化，RMSE為3.18 m，冠頂高估算效果較好。本文在基于去噪后的GLAS波形建立模型估算森林冠頂高時，去噪后的各模型精度顯著提高。以上研究表明，在基于GLAS數(shù)據(jù)估算冠頂高的過程中，對波形去噪方法的優(yōu)選很有必要。

在以地形指數(shù)和波形長度為因變量建立的冠頂高估算模型研究中，Lefsky 等(2005)所建估算模型能解釋59% ～68%的森林最大冠頂高變化，RMSE范圍為4.85～12.66 m； Rosette等(2008)所建冠頂高估算模型R2為0.89，RMSE=2.99 m； 僅利用波形自身參數(shù)所建冠頂高估算模型中，Lefsky等(2007)在7個研究區(qū)基礎(chǔ)上利用3 個波形指數(shù)(波形寬度、下端邊緣糾正系數(shù)和上端邊緣糾正系數(shù))所建模型能解釋83%的森林冠頂高變化，RMSE為5 m。相比于湯旭光(2013)利用地形指數(shù)和波形長度以及波形半能量高分林型建立多元線性回歸模型，各林型調(diào)整后R2最大為0.904，RMSE最小為2.021 m，本研究在此基礎(chǔ)上還引入了波形前緣長度和波形后緣長度，各林型所建模型R2最大為0.972，RMSE最小為1.001 4 m，模型的模擬效果相比前人研究有所改善。各林型中，針葉林和針闊混交林冠頂高模擬效果最好，這與湯旭光(2013)的結(jié)論一致。盡管在缺少林型數(shù)據(jù)時，本研究所建模型對冠頂高的反演能力相對較差，R2最大為0.79，RMSE最小為2.559 7 m，但總體來說，也能獲取相對準確的森林冠頂高數(shù)據(jù)。以上分析表明，分林型引入多種波形參數(shù)建立冠頂高反演模型在地形復(fù)雜地區(qū)具有很強的適用性。本研究利用地形指數(shù)和波形參數(shù)雖在一定程度上修正了坡度對冠頂高提取的影響，但研究中影響冠頂高反演模型對森林冠頂高模擬精度的因素還包括時間差異、樹種差異、樣地地表狀況等不確定性因素，因此在今后的研究中，基于地形數(shù)據(jù)和波形參數(shù)聯(lián)合光學(xué)遙感數(shù)據(jù)，修正影響森林冠頂高反演精度的多種因素，對進一步提高森林垂直結(jié)構(gòu)參數(shù)的提取精度有重大意義。

由于GLAS光斑呈條帶形分布，本研究中樣地布設(shè)原則是與GLAS光斑的位置大小相匹配，來估算GLAS光斑內(nèi)森林冠頂高并顧及樣地的可達性，所以本研究樣地也呈條帶形分布，這在一定程度上會影響模型的推廣，為此，本研究在樣地設(shè)置時，已綜合考慮了研究區(qū)主要林型及地形情況，樣地盡可能涵蓋研究區(qū)森林類型及不同地形條件。為了進一步增強冠頂高反演模型的適用性，今后將通過補充野外調(diào)查，獲取更多的樣地數(shù)據(jù)，進一步優(yōu)化模型，并結(jié)合光學(xué)遙感數(shù)據(jù)，綜合考慮森林冠頂高空間擴展的可能有效因子，確定光斑內(nèi)冠頂高空間外推的最佳尺度和模型，進行森林冠頂高的空間擴展研究。

5 結(jié)論

本文以吉林省露水河林區(qū)為研究區(qū)，基于大光斑激光雷達數(shù)據(jù)，采用布萊克曼窗函數(shù)和高斯窗函數(shù)對波形進行去噪，并定量比較2種方法的去噪效果。分別基于去噪后的GLAS波形，引入地形指數(shù)和波形參數(shù)建立不同森林類型冠頂高反演模型，通過對比2類冠頂高反演模型比較復(fù)雜地形中各反演參數(shù)對不同林型的冠頂高反演能力。結(jié)果表明：

1) 在對GLAS數(shù)據(jù)進行去噪處理中，與布萊克曼窗函數(shù)去噪相比，高斯窗函數(shù)去噪得到的RMSE較低、SNR較高，信噪分離效果較好，適應(yīng)性較強?；贕LAS波形參數(shù)建立冠頂高反演模型時，對波形去噪后，所建模型顯著提高了模型對森林冠頂高的模擬效果。因此在基于GLAS數(shù)據(jù)提取冠頂高時，對波形進行去噪處理很有必要。

2) 在地形復(fù)雜地區(qū)，基于GLAS數(shù)據(jù)分林型對冠頂高的反演能力比不分林型的反演能力強，對針葉林和針闊混交林的反演效果較好，所建模型的R2最高為0.972，RMSE最小為1.001 4 m； 闊葉林反演效果較差，這可能是由于樣地內(nèi)最大冠頂高林木為落葉闊葉樹種，GLAS獲取時間為10月底，由于林木凋零，導(dǎo)致冠頂回波信號較弱而難以識別，從而導(dǎo)致估測值偏低。

3) 在所建森林冠頂高反演模型中，各林型的冠頂高估算效果有所不同。其中地形因子模型中針葉林效果最好：R2=0.853， RMSE=2.519 7 m； 全模型中針闊混交林效果最好：R2=0.972， RMSE=1.001 4 m?？傮w來說，各模型中，針葉林和針闊混交林冠頂高模擬效果最好。

4) GLAS數(shù)據(jù)反演森林冠頂高時，引入多種波形參數(shù)顯著改善了由地形指數(shù)建立的校正模型，其中效果最好的反演參數(shù)為波形長度、地形指數(shù)、波形半能量高和波形后緣長度，說明結(jié)合波形自身參數(shù)建立模型在森林冠頂高反演方面具有很大的應(yīng)用潛力。

董立新. 2008.基于多源遙感數(shù)據(jù)的三峽庫區(qū)森林冠層高度與生物量估算方法研究.北京： 中國科學(xué)院研究生院博士學(xué)位論文.

(Dong L X. 2008.Study on canopy height and biomass estimation method of Three Gorges Reservoir Area based on multi-source remote sensing data. Beijing： PhD thesis of Graduate University of Chinese Academy of Sciences. [in Chinese])

李立存, 張淑芬, 邢艷秋. 2011.全站儀和測高儀在樹高測定上的比較分析. 森林工程, 27(4): 38-41.

(Li L C, Zhang S F, Xing Y Q. 2011.Comparative analysis on tree height measured by total station and vertex IV altimeter. Forest Engineering, 27(4): 38-41. [in Chinese])

劉經(jīng)南,張小紅. 2005. 利用激光強度信息分類激光掃描測高數(shù)據(jù). 武漢大學(xué)學(xué)報：信息科學(xué)版, 30(3): 189-193.

(Liu J N,Zhang X H. 2005.Classification of laser scanning altimetry data using laser intensity. Geomatics and Information Science of Wuhan University, 30(3): 189-193. [in Chinese])

婁雪婷, 曾 源, 吳炳方. 2011.森林地上生物量遙感估測研究進展.國土資源遙感,23(1):1-8.

(Lou X T, Zeng Y, Wu B F. Advances in the estimation of above-ground biomass of forest using remote sensing.Remote Sensing For Land & Resources,23(1): 1-8. [in Chinese])

龐 勇, 孫國清, 李增元.2006.林木空間格局對大光斑激光雷達波形的影響模擬. 遙感學(xué)報, 10(1): 97-103.

(Pang Y, Sun G Q, Li Z Y. 2006.Large footprint lidar waveform modelling of forest spatial patterns. Journal of Remote Sensing, 10(1): 97-103. [in Chinese])

邱 賽, 邢艷秋, 李立存, 等. 2012.基于小波變換的ICESat-GLAS波形處理. 森林工程, 28(5): 33-35,59.

(Qiu S, Xing Y Q, Li L C,etal. 2012. ICESat-GLAS data processing based on wavelet transform. Forest Engineering, 28(5): 33-35,59. [in Chinese])

孫國清, Ranson K J, 張鐘軍. 2006.利用激光雷達和多角度頻譜成像儀數(shù)據(jù)估測森林垂直參數(shù). 遙感學(xué)報, 10(4): 523-530.

(Sun G Q, Ranson K J, Zhang Z J. 2006.Forest vertical parameters from lidar and multi-angle imaging spectrometer data. Journal of Remote Sensing, 10(4): 523-530. [in Chinese])

湯旭光. 2013.基于激光雷達與多光譜遙感數(shù)據(jù)的森林地上生物量反演研究. 北京: 中國科學(xué)院大學(xué)博士學(xué)位論文.

(Tang X G.2013.Inversion of forest ground biomass based on LiDAR and multispectral remote sensing data. Beijing: PhD thesis of University of Chinese Academy of Sciences. [in Chinese])

王愛娟, 邢艷秋, 邱 賽, 等. 2016.基于窗函數(shù)的林區(qū)ICESAT-GLAS波形數(shù)據(jù)消噪研究. 西北林學(xué)院學(xué)報, 31(1): 214-220.

(Wang A J, Xing Y Q,Qiu S,etal. 2016.Denoising of forest ICESat-GLAS waveform data based on window function. Journal of Northwest Forestry University, 31(1): 214-220. [in Chinese])

邢艷秋, 王立海. 2008.基于森林生物量相容性模型長白山天然林生物量估測. 森林工程, 24(2): 1-4.

(Xing Y Q, Wang L H. 2008.Biomass estimation of natural forest in Changbai Mountains based on forest biomass compatibility model. Forest Engineering, 24(2): 1-4. [in Chinese])

邢艷秋, 王立海. 2009.基于ICESat-GLAS完整波形的坡地森林冠層高度反演研究——以吉林長白山林區(qū)為例. 武漢大學(xué)學(xué)報：信息科學(xué)版, 34(6): 696-700.

(Xing Y Q, Wang L H.2009. ICESat-GLAS full waveform-based study on forest canopy height retrieval in sloped area—a case study of forests in Changbai Mountains, Jilin. Geomatics and Information Science of Wuhan University, 34(6): 696-700. [in Chinese])

Chen Q. 2010. Retrieving vegetation height of forests and woodlands over mountainous areas in the pacific coast region using satellite laser altimetry. Remote Sensing of Environment, 114(7): 1610-1627.

Hayashi M, Saigusa N, Oguma H,etal. 2013.Forest canopy height estimation using ICESat/GLAS data and error factor analysis in Hokkaido, Japan. ISPRS Journal of Photogrammetry and Remote Sensing, 81(7): 12-18.

Iqbal I A, Dash J, Ullah S,etal. 2013.A novel approach to estimate canopy height using ICESat/GLAS data: a case study in the New Forest National Park, UK. International Journal of Applied Earth Observation and Geoinformation, 23(1): 109-118.

Lefsky M A, Harding D J, Keller M,etal. 2005.Estimates of forest canopy height and aboveground biomass using ICESat. Geophysical Research Letters, 32(22): 441.

Lefsky M A, Keller M, Pang Y. 2007.Revised method for forest canopy height estimation from geoscience laser altimeter system waveforms. Journal of Applied Remote Sensing, 1(1): 1-18.

Nie S, Wang C, Zeng H,etal. 2015.A revised terrain correction method for forest canopy height estimation using ICESat/GLAS data. ISPRS Journal of Photogrammetry and Remote Sensing, 108: 183-190.

Pang Y, Li Z Y, Lefsky M,etal.2006. Model based terrain effect analyses on ICESAT GLAS waveforms. IEEE International Conference on Geoscience and Remote Sensing Symposium, 3232-3235.

Popescu S C, Zhao K, Neuenschwander A,etal.2011. Satellite lidar vs.small footprint airborne LiDAR:comparing the accuracy of aboveground biomass estimates and forest structure metrics at footprint level. Remote Sensing of Environment,115(11):2786-2797.

Rosette J A B, North P R J, Suarez J C. 2008.Vegetation height estimates for a mixed temperate forest using satellite laser altimetry.International Journal of Remote Sensing,29(5):1475-1493.

Xing Y, de Gier A, Zhang J,etal.2010. An improved method for estimating forest canopy height using ICESat/GLAS full waveform data over sloping terrain: a case study in Changbai Mountains, China. International Journal of Applied Earth Observation and Geoinformation, 12(5): 385-392.