褚 喆

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

1 引 言

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

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

2 研究理論與方法

2.1 波形數(shù)據(jù)預(yù)處理

2.1.1 自適應(yīng)波形背景噪聲濾除

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

(1)

2.1.2 高斯平滑濾波

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

(2)

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

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

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

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

(3)

(4)

式中,u為波形采樣數(shù)據(jù)均值；σ2為波形采樣數(shù)據(jù)標(biāo)準(zhǔn)差。

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

(5)

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

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

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

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

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

(6)

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

2.3 植被高度估算

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

(7)

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

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

htree=(tend-tbegin)×BinRange

(8)

式中,tend為篩選后的最后一個(gè)波峰對(duì)應(yīng)幀時(shí)刻；tbegin為篩選后的第一個(gè)波峰對(duì)應(yīng)幀時(shí)刻；BinRange波形每幀數(shù)據(jù)對(duì)應(yīng)測(cè)距值,對(duì)于GLAS為0.15m,對(duì)于我國GF-7激光為0.075 m。

3 試驗(yàn)區(qū)域與數(shù)據(jù)

3.1 試驗(yàn)區(qū)域

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

3.2 試驗(yàn)數(shù)據(jù)

3.2.1 林高實(shí)測(cè)數(shù)據(jù)

本文選用文獻(xiàn)[13]已經(jīng)公開的林高數(shù)據(jù)作為實(shí)驗(yàn)驗(yàn)證數(shù)據(jù),該數(shù)據(jù)由東北林業(yè)大學(xué)邢艷秋教授團(tuán)隊(duì)在2006年進(jìn)行現(xiàn)場(chǎng)實(shí)測(cè)獲取。試驗(yàn)區(qū)內(nèi)共有20塊實(shí)測(cè)樣地林高數(shù)據(jù),如圖1(b)同心圓圖標(biāo)(實(shí)測(cè)點(diǎn))所示。這20個(gè)點(diǎn)基本信息如表1所示。

圖1 長(zhǎng)白山林地試驗(yàn)區(qū)域圖

表1 實(shí)測(cè)數(shù)據(jù)及其鄰近GLAS激光點(diǎn)坐標(biāo)

3.2.2 GLAS激光點(diǎn)數(shù)據(jù)

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

表2 試驗(yàn)區(qū)內(nèi)GLAS 01級(jí)數(shù)據(jù)包信息表

4 試驗(yàn)結(jié)果與分析

4.1 波形分解結(jié)果

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

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

圖2 激光點(diǎn)120954008.76波形分解過程圖

圖3 激光點(diǎn)184546382.17波形分解過程圖

圖4 激光點(diǎn)184936425.75波形分解過程圖

圖5 激光點(diǎn)216385708.81波形分解過程圖

圖6 激光點(diǎn)216775745.77波形分解過程圖

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

4.2 植被高度估算

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

表3 樣地實(shí)測(cè)樹高與激光估算的樹高對(duì)比表

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

5 結(jié) 論

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

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