基于星載激光測高儀的森林植被高度估算

褚 喆

(河南測繪職業(yè)學院測繪工程系,河南 鄭州 451460)

1 引 言

森林是以樹木等植被為主體的生物群落,其作為“地球之肺”,越來越受到人類的關注。尤其是在近幾十年來全球氣候變暖的大環(huán)境下,人類越發(fā)重視森林的覆蓋率、森林碳儲量,而森林內部植被高度直接決定著森林的碳儲量[1],因此,針對森林植被高度的調查研究不斷被提出。由于光學影像無法透視植被的垂直結構,因此常用的研究方法有兩種,一種是人工對森林樣地進行實地測量[2],該方法精度高,但是測量范圍十分有限。另一種是利用機載激光雷達進行航飛測量[3],雖然這種方法相對便捷且具備較高精度,但該方法花費巨大,對于全球或者國內的森林地區(qū)仍然無法實現(xiàn)大面積的普查。

星載激光測高儀憑借其精度高、穿透能力強、時效性短、觀測范圍廣等特點[4],廣泛用于了全球極地冰蓋監(jiān)測、陸地高程測量、以及森林植被高度估算等(如2003年美國發(fā)射的GLAS(Geoscience Laser Altimeter System)激光測高系統(tǒng)[5])。尤其是近幾年來隨著我國星載激光測高儀的快速發(fā)展,如2016年我國發(fā)射的資源三號02星(ZY3-02)激光測高儀[6],2019年我國發(fā)射的高分七號(GF-7)星載激光測高系統(tǒng)[7-8],以及我國即將發(fā)射的專門用于森林普查的陸地生態(tài)碳衛(wèi)星激光測高系統(tǒng),這使得利用星載激光測高儀進行國內以及全球森林植被高度估算成為了新的手段。利用星載激光測高儀研究森林植被高度,主要是對經植被返回的激光波形進行研究,估算出激光光斑內部植被的高度。目前針對激光雷達回波波形,已有一些學者開展了相關的研究。Yuchu Qin等人提出了一種通過找到激光雷達回波峰值與波形拐點特殊關系的波形分解方法,該方法針對的是小光斑的激光雷達[9]。唐福鑫等人[10]和李芬芬等人[11]采用小波變化方法對激光雷達回波波形進行分解,并提取出波形有效長度。趙欣等人采用常用的高斯分解方法對波形進行分解,利用地面激光器估算出短距離的樹木高度[12]。崔成玲提出了一種改進的基于波峰自動識別的全波形高斯分解方法,進行植被高度估算[13]。劉峰等人采用廣義高斯函數模型對激光波形進行擬合,提取出了激光點的振幅、距離等參數[14]。Michael 等人利用GLAS數據與SRTM(Shuttle Radar Topograpphy Mission)數據聯(lián)合估算了植被高度與森林蓄積量[15]。以上研究均取得了一定的成果,并對本文有著重要啟發(fā),但這些研究大多注重于對簡單波形的分解,并計算光斑內地表高差；或是采用地面激光器進行測量實驗,驗證方法精度,其估算植被高度的精度還需進一步提高,同時對于實際在軌的星載激光測高儀進行森林植被高度估算仍需進一步研究。因此,本文在已有研究基礎上,提出了基于一種基于高斯混合函數模型對星載激光測高儀森林植被波形進行分解,從而估算森林植被高度的方法。實驗以GLAS激光測高系統(tǒng)為例,先后開展自適應波形背景噪聲濾除、高斯波形平滑濾波波形預處理工作,對平滑后的波形采用基于高斯混合函數模型的波形分解方法實現(xiàn)對有效波形的提取。最后通過設置一定閾值找到森林地區(qū)回波波形中首末高斯峰值,從而估算出該激光光斑內植被的最大高度。

2 研究理論與方法

2.1 波形數據預處理

2.1.1 自適應波形背景噪聲濾除

受大氣、云層、太陽背景噪聲,以及激光雷達自身的因素,星載激光回波波形中存在大量噪聲。背景噪聲的剔除有利于判定有效信號的起始位置,對于提高后續(xù)的波形特征參數提取精度。然而,激光落在不同地物表面時,波形不一致,噪聲量也不同。因此,針對森林地球回波波形,本文提出一種自適應波形背景噪聲濾除方法。首先選取1×10的一維窗口,從波形起點逐步向后檢索,當窗口內波形能量強度開始呈現(xiàn)遞增趨勢,停止檢索,獲取波形起點至遞增前一個位置處所有波形前端數據Nb。隨后,以相同方式從波形末端往前端搜索獲取波形后端數據Na。最后,將波形前端數據與波形后端數據合并,認為時該波形噪聲數據ν,ν=Nb+Na。則單個波形背景噪聲閾值Noise的計算公式為:

(1)

2.1.2 高斯平滑濾波

由于星載激光雷達回波波形由各個采樣點連接形成,前后采樣點之間可能存在較大采樣間隙,同時回波波形中存在較多的小起伏現(xiàn)象(“波形毛刺”)。波形平滑后可消除該類波形中小“毛刺”,從而減小后續(xù)波形分解難度、速度,并提高其精度。平滑后的波形更加容易反映地表垂直信息,針對這一現(xiàn)象本文采用高斯平滑濾波的方式對回波波形進行平滑處理。高斯平滑濾波公式如下所示。

(2)

式中,a為平滑的曲線幅值;b為平滑的曲線在x軸的中心;c指高斯曲線半峰全寬。在對激光回波波形進行高斯平滑處理時,可采用一定窗口進行遍歷處理,窗口大小1×10。

2.2 基于高斯混合函數模型的波形分解

星載激光在復雜地形回波波形呈現(xiàn)多波峰現(xiàn)象,特別是森林地區(qū),星載激光回波波形更為復雜。傳統(tǒng)單一高斯函數模型通常用于描述單個波峰的波形數據,而當波形數據出現(xiàn)多峰值、多拐點、不規(guī)則現(xiàn)象時,該函數模型無法表達該波形。故,本文針對森林地星載激光多峰值回波波形,提出采用高斯混合函數模型(Gaussian mixture model,GMM)進行波形分解。

高斯混合函數模型采用聚類算法思想,利用多個單一高斯函數模型作為一個離散值數值的參數模型,并根據期望最大算法進行訓練,從而獲得復雜的森林地區(qū)回波波形。高斯混合函數模型的概率密度函數為:

(3)

(4)

式中,u為波形采樣數據均值；σ2為波形采樣數據標準差。

實際操作中,由于不知道xi屬于哪個高斯函數以及每個高斯函數生成器在混合模型中所占的比例是未知的?，F(xiàn)采用最大似然法估計波形樣本的高斯混合函數概率密度函數的概率最大值,從而選擇最佳的函數模型。為防止在計算過程中產生溢出現(xiàn)象,可將目標函數取對數,得到最大化對數似然函數,如下式所示:

(5)

基于高斯混合函數模型進行復雜森林地區(qū)波形分解的基本流程如下:

(1)計算星載激光測高儀回波波形數據中單高斯函數個數,以及每個高斯函數對整個波形生成的概率；

(2)估計星載激光測高儀回波波形數據中每個高斯函數參數uk與σk的值；

(3)重復步驟(1)與(2),直到似然函數收斂,停止計算,即可得到k初始參數(xk|uk,σk)；

(4)將步驟(3)得到的初始參數代入公式(6),可擬合波形,即可完成回波波形分解。

(6)

式中,Ak為第k可高斯函數的峰值,可直接從回波波形中提??；ε為常量,通常為噪聲。

2.3 植被高度估算

通常情況下認為波形分解得到的第一個波形為植被冠層回波,最后一個波形為地面回波,最后一個波峰與第一個波峰時間差轉換為距離,即為激光光斑內植被高度。但實際情況下,第一個回波波形幅值可能非常小,其主要由個別樹木冠層作用生成,利用該波形計算會導致光斑內植被高度偏大。本文采樣波峰能量平均閾值法篩選出有效波形,首先計算波形分解得到的所有波峰均值,將其的一定倍數設置為波峰閾值,剔除小于閾值的波形。波峰閾值公式如下所示:

(7)

式中,Ai為第i個高斯波形峰值；k為回波波形分解得到高斯波形個數；λ為比例因子,森林地區(qū)可設置為0.1～0.5。

基于上述波峰閾值篩選得到高斯波形后,根據時間順序上首末高斯波形即可計算出該光斑內植被高度htree,如下式所示:

htree=(tend-tbegin)×BinRange

(8)

式中,tend為篩選后的最后一個波峰對應幀時刻；tbegin為篩選后的第一個波峰對應幀時刻；BinRange波形每幀數據對應測距值,對于GLAS為0.15m,對于我國GF-7激光為0.075 m。

3 試驗區(qū)域與數據

3.1 試驗區(qū)域

本文試驗區(qū)域選取我國東北長白山林區(qū),位于吉林省汪清縣境內,選取的試驗區(qū)經緯度為43°6′N～43°42′N,130°18′～130°49′E,如圖1所示。該地區(qū)處于寒溫帶森林生態(tài)系統(tǒng),地形為丘陵。受氣候因素影響,該地區(qū)植被為針葉林和落葉林混交分布,以針葉林為主。其中主要植株為紅松、云杉等針葉林木,以及段樹、楓華闊葉林木。

3.2 試驗數據

3.2.1 林高實測數據

本文選用文獻[13]已經公開的林高數據作為實驗驗證數據,該數據由東北林業(yè)大學邢艷秋教授團隊在2006年進行現(xiàn)場實測獲取。試驗區(qū)內共有20塊實測樣地林高數據,如圖1(b)同心圓圖標(實測點)所示。這20個點基本信息如表1所示。

圖1 長白山林地試驗區(qū)域圖

表1 實測數據及其鄰近GLAS激光點坐標

3.2.2 GLAS激光點數據

依據區(qū)域范圍,檢索2003—2009年間所有GLAS激光數據(GLAS數據下載地址為:https://search.earthdata.nasa.gov/search?q=GLAH),試驗區(qū)內經過實測樣地的GLAS 01級數據如表2所示(GLAS激光點如如圖1(a)所示)。其中上述20塊樣地最鄰近的GLAS點詳細經緯度如表 1所示,這些GLAS數據點與實測樣地分布如圖1(b)所示。

表2 試驗區(qū)內GLAS 01級數據包信息表

4 試驗結果與分析

4.1 波形分解結果

根據本文所述對選取的20個GLAS激光點先后進行波形背景噪聲濾除、高斯平滑濾波以及基于高斯混合模型的波形分解。由于這20個激光點大致分布于五個區(qū)域,圖2展示出每個區(qū)域的一個GLAS激光點處理結果。

圖2～圖6給出了基于本文方法對GLAS波形數據進行波形分解的實驗過程圖,每幅圖的左圖為對原始波形數據進行背景噪聲濾除結果,中間圖為對去噪后的波形進行高斯平滑后結果,右圖為平滑后波形進行基于高斯混合函數模型的波形分解結果。

圖2 激光點120954008.76波形分解過程圖

圖3 激光點184546382.17波形分解過程圖

圖4 激光點184936425.75波形分解過程圖

圖5 激光點216385708.81波形分解過程圖

圖6 激光點216775745.77波形分解過程圖

上述結果可看出,對于去噪后的波形,通過高斯平滑濾波方法能夠很好消除波形中細微的波形“毛刺”,且能保留下波形中主要的峰值,使得整個波形更加平滑,易于后續(xù)開展波形分解。根據上述本文基于高斯混合函數模型的波形分解結果可看出,該方法與已有方法不同,其得到的子波形個數不局限于6個[13],如圖2和圖3所示,本文方法分解后得到9個子波形,若按照已有方法僅保留6個波形,舍去峰值最小的3個子波形,會導致根據首末波峰提取時間差變小,直接導致根據波形分解得到的植被高度變小,從而降低星載激光波形提取植被高度精度。而且,通過設置波峰閾值可剔除波峰能量非常小的子波形,該波形一般由植被局部少量枝葉遮擋形成,該波形信息將會影響基于波形提取的植被高度。綜上,經本文方法能夠有效的將復雜波形分解開來,且便于后續(xù)進行植被高度計算。

4.2 植被高度估算

經本文方法對GLAS原始波形進行波形分解后,利用分解出來的最后一個子波形峰值對應時刻減去第一個子波形峰值對應時刻,作為激光在植被冠層至地面間的渡越時間。激光在空氣中傳播速度近似為0.3 m/ns,但由于激光從植被冠層傳輸至地面后還需返回,故計算植被高度時,將該渡越時間乘以0.15 m/ns作為是植被的高度。表1中GLAS數據估算的植被高度,以及樣地實測樹高與估算樹高差值如表3所示。

表3 樣地實測樹高與激光估算的樹高對比表

由表3結果可以看出,整體波形估算的樹高比實測樹高要低,其原因是,實測數據可直接取到冠層頂部高度,而波形中冠層頂部少量的葉片和頂端無法形成獨立有效峰值,從而導致激光波形估算的植被高度要比實測的低。且前10號點激光估算的較實測樹高相差更大,分析其原因是前10號點為GLAS在2003年和2005年測得的,實測數據為2006年,樹木存在一定生長。后10號點雖是GLAS于2006年獲取的,但其與實測樹高仍存在一定差距,分析認為除了冠層和樹頂無法形成有效波峰外,本文尚未考慮光斑內部地表坡度,坡度會影響提取的植被高度,當坡度較大時,最后波形分解得到最后波形時間延長,從而導致估算植被高度增加(如表3中13、15、20號點估算樹高大于實測樹高)。但是,從表3仍然可以看出,整體上本文方法提取的植被高度與實際樣地測量的植被高度相差不大,此20個點激光波形計算的植被與實測樹高偏差為(0.3±1.4)m,相對于文獻[13]精度(0.5±2.9)m有明顯提升。綜上所述,說明本文基于高斯混合函數模型的波形分解方法精度較高,且基于該方法用于提取植被高度是有效可行的。

5 結 論

星載激光測高儀主要利用其回波波形對森林植被垂直高度進行估算。本文針對森林地區(qū)的星載激光回波波形,提出了一種基于高斯混合函數模型的波形分解方法,對預處理后的激光波形進行分解,破除了已有研究中波形分解最多分解出6個子波形的限制,從而提升了波形分解估算森林植被高度的精度。

本文以長白山林區(qū)為實驗區(qū)域,選取了20個GLAS激光點作為實驗對象,根據本文方法估算出各個光斑內植被最大高度,并利用該光斑內人工樣地實測的森林植被高度進行驗證。結果表明本文方法估算的植被高度精度為0.3±1.4m,相對已有的方法精度更高。從而證明該方法可用于星載激光測高儀進行森林植被高度估算,為后續(xù)我國高分七號、陸地生態(tài)碳衛(wèi)星等國產星載激光測高儀進行森林植被高度估算提供可行的技術方法。