全波形星載激光測距誤差抑制的滑動窗口高斯擬合算法

謝俊峰,劉 仁

1. 自然資源部國土衛星遙感應用中心, 北京 100048； 2. 河海大學地球科學與工程學院, 江蘇 南京 211100

星載激光測高儀作為主動式測量設備，因其較高的測距精度，已廣泛被用于深空探測和對地觀測[1]。在對地觀測方面，2003年1月美國首次發射了GLAS(geoscience laser altimeter system)激光測高系統[2-3]，而后我國于2016年5月發射資源三號02星，該星搭載了一臺激光測高儀作為試驗性載荷用于對地觀測[4-6]。隨著我國激光技術的發展，我國于2019年11月發射了高分七號(GF-7)星載激光測高系統[7-8]，2020年7月發射了資源三號03星激光測高儀[9]，同年12月我國發射的高分十四號衛星(天繪三號)也搭載了一套激光測高系統，2021年還將發射陸地生態碳衛星激光測高系統。

星載激光測高儀在軌后，研究人員主要關注于修正衛星在軌后激光指向角及測距系統誤差[1,4,10-14]，卻鮮有對由激光時間測量引起的測距隨機誤差進行改正的研究。基于閾值時刻鑒別體制的星載激光測高儀，一般利用前后緣閾值時刻估算波形重心位置作為激光出光與返回時刻，進行激光渡越時間測量[15-17]，計算激光測距值。這種模式下測距隨機誤差難以被發現，測距精度相對較低。對于全波形星載激光測高儀(如GLAS，GF-7、陸地生態碳衛星、天繪三號衛星)，其測距方式通過發射與返回波形峰值時間差計算得到[18]，其測時原理如圖1所示。由于激光模擬回波信號經數字化處理后的量化誤差存在，導致基于波形獲取的時間存在一個微小的隨機誤差，這將產生厘米級的測距隨機誤差。而該測距誤差，對于如高分七號搭載的厘米級別測距精度的全波形激光測高儀而言，往往不應被忽略。

圖1 全波形星載激光測高儀測時原理示意Fig.1 Schematic of time measurement of full-waveform spaceborne laser

為了減小全波形星載激光測高儀測距隨機誤差，本文針對其發射與返回波形，提出了一種依賴峰值初始位置的全波形星載激光測距誤差抑制的滑動窗口高斯擬合算法。該方法在確定峰值位置后，應用直線判斷原理，抑制了波形中噪聲，隨后利用有效波形點開展波形擬合，確定最優峰值時刻，從而提取精確的激光測距值。由于GF-7衛星在軌運行期間激光指向整體精度優于1.5″，使得由激光指向抖動引起激光在平坦地形的測距誤差可以忽略不計。故選取平靜的內陸湖面、冰面及江蘇平地為試驗區域，以GF-7星載激光數據為試驗對象，利用一般峰值法與本文方法提取的測距值，通過計算湖面與冰面的激光相對高程精度，以及利用高精度LiDAR[19]計算平地的絕對高程精度，從而驗證了本文方法對測距隨機誤差的改進程度。

1 全波形星載激光測距提取算法

1.1 基于一般峰值法提取的測距誤差分析

通常情況下，波形采樣間隔非常小，大多數大氣和地形條件下激光獲取的回波近乎平滑。一般全波形激光發射和返回波形的最大值對應時刻即為激光發射或者返回時刻，根據它們的時間差計算渡越時間，稱為一般峰值法。然而，實際在軌工作中，星載激光模擬回波信號經數字化處理后產生量化誤差，從而導致了波形峰值處離散采樣點存在似噪聲的波動，使得部分情況下發射與返回波形出現多個最大值，或是峰值點并不是真實的波形峰值點(圖2、圖3)等情況，最終導致基于一般峰值法計算的激光測距存在明顯的隨機誤差。

圖2 GF-7星載激光測高儀3種典型非標準高斯發射波形峰值Fig.2 Three typical non-standard Gaussian peak values of GF-7 spaceborne laser emission waveform

圖3 GF-7星載激光測高儀三種典型非標準高斯返回波形峰值Fig.3 Three typical non-standard Gaussian peak graphs of GF-7 spaceborne laser echo

如圖3所示，一般峰值法依據發射與返回波形提取的激光渡越時間將產生1～5個采樣間隔的隨機誤差，對于目前較高的波形0.5 ns采樣頻率(2 GHz)的GF-7星載激光會導致0.05～2.5 ns時間測量誤差，即0.075～0.375 m的測距誤差。對于其他常用的低采樣頻率星載激光測高儀來說，一般峰值法帶來的隨機測距誤差將翻倍。綜上，對于全波形星載激光測高儀，一般峰值法提取的激光渡越時間直接影響了激光的測距精度，繼而影響星載激光測高精度。

1.2 基于滑動窗口的高斯曲線擬合峰值的測距誤差抑制

1.2.1 波形噪聲抑制與激光渡越時間提取

針對全波形激光模擬回波信號經數字化處理后產生量化誤差，引起的激光測距隨機誤差，本文提出了一種全波形星載激光測距誤差抑制的滑動窗口高斯擬合方法，適用于星載激光各種情況下的發射與返回波形峰值擬合，方法基本流程如圖 4所示。

圖4 基于滑動窗口的高斯曲線擬合方法流程Fig.4 Flow chart of Gaussian curve fitting method based on sliding window

依據上述流程，本文方法核心步驟如下：

(1) 波形初始峰值檢索。

根據輸入波形查找到波形最大值作為初始峰值，以它為中心將波形分為左側(上升沿)與右側(下降沿)波形進行算法實現。

(2) 基于滑動窗口波峰左側(右側)波形有效點提取。

算法設置1×3窗口，由初始峰值點沿波形下降方向(作為正向)進行滑動檢索，沿正向方向單個窗口內第1個點定義為P1，第2、3個點分別定義為P2、P3。若P2=P1，剔除P1，窗口起點滑動至P2重新獲取3個波形采樣點計算，直至P2≠P1。此時利用窗口內前兩個點P1(x1,y1)和P2(x2,y2)構建的正向直線L為

L12(x,y)=(x2-x1)(y-y1)-

(y2-y1)(x-x1)

(1)

對于已構建的正向直線，窗口內第3個波形點所在位置可定義為3類。①直線上點；②內點：沿正向直線L前進方向右側點；③外點：沿正向直線L前進方向左側點。將窗口內波形第3個點P3(x3,y3)，代入式(1)，進行內外點判斷[20]，計算公式為

(2)

對于標準高斯波形，波峰右側波形沿下降方向窗口內第3個點始終為內點；波峰左側恰好相反，窗口內第3個點始終為外點，如圖5(a)所示。

圖5 高斯波形采樣點分布與波峰右側波形噪聲剔除Fig.5 Gaussian waveform sampling point and noise eliminating on right waveform

實際在軌星載激光波形，1×3的窗口內第3個點P3僅存在兩類情況：①波峰左側波形P3為外點或者右側波形P3為內點；②P3為內點(波峰左側波形)或為外點(波峰右側波形)。當在第1種情況下，P1、P2、P3滿足標準高斯曲線分布，直接保留P1、P2、P3，窗口起點滑動至P2重新計算。當在第2種情況下，P3為波形凸點；其中，若P3y=P1y，P1,P2為噪聲點，直接剔除且保留P3，窗口起點滑動至P3再次計算；若P3y≠P1y，P2為噪聲點直接剔除，保留P1,P3且窗口起點滑動至P3再次計算(以波峰右側波形為例，示意圖如圖5(b)所示)。上述兩種情況下，直至波峰左側和右側波形同時各保留N個凸點(N≥2)時停止試驗，對于GF-7星載激光，本文設置N為6。

(3) 高斯曲線擬合。

星載激光測高儀發射與返回波形均滿足標準高斯分布，利用高斯曲線方程對選定的N個波形擬合點進行波形擬合，曲線擬合方程為[21]

(3)

式中，y為波形擬合點幅值；x為波形擬合點時刻；x0為擬合波形峰值時刻；δ為擬合波形脈寬；A為擬合波形幅值。

1.2.2 激光測距值計算

根據激光計時系統記錄的發射與返回波形起始時間，加上它們各自到波峰位置的時間差即可計算出激光渡越時間，并轉換為激光測距值。以GF-7全波形星載激光測高儀為例，其測距公式為[15]

τR=a*((T2+Techo)-(T1+Twf))+b

(4)

(5)

式中，Rlaser為星載激光精確測距值；τR為激光渡越時間；T1為發射波形起點時間；T2為返回波形起點時間；Twf為發射波形峰值至發射波形起點時間；Techo為返回波形峰值至返回起點時間；a為激光計時系統常數因子；b為激光計時偏移改正量；c為光速，c=299 792 458 m/s。

1.3 基于激光腳點高程的測距精度驗證

本文采用控制變量法，在僅改變激光測距值條件下進行激光腳點高程計算，對比在不同測距值下的高程精度。首先，根據星載激光在軌幾何定位原理，考慮大氣[22-23]、潮汐[24-25]引起的測距誤差，構建嚴密幾何定位模型[4]；然后，利用一般峰值法與本文方法提取的測距值分別計算激光點高程；最后，利用平靜湖面、冰面及高精度LiDAR數據，進行兩類測距值的激光點高程相對精度與絕對精度對比。其流程如圖6所示。

圖6 基于激光腳點高程的測距精度驗證流程Fig.6 The flow of ranging accuracy verification based on laser footprint elevation

(1) 基于激光測高相對精度的測距驗證。

借助內陸平靜的湖面和冰面，計算激光在湖面或者冰面上各點的高程，隨后統計該湖面或冰面上的激光高程標準偏差

(6)

理想條件下，在同一平靜的湖面或者冰面激光腳點高程基本一致。由于激光測距隨機誤差存在，同一平靜湖面和冰面激光高程會存在一定波動，可采用激光點高程標準偏差衡量測距精度，即同一湖面或冰面激光高程標差偏差越小，測距精度越高[4]。

(2) 基于激光測高絕對精度的測距驗證。

將激光測距值作為影響激光腳點高程唯一變量，選取高精度的地面控制數據，如：RTK控制點，高精度LiDAR點云，驗證激光腳點高程絕對精度。通過計算激光點與地面控制點高程差值(絕對精度)的均值hmean與中誤差hRMSE，用于對比不同測距值下的激光高程絕對精度，其均值與中誤差計算公式為[4]

(7)

(8)

式中，hli為第i激光點高程；hgi為第i激光點地面實際高程；n為試驗中的激光點個數。

2 試驗與驗證

2.1 試驗數據

本文以GF-7星載激光測高儀為試驗對象，其采用雙波束激光同時對地觀測，每條波束發射能量均為100～180 mJ，脈寬為4～8 ns，發散角為30 μrad[8]。激光接收系統采用數字化回波采樣設備，同時記錄下發射與返回高增益和低增益波形，波形采樣間隔為0.5 ns，回波最大采樣長度為1200個采樣間隔[7]。

由于GF-7衛星過頂時間短，激光地面點數據較少，地面密度小，故而在基于激光測高相對精度的測距驗證試驗中，選取湖泊較大的瑞典維納恩湖和我國江蘇省太湖作為試驗區域。維納恩湖位于瑞典南部，最北邊緯度約為北緯59.3°，GF-7星載激光經過該湖泊的數據為2019年12月19日的第691軌，此時維納恩湖已經進入寒冷冰期，湖面已變成厚厚的冰面。激光在該湖泊分布如圖7所示，由于天氣因素，波束1點有效點8個，波束2有有效點5個。太湖常年無冰期，風浪較小，適于進行試驗，其中GF-7過太湖的激光數據為2020年5月8日第2844軌，如圖8中東南角所示，但由于激光測量時地表云層較厚導致波束1部分數據無返回信號，波束2湖面無激光點。故采用波束1位于太湖湖面8個有效激光點(南北跨越50 km)。

圖7 維納恩湖冰面試驗區域及其試驗數據分布Fig.7 The experimental area and its experimental data map of the ice surface of Lake Vanern

圖8 江蘇省試驗區域與試驗數據分布Fig.8 The experimental area and its experimental data map of Jiangsu province

基于激光測高絕對精度的測距驗證試驗，本文選取GF-7星載激光過江蘇省西北部宿遷地區的激光點，如圖8中西北角所示，該數據為2019年12月12日GF-7衛星第595軌。該區域范圍為：33.152 46°N—34.683 58°N，118.104 62°E—118.619 99°E，其中波束1共55個激光點，波束2共67個激光點。用于進行激光點高程絕對精度驗證數據，為采用徠卡ALS70機載激光掃描系統獲取區域內的LiDAR點云數據。該設備最大脈沖頻率為500 kHz，最大掃描頻率為200 Hz，最大視場角為75°，可接收無限次回波記錄。現場測量時飛行航高為2400 m，地面光斑大小約為60 cm，獲取的點云密度約為1.53/m2，高程精度中誤差為0.12 m，數據采集時間為2018年，獲取的數據面積約3432 km2。

2.2 基于平靜湖面的激光高程相對精度驗證

2.2.1 波形噪聲抑制試驗與分析

采用本文基于滑動窗口高斯擬合方法，分別對GF-7星載激光第2844軌波束1湖面上8個激光點，以及第691軌冰面上波束1的8個激光點和波束2的5個激光點，共計21個激光點的發射與返回波形進行試驗。剔除波形中噪聲點，利用自動挑選出波峰周圍12個波形點進行波形峰值擬合，3類典型的發射與返回波形擬合結果如圖9和圖10所示。衛星實際在軌后，激光波形峰值形狀主要為6種情況，其中發射波形3種，返回波形3種，以下利用本文方法對這6種情況下的波形擬合結果進行詳細分析。

圖9 發射波形擬合結果Fig.9 Fitting results of emission waveform

圖10 返回波形擬合結果Fig.10 Fitting results of echo waveform

針對以上6種不同的情況下的波形，由整體結果可以看出，本文算法均能準確地擬合出回波波形，其中每類波形脈寬吻合度極高，擬合波形峰值與實際回波最大值點幅值相當。從擬合波形細節觀察，本文算法避開了各種情況下的原始波形中噪聲點；且在上述6種典型波形情況下，擬合后波形峰值均能位于原始波形幅值最大的兩個采樣點之間。充分說明本文算法能夠有效地擬合出實際波形的峰值，抑制激光發射與回波波形中噪聲，從而可根據擬合后的波峰即可精確獲取激光發射或返回波形時刻。

根據一般峰值法與本文方法，計算第691軌冰面上波束1的8個激光點和波束2的5個激光點的發射脈沖時刻(Twf)與返回時刻(Techo)，結果見表1。兩種方法計算的發射和返回波形峰值時間差分別為ΔTwf、ΔTecho；由兩種方法計算的測距差值為ΔTrange。

表1 一般峰值法與本文方法提取的激光渡越時間差值Tab.1 The time difference of laser transit time between general peak method and this paper method

由表1結果可分析，無論是星載激光測高儀發射波形還是返回波形，經本文方法擬合后波形峰值時刻發生了一定偏移。其中，相對于一般峰值法，本文方法提取的發射波形峰值點位置偏移了0.235 ns；本文方法提取的返回波形峰值點位置偏移了0.49 ns。將發射波形與返回波形峰值偏移量應用于GF-7星載激光測高儀測距計算式(4)之中，本文方法提取的激光渡越時間與一般峰值法提取的渡越時間偏差為0.50 ns，相當于7.5 cm的測距誤差。根據上述分析可知，本文方法能夠明顯減小激光測距隨機誤差。

2.2.2 平靜湖面激光高程相對精度驗證

為分析本文方法提取的激光測距值，能否有效提升星載激光測高儀相對測高精度。根據在軌檢校后激光指向角，分別利用一般峰值法與本文方法提取的測距值，對第2844軌波束1太湖湖面上8個激光點，以及第691軌冰面上波束1的8個激光點和波束2的5個激光點進行激光腳點高程解算，并轉換至正常高，結果見表2、表3。試驗中，僅有激光測距唯一變量，其余參數，如激光指向、大氣改正/潮汐改正值、高程異常等完全一致。

根據表2、表3結果可以看出，在維納恩湖冰面上，GF-7星載激光波束1與波束2高程相對精度基本一致，其中利用一般峰值法提取的激光測距值計算的激光高程標準偏差為0.125 m，基于本文方法提取的激光測距值計算相同激光點高程的標準偏差為0.082 m。對于內陸太湖，利用一般峰值法提取的激光測距值計算太湖湖面8個激光點正常高的標準偏差為0.123 m；而基于本文方法提取的激光測距值計算相同激光點正常高的偏差為0.065 m。其中本文方法對于太湖湖面的激光點精度提升更高，然而實際情況下維納恩湖冰面結果更加準確、可信，因為太湖湖面或多或少存在一定風浪，產生了該現象。綜上，本文方法提取的激光測距值對星載激光測高儀相對高程精度提升了近4.2 cm，由原始的12.5 cm提升至8.3 cm。

表3 GF-7激光點在瑞典維納恩湖冰面高程相對精度Tab.3 The relative elevation accuracy of GF-7 laser footprint on the ice surface of Vanern Lake in Sweden

表2 GF-7激光點在太湖湖面高程相對精度

假設理想條件下，即激光測距無誤差、同一無風浪湖面或冰面上的激光測高應該完全一致。因此，激光在冰面或無風浪湖面的高程相對偏差直接大致可以反映出星載激光測高儀測距精度。結合表3冰面的試驗結果分析，由表1中可以看出，相對于一般峰值法，基于本文方法提取的激光測距精度提升了7.5 cm。

2.3 基于高精度機載LiDAR的激光高程絕對精度驗證

分別對GF-7第595軌江蘇境內的波束1的55個激光點和波束2的67個激光點進行波形擬合，并利用一般峰值法與本文方法提取出兩套測距值。在其他參數完全一致的條件下，根據激光幾何定位模型計算每個激光腳點大地坐標。從激光足印落點位置，發現GF-7衛星第595軌激光過宿遷市等多個城區，如圖8所示，部分星載激光點落在房屋等建筑與林木上。剔除這些異常激光點后，波束1和波束2分別剩余48和61個激光點。利用機載LiDAR點云內插激光點高程，將激光高程與內插的LiDAR的高程差的均值和中誤差作為GF-7星載激光高程絕對精度。GF-7星載激光平地地區兩類測距計算的高程絕對精度如圖11和圖12所示。

由圖11和圖12結果顯示，無論是GF-7激光波束1或波束2，本文方法提取的測距值計算的激光高程差曲線向0值附近收縮，可判斷出本文方法能夠明顯提升星載激光高程絕對精度。為定量分析其對激光高程絕對精度提升空間，統計上述波束1與波束2所有點在兩類測距值下的均值與中誤差，結果見表4。

圖11 兩類測距值下第595軌波束1激光高程絕對精度Fig.11 The 595th track beam 1 laser elevation absolute precision under two sets of ranging value

圖12 兩類測距值下第595軌波束2激光高程絕對精度Fig.12 The 595th track beam 2 laser elevation absolute precision under two sets of ranging value

由表4可以看出，在平原地區，基于一般峰值法提取的測距值計算的GF-7星載激光波束1激光高程絕對精度為-0.056±0.192 m；基于本文方法提取的測距值計算的波束1激光高程絕對精度為-0.040±0.177 m。同理，GF-7星載激光波束2在兩類測距值下的高程絕對精度分別為-0.038±0.238 m、-0.029±0.189 m。總體來看，本文方法可將GF-7星載激光高程絕對精度由初始-0.047±0.215 m提升至-0.035±0.183 m，提升了4.5 cm。從量級上看，高程絕對精度提升不大，但對于絕對高程精度本已非常高的GF-7星載激光測高儀而言，再提升4.5 cm精度具有十分重要的現實意義。

表4 兩類測距值下平地地區激光測高絕對精度

3 結 論

本文在分析星載激光測高儀全波形數據過程中，發現波峰附近存在明顯類似噪聲的現象，采用一般峰值法提取的測距值存在較大隨機誤差的問題，故而提出了一種全波形星載激光測距誤差抑制的滑動窗口高斯擬合算法，用于提取激光測距值，該方法明顯提高了星載激光測距精度。利用一般峰值法與本文方法提取的測距值，對瑞典維納恩湖冰面、江蘇太湖湖面，以及江蘇平地地區的GF-7星載激光數據，分別進行激光高程相對與絕對精度驗證，得到相關結論如下：

(1) 本文方法對噪聲不敏感，可適用于星載激光各類情況下波形，正確擬合出波峰位置。本文方法提取的測距值，較一般峰值法，測距精度提升了7.5 cm。

(2) 利用瑞典維納恩湖冰面的GF-7激光點為試驗對象，經驗證結果表明，經本文方法提取的激光測距值，計算的高程相對精度提升了4.2 cm。

(3) 以平地地區高精度機載LiDAR點云為高程驗證數據，本文方法相對一般峰值法提取的測距值，計算的激光高程絕對精度提升了4.5 cm。

綜上，本文提出的方法能夠抑制激光發射與回波波形噪聲，有效地減小全波形星載激光測距隨機誤差，提升了激光腳點高程精度，目前已用于GF-7星載激光數據業務化處理中。針對不同地物區域的星載激光返回波形引起的測距誤差，還需進一步開展深入的研究分析。