譚建偉，程春泉，王志勇，徐志達

(1.山東科技大學 測繪科學與工程學院，山東 青島 266590；2.中國測繪科學研究院，北京 100036)

0 引言

冰、云和陸地高程測量衛(wèi)星(ice，cloud and land elevation satellite，ICESat)搭載地球科學激光測高系統(tǒng)(the geoscience laser altimeter system，GLAS)，是全球首個具備全波形記錄和采樣功能的大光斑激光雷達測高衛(wèi)星[1]。其主要科學任務是監(jiān)測極地冰蓋變化情況，另外也用于準確估計全球地區(qū)的森林結構特征，如林區(qū)冠高、生物量和碳含量等[2]。美國航空航天局(National Aeronautics and Space Administration，NASA)于2018年發(fā)射的ICESat-2衛(wèi)星進一步監(jiān)測海冰厚度變化，并收集有關森林生長和地形高度的數(shù)據(jù)[3]。此外，作為我國高分辨率對地觀測系統(tǒng)重大專項之一的高分七號已經(jīng)發(fā)射升空。該星具備優(yōu)于1 m分辨率的立體觀測能力以及精度優(yōu)于1 m的激光測高能力，主要用于1∶50 000、1∶10 000比例尺地圖測圖和更新，以及海冰和陸地植被信息等資源的調(diào)查工作。為更好了解并利用衛(wèi)星測高數(shù)據(jù)產(chǎn)品進行極地冰蓋變化監(jiān)測、陸地資源調(diào)查以及提升自主處理星載激光雷達數(shù)據(jù)的能力，對星載全波形數(shù)據(jù)進行數(shù)據(jù)處理研究工作顯得非常必要。

由于激光雷達具有穿透植被的特性，回波樣本中記錄的是多個目標反射波形相互疊加的結果[4]。為更好獲取地物特征參數(shù)，需對回波波形進行分解，中外學者對波形分解的方法進行了大量而深入的探討。Hofton等[5]提出經(jīng)典的高斯分解算法，通過相鄰奇偶拐點確定高斯分量，該方法簡單有效，但受噪聲影響較大，會產(chǎn)生許多“偽”拐點。王成等[6]采用二階偏導數(shù)求取拐點來解算高斯分量的中心位置、振幅和半寬，并對最大振幅處的高斯分量位置和半寬予以調(diào)整，用來削弱由于飽和前向散射使波形曲線扭曲的影響，最后采用非線性最小二乘來擬合優(yōu)化，其分解精度的高低依賴于高斯分量個數(shù)初始估計的準確度。賴旭東等[7]采用逐層剝離的策略不斷分解出波形分量，直到最大峰值小于設定的閾值，該方法分解速度快，對弱波的探測能力強，但其利用1/2峰值位置確定半寬值的方法，導致在連續(xù)波形震蕩處往往得不到正確結果。李國元等[8]提出一種基于波峰自動識別的全波形高斯分解算法，將峰值檢測法和奇偶拐點法進行結合，大大提高了分解的速度和精度，但快速獲取有效拐點依然受到限制。段乙好等[9]提出一種高斯拐點匹配法，利用平面曲線離散點、集拐點的快速查找算法檢測拐點，通過相鄰左、右拐點來確定高斯分量，該方法能夠刪減大量的“偽拐點”，對弱波的檢測能力較強，但其分量峰值位置會產(chǎn)生一定的偏移。

針對目前全波形分解方法存在的“偽”拐點數(shù)量多、弱子波形提取困難以及擬合精度不高等問題，本文提出一種在總波形約束下的子波形漸進剝離分解方法。該方法能快速、準確探測分量參數(shù)初值，有效削弱次波對主波的疊加效應影響，對弱波也有一定的探測能力，同時在總波形的約束下進行所有提取子波形的整體擬合，擬合精度較高。

1 次波對主波的影響分析

復雜地物的反射信號通常是一個連續(xù)震蕩的波形，子波形之間的相互疊加效應使真實地物回波信號產(chǎn)生展寬及峰值位置偏移等現(xiàn)象。子波對主波強度的瞬時影響絕對量和相對量與該時刻次波的強度及次波與主波的相對強度直接相關，離子波交點越近的時刻，主波波形的相對和絕對影響均越大。也可認為，波形拐點和峰值位置的偏移程度與主、次波之間的疊加面積大小成正比，次波與主波之間的重疊面積越大，對主波的拐點和峰值位置影響也就越大；反之，則越小。

如圖1所示，當次波位于主波的左側(cè)時，峰值位置左側(cè)的拐點位置受影響程度較大，往左偏離了主波波形，而峰值位置右側(cè)拐點偏移程度很??；同樣當次波位于主波的右側(cè)時，對峰值右側(cè)拐點位置影響較大，因此根據(jù)整體波形的奇偶拐點測定各波形分量的實際峰值位置和半寬值誤差較大，會產(chǎn)生偏移和波形展寬等現(xiàn)象。本文透過次波對主波的影響探討子波形峰值位置偏移和波形展寬的原因，從提升波形參數(shù)的初始精度和減少“偽”拐點的錯誤數(shù)2個方面入手，利用受次波影響較小的主波峰值和近峰值拐點代替左、右拐點提取波形參數(shù)初值，并與峰值配合多種策略，減少“偽”拐點的數(shù)目，穩(wěn)健提取各級子波形，實現(xiàn)高精度整體擬合。

圖1 次波對主波的影響

2 子波形的分解與擬合

2.1 高斯波形及其預處理

一般認為，復雜地物的接收激光脈沖信號是若干個單一地物反射的高斯信號相互疊加的結果，接收器接收的返回激光脈沖可認為是發(fā)射脈沖與探測目標在時間域上的卷積[10]。經(jīng)典的高斯函數(shù)模型和疊加情況如式(1)、式(2)所示。

(1)

(2)

式中：φk(xi)為第k個高斯分量函數(shù)；Ak、Tk和σk分別為第k個高斯分量的振幅、中心位置和半寬值；N為高斯分量個數(shù)，一般不大于6；noise為回波信號的背景噪聲值。

波形復原是將原始回波記錄樣本中反向記錄的數(shù)據(jù)通過波形倒置并進行電壓值轉(zhuǎn)換，將其以真實強度值正向輸出來。本文采用5×5的高斯濾波模板從左往右逐步掃描，將鄰域內(nèi)的加權平均值作為模板中心點的幅值來平滑濾波整個回波信號，以提高波形的信噪比[11]，利用首尾各20幀波形數(shù)據(jù)的均值作為背景噪聲均值并統(tǒng)計其標準差，以2、3倍標準差與噪聲均值之和作為整個波形的振幅閾值[12]，這種方法對符合高斯分布的波形尤為適用?；夭ㄐ盘柕娜ピ胄Ч鐖D2所示，其中，橫坐標序為等間距1 ns的時間序列；縱坐標為波形脈沖的振幅大??；綠線表示原始回波信號；藍線表示去噪后波形；紅線表示振幅閾值線。

圖2 GLA01波形去噪圖

2.2 子波形的漸進剝離

由于GLAS回波信號是以橫坐標為等間距1 ns的離散點集構成，本文采用平面曲線離散點集拐點的快速查找方法[13]來探測拐點的位置及數(shù)量。在拐點處作一條正向切線，其左右凹凸曲線上的離散點集都分別位于切線的兩側(cè)，因此需要4個點就能判斷拐點的位置。設4個點的坐標分別為(xk-2，yk-2)、(xk-1，yk-1)、(xk，yk)、(xk+1，yk+1)，當回波信號相鄰4點坐標滿足式(3)時，(xk，yk)則為拐點。表達式如式(3)所示。

(yk-2+yk-2yk-1)*(yk-1+yk+1-2yk)<0

(3)

通過以上算法可探測整個回波信號中的拐點。由于復雜地物背景噪聲值的干擾會產(chǎn)生許多的“偽”拐點，需要予以刪減。首先，將處于振幅閾值以下的拐點全部刪除；其次，將檢測出來的拐點位于波峰位置左側(cè)稱為左拐點，位于右側(cè)的稱為右拐點，通過判斷波形分量峰值位置，兩側(cè)拐點跟其后一個波形點的單調(diào)性，將左、右拐點之間的異常點進行剔除，當存在多個左、右拐點時，取其平均值作為最終的左、右有效拐點值。由于噪聲和波形疊加的影響，左右拐點經(jīng)常不對稱，通常取左、右拐點分別到波峰位置處的距離較小值作為半寬值。

本文通過計算回波局部最大值來確定高斯分量的峰值和中心位置[14]，采用高斯拐點匹配法求取左、右有效拐點，取左、右拐點到分量中心位置的較小值為半寬值，漸進剝離波形直到剩余波形低于振幅閾值線。具體步驟如下。

1)計算濾波后回波的局部最大值，即為第1個高斯分量的振幅值A1，其對應的橫坐標值為高斯分量的中心位置T1。

2)根據(jù)高斯拐點匹配法，判斷波峰兩側(cè)的左、右有效拐點Tm-1和Tm+1，半寬值σ取|Tm-1-T|和|Tm+1-T|的較小值，如式(4)所示。

σ=min(|Tm-1-T|，|Tm+1-T|)

(4)

3)根據(jù)3個基本高斯參數(shù)即可確定第1個高斯分量，并從波形中剝離出第1個高斯分量，將剩余波形重復步驟1)、步驟2)得到第2個高斯分量，直到剝離的剩余波形局部最大值小于設定的振幅閾值。

4)一般認為，一個復雜地物的回波信號最多分離出6個高斯分量，當探測的高斯分量多余6個時，根據(jù)脈沖和面積合并原則進行高斯分量合并[15]，一般將半寬小于發(fā)射脈寬的一半或者面積較小的分量合并到其鄰近面積最大的高斯分量上，合并后高斯分量振幅取其中較大的高斯分量振幅值，中心位置和半寬值皆取二者合并前的平均值。

5)將剝離出的高斯分量進行整體擬合，求取擬合后的波形與原始波形之間的均方根誤差，判斷其是否滿足擬合精度要求。

2.3 子波形的整體擬合

當激光足印內(nèi)存在坡地或建筑物等具有復雜地物類型時，其回波波形通常較為復雜且伴有波形展寬現(xiàn)象，以及采用高斯濾波去噪會使各子波形振幅變小，因此需要對分解的子波形進行整體優(yōu)化擬合，以使擬合波形最大程度逼近原始回波信號。最小二乘擬合法主要是通過最小化誤差的平方和找尋數(shù)據(jù)的最佳函數(shù)匹配，通過判斷原始數(shù)據(jù)和擬合結果的差值來調(diào)節(jié)迭代系數(shù)的大小，以達到最優(yōu)的擬合結果。其在一定程度上具有彈性擬合的優(yōu)勢，被廣泛應用于非線性方程的求解過程中。

本文采取最小二乘擬合算法對高斯分量參數(shù)進行整體擬合優(yōu)化工作，使優(yōu)化后的擬合波形盡可能逼近原始回波。首先將疊加式(2)展開如式(5)所示。

(5)

式中：An、Tn和σn分別代表第n個高斯分量的振幅值、峰值中心位置和半寬值。

對其進行全微分，表達式如式(6)所示。

f(x)=f(x0)+a1·dA1+b1·dT1+c1·dσ1+

a2·dA2+b2·dT2+c2·dσ2+…

(6)

誤差方程如式(7)所示。

V=AX-L

(7)

A=[a1b1c1a2b2c2…]

(8)

式中：[anbncn]分別為回波函數(shù)f(x)對各子波的振幅A、中心位置T和半寬σ的一階導數(shù)值。

L=[f(x)-f(x0)]

(9)

式中：f(x0)為子波形擬合后的回波信號。

求解誤差方程如式(10)所示。

X=(ATA)-1ATL

(10)

在擬合前后差值L達到最小的情況下，通過解算式(10)可獲得各高斯分量的最佳擬合參數(shù)。漸進剝離與整體擬合方法的流程如圖3所示。

圖3 漸進剝離與整體擬合流程圖

3 實驗結果與分析

3.1 實驗數(shù)據(jù)

GLA01全球測高數(shù)據(jù)產(chǎn)品由美國冰雪數(shù)據(jù)中心(National Snow and Ice Data Center，NSIDC)發(fā)布，記錄回波樣本中包括激光足印的位置、脈沖強度及接收時間等信息。其1 s內(nèi)的40個激光腳點共享一個地理坐標，地表每個激光足印的直徑為72 m，相鄰足印中心間距170 m[16]，從2017年8月份開始，數(shù)據(jù)產(chǎn)品由二進制格式改為HDF5格式進行分發(fā)。本文實驗采用的是1A級GLA01數(shù)據(jù)產(chǎn)品，其主要屬性如表1所示。

表1 實驗數(shù)據(jù)主要屬性

3.2 拐點提取實驗

本文拐點提取實驗以子號為20的波形數(shù)據(jù)為例，利用平面曲線離散點集拐點的快速查找算法對去噪后波形進行拐點提取，然后根據(jù)高斯拐點匹配法的相關原則，將低于振幅閾值線以及不符合左、右拐點條件的“偽”拐點予以刪減，具體方法步驟參考上文。提取的整個波形拐點和經(jīng)刪減后獲得的波形有效拐點分別如圖4、圖5所示，圖中紅色星號代表提取的拐點位置，藍線代表去噪后波形，紅線是根據(jù)背景噪聲值設置的振幅閾值線。

圖4 波形拐點圖

圖5 波形有效拐點圖

3.3 子波形的漸進剝離實驗

首先，進行波形峰值最大值檢測提取第1個高斯分量，將波峰最大值及對應的波形中心位置作為分解的第1個高斯分量的振幅和中心位置，半寬則依據(jù)高斯拐點匹配法確定的左、右拐點來確定，分解的第1個高斯分量如圖6(a)所示。接著，從完整去噪波形中剝離出第1個高斯分量，將剩余波形進行新一輪的波峰最大值檢測和高斯拐點匹配，以提取第2個高斯分量，在圖6(a)的基礎上剝離的第2個高斯分量圖如圖6(b)所示。依此類推，漸進剝離出的高斯分量分別如圖6(c)、圖6(d)、圖6(e)所示，直至剩余波形的最大振幅值低于設定的振幅閾值。圖6(f)是漸進剝離分解后的剩余波形全部位于振幅閾值線以下的結果。圖6(a)至圖6(f)中的橫坐標和縱坐標分別代表回波數(shù)據(jù)的相對時間和電壓幅值，藍線代表原始去噪后波形以及漸進剝離后的剩余波形，綠線代表分解出的高斯分量，紅線代表振幅閾值線。

圖6 子波形漸進剝離分解圖

3.4 子波形的整體擬合實驗

子波形的整體擬合實驗選取主峰在中間、左側(cè)以及右側(cè)3個典型位置的復雜激光測高波形數(shù)據(jù)，分別命名為波形1、波形2和波形3，其索引子號分別為20、15和22。復雜原始回波通過波形分解可獲得若干個子波形，將各子波形的3個高斯基本參數(shù)進行最小二乘擬合，獲取精化后的高斯參數(shù)，以盡可能使擬合后波形最大限度地逼近原始回波信號，從而達到整體擬合的目的。圖7為主峰在中間位置的波形1經(jīng)本文方法分解后，獲得的各個高斯分量及擬合波形與原始去噪后回波信號的對比情況。圖8和圖9分別為主峰位于左側(cè)的波形2和位于右側(cè)的波形3經(jīng)漸進剝離分解并擬合后的結果。圖7至圖9中的藍線代表原始去噪后波形，綠線代表經(jīng)漸進剝離后獲得的高斯分量，紅線代表整體擬合后波形。

圖7 波形1高斯分量擬合圖

圖8 波形2高斯分量擬合圖

圖9 波形3高斯分量擬合圖

3.5 精度對比

為了驗證本文方法在ICESat-GLAS回波信號分解的可靠性和適用性，通過高斯分量數(shù)目以及擬合波形與去噪回波間的均方根誤差、相關系數(shù)和擬合優(yōu)度等評判指標進行定量綜合對比。均方根誤差指的是實際觀測值與擬合值之間的偏差，常用來分析數(shù)據(jù)的可靠性[17]，其計算如式(11)所示。

(11)

式中：N為波形的采樣點個數(shù)；fest(x)為擬合波形；f(x)為原始去噪后波形。

相關系數(shù)反映的是擬合波形和原始去噪后回波之間的密切程度，取值位于[0，1]之間，越逼近1代表二者相關性越密切[18]，其計算如式(12)所示。

(12)

式中：Cov()和Var[]分別指f(x)和fest(x)的協(xié)方差和方差計算。

擬合優(yōu)度也稱為決定系數(shù)，其值是相關系數(shù)的平方，擬合優(yōu)度值越接近1，代表擬合效果越好。

文中選取3種典型的激光測高數(shù)據(jù)波形1、波形2和波形3為研究對象進行測試，以經(jīng)典的奇偶高斯拐點法、峰值檢測法、波峰自動識別的高斯分解法作對比，以驗證本文方法的可靠性和實用性。其中，經(jīng)典的奇偶高斯拐點法由相鄰的奇偶拐點來確定高斯參數(shù)；峰值檢測法是通過查找波形中比左右相鄰元素值都大的峰值點以確定其位置和幅值，再人為定義半寬取值區(qū)間為[3，6]；波峰自動識別法是結合上述2種方法的優(yōu)勢，分解速度提高且精度理論上更為可靠，各方法得到的結果統(tǒng)計如表2所示。同時以波形1為例，通過計算提取的高斯參數(shù)初值與精化后參數(shù)值間的誤差，定量對比分析本文方法分解子波形參數(shù)初值的準確度，其中，ΔA、ΔT和Δσ分別代表分量振幅、中心位置和半寬精化前后的差值，具體結果如表3所示。

表2 不同波形和分解方法的結果統(tǒng)計

表3 高斯分量參數(shù)初值與精化值對比

3.6 實驗分析

提取數(shù)目上，以波形1為例，奇偶高斯拐點法得到的高斯分量個數(shù)最多，達到8個，可見該方法受噪聲影響較大，因噪聲產(chǎn)生的“偽”拐點沒有得到有效去除；峰值檢測法和波峰自動識別法確定的高斯分量相對較少，均為4個，“偽”拐點問題得到了有效解決，但忽略了波形首尾的部分弱子波形；本文方法確定的高斯分量個數(shù)為5個，包括漸進剝離分解出4個較為明顯的子波形和1個隱含的弱子波形，剩余波形的噪聲屬性明顯，表明分解結果的合理性與可靠性。

提取精度上，與其他3種方法相比，均方根誤差大小分別下降75%、66.85%和64.1%；相關系數(shù)分別提升11.16%、1.53%和0.81%；擬合優(yōu)度分別提升23.54%、3.12%和1.64%。由此可見，通過本文方法得到的擬合波形和原始去噪后回波信號偏差較小，二者間的相關性也較為密切。

通過本文方法在典型位置波形1中提取的高斯參數(shù)初值與精化后值對比，各子波形電壓幅值初值與精化值的最大差值為0.186 V，最小差值為0.010 V；峰值中心位置誤差最大差值為1.514 ns，最小差值為0.016 ns；半寬誤差最大差值為1.479 ns，最小差值為0.366 ns?？梢钥闯觯烙嫷母咚狗至繀?shù)初值與精化值誤差較小，表明本文方法能較為準確地測定高斯分量參數(shù)初值，同時經(jīng)過參數(shù)精化和整體擬合后的波形也能近似原始回波信號，從而達到在總波形的約束下精確分解高斯分量的目的。

4 結束語

每個波形分量都對應地物的一次反射特征，波形分解的效果直接影響反演地物特征參數(shù)的精度。本文在分析經(jīng)典波形分解的基本原理以及改進方法的基礎上，利用漸進剝離波形與整體擬合的方法對ICESat-GLAS原始全波形信號進行波形分解，通過漸進剝離波形尋找局部最大值來確定分量峰值和位置，并結合高斯拐點匹配法左、右拐點原則來確定半寬值，最后對3種典型位置的波形利用4種分解方法實驗并進行精度定量對比。實驗結果表明，本文方法能準確測定分量參數(shù)初值和部分弱波，有效避免了傳統(tǒng)分解算法中因噪聲拐點產(chǎn)生的錯誤分量過多以及不合理給定半寬值等問題，同時對高斯分量參數(shù)初值進行精化和整體擬合，使分解的高斯分量受到總波形的約束，從而提高波形分解的精度。

此外，由于次波對主波的疊加效應影響,導致擬合波形主峰峰值變大,以及噪聲去除不完全都會對高斯分量提取和擬合的精度產(chǎn)生一定影響。因此，接下來的研究工作將繼續(xù)深入探討次波對主波的影響因素，如采用最優(yōu)的擬牛頓算法擬合子波形來降低主波峰值偏大的問題，以及結合更高效的波形去噪方法，如用改進的小波變換去噪方法,來降低因高斯濾波去噪引起的各子波振幅變小誤差。