諶炳漢, 左崢嶸,*, 夏魯瑞
(1. 華中科技大學自動化學院多譜信息處理技術國家級重點實驗室, 武漢 430074; 2. 航天工程大學太空安全研究中心, 北京 101416)
目標的遙感測量可有效獲取目標運動特性,有助于對目標飛行試驗效果進行評估。目標的軌跡重建是遙感測量系統獲取目標信息的關鍵技術,目標的空間軌跡及運動參數通常需要雙平臺觀察,并采用雙目視覺方法進行重建,這提高了觀測系統的復雜度,例如,天基紅外預警系統的高軌部分包括4顆靜止軌道衛星和2顆大橢圓軌道衛星,低軌部分由21~28顆小型衛星組成[1],其部署和維護成本極大,因此,在實際觀測系統中,常常只進行單平臺的觀測[2]。這對單平臺觀測條件下的目標運動特性反演提出了迫切的需求。研究人員對單星條件下的目標運動參數估計方法進行了較多的研究,參數估計方法主要包括基于彈道模板的彈道參數估計方法[3-4]和基于目標運動特性建模的參數估計方法[5]。為提高估計精度,研究人員進一步將目標動力特征的先驗知識應用于單星彈道導彈參數估計中[6],但目前基于單星的目標運動特性反演方法仍存在反演精度較低的問題,該類方法還不足以支撐目標參數測量的需求。因此,本文結合實際應用中所能獲得的目標先驗知識,研究具有更高估計精度的目標特性反演方法。
在單平臺下目標軌跡重建思路為:對獲取的序列圖像進行目標的提取,確定目標在圖像上的位置,獲得不同成像時刻的目標圖像位置。對每一個目標位置,根據對應的成像參數,計算目標視線的方向矢量,從成像器所在位置出發,考慮大氣折射,對目標視線的光學路徑進行反向追蹤,追蹤路徑與發射面的交點即為目標的空間位置,并通過多幀圖像反演獲得的序列目標空間位置進行擬合處理,修正目標空間位置,反演目標運動參數。計算過程中所需的發射面參數可根據經驗設定,也可以由其他計算方式獲得。單平臺下成像目標運動特性反演的總流程如圖1所示。
通過目標的提取,獲得成像時刻在相機成像的第k幀圖像上目標的位置p=(xk,yk)T,圖像行數r及列數c,以及像元視場角φpixel,進而可求出該位置對應于目標相對于光軸的偏角(αx,αy)為
目標視線方向矢量示意圖如圖2所示,Ose-XseYseZse為傳感器坐標系。根據光軸的偏角,可計算得到由相機發出并最終到達目標的光線在傳感器成像坐標系下的初始方向向量[7-8]Tc=(sinαxcosαy,sinαy,cosαxcosαy)。
根據如圖3所示的成像投影相關坐標系的關系,可將先驗獲得的發射點位置及發射面轉換到地心坐標系下,同時將目標視線矢量和發射面方程轉換到統一的地心坐標系,坐標系變換涉及傳感器坐標系、平臺坐標系、軌道坐標系、地心坐標系等坐標系之間的轉換[9-10]。
參見圖3,成像時刻tk時的坐標系轉換矩陣為

圖1 實現方案總流程圖Fig.1 General flowchart of implementation scheme

圖2 目標視線方向矢量示意圖Fig.2 Schematic diagram of target line-of-sight direction vector

圖3 成像投影相關坐標系的關系示意圖Fig.3 Schematic diagram of relationship between image projection related coordinate systems

(1)
式中:Mearth2orbit為地心坐標系Oe-XeYeZe到軌道坐標系Oo-XoYoZo的轉換矩陣;Morbit2sat為軌道坐標系Oo-XoYoZo到平臺坐標系Os-XsYsZs的轉換矩陣;Msat2sensor為平臺坐標系Os-XsYsZs到傳感器坐標系Ose-XseYseZse的轉換矩陣。圖3中的ept和αpt分別為傳感器指向在Os-XsYsZs下的俯仰角和方位角。
根據目標視線方向及平臺位置可得到過目標的射線方程為
(2)
式中:(xs,ys,zs)為傳感器的空間位置;(tx,ty,tz)為射線的法向量;(x,y,z)為目標點P的空間位置,可由以下聯立方程解得
(3)
式中:(qx,qy,qz)為平面法向量。
由于受到地心引力的作用,地表大氣密度分布不均。大量的實驗研究表明,大氣垂直方向比水平方向的變化要大1~3個數量級[11]。因此在實際研究大氣折射對光線傳播的影響時,常常忽略大氣水平方向的變化,將大氣視為水平方向均勻分布。將大氣視為球面分層,每層大氣均勻分布并且折射率相同,從而大氣折射率可簡化成僅隨海拔高度h而變化的量,即N=N(h)[11-14]。
如圖4所示,設相機C在地心坐標系中坐標為(xc,yc,zc),目標P在地心坐標系中坐標為(xp,yp,zp),地心Oe坐標為(0,0,0)。相機、目標、地心3 點組成的平面坐標系Oe-XpYpZp定義如下:平面坐標系的原點與地心直角坐標系的原點重合,平面OeCP與平面OeXeYe的交線OeXp為平面坐標系的Xp軸,Yp軸垂直平面OeCP,則由右手螺旋定則確定Zp軸。平面坐標系Oe-XpYpZp到地心坐標系Oe-XeYeZe的轉換矩陣為
(4)
式中:Mearth2plan為地心坐標系Oe-XeYeZe到平面坐標系Oe-XpYpZp的轉換矩陣。
為了方便計算光在每層中傳輸路徑,相機、目標、地心平面內的點坐標和向量可轉化到平面坐標系下,因為Yp軸坐標始終為0,所以將坐標記為二維平面坐標(xp,zp)的形式[15]。

當光線由第i層射向第i+1層時,根據第i+1層入射角,第i層大氣折射率,第i+1層大氣折射率,計算第i+1層出射角。如圖5所示,當光線由第i層射向第i+1層時,第i+1層光線的入射方向向量vi+1,入射點坐標Ti+1=(xplan(i+1),zplan(i+1)),第i+1層入射角λi+1,第i層大氣折射率ni,第i+1層大氣折射率ni+1,每層的地心夾角dθ已知。光線在每層中傳輸路徑的計算原理如下:

圖4 地心坐標系與平面坐標系之間的關系圖Fig.4 Relationship between geocentric coordinate system and plane coordinate system

圖5 相機、目標和地心組成的平面內,光在每層大氣中的傳輸路徑圖Fig.5 In plane composed of camera, target and center of the earth, transmission path of light in each layer of atmosphere
根據大氣球面分層Snell定律,由入射角λi+1、大氣折射率ni、ni+1,求得第i+1層出射角φi+1。
大氣折射率采用Hopfield模型進行計算,利用該模型可以求取任意海拔高度h的大氣折射率n=N(h),Hopfield模型考慮了大氣含有水汽的情況,將大氣折射率分為干、濕2項[16-17]。
根據第i+1層入射點法線向量和出射角計算出射方向的單位向量。由第i+1層入射點法線向量Qi+1=[xplan(i+1),zplan(i+1)]和出射角φi+1,利用向量內積求得出射方向的單位向量ue=[ux,uz],具體求解方法為
(5)
根據第i+1層的地心夾角、出射角、其距地心的距離以及出射方向的單位向量,計算該層出射光線的長度,進而得到其出射光線方向向量,即第i+2層入射光線方向向量。在△OTi+1Ti+2中,根據正弦定理,利用公式:
(6)
求得第i+1層出射光線的長度Li+1,進而求得第i+1層出射光線方向向量Li+1ue,該方向向量即為第i+2層的入射光線方向向量vi+2,其中第i+1層的地心夾角為dθ、第i+1層距地心的距離為Hi+l。
由第i+1層入射點法線向量和第i+2層入射光線方向向量,計算得到第i+2層光線入射點法線向量。進而根據向量間的關系Qi+2-Qi+1=Li+1ue,求得第i+2層光線入射點Ti+2法線向量為
Qi+2=Qi+1+Li+1ue=(xplan(i+2),zplan(i+2))
(7)
根據第i+1層光線的入射點坐標和入射點向量計算第i+1層光線的交點;判斷是否達到最后一層。在每層中求解射線與發射面的交點,已知發射點F(Xf,Yf,Zf)和發射面法向量fαj=[μαj,ναj,γαj],可使用點法式求解交點K(xk,yk,zk),當第i+1層光線的入射點Ti+1和第i+1層中交點Ki+l之間距離小于第i+1層光線的入射點Ti+1和第i+2層光線的入射點Ti+2之間距離時,交點Ki+l(xk,yk,zk)即為目標點。

(8)


(9)

在用二階多項式擬合時,為了保證足夠高的擬合精度,采用分段擬合的方式來減小擬合誤差。圖6為對實際軌跡采用上述擬合方式時,位置平均誤差和徑向速度平均誤差曲線。
考慮到實際檢測數據點數,僅僅比對兩者擬合點數在10~30范圍內的情況。可以看出,在擬合點數為10~30范圍內,位置平均誤差隨擬合點數增加而增加,而徑向速度平均誤差隨擬合點數增加而減小。發生上述情況的原因是擬合時3個分量各自分別擬合,而實際中軌跡的3個分量存在一定的耦合,隨著點數增多位置擬合誤差會增大,而擬合點數增多時擬合出的曲線形狀與實際更相似,所以速度誤差會降低。取折中的情況,采用擬合點數為20的分段擬合。

圖6 位置及徑向速度平均誤差隨擬合點數變化曲線Fig.6 Curves of location and radial velocity average error changing with number of fitting points
先驗給出的發射面參數經常存在誤差,且該誤差對重建精度有極大影響,因此需要在重建過程中同時優化發射面參數,以提高重建精度。目標運動在短時間內可近似為勻加速直線運動,發射面準確時即理想情況下,所重建的運動加速度應為一常量。若發射面參數不準確即發射面發生偏轉,目標將變為變加速運動。故在對重建目標點進行分段擬合時,理論上每一段的加速度較之于前一段的加速度差值其絕對值和應接近于零,而對于變加速運動,發射面偏轉越大,差值越大。

(10)
對應于發射面偏角αj,可反演出相應的目標運動參數(x,y,z,v,a,da)aj,分別為位置坐標(x,y,z),速度v,加速度a和加速度差da,則最優目標運動參數為
(11)
圖7為發射點經緯高為(E97.306°,N41.090°,6 369.853 km),發射面法向量為(-0.137,0.604,0.784)的軌跡在不同偏角情況下的反演加速度差值和的變化曲線,記錄加速度變化平均值的最小值以及相對應的發射面偏角、發射點位置、發射面法向量。將此條件下得到的軌跡點作為最后的反演軌跡。

圖7 發射面優化分析圖Fig.7 Diagram of launch plane optimization analysis
在發射點和發射面等先驗信息準確情形下,不考慮成像和檢測的影響因素,單平臺定位的位置平均誤差標準差和徑向速度平均誤差與指向誤差標準差的關系如圖8、圖9所示。
從圖8可以看出,單平臺定位時,分段擬合的方法在指向誤差極小的情形下會帶來誤差增大,但是指向誤差標準差在5″內時能夠有效地抑制指向誤差引起的定位誤差,并且誤差在500 m以下,有足夠空間考慮成像和檢測的誤差。圖9顯示指向誤差標準差變化范圍為(0″,50″)時,徑向速度平均誤差與指向誤差的標準差成正相關,徑向速度平均誤差小于60 m/s。

圖8 單平臺定位的位置平均誤差隨指向誤差標準差變化曲線Fig.8 Curves of single platform positioning’s location average error changing with pointing error standard deviation

圖9 單平臺定位的徑向速度平均誤差隨指向誤差標準差變化曲線Fig.9 Curve of single platform positioning’s radial velocity average error changing with pointing error standard deviation
由單平臺下檢測到的目標像面位置重建三維軌跡是典型的病態問題,先驗知識發射點位置和發射面指向顯的尤為重要,其中發射點位置能夠精確獲得,故發射面指向是影響單平臺定位誤差的關鍵因素。
從圖10可以看出,即使是根據已知的準確軌跡擬合出來的發射面都不是十分準確,而在實際目標飛行過程中,受到環境影響因素太多,發射面不能保持準確,并會產生很大的誤差。從圖中可以看出在發射面偏轉500″,也就是0.139°,在實際中可看做是較小的誤差,但是目標重建位置誤差已經在1 km以上。故單平臺定位在實際應用中存在很大的難題也就是發射面準確度的問題,其定位誤差遠大于雙平臺定位誤差。
在空間目標點位置確定的情況下,影響定位精度的另一個重要因素就是平臺位置誤差,在高空中平臺位置的變化給定位帶來很大的影響,分別對平臺位置誤差添加均值為0,不同標準差的隨機誤差,統計分析平臺位置誤差給定位精度帶來的影響,圖11為平臺位置誤差給目標重建位置和速度帶來的誤差。
從圖11可以看出平臺位置相較于其他影響因素,對目標重建速度誤差帶來的誤差影響是很大的。在平臺位置誤差為420 m時,重建速度誤差已經比較大,而重建位置誤差符合定位要求。


圖10 單平臺定位的位置誤差隨發射面誤差變化曲線Fig.10 Curve of single platform positioning’s location error changing with launch plane error

圖11 平臺位置誤差對目標重建位置誤差及速度誤差的影響Fig.11 Effect of platform location error on target reconstruction location error and speed error
本文針對單平臺下的觀測數據進行空間軌跡重建是典型的病態問題,提出單平臺成像目標運動特性反演方法,得到:
1) 本文提出從平臺位置出發,反向追蹤,將大氣劃分多層,計算在大氣傳輸路徑,最終與發射面相交得到目標位置,并通過迭代優化發射面和分段擬合等方法抑制重建誤差,此方法能夠有效反演單平臺目標運動特性。
2) 通過本文提出的單平臺目標運動特性反演方法,目標反演位置誤差在200 m以下,速度誤差在60 m/s以下,能夠用作實際場景。
3) 對重建出的原始數據進行分段擬合,能夠有效的抑制重建的位置誤差及速度誤差。
4) 在各影響因素中,發射面的精確度的影響最大,較小的偏轉就會帶來很大的反演誤差,而通過發射面的迭代優化,可以較為有效地抑制小范圍發射面的偏轉帶來誤差。