劉迎春,高顯連,賀 巖,崔晨彥,高劍新,俞家勇,蔡龍濤,吳發(fā)云

(1.國家林業(yè)和草原局調(diào)查規(guī)劃設(shè)計(jì)院,北京 100714；2.中國科學(xué)院上海光學(xué)精密機(jī)械研究所,上海 201800；3.安徽建筑大學(xué),合肥 230601；4.東北林業(yè)大學(xué),哈爾濱 150040)

森林是全球重要的陸地生態(tài)系統(tǒng),可通過樣地調(diào)查的方法評估其資源量和生物量[1-2]。隨著激光雷達(dá)技術(shù)的發(fā)展,采用星載大光斑激光雷達(dá)估算大區(qū)域森林地上生物量將成為另一種選擇。由于星載激光雷達(dá)數(shù)據(jù)少,近20年基本都圍繞著美國的GLAS(Geoscience Laser Altimeter System)數(shù)據(jù)開展激光雷達(dá)林業(yè)應(yīng)用研究[3-6],或者采用仿真大光斑激光雷達(dá)的方法進(jìn)行前瞻性研究[7]。直到近年,GEDI(Global Ecosystem Dynamics Investigation)和高分7號發(fā)射,以及專門用于森林觀測的陸地生態(tài)系統(tǒng)碳監(jiān)測衛(wèi)星(陸地碳衛(wèi)星)預(yù)計(jì)于2022年發(fā)射,彌補(bǔ)激光雷達(dá)數(shù)據(jù)的不足。陸地碳衛(wèi)星搭載了5個波長為1 064nm的激光器,能夠在大光斑足印處獲取高精度森林高度。

為發(fā)展和驗(yàn)證陸地碳衛(wèi)星光斑尺度森林高度和生物量模型,國家林業(yè)和草原局調(diào)查規(guī)劃設(shè)計(jì)院按照陸地碳衛(wèi)星指標(biāo)研制的全國林草資源調(diào)查機(jī)載大光斑激光雷達(dá)系統(tǒng)(National Forest and Grassland Inventory Airborne Large-footprint LiDAR,NFGI-LIDAR-L),可以在2 500m航高獲取水平方向直徑25m、垂直方向間隔0.15m的激光雷達(dá)波形數(shù)據(jù)。由于森林的空間異質(zhì)性高,激光雷達(dá)定位精度直接影響了波形與樹高調(diào)查數(shù)據(jù)的匹配精度,進(jìn)而影響到大光斑激光雷達(dá)估算森林高度模型的精度。因此,在大光斑激光雷達(dá)應(yīng)用到森林資源調(diào)查之前,需要實(shí)現(xiàn)高精度幾何檢校。

大光斑激光雷達(dá)幾何檢校方法包括地面探測器法、機(jī)載紅外相機(jī)成像法、角棱鏡反射法、平坦地形檢校法、傾斜地形檢校法、森林高度匹配法等[8-11]。目前來看,地面探測器法的精度高、適用范圍廣、使用較多,已用于GLAS、資源三號02星、高分7號檢校。星載大光斑激光雷達(dá)檢校水平精度約為10m、垂直精度約為1m[9,11],這對于直徑20～30m、平均樹高10m的森林調(diào)查樣地來說誤差偏高。因此,高精度的幾何檢校方法對陸地碳衛(wèi)星數(shù)據(jù)能否用于林業(yè)至關(guān)重要。

1 機(jī)載大光斑激光雷達(dá)測樹原理

機(jī)載大光斑激光雷達(dá)以一定高度飛行,激光器發(fā)射激光脈沖,并記錄發(fā)射波形能量和時間；激光脈沖到達(dá)地面并返回,記錄接收波形能量和時間。接收波形包括了植被、地面、背景等信息的混合波形。其中:背景信息通過濾波去除；植被和地面波形及其高程通過波形分解提取,并計(jì)算森林高度(圖1)。按照探測器采樣頻率1GHz,即采樣間隔1ns計(jì)算,以光速為X299 792 458m/s換算為測距分辨率,其值為0.149896m/ns≈299792458/(2×109)m/ns,這樣可將波形時間換算為距離。

圖1 機(jī)載大光斑激光雷達(dá)接收波形特征Fig.1 Characteristics of received waveform for airborne large-footprint LiDAR

2 激光雷達(dá)掛飛檢校方法

2.1 激光雷達(dá)檢校試驗(yàn)

大光斑激光雷達(dá)檢校場地位于海南省儋州市西慶機(jī)場(北緯19.55°,東經(jīng)109.43°)。機(jī)場周邊以三葉橡膠(Heveabrasiliensis)人工林為主,零散分布有居住地、道路、農(nóng)地等。機(jī)場周邊地勢平坦,機(jī)場內(nèi)以草地、跑道和裸地為主。在機(jī)場內(nèi)平坦空地上,選取以砂石為主的地面,并清理了地表較高的雜草和磚塊,用于布設(shè)地面激光探測器(圖2)。機(jī)載大光斑激光雷達(dá)數(shù)據(jù)采集時間為2020年12月2日16:00—17:30,云層高度為2 400m。用Cessna 208飛機(jī)搭載NFGI-LIDAR-L,以2 000m和1 500m兩個相對航高、240km/h航速,采用交叉航線飛行,共獲取6條航線數(shù)據(jù)。

圖2 全國林草調(diào)查機(jī)載大光斑激光雷達(dá)檢校場和地面激光探測器布設(shè)點(diǎn)位Fig.2 Calibration site of NFGI-LIDAR-L and distribution of ground Laser detectors

NFGI-LIDAR-L主要由大光斑激光雷達(dá)系統(tǒng)、位置和姿態(tài)系統(tǒng)、穩(wěn)定平臺系統(tǒng)、飛行管理系統(tǒng)、航空攝影系統(tǒng)、數(shù)據(jù)后處理系統(tǒng)組成。相較于上一套機(jī)載大光斑激光雷達(dá)系統(tǒng)[12],NFGI-LIDAR-L的激光器發(fā)射頻率由40Hz改為40～200Hz可調(diào),脈沖寬度由1.57ns改為3.037ns,掃描角由0°改為±15°,激光發(fā)散角由5mrad改為10mrad,波形采樣間隔均為1ns(表1)。為保障高精度快速解算數(shù)據(jù),在機(jī)場架設(shè)了一臺采用實(shí)時動態(tài)載波相位差分定位技術(shù)(Real-Time Kinematic,RTK)的基站。

表1 全國林草資源調(diào)查機(jī)載大光斑激光雷達(dá) 地面激光探測器參數(shù)及航飛和檢校場情況Tab.1 Parameters of NFGI-LIDAR-L,ground Laser detector,flight and calibration site

用RTK流動站從基站引點(diǎn)到檢校場,按照5m×5m間隔,順序布設(shè)了100個直徑11cm的地面激光探測器(圖2)。在航飛前后,用全站儀多次測定探測器的位置,精確到1mm。地面激光探測器專門用于NFGI-LIDAR-L的幾何標(biāo)定研究,探測激光的中心波長為1 064nm,頻率高于200Hz,記錄接收到激光的能量和時間。探測激光能量密度的范圍為0.05～1.00 nJ/cm2,按照能量密度線性劃分4 095個等級。探測器接收GPS時間,并與之時間同步；探測器之間時間同步誤差為6μs。

航飛結(jié)束后,篩選探測器獲得的完整大光斑激光雷達(dá)足印,共8個。其中:2 000m航高由18個、16個和15個探測器探測到的足印分別為3個、2個和1個,共6個；1 500m航高由12個探測器探測到的足印2個。

2.2 基于激光能量分布特征計(jì)算大光斑中心位置方法

由于地面激光探測器間隔5m,用探測器難以直接捕獲到激光雷達(dá)的中心位置。本文基于激光能量自中心向邊緣逐漸降低且能量梯度成比例分布的特征,計(jì)算地面光斑理論中心和能量密度。假定在相對航高2 000m時,同一時間(1ms內(nèi))激光雷達(dá)探測器aij(i和j分別為行號和列號)的數(shù)值最大,計(jì)算以aij為中心的8個方向的能量梯度G(nJ/(m2·m)),并以能量梯度和探測器間距計(jì)算大光斑激光雷達(dá)中心能量密度(Wmax)和坐標(biāo)(Xd,Yd,Zd)。通過實(shí)驗(yàn)室測定得到的激光能量密度W(nJ/m2)與探測器數(shù)值a之間的關(guān)系方程為:W=2.45638a+124.8。

(1)

以探測器a44(i=4,j=4)測得的能量密度最大(W44；i=4,j=4)為例((1)式),x方向最小的能量梯度為G44-43=(W44-W43)/5,y方向最小的能量梯度為G44-34=(W44-W34)/5。假定激光雷達(dá)的能量密度梯度在一個方向上是等比例變化,則

(2)

(2)式中:dx和dy分別為探測器坐標(biāo)x和y方向W44到Wmax的距離。首先,通過解方程得到各參數(shù)值,并通過坐標(biāo)變換得到激光雷達(dá)中心坐標(biāo)；然后,通過克里金插值法得到大光斑激光雷達(dá)水平方向的能量分布。

2.3 機(jī)載大光斑激光雷達(dá)幾何檢校方法

大光斑激光雷達(dá)通過發(fā)射波形和接收波形的時間差測定載荷平臺到地面足印的距離(d),根據(jù)平臺的位置和姿態(tài)、激光發(fā)射角等計(jì)算地面足印的位置((3)式)。

(3)

(3)式中:

1)[XsYsZs]′為激光雷達(dá)地面點(diǎn)在WGS84地固直角坐標(biāo)系的坐標(biāo)。

2)[XgYgZg]′為Inertial Explorer解算后定位系統(tǒng)(GNSS/IMU)在WGS84地固直角坐標(biāo)系的坐標(biāo)。

3)C為載荷平臺坐標(biāo)系向當(dāng)?shù)厮阶鴺?biāo)系的轉(zhuǎn)換矩陣,計(jì)算公式如(4)式所示。(4)式中:L和B分別為平臺在WGS84地固直角坐標(biāo)系的經(jīng)度和緯度(rad)。

4)R為載荷平臺坐標(biāo)系向當(dāng)?shù)厮阶鴺?biāo)系的轉(zhuǎn)換矩陣((5)式),涉及到由右手坐標(biāo)系轉(zhuǎn)為左手坐標(biāo)系和坐標(biāo)旋轉(zhuǎn)；當(dāng)飛機(jī)從南往北飛時,飛機(jī)坐標(biāo)系跟地固坐標(biāo)系的投影坐標(biāo)系相一致。h,p,r分別為航向角、俯仰角、翻滾角。

(4)

(5)

5)CI-G為IMU坐標(biāo)系到載荷平臺坐標(biāo)系的旋轉(zhuǎn)矩陣。

6)[xyz]′為激光雷達(dá)出光點(diǎn)與載荷平臺原點(diǎn)在載荷平臺坐標(biāo)系中的位置偏移量(m)。

7)d為激光雷達(dá)測距(m),由激光雷達(dá)數(shù)據(jù)采集系統(tǒng)獲得。

8)CL-I為激光雷達(dá)坐標(biāo)系轉(zhuǎn)換為IMU坐標(biāo)系的轉(zhuǎn)換矩陣。本文采用了二參數(shù)法((6)式)和三參數(shù)法兩種轉(zhuǎn)換矩陣。

二參數(shù)法指通過確定激光雷達(dá)出光軸在激光雷達(dá)坐標(biāo)系YOZ面投影與Z軸正方向的夾角(α)和激光雷達(dá)出光軸與激光雷達(dá)坐標(biāo)系YOZ面投影的夾角(β)對設(shè)備進(jìn)行幾何檢校[11]。

(6)

(6)式中:α=θ+ρ。θ由激光雷達(dá)測定的碼盤讀數(shù)(Ag)換算得到,θ=(Ag-93594)×360/163840(°)；ρ為設(shè)備左右安裝角,即激光出光軸在激光雷達(dá)坐標(biāo)系YOZ面投影與碼盤角差值(°)。

三參數(shù)法指通過確定激光雷達(dá)坐標(biāo)系與慣導(dǎo)坐標(biāo)系X,Y,Z三個方向的安置角(h0,p0,γ0)對設(shè)備進(jìn)行幾何檢校((7)式)。

由(3)式得到的Xs,Ys,Zs,需要從WGS84地固直角坐標(biāo)系通過國際地球參考框架(ITRF)的7參數(shù)換算到CGCS2000地固直角坐標(biāo)系[13-14],再換算到CGCS2000大地坐標(biāo)系,最后換算為CGCS2000高斯投影3度帶[15](中央子午線108°E,記為X2000,Y2000,Z2000)。

(7)

由于NFGI-LIDAR-L為±15°掃擺式200Hz激光雷達(dá),所以在做幾何檢校之前,需要先使探測器數(shù)據(jù)與機(jī)載數(shù)據(jù)匹配。本文通過波形特征、碼盤角及試算后檢校參數(shù)是否能夠適用到其他大光斑激光雷達(dá)點(diǎn)位等方法確定地面激光探測器探測到的大光斑激光雷達(dá)。由8個檢校點(diǎn)得到探測器記錄脈沖時間比NFGI-LIDAR-L記錄時間晚1004±2μs。

以遍歷法確定檢校參數(shù)ρ和β,使得由激光發(fā)射計(jì)算得到的大光斑中心位置(X2000,Y2000,Z2000)與由探測器計(jì)算的理論大光斑中心位置(Xd,Yd,Zd)的殘差(Δp)最小((8)式)。每個大光斑對應(yīng)一套檢校參數(shù)。

(8)

以遍歷法確定檢校參數(shù)ρmin和βmin,以及h0min,p0min,γ0min,使獲得的所有大光斑位置殘差的均方根差RMSE(εmin)最小((9)式),則認(rèn)為該參數(shù)是最優(yōu)檢校參數(shù)。

(9)

本試驗(yàn),n=8；最優(yōu)參數(shù)計(jì)算的大光斑中心位置記為X2000_min,Y2000_min,Z2000_min。

2.4 精度分析和驗(yàn)證

由單個檢校點(diǎn)最優(yōu)參數(shù)得到的大光斑中心位置(X2000_i,Y2000_i,Z2000_i)與最優(yōu)參數(shù)計(jì)算的大光斑中心位置(X2000_min,Y2000_min,Z2000_min)差(Δpi),可以視為探測器單次檢校的大光斑中心位置與大光斑理論中心位置的誤差((10)式)。

因檢校參數(shù)導(dǎo)致的誤差,可以通過比較三參數(shù)法與二參數(shù)法計(jì)算大光斑位置的均方根誤差(RMSE)來評估計(jì)算精度((11)式)。

(10)

(11)

絕對高程精度的驗(yàn)證:2021年2月3日,利用RTK測量了12個高程精度驗(yàn)證點(diǎn)(Zki),測量精度為2cm,用于評價(jià)大光斑激光雷達(dá)絕對高程精度(ΔZ)。評價(jià)指標(biāo)為均方根誤差((12)式)、誤差平均值和標(biāo)準(zhǔn)差。

(12)

相對高程精度的驗(yàn)證:采用與由機(jī)載小光斑激光雷達(dá)生產(chǎn)的高精度數(shù)字高程模型(DEM)進(jìn)行比較的方法。2020年5月31日,通過搭載的Riegl VQ-1560i,獲取了西慶機(jī)場及周邊的小光斑激光雷達(dá)數(shù)據(jù),點(diǎn)云密度10個點(diǎn)/m2。經(jīng)過數(shù)據(jù)解算、航帶匹配、去噪和精分類,以及提取地面點(diǎn)等處理,生成了1m分辨率DEM。通過地面控制點(diǎn)對DEM數(shù)據(jù)進(jìn)行修正,得到高精度DEM數(shù)據(jù)。在試驗(yàn)區(qū)選擇落在裸地、草地、森林、農(nóng)田的大光斑激光雷達(dá),用其地面回波的高程與DEM進(jìn)行比較,評估其相對高程精度。

3 結(jié)果

3.1 大光斑中心位置和能量密度

與通過直接插值的結(jié)果相比(圖3),根據(jù)激光能量分布特征解算的大光斑理論中心位置和能量密度分布結(jié)果(圖4),大光斑能量密度分布更完整,跟發(fā)射波形(圖5)的能量分布更一致,中心位置也更明確。

圖3 探測器測定的大光斑能量密度用克里金插值法Fig.3 Energy density from ground Laser detector and interpolated by Kriging

圖4 根據(jù)激光能量空間分布解算后大光斑理論中心和能量密度分布Fig.4 The geocentric and energy density of large-footprint LiDAR calculated based on spatial pattern of laser from ground Laser detector

圖5 NSGI-LIDAR-L記錄的發(fā)射波形Fig.5 Transmitted waveform recorded by NSGI-LIDAR-L

在2 000m航高處,大光斑足印中心相較最大觀測能量位置移動了1.27～3.47m,平均移動2.23m,移動方向沒有明確規(guī)律；在1 500m航高處,2個大光斑足印中心相較最大觀測能量位置分別移動1.36m和3.04m。在1 500m和2 000m航高,大光斑足印中心能量密度為2 197～2 221nJ/m2和1 554～2 060nJ/m2,平均為2 209nJ/m2和1 783nJ/m2。

3.2 激光雷達(dá)檢校參數(shù)

激光雷達(dá)系統(tǒng)檢校前(ρ=0,β=0),x,y,z方向的位置差分別為60.51,52.76,1.88m,總RMSE為80.30m。單個檢校點(diǎn)使大光斑幾何精度提高2個數(shù)量級,使總位置差降低到0.85～1.41m之間,平均為1.14m。其中,x,y,z方向的RMSE分別為0.82,0.78,0.05m。8個檢校點(diǎn)得到的最優(yōu)檢校參數(shù)為ρ=2.2345°,β=-0.7732°；對應(yīng)x,y,z方向的位置差分別為0.61,0.54,0.05m,總RMSE為0.82m(表2)。8個檢校點(diǎn)得到的平均檢校參數(shù)已經(jīng)接近最優(yōu)參數(shù),說明通過8個檢校點(diǎn)，能夠獲得比較精確的大光斑激光雷達(dá)中心位置。

表2 激光雷達(dá)系統(tǒng)二參數(shù)檢校法的參數(shù)和檢校前后位置差Tab.2 Constants of Two-parameter Calibration,and x,y,z and total RMSE pre-and post-calibration

從三參數(shù)檢校結(jié)果看(表3),最優(yōu)檢校參數(shù)為h0=0.02,p0=0.7734,r0=2.231；x,y,z方向的RMSE分別為0.62,0.54,0.05m,總RMSE為0.82m。跟二參數(shù)檢校精度基本一致。單檢校點(diǎn)的位置差范圍,x,y,z方向分別為0.64～1.12m,0.56～0.93m,0.05～0.06m,RMSE分別為0.82,0.78,0.05m。

表3 激光雷達(dá)系統(tǒng)三參數(shù)檢校和檢校前后位置差Tab.3 Constants of Three-parameter Calibration,and x,y,z and total RMSE pre- and post-calibration

從三參數(shù)法與二參數(shù)法最優(yōu)參數(shù)計(jì)算的大光斑中心位置差(表4)來看,x,y,z方向的平均位置差分別為0.083,0.094,0.002m,總平均位置差為0.125m。從影響因素來看,兩種方法受到航高和碼盤角影響。航高越高,誤差越大；1 500m和2 000m航高的平均位置差0.101m和0.134m,并且主要集中在水平誤差上。同一航高時,碼盤角越大位置差越大。比如2 000m航高,碼盤角為-0.74°,-1.67°,-2.38°,-3.48°時,位置差分別為0.129,0.132,0.134,0.137m。

表4 三參數(shù)與二參數(shù)法最優(yōu)參數(shù)計(jì)算的大光斑中心位置差 相對航高和碼盤角Tab.4 RMSE between Three-parameter Calibration and Two-parameter Calibration,relative flight heights and angles of code wheel

3.3 激光雷達(dá)高程檢校

以RTK檢校點(diǎn)評估大光斑激光雷達(dá)絕對高程精度,高程差的變異范圍為0.01～0.18m,均方根誤差為0.12m,平均誤差為0.10m,標(biāo)準(zhǔn)差為0.06m,檢校點(diǎn)12個。大光斑激光雷達(dá)與檢校點(diǎn)的高程差未出現(xiàn)隨距離或激光掃描角度的增大而增大的現(xiàn)象。因此,NFGI-LIDAR-L檢校后絕對高程精度在0.15m以內(nèi)。

與小光斑激光雷達(dá)生成的高精度DEM相比,大光斑激光雷達(dá)得到的裸地和跑道、草地、農(nóng)田、森林等地面高程均方根誤差分別為0.13,0.15,0.38,0.58m,樣本量分別為8,10,3,8。

4 討論

4.1 激光雷達(dá)檢校試驗(yàn)

2017年采用地形法得到的機(jī)載大光斑激光雷達(dá)幾何檢校的精度約為5m。衛(wèi)星的地面激光探測器間距為10～25m,水平方向檢校精度在10m左右[9,11],說明探測器間隔對水平檢校精度影響較大。本文基于激光能量梯度計(jì)算大光斑中心位置的方法,使檢校精度突破了0.5倍探測器間距限制,達(dá)到0.82m。

由于本次檢校飛行的掃描角比較小(碼盤角在-0.74°～-3.48°),采用三參數(shù)和二參數(shù)檢校法的差異在0.1m級別,說明當(dāng)掃描角小時(±1.5°),兩種方法都能夠得到比較精確的位置。隨著掃描角的增加,三參數(shù)法可能會比二參數(shù)法檢校精度高。相對航高也影響檢校參數(shù)的精度,未來需要進(jìn)一步通過試驗(yàn)驗(yàn)證。

4.2 激光雷達(dá)檢校試驗(yàn)改進(jìn)

由于地面探測器靈敏度高,可以支持更大角度機(jī)載大光斑激光雷達(dá)幾何檢校。未來可以通過增加探測器數(shù)量、提高探測器間距、變換不同航高(1 000m～3 000m)、增大掃擺角度(±5°,±10°)等試驗(yàn),進(jìn)一步獲得更準(zhǔn)確的大光斑激光雷達(dá)能量分布和中心位置,以提高大角度激光雷達(dá)的檢校精度。

4.3 星載激光雷達(dá)檢校展望

本文驗(yàn)證了采用地面探測器檢校機(jī)載大光斑激光雷達(dá)的可行性。隨著激光器數(shù)量增加,地面探測器檢校法的工作量將大幅增加。如何對陸地生態(tài)碳衛(wèi)星搭載的5個激光雷達(dá)進(jìn)行幾何檢校,需要繼續(xù)探討。比如,探討改變衛(wèi)星姿態(tài)角,將衛(wèi)星平臺旋轉(zhuǎn)90°,使5個激光器由垂軌分布改為沿軌分布,5個激光器發(fā)射的大光斑依次通過同一個地面探測器場,獲得5個激光雷達(dá)相對位置；再將衛(wèi)星平臺轉(zhuǎn)-90°,使激光器正常沿軌運(yùn)行,之后多次觀測其中1個激光器的光斑位置,從而實(shí)現(xiàn)5個激光器的精確檢校。

5 結(jié)論

針對大光斑激光雷達(dá)檢校精度不能滿足林業(yè)應(yīng)用的問題,本研究通過新研制的NFGI-LIDAR-L和地面激光探測器,用編寫的大光斑激光雷達(dá)檢校算法和基于激光空間分布特征計(jì)算大光斑激光雷達(dá)中心位置算法,實(shí)現(xiàn)了NFGI-LIDAR-L的幾何檢校。大光斑激光雷達(dá)幾何精度從檢校前的80.3m提高到0.82m,絕對高程精度達(dá)到0.12m。因此,通過地面激光探測器法對大光斑激光雷達(dá)進(jìn)行幾何檢校的方法能夠基本達(dá)到林業(yè)應(yīng)用要求。

致謝:

野外工作期間,西慶機(jī)場提供試驗(yàn)場地,上海大恒精密機(jī)械有限公司吳致昊工程師和曹健東工程師獲取飛行數(shù)據(jù),山東科技大學(xué)研究生徐飛、劉冠杰、林景峰、宋成航參加檢校場布設(shè)。他們?yōu)楸狙芯刻峁┝舜罅χС趾椭T多幫助,在此一并表示感謝。