韓曉爽，潘 超，劉宇哲，王麗東，戴學(xué)兵，趙一鳴

（北京遙測技術(shù)研究所 北京 100094）

引 言

激光高度計衛(wèi)星的研制已有近30 年的發(fā)展歷史。從1994 年美國NASA 首次將激光測高儀搭載在航天飛機(jī)進(jìn)行對地觀測探索開始，火星軌道激光高度計MOLA（Mars Orbiter Laser Altimeter）、搭載于ICESat-1（Ice，Cloud，and land Elevation Satellite）衛(wèi)星的地球科學(xué)激光測高儀系統(tǒng)GLAS（Geoscience Laser Altimeter System）、水星激光高度計MLA（Mercury Laser Altimeter）及月球軌道激光高度計LOLA（Lunar Orbiter Laser Altimeter）分別在陸地高程、森林、冰蓋等探測中取得了豐富的成果。2018 年發(fā)射的ICESat-2 衛(wèi)星搭載更為先進(jìn)的地形激光高度計載荷ATLAS（Advanced Topographic Laser Altimeter System），率先采用6 波束激光、微脈沖、單光子探測體制，憑借其精確的測高能力，在山岳、冰川及冰帽、地表地貌、河流及湖泊深度、近海水深、海洋潮汐、高緯地區(qū)大地水準(zhǔn)面、海洋重力場等地球科學(xué)的多個領(lǐng)域均廣泛應(yīng)用。NASA 規(guī)劃的最新一代星載全球高分辨率成像測繪激光雷達(dá)系統(tǒng)LIST（Lidar Surface Topography）采用線陣推掃模式，具備多達(dá)1000 束激光的探測能力，可提供全球三維地形數(shù)據(jù)，預(yù)計2025 年發(fā)射。

多年來，在國家各類項目的支持下，我國星載激光雷達(dá)技術(shù)也得到了長足發(fā)展，已經(jīng)具有較深厚的技術(shù)儲備。從“嫦娥”探月系列衛(wèi)星上搭載的激光測高儀開始，搭載了激光雷達(dá)載荷的“資源三號”02星、“高分七號”已經(jīng)分別于2016 年和2019 年成功發(fā)射并在軌應(yīng)用，搭載了星載多波束激光雷達(dá)載荷的“陸地生態(tài)系統(tǒng)碳監(jiān)測衛(wèi)星”也將于2021 年發(fā)射，實現(xiàn)對全球植被高度與大氣分布的協(xié)同探測。多家單位研制的機(jī)載多波束三維成像激光雷達(dá)類載荷也開展了飛行試驗。其中，由北京遙測技術(shù)研究所自主研制的機(jī)載大光斑激光雷達(dá)于2017 年分別在石家莊、張家界開展了多架次飛行試驗，實現(xiàn)了星載多波束激光雷達(dá)載荷的機(jī)載驗證與反演算法優(yōu)化。

在星載激光測高儀激光地面腳點坐標(biāo)解算方面，文獻(xiàn)[1-4]介紹了ICESat-1/GLAS 的腳點定位模型和誤差。ICESat-1/GLAS 和ICESat-2/ATLAS 的ATBD（Algorithm Theoretical Basis Document）報告分別對腳點定位流程及誤差有較為詳細(xì)的介紹。文獻(xiàn)[5-7]也較為系統(tǒng)地闡述了衛(wèi)星激光測高嚴(yán)密幾何模型及精度分析，同時，文獻(xiàn)采用“資源三號”02 星的激光測高數(shù)據(jù)對其進(jìn)行了處理和驗證。

本文對單波束、全波形和多波束、單光子兩種探測體制激光載荷的坐標(biāo)解算方法進(jìn)行分析，在上述方法的指引下，構(gòu)建機(jī)載平臺坐標(biāo)轉(zhuǎn)換模型及方法，并通過機(jī)載飛行試驗進(jìn)行了驗證。

1 星載激光雷達(dá)坐標(biāo)解算方法

星載激光雷達(dá)地面腳點坐標(biāo)的解算方法可基本概括為：通過激光發(fā)射和接收的時間差t，以及激光傳播速度c，求得測距值，再結(jié)合衛(wèi)星搭載的GPS 和星敏感器獲得衛(wèi)星位置和姿態(tài)信息，求得激光地面點在目標(biāo)坐標(biāo)系下的三維坐標(biāo)。

根據(jù)衛(wèi)星實際情況，星載激光雷達(dá)地面腳點坐標(biāo)解算涉及瞬時激光束坐標(biāo)系、測量參考坐標(biāo)系（光學(xué)平臺坐標(biāo)系）、衛(wèi)星本體坐標(biāo)系、國際天球參考框架（ICRF）下的J2000 坐標(biāo)系和國際地球參考框架（ITRF）下的WGS84 坐標(biāo)系五個坐標(biāo)系[8]。

轉(zhuǎn)換流程如圖1 所示[9]，其中下標(biāo)SC、ICRF 表示衛(wèi)星本體坐標(biāo)系，COM 為衛(wèi)星質(zhì)心，Ref 為激光參考點，Laser Spot 為激光地面腳點，Δref為激光參考點和衛(wèi)星質(zhì)心的位置偏移量，h為激光傳輸距離。由RSPOT=RSC+Δref+h可求得激光地面腳點在ICRF 框架下的坐標(biāo)，再將其轉(zhuǎn)換至ITRF 參考框架下的坐標(biāo)。獲得足印點的三維坐標(biāo)（X,Y,Z）后，也可用大地橢球坐標(biāo)系（B,L,H）來表示。

在星載激光雷達(dá)地面坐標(biāo)點的解算過程中，基本沿用以上原理進(jìn)行坐標(biāo)轉(zhuǎn)換。根據(jù)測高體制、波束分布的不同，在具體的轉(zhuǎn)換流程中也存在區(qū)別。下文分別對ICESat-1、ICESat-2的坐標(biāo)解算方法進(jìn)行分析，探尋不同載荷體制下坐標(biāo)解算過程中的異同。

1.1 ICESat-1/GLAS

ICESat-1 于2003 年1 月發(fā)射成功，其搭載的GLAS 是目前較為成功的對地星載激光測高系統(tǒng)。GLAS為單波束、全波形測高系統(tǒng)，脈沖頻率40Hz，沿軌方向足印間隔170m，地面足印大小55m～110m，其科學(xué)探測目標(biāo)為南北極冰蓋監(jiān)測、陸地及水文地貌監(jiān)測、植被高度監(jiān)測和云-氣溶膠高度分布監(jiān)測。公開發(fā)布的數(shù)據(jù)產(chǎn)品包括GLA0-GLA15 在內(nèi)的0 級、1 級和2 級數(shù)據(jù)產(chǎn)品。

GLAS 的激光腳點解算方法符合星載激光雷達(dá)坐標(biāo)解算流程，但過程有所簡化。采用ANC04 慣性坐標(biāo)系旋轉(zhuǎn)矩陣、ANC08 衛(wèi)星精密定軌（POD）數(shù)據(jù)、ANC09 激光指向和精密定姿（PAD）數(shù)據(jù)、ANC25時間轉(zhuǎn)換數(shù)據(jù)，與計算得到的激光測距值聯(lián)合解算，得到激光地面點在ITRF 框架下的坐標(biāo)。其計算流程如圖2 所示。

圖1 星載激光雷達(dá)地面腳點解算流程Fig.1 Spaceborne laser altimeter footprint coordinate calculation process

圖2 ICESat-1/GLAS 地面坐標(biāo)解算流程Fig.2 Footprint coordinate calculation process of ICESat-1/GLAS

由于GLAS 發(fā)射波形可近似為高斯脈沖，激光經(jīng)地表反射后的回波波形可視作一個或多個高斯脈沖的疊加，采用波形分解的方法對回波波形進(jìn)行擬合，從而求得激光單向傳輸?shù)木嚯x。

1.2 ICESat-2/ATLAS

已于2018 年9 月15 日發(fā)射的ICESat-2 衛(wèi)星是ICESat-1 衛(wèi)星的延續(xù)，搭載了更先進(jìn)的地形激光測高系統(tǒng)載荷ATLAS。ATLAS 為非掃描推帚式、高重頻、微脈沖、光子計數(shù)式的多光束激光測高載荷，軌道高度500km，軌道傾角92°，重訪周期91 天。ATLAS 激光脈沖頻率達(dá)10kHz，激光腳點沿軌方向間距0.7m，地面足印直徑為17m；采用6 波束探測，其中3 波束能量較強(qiáng)，3 波束能量較弱，能量比為4:1；每組強(qiáng)激光腳點之間的垂軌距離3.3km，沿軌距離2.5km，強(qiáng)-弱激光腳點之間的垂軌距離為90m，如圖3 所示[10,11]。相較GLAS 載荷，ATLAS 載荷可提供空間分辨率和平面定位精度更高的數(shù)據(jù)產(chǎn)品。

ATLAS 的激光腳點解算方法總體上與GLAS 相似，采用ATL02 數(shù)據(jù)、ANC03 慣性坐標(biāo)系旋轉(zhuǎn)矩陣、ANC04 衛(wèi)星精密定軌（POD）數(shù)據(jù)、ANC05 激光指向和精密定姿（PAD）數(shù)據(jù)聯(lián)合解算，得到激光地面點在ECF 坐標(biāo)系下的坐標(biāo)，其實現(xiàn)流程如圖4 所示。

圖3 ICESat-2 激光地面足印示意圖Fig.3 Laser footprints of ICESat-2

圖4 ICESat-2/ATLAS 地面坐標(biāo)解算流程Fig.4 Footprint coordinate calculation process of ICESat-2/ATLAS

基于脈沖時間飛行法TOF（time of flight）的光子計數(shù)激光雷達(dá)通過記錄起始脈沖信號與目標(biāo)返回信號的時間差來測量距離。設(shè)激光光源發(fā)射脈沖信號的時刻為tT，脈沖信號傳播至目標(biāo)表面經(jīng)漫反射、接收光學(xué)系統(tǒng)收集至探測器的時刻為tR，介質(zhì)中的光速為c，則理想情況下探測距離為

針對光子計數(shù)體制的特點，載荷ATLAS 的地面腳點坐標(biāo)解算方法也進(jìn)行了相應(yīng)的變化。在ATL02數(shù)據(jù)中，將N個光子的數(shù)據(jù)合并為一組進(jìn)行處理，每組光子對應(yīng)時間約為5ms，對應(yīng)沿軌距離約為40m，通過解算獲得所有N個光子的地面腳點坐標(biāo)。后續(xù)在進(jìn)行大氣參數(shù)和大氣延遲校正時，為節(jié)省存儲與運算資源，從N個光子中選取第j個光子，其指向向量作為后續(xù)校正所用向量[12,13]。

具體解算步驟如下：

⑥利用ANC03 數(shù)據(jù)得到光子到達(dá)地面點時刻tB(i)下ECI 坐標(biāo)系向ECF 坐標(biāo)系的旋轉(zhuǎn)矩陣，并計算激光地面點在ECF 坐標(biāo)系下的坐標(biāo)向量。

⑦將激光地面點在ECF 坐標(biāo)系空間直角坐標(biāo)系下的坐標(biāo)轉(zhuǎn)化為大地坐標(biāo)系下的坐標(biāo)形式。

2 機(jī)載激光雷達(dá)坐標(biāo)解算方法

機(jī)載激光雷達(dá)坐標(biāo)解算模型與星載模型較為相似，通過發(fā)射點到地面反射點距離、慣性導(dǎo)航系統(tǒng)（INS）測定的飛機(jī)姿態(tài)參數(shù)及GPS 測定的飛機(jī)精確位置聯(lián)合解算獲取地面點在目標(biāo)坐標(biāo)系下的三維坐標(biāo)。星載平臺通過星敏感器或陀螺儀定姿，機(jī)載激光雷達(dá)通過慣導(dǎo)定姿，在坐標(biāo)轉(zhuǎn)換過程中不涉及天球坐標(biāo)系，不需考慮歲差、章動、極移等現(xiàn)象。參考激光測高衛(wèi)星定位模型和ICESat-1/GLAS 公布的處理流程，可構(gòu)建機(jī)載平臺坐標(biāo)轉(zhuǎn)換模型，進(jìn)行激光地面腳點坐標(biāo)解算，以下分別對機(jī)載激光雷達(dá)涉及的坐標(biāo)系統(tǒng)和解算方法進(jìn)行介紹。

2.1 坐標(biāo)系統(tǒng)

坐標(biāo)系統(tǒng)包括瞬時激光束坐標(biāo)系、激光掃描參考坐標(biāo)系、慣性平臺參考坐標(biāo)系、導(dǎo)航坐標(biāo)系和大地坐標(biāo)系。其定義分別如下[14]：

①瞬時激光束坐標(biāo)系：定義激光發(fā)射參考點為原點O，X軸為飛行方向，Z軸為瞬時光束方向，O-XYZ構(gòu)成右手系。

②激光掃描參考坐標(biāo)系（L 系）：定義激光發(fā)射參考點為原點O，XL軸為飛行方向，ZL軸為掃描系統(tǒng)的零點方向，O-XLYLZL構(gòu)成右手系。

③慣性平臺參考坐標(biāo)系（I 系）：定義慣性平臺中心為原點O，XI軸、YI軸ZI軸方向按照IMU 參考標(biāo)架定義，O-XIYIZI構(gòu)成右手系。

④導(dǎo)航坐標(biāo)系（N 系）：又稱為當(dāng)?shù)厍衅矫孀鴺?biāo)系，GPS 天線相位中心為原點O，XN軸指向真北，YN軸指向東，ZN軸為鉛直向下方向，O-XNYNZN構(gòu)成右手系。

⑤大地坐標(biāo)系（E 系）：采用WGS84 坐標(biāo)系，地球質(zhì)心為原點O，Z軸指向國際時間局（BIH）1984.0定義的協(xié)議地級（CTP）方向，XE軸指向BIH1984.0 的協(xié)議子午面和CTP 赤道的交點，O-XEYEZE構(gòu)成右手系。

2.2 坐標(biāo)轉(zhuǎn)換過程

機(jī)載激光雷達(dá)對地定位中的坐標(biāo)轉(zhuǎn)換順序為瞬時激光束坐標(biāo)系→激光掃描參考坐標(biāo)系（L 系）→慣性平臺參考坐標(biāo)系（I 系）→導(dǎo)航坐標(biāo)系（N 系）→大地坐標(biāo)系（E 系）[15,16]，各坐標(biāo)系的定義2.1 節(jié)已進(jìn)行闡述，具體轉(zhuǎn)換過程如下。

①瞬時激光束坐標(biāo)系轉(zhuǎn)L 系

在激光掃描參考坐標(biāo)系下定義掃描角θ為激光束與ZL的夾角，逆時針轉(zhuǎn)動方向為正，如圖5 所示。

當(dāng)激光測距值為ρ時，地面點在L 系下的坐標(biāo)可表示為

圖5 激光掃描參考坐標(biāo)系Fig.5 Laser scanning reference system

② L系轉(zhuǎn)I 系

可通過測量得到L 系和I 系坐標(biāo)原點間的距離，即慣性平臺中心與激光發(fā)射參考點的偏移量為[ΔXLΔYLΔZL]T；坐標(biāo)軸之間的夾角為偏心角α、β、γ，即慣導(dǎo)與激光掃描儀在安裝過程中三個坐標(biāo)軸產(chǎn)生的固定軸線偏差轉(zhuǎn)角。地面點在I 系下的坐標(biāo)可表示為

③I 系轉(zhuǎn)N 系

N 系和I 系坐標(biāo)原點間的距離為偏移量[ΔXNΔYNΔZN]T，坐標(biāo)軸之間的夾角可由IMU 記錄的三個姿態(tài)角航向角H(heading)、俯仰角P(pitch)和側(cè)滾角R(roll)表示，地面點在N 系下的坐標(biāo)可表示為

④N 系轉(zhuǎn)E 系

GPS 系統(tǒng)可以實時獲得天線相位中心即導(dǎo)航坐標(biāo)系原點的參心大地坐標(biāo)（B,L,H）（緯度，經(jīng)度和橢球高），[XgpsYgpsZgps]T為GPS 在WGS84 坐標(biāo)系下的坐標(biāo)，可利用下面的公式進(jìn)行計算：

其中，e 為第一偏心率，N為卯酉圈曲率半徑，，a 為WGS84 橢球的長半軸。

地面點在WGS84 下的坐標(biāo)可表示為

根據(jù)上述的坐標(biāo)轉(zhuǎn)換公式，最終可求得激光地面點在WGS84 坐標(biāo)系下的坐標(biāo)表示為

其中，[XoffsetYoffsetZoffset]T為激光參考點到GPS 天線中心的偏移量。

通過上述坐標(biāo)系統(tǒng)的轉(zhuǎn)換，可以得到每個激光地面點在地心坐標(biāo)系WGS84 下的三維坐標(biāo)值?？梢愿鶕?jù)需要將激光地面點的空間直角坐標(biāo)轉(zhuǎn)換為大地坐標(biāo)。

3 機(jī)載驗證系統(tǒng)——大光斑激光雷達(dá)

機(jī)載大光斑激光雷達(dá)是為滿足碳監(jiān)測衛(wèi)星多波束激光雷達(dá)載荷的算法驗證需求而研制的機(jī)載驗證載荷，采用單波束、全波形體制，通過將機(jī)載大光斑激光雷達(dá)與航空相機(jī)、航空慣導(dǎo)系統(tǒng)、高穩(wěn)平臺、控制及存儲系統(tǒng)等設(shè)備的集成和優(yōu)化，形成機(jī)載大光斑激光雷達(dá)系統(tǒng)，以40Hz 的重頻發(fā)射激光脈沖進(jìn)行地表探測[17,18]。系統(tǒng)參數(shù)指標(biāo)如表1 所示，系統(tǒng)組成如圖6 所示。

表1 機(jī)載大光斑激光雷達(dá)系統(tǒng)參數(shù)指標(biāo)Table 1 Parameters of airborne large spot LiDAR system

4 機(jī)載飛行驗證試驗數(shù)據(jù)結(jié)果

機(jī)載大光斑激光雷達(dá)系統(tǒng)于2017年4 月上旬，在石家莊欒城機(jī)場附近區(qū)域開展了首次機(jī)載飛行試驗，并開展了點云激光雷達(dá)同步試驗比對、地面外業(yè)試驗交叉比對與精度評估。此次飛行試驗包括兩個架次的校飛試驗，獲取了檢校區(qū)與測區(qū)共計10 條航線，約300GB 的有效數(shù)據(jù)。系統(tǒng)再次于2017 年12 月下旬，在張家界荷花機(jī)場附近區(qū)域開展了兩個架次的飛行試驗。此次飛行試驗航線總長超過980km，航高超過3000m，航時總計11 小時，獲取了檢校區(qū)與測區(qū)共計四個測區(qū)的有效數(shù)據(jù)。

圖6 機(jī)載大光斑激光雷達(dá)系統(tǒng)構(gòu)成Fig.6 System configuration of airborne large spot LiDAR system

對機(jī)載大光斑激光雷達(dá)系統(tǒng)于2017 年12 月25 日在張家界附近采集到的一段數(shù)據(jù)進(jìn)行解算，可得到地面激光點的軌跡如圖7 所示，部分?jǐn)?shù)據(jù)如表2 所示。

圖7 部分地面激光點軌跡Fig.7 Partial laser footprint track

表2 經(jīng)解算得到的部分激光地面點坐標(biāo)Table 2 Partial calculated laser footprint coordinates

獲得回波信號數(shù)據(jù)后，通過高斯濾波方法對信號進(jìn)行去噪處理，再用最小二乘曲線擬合法或最大值檢測法可以得到該點位測距值[19,20]。本論文通過計算回波信號的二階差分獲取峰值，采用公式如下：

其中，w為波形數(shù)據(jù)，有：

利用find 函數(shù)獲得滿足式的元素的序號，k值所在位置即為回波峰值所在的點位，當(dāng)回波信號有多個峰值時，通常選取回波強(qiáng)度最大的位置作為峰值所在點。利用峰值位置計算出激光測距值后，根據(jù)第

2.2 節(jié)介紹的坐標(biāo)轉(zhuǎn)換方法可解算出對應(yīng)的激光地面點坐標(biāo)。

以表1 第8 行數(shù)據(jù)為例，此時UTC 為12987.5，采集到的回波信號波形如圖8 所示，通過地圖定位可知該位置為樹木，坐標(biāo)解算結(jié)果與信號符合。表1 第25 行數(shù)據(jù)對應(yīng)UTC 為13017.85，采集到的回波信號波形如圖9 所示，通過地圖定位可知該位置為裸地，信號強(qiáng)度過大，存在飽和現(xiàn)象，坐標(biāo)解算結(jié)果與信號符合。

圖8 激光回波信號數(shù)據(jù)與對應(yīng)地面腳點坐標(biāo)（UTC=12987.5）Fig.8 Reflected laser echo signal data and footprint coordinate（UTC=12987.5）

圖9 激光回波信號數(shù)據(jù)與對應(yīng)地面腳點坐標(biāo)（UTC=13017.85）Fig.9 Reflected laser echo signal data and footprint coordinate（UTC=13017.85）

5 結(jié)束語

本文通過分析ICESat-1 和ICESat-2 星載激光測高載荷地面腳點坐標(biāo)解算方法，比較了多波束、單光子體制對地激光三維測繪雷達(dá)與傳統(tǒng)的單波束、全波形體制激光測高雷達(dá)在地面坐標(biāo)解算過程的異同，并建立了機(jī)載激光雷達(dá)坐標(biāo)解算方法與流程，通過自主研制的機(jī)載大光斑激光雷達(dá)與飛行試驗，對坐標(biāo)解算方法進(jìn)行驗證與分析。隨著單光子探測技術(shù)的發(fā)展，光子計數(shù)激光雷達(dá)正在逐步應(yīng)用于機(jī)載、星載對地觀測領(lǐng)域，北京遙測技術(shù)研究所也正在進(jìn)行機(jī)載光子計數(shù)激光雷達(dá)研制。通過對ICESat-2 地面坐標(biāo)解算流程作進(jìn)一步理論分析，改進(jìn)并構(gòu)建適用于機(jī)載光子計數(shù)激光雷達(dá)的坐標(biāo)解算方法，并應(yīng)用于機(jī)載飛行驗證試驗，是下一步工作的重點。