李榮華,王振宇,陳 鳳,肖余之,薛豪鵬
(1. 大連交通大學機械工程學院,大連 116028;2. 上海宇航系統(tǒng)工程研究所,上海 201109)
太空垃圾又稱空間碎片或軌道碎片,是指圍繞地球軌道的無用人造物體。太空垃圾小到由人造衛(wèi)星碎片、漆片、粉塵,大到整個飛船殘骸構(gòu)成。其中大型太空垃圾(衛(wèi)星由于故障或任務結(jié)束而被放棄)不但占用其原在的有限軌道資源,還對其他航天器的安全造成不可估計的危害。為解決這一問題,各國陸續(xù)開展了對失效衛(wèi)星的環(huán)繞檢測、運動跟蹤、交會對接、在軌維修等操作[1-4]。而實現(xiàn)這些目的,就必需獲取非合作目標的三維空間信息,進行三維模型重建。激光雷達由于其受光照和距離等因素影響小,抗干擾能力強等優(yōu)點,是目前監(jiān)測空間目標較為理想的工具,常見的成像方式為線陣式與面陣式[5]。
面陣式激光雷達[6]采用多束密集的激光束直接向各個方向漫射,最大優(yōu)勢是能夠快速記錄整個場景,避免了掃描過程中目標或激光雷達移動帶來的各種問題。但對遠距離空間運動,面陣式成像分辨率低,對目標的測量精度仍然不能滿足測量要求。線陣式激光雷達在空間作業(yè)中,相較于面陣式負載較小,成像分辨率高,獲取的線陣數(shù)據(jù)同樣能夠用于運動參數(shù)估計及三維模型重建。對于靜態(tài)目標,已知目標與雷達相對運動信息,運用線陣式激光雷達獲取不同時刻線陣數(shù)據(jù),可直接組合線陣信息得到測量目標的三維模型;但對于運動目標,目標位姿一直處于變化狀態(tài),不同時刻獲取的是當前時刻目標掃描位置線陣數(shù)據(jù),掃描數(shù)據(jù)存在畸變,所測量的運動目標輪廓與目標的實際外形相比一般會存在形變,直接組合線陣信息并不能準確地對被測目標進行三維重建。對于非平衡(失穩(wěn))狀態(tài)下,蔡銀橋等[7]對星載三維激光雷達圖像畸變進行研究,研究了在不同掃描周期內(nèi)三維激光雷達圖像變形,得出了控制成像時間是降低圖像畸變的關(guān)鍵因素。針對線陣畸變對目標方位估計性能的影響,李建等[8]通過理論、仿真以及試驗數(shù)據(jù)處理等方面進行分析,并提出了一種初步的陣形校正方法。
解決運動目標線陣掃描中存在的畸變問題,主要分為合作目標與非合作目標兩種情況下進行分析[9-10]。對于合作目標,即目標運動參數(shù)已知或目標可添加標記,趙竹新等[11-12]給出了一種利用線陣相機立體交會測量彈丸設計,根據(jù)彈丸目標的外形特點設計了一種新的標記線方法,避免了線陣相機掃描速度不足而難以實現(xiàn)彈道同步測量的難題。高俊釵等[13]根據(jù)線陣相機圖像采集特點,以標準圓形作為標定目標確定了物像關(guān)系和校正插值系數(shù),并設計線性插值算法進而能夠?qū)θ我庑螤畹哪繕诉M行成像校正。在合作模式下,被測目標運動估計精度通常能夠達到要求,但是對于無法進行標定的測量目標或者無法增加輔助設備的測量情況并不適用。對于非合作目標,王盈等[14]提出了一種適用于復雜空間目標的在軌激光雷達成像仿真流程,并強調(diào)了面積法在復雜目標成像可見判斷方面的適用性,適用于大表面目標測量。本課題組前期通過相鄰兩次掃描的畸變點云配準來獲取測量目標的瞬時運動參數(shù),該方法僅使用單一載荷實現(xiàn)了非合作空間目標相對位姿的測量,但是僅能估算自旋狀態(tài)的空間失穩(wěn)目標,對于復雜運動的空間失穩(wěn)目標或者兩次獲取點云重合量達不到配準要求的情況有待進一步研究。
綜上所述,本研究以實現(xiàn)空間失穩(wěn)目標三維重建為目的,針對空間失穩(wěn)目標,因其非合作性、形狀多樣性、運動狀態(tài)不確定性,且線陣掃描式激光雷達測量動態(tài)目標發(fā)生畸變而無法精確恢復目標三維模型問題,首先對掃描畸變從測量角精度與測距精度兩大判定條件進行分析,然后根據(jù)測量工況對目標點云進行帶畸變配準。該方法不僅能夠有效解決空間失穩(wěn)目標線陣成像畸變問題,而且可以實現(xiàn)較高精度的動態(tài)非合作目標三維重建,為線陣式測量技術(shù)進一步發(fā)展提供新的研究方向。
線陣掃描型激光雷達在探測失穩(wěn)目標過程中,會發(fā)生時空非定常點云數(shù)據(jù)畸變,所掃描的運動目標數(shù)據(jù)與目標的實際模型相比會存在差異。這種差異是由于線陣雷達的掃描等效速度與目標運動速度不一致所導致的。所謂的線陣雷達的掃描等效速度是指線陣雷達單個像素的寬度與掃描頻率的乘積,比如傳感器的單個像素寬度為0.001 mm,掃描頻率為100k 線/s,則掃描等效速度為1 m/s。當雷達的掃描等效速度與目標的運動速度相等時,目標的線陣成像與目標的實際模型一致;當雷達的掃描等效速度大于目標的實際運動速度時,目標的線陣成像與其實際目標模型相比顯得被拉伸了,相反,當雷達的掃描等效速度小于目標的運動速度時,目標的線陣成像與其實際目標模型相比顯得被壓縮了。運動目標線陣成像的這種特點,從成像效果看,是一種成像失真。如圖1所示,假設目標總共需要線陣雷達單程掃描8次,奇數(shù)次雷達掃描方向與目標運動方向一致,雷達的掃描等效速度小于目標的實際運動速度,目標表面同向壓縮;偶數(shù)次雷達掃描方向與目標運動方向相反,雷達的掃描等效速度大于目標的實際運動速度,目標表面反向拉伸。最終導致線陣雷達掃描運動目標后產(chǎn)生不容忽視的畸變誤差。
圖1 目標表面掃描畸變產(chǎn)生示意圖
因為線陣掃描雷達受掃描棱鏡加工精度影響,角分辨率大小會出現(xiàn)周期性變化,所以即使目標是靜止狀態(tài),線陣雷達測量的點云也會出現(xiàn)畸變。由此可見,當目標受運動影響出現(xiàn)的畸變大小如果不大于雷達自身測量精度,則恢復畸變的過程將失去意義;反之,如果運動畸變大于雷達測量畸變,則有理由認定為目標失穩(wěn),需要畸變補償。此外,由于受雷達角分辨率精度影響的畸變程度無法通過疊加來判斷,而且考慮到為了提高畸變恢復效果及重建精度,以雷達測角精度為參考的失穩(wěn)臨界值向下降低一個數(shù)量級。
以雷達測量精度為衡量基準計算目標畸變?nèi)萑潭扰R界值,即當目標由于運動導致的數(shù)據(jù)畸變程度超過了雷達自身測量精度的10%(降低一個數(shù)量級)時的目標運動速率認為是畸變?nèi)萑潭扰R界速率。現(xiàn)從雷達測角精度與測距精度兩方面對畸變?nèi)萑潭扰R界值進行判定。
1.2.1角精度判定條件
當畸變?nèi)萑探嵌葹閞時,根據(jù)雷達指標:雷達環(huán)繞測量分辨率為m×n=400×400,線掃描寬度為Ls=50,單方向掃描次數(shù)為m/Ls=8,單次掃描時間為t=100 ms=0.1 s??梢酝扑愠鱿噜忺c掃描所花時間t′=3.125×10-5s:
(1)
重建所能容忍的相鄰點間畸變角度臨界值r為激光雷達測角精度(30″),即r=30″,則目標失穩(wěn)臨界自旋速率ωr=26.6667(°)/s:
(2)
目標畸變?nèi)萑潭扰R界自旋速率為ω=26.6667(°)/s;即當目標自旋速率大于該值時,即使忽略雷達測角精度,章動的復合運動的角重建誤差也會超過指標。
1.2.2測距精度判定條件
假定已知信息:目標自旋角速率為ωm=12(°)/s,章動角速率為ωz=2.5(°)/s,雷達視場角為20°×20°,目標與雷達工作距離為L=80 m,成像分辨率為m×n=400×400,目標主體參數(shù)為2.000 m×3.500 m×2.500 m,雷達測距精度為dm=0.05 m。目標章動方向和自旋運動方向在某一時刻同向時,目標表面一點達到最大的運動角速率ωmax=14.5(°)/s:
ωmax=ωm+ωz
(3)
考慮最大畸變情況,對下面兩種情況分別計算。
1)目標運動方向與掃描方向平行時
目標運動方向與掃描方向平行時,雷達掃描運動示意圖如圖2所示,根據(jù)雷達掃描距離L=80 m,雷達視場角為20°×20°,可得X軸單方向雷達最大掃描距離為LX=28.212 m(Y軸方向同理):
圖2 平行時雷達掃描運動示意圖
(4)
X軸單方向掃描次數(shù)為m/Ls=8,則單次掃描距離L′X=3.527 m:
L′X=LX/(m/Ls)
(5)
目標主體參數(shù)中寬度3.500 m與高度2.500 m均小于3.527 m,因此雷達掃面目標主體掃描一次或兩次。當雷達掃面目標主體掃描兩次時,用時最長,且畸變最大。最大用時tmax=0.014 s:
tmax=(1+3.500/28.212)t/8
(6)
圖3 目標運動半徑
如圖3所示,目標運動半徑r為:
(7)
圖4 目標表面點畸變量示意圖
圖4中,目標表面P點為真實位置,P′為畸變位置?;兞縇dis=0.0071 m:
Ldis=2rsin(tmaxωmax/2)
(8)
計算顯示:在畸變最大的情況下,表面一點P的畸變量為0.0071 m,而雷達測距精度小于等于0.05 m,因此畸變可忽略。
2)當目標運動方向與掃描方向垂直時
當目標運動方向與掃描方向垂直時雷達掃描運動示意圖如圖5所示。
圖5 垂直時雷達掃描運動示意圖
同理當目標運動方向與掃描方向平行時,可得目標運動方向與掃描方向垂直狀態(tài)畸變量Ldis=0.0071 m:
Ldis=2rsin(tmaxωmax/2)
(9)
計算顯示:在畸變最大的情況下,表面一點P的畸變量為0.0071 m,而雷達測距精度小于等于0.05 m,因此畸變可忽略。
3)根據(jù)測距精度計算最大運動速度容忍度ωr
根據(jù)測距精度0.05 m,畸變影響在10%獲得失穩(wěn)速率臨界:
ωr=2arcsin(0.05×10%/2/r)/tmax=10.15(°)/s
(10)
若章動為2.5(°)/s,自旋速率最大容忍度為7.65(°)/s。
綜合兩種臨界值判定條件:根據(jù)角精度判定條件,目標畸變?nèi)萑潭扰R界自旋速率ω=26.6667(°)/s;根據(jù)測距精度判定條件,當目標運動方向與掃描方向平行與垂直兩種狀態(tài)下,畸變量均低于測距精度一個數(shù)量級,且根據(jù)測距精度,目標畸變?nèi)萑潭扰R界自旋速率ωr=10.15(°)/s。在本測量系統(tǒng)已知數(shù)據(jù)下,結(jié)合費馬引理,目標畸變?nèi)萑潭扰R界自旋速率取極大值,即ωr=26.6667(°)/s。當目標自旋速率大于等于ωr時,初次配準時需要對點云數(shù)據(jù)進行畸變恢復;當目標自旋速率小于ωr時,初次配準時對點云的畸變恢復將失去意義。
目標坐標系坐標原點設為待測目標幾何中心位置,ZO軸方向設為待測目標自旋軸方向,YO軸方向設為太陽翼中軸線方向,基于右手法則確定XO軸方向。雷達坐標系以雷達中心為坐標原點,環(huán)繞檢測半徑設為80 m,三軸方向與初始獲取目標反饋信息的方向一致。世界坐標系與雷達坐標系重合,掃描工況笛卡爾坐標系如圖6所示。
圖6 掃描工況坐標系
點云配準是通過匹配具有重疊部分的數(shù)據(jù)集,將不同視場下得到的三維數(shù)據(jù)變換到同一坐標系下,得到信息更加豐富完整的目標三維點云數(shù)據(jù),計算出點云間的歐拉變換關(guān)系,即旋轉(zhuǎn)變換矩陣R和平移向量T。
現(xiàn)采用線陣掃描型激光雷達進行N次掃描,得到N個視場點云數(shù)據(jù),假設兩組相鄰點云數(shù)據(jù)Ni和Ni+1是對同一目標在不同成像視角得到的,配準目的便是通過平移和旋轉(zhuǎn)變換將Ni+1變換到Ni坐標系下,定義為N′i+1:
N′i+1=R(α,β,γ)×Ni+1+T(tX,tY,tZ)
(11)
其中,R(α,β,γ)為旋轉(zhuǎn)矩陣;T(tX,tY,tZ)為平移矢量;α,β,γ分別為繞X軸、Y軸、Z軸逆時針旋轉(zhuǎn)的歐拉角;tX,tY,tZ為三個坐標軸正方向的平移量。旋轉(zhuǎn)矩陣R(α,β,γ)可以表示為:
R(α,β,γ)=
(12)
通過選取合適的變換角度和平移量,將Ni+1變換到Ni坐標系下,進行坐標系的統(tǒng)一處理,最后得到的N′i+1與Ni即為所需要的配準結(jié)果。
1)統(tǒng)計配準參數(shù)與對應時間數(shù)組
通過上一步配準,獲取目標相鄰掃描時間內(nèi)運動的旋轉(zhuǎn)矩陣與平移向量,以初始時刻目標所在位置為位姿起始點,對應位姿增量疊加可以獲得目標位姿-時間函數(shù)。通過配準獲得參數(shù)如式(13)所示:
(13)
2)計算第i次旋轉(zhuǎn)平移運動增量
旋轉(zhuǎn)平移矩陣形式為(旋轉(zhuǎn)次序依次為:Z-Y-X):
(14)
(15)
旋轉(zhuǎn)運動增量:
(16)
平移運動增量:
(17)
式中:sX與cX代表繞X軸旋轉(zhuǎn)的正弦與余弦值;Y軸、Z軸同理。
實現(xiàn)點云精配準現(xiàn)今主流的核心算法是ICP算法,該算法精度很大程度上取決于預先給定的位置姿態(tài)初值,即需要先進行點云粗配準,航天器一般具有較為明顯的平面幾何特征,現(xiàn)采用RANSAC算法提取目標點云特征點進行粗配準[15]。
1)點云粗配準
設由線陣雷達獲取的點云數(shù)據(jù)集為P,從中隨機選取三個點:P1(x1,y1,z1);P2(x2,y2,z2);P3(x3,y3,z3)。然后根據(jù)三點位置確定平面S:z=Ax+By+C,參數(shù)A,B,C由式(18)確定:
(18)
統(tǒng)計平面S上的點的個數(shù),設定平面厚度為2σ,計算P中任何一點Pi到平面S的距離di,di由下式計算得到:
(19)
統(tǒng)計di<σ的點的個數(shù),記為S的得分數(shù)。重復上述步驟K次,選擇出得分最高的平面S*,其中K的計算公式為:
(20)
式(20)中m為點云中總的點數(shù),n為點云中特征點點數(shù),由于m和n值都很大,采用近似計算,由下式計算K值:
(21)
式中:ε為位于平面S*之外點所占比例的值(估計值),Ф為經(jīng)過K次采樣之后最后被選中的概率。
如圖7所示,運用RANSAC算法提取相鄰點云平面幾何特征,圖7(a)為相鄰點云數(shù)據(jù)顯示,圖7(b)為兩組點云相似特征點連線,據(jù)此進行點云粗配準。
圖7 RANSAC算法
2)點云精配準
精細配準是一種更深層次的配準方法。經(jīng)過前一步粗配準,得到了變換估計值。將此值作為初始值,在經(jīng)過不斷收斂與迭代的精細配準后,達到更加精準的效果。ICP算法(最鄰近點迭代算法)最早由文獻[16-17]提出,文獻[16]采用的是待配準表面上點的配準方法,文獻[17]采用曲面點集與切平面的配準方式。本文處理的點云數(shù)據(jù)形式為點集,而后者所采用的ICP算法對每個點的法向量進行計算,計算量較大,因此本文主要采用前者的待配準表面上點的配準方法。
算法基本原理:兩點云P與Q,對于點云P中的每個點pi,搜索其在點云Q上的最近點作為與之相對應的點qi,然后依據(jù)對應關(guān)系求解使得式(22)中目標函數(shù)達到最小的剛體變換,即旋轉(zhuǎn)矩陣R和平移向量T,并將該變換作用于源點云,迭代進行這一過程直到滿足某一設定的收斂準則。
(22)
迭代終止條件:當目標相對于真值位姿精度不大于1°(3σ,三軸)時迭代終止。
假定失穩(wěn)衛(wèi)星最大自旋速度約為12(°)/s,定義旋轉(zhuǎn)軸過模型質(zhì)心,方向沿Z軸方向,并且目標在短時間內(nèi)為勻速運動?,F(xiàn)對Y軸方向采用最大自旋速度進行高壓仿真,即三個方向速率分別為:vX=0(°)/s,vY=12(°)/s,vZ=0(°)/s。環(huán)繞目標沿Y軸按照每12度1幀提取點云可視化數(shù)據(jù)信息,如圖8所示。
圖8 多視場點云數(shù)據(jù)獲取
以雷達第一次掃描的初始時刻,目標所在空間位置建立目標坐標系,如果想把每一條畸變線信息恢復到目標形貌真實位置(x0,y0,z0),需要計算該條畸變點云線相對初始時刻所在坐標系的平移量和旋轉(zhuǎn)分量,經(jīng)過前面的步驟已經(jīng)獲取了目標天體的運動參數(shù)模型,可以通過對不同時刻采集的線信息進行逆運算以把多組線信息整合到同一時刻的空間位置,ti代表統(tǒng)計雷達在獲取每一條線信息相對初始時刻所經(jīng)過的掃描時間。則每一條雷達掃描線點云的逆向數(shù)據(jù)恢復模型為:
(23)
式中:
(24)
(25)
通過對多組相鄰畸變點云數(shù)據(jù)的配準,可以獲取沿時間延長的多組有限個相鄰掃描時間內(nèi)的相對位姿增量,如圖9所示。
圖9 目標瞬時運動狀態(tài)
本研究所采用的目標模型由衛(wèi)星主體及太陽能帆板構(gòu)成。在激光雷達環(huán)繞檢測過程中,兩側(cè)的太陽能帆板會因主體及自身遮擋而導致點云關(guān)鍵部位特征突變,從而導致相鄰點云配準時產(chǎn)生極大躍性偏差,對重建結(jié)果造成極大干擾。運用自修正原理,通過去除運動參數(shù)中存在的跳躍點,對環(huán)繞檢測一周數(shù)據(jù)進行擬合,恢復躍性點處數(shù)據(jù),并據(jù)此進行三維重建。
失穩(wěn)目標在雷達相鄰的兩次掃描時間內(nèi)能夠保證有效重合掃描區(qū)域,基于此,對目標相鄰掃描點云進行拼接;從而在目標旋轉(zhuǎn)周期內(nèi),重復掃描、配準過程完成目標表面點云的三維重建。過程如下:1)獲取第j次掃描點云,并作為重建起始,經(jīng)上一步實現(xiàn)點云畸變信息恢復;2)獲取第j+1次掃描點云,同樣恢復畸變信息,與第i次點云進行配準,舍去重合區(qū)域;3)繼續(xù)獲取下一時刻掃描點云,重復上述過程;4)以當前時刻掃描點云與第j次點云配準重合度作為終止重建條件:當兩處點云存在重合區(qū)域或者重合度大于設定閾值,認為目標完成一次周期自轉(zhuǎn),即完成重建工作。依據(jù)配準擬合結(jié)果,將各旋轉(zhuǎn)角度對應結(jié)果還原到目標初始坐標系下,點云數(shù)據(jù)依次疊加,得到目標初始三維重建模型。
為了證明算法的有效性,進而驗證本文提出的帶畸變自修正式三維重建方法的可行性,現(xiàn)應用VS2013結(jié)合點云庫(Point cloud library,PCL)對上述思想逐一進行仿真試驗與校驗。(文中算法結(jié)果均在CPU為Intel(R)Core(TM)i5-3210M CPU,2.50 GHz,內(nèi)存為4 GB的PC機上運行所得,操作系統(tǒng)為Windows 7)
現(xiàn)對目標進行仿真重建,如圖10所示,外形尺寸為16000 mm×6000 mm×3000 mm,本體尺寸為3500 mm×2000 mm×2500 mm。
圖10 目標點云模型
按照每12度1幀,逆時針環(huán)繞目標沿Y軸提取點云數(shù)據(jù),可獲得N個視場點云數(shù)據(jù):
(26)
其中第1組0°為基準初始位置,其后每隔12°取樣一次,最后一組為環(huán)繞檢測360°后得到的可視化點云數(shù)據(jù)。32組目標環(huán)繞點云數(shù)據(jù)組如圖11所示。
圖11 32組目標環(huán)繞點云
按照上面提出的帶畸變自修正式三維重建方法,針對獲取的32組點云數(shù)據(jù)進行校驗。
1)第i組與第i+1組相鄰視角點云數(shù)據(jù)ICP配準。
2)獲取該組旋轉(zhuǎn)矩陣Ri與平移向量Ti。
3)將第i+1組點云根據(jù)旋轉(zhuǎn)矩陣Ri與平移向量Ti,配準恢復到第i組點云狀態(tài)下,相鄰組點云配準如圖12所示。
圖12 相鄰點云配準
4)將32組點云數(shù)據(jù)相鄰兩兩進行配準,獲取目標運動增量、解算目標瞬時運動參數(shù)。
5)然后對目標運動參數(shù)進行自修正式處理,去除特征突變干擾點。
6)依次恢復點云到初始幾何位置,完成初始目標三維重建。
7)通過系統(tǒng)進一步迭代,位姿增量自收斂,直至最終滿足重建精度要求。
將32組點云數(shù)據(jù),按照相鄰組旋轉(zhuǎn)配準,按照每12度1幀取樣,每組旋轉(zhuǎn)角度真值為12°,結(jié)果如表1所示。
表1 相鄰組旋轉(zhuǎn)配準結(jié)果
如圖13所示,實線為32組點云數(shù)據(jù)第一次配準結(jié)果,虛線為真值,點劃線為擬合結(jié)果。其中當旋轉(zhuǎn)角度為72°及252°時,模型單側(cè)帆板在視角可視范圍內(nèi)消失,使得ICP配準時產(chǎn)生明顯偏差。根據(jù)第2節(jié)闡述的方法,剔除兩組較大偏差點,將剩余30組配準數(shù)據(jù)進行擬合,得到擬合結(jié)果。
圖13 相鄰點云第1次配準
依據(jù)配準擬合結(jié)果,將各旋轉(zhuǎn)角度對應結(jié)果還原到目標初始坐標系下,點云數(shù)據(jù)依次疊加,得到目標初始三維重建模型,第1次配準耗時為18337 ms。
本文獲取的目標模型,在進行第6次迭代后,單次耗時為6779 ms,目標單位時間內(nèi)運動增量均位于12°±1°內(nèi);在現(xiàn)有硬件配置環(huán)境下,配準總耗時為72467 ms,符合三維重建技術(shù)指標。第6次配準結(jié)果如表2所示,擬合結(jié)果及目標最終三維重建模型如圖14~圖15所示。
本研究以實現(xiàn)空間失穩(wěn)目標三維重建為目的,解決目標失穩(wěn)及線陣雷達掃描動態(tài)目標等多種因素所導致的目標三維模型無法精確恢復問題。首先根據(jù)角精度與測距精度兩大判定條件計算畸變?nèi)萑潭茸畲笈R界值,其次對目標點云進行自修正式迭代配準,直至滿足三維重建要求。針對本文獲取的目標模型,三維重建精度可達到±1°以內(nèi),重建總耗時在現(xiàn)有配置環(huán)境下為72467 ms。該方法不僅能夠在理論上指導三維目標的點云數(shù)據(jù)配準、非合作目標的三維姿態(tài)測量和目標三維重建等工作,還能在實際應用中為線陣式測量技術(shù)的進一步發(fā)展提供新的研究方向。
表2 相鄰組旋轉(zhuǎn)配準結(jié)果
圖14 相鄰點云第6次配準
圖15 目標最終三維重建模型