一種漸進剝離與整體擬合的子波形提取方法

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

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

0 引言

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

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

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

1 次波對主波的影響分析

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

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

圖1 次波對主波的影響

2 子波形的分解與擬合

2.1 高斯波形及其預處理

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

(1)

(2)

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

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

圖2 GLA01波形去噪圖

2.2 子波形的漸進剝離

由于GLAS回波信號是以橫坐標為等間距1 ns的離散點集構成，本文采用平面曲線離散點集拐點的快速查找方法[13]來探測拐點的位置及數量。在拐點處作一條正向切線，其左右凹凸曲線上的離散點集都分別位于切線的兩側，因此需要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)

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

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

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

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

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

(4)

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

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

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

2.3 子波形的整體擬合

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

本文采取最小二乘擬合算法對高斯分量參數進行整體擬合優化工作，使優化后的擬合波形盡可能逼近原始回波。首先將疊加式(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]分別為回波函數f(x)對各子波的振幅A、中心位置T和半寬σ的一階導數值。

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

(9)

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

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

X=(ATA)-1ATL

(10)

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

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

3 實驗結果與分析

3.1 實驗數據

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

表1 實驗數據主要屬性

3.2 拐點提取實驗

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

圖4 波形拐點圖

圖5 波形有效拐點圖

3.3 子波形的漸進剝離實驗

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

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

3.4 子波形的整體擬合實驗

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

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

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

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

3.5 精度對比

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

(11)

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

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

(12)

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

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

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

表2 不同波形和分解方法的結果統計

表3 高斯分量參數初值與精化值對比

3.6 實驗分析

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

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

通過本文方法在典型位置波形1中提取的高斯參數初值與精化后值對比，各子波形電壓幅值初值與精化值的最大差值為0.186 V，最小差值為0.010 V；峰值中心位置誤差最大差值為1.514 ns，最小差值為0.016 ns；半寬誤差最大差值為1.479 ns，最小差值為0.366 ns。可以看出，估計的高斯分量參數初值與精化值誤差較小，表明本文方法能較為準確地測定高斯分量參數初值，同時經過參數精化和整體擬合后的波形也能近似原始回波信號，從而達到在總波形的約束下精確分解高斯分量的目的。

4 結束語

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

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