一種基于GPS單差模型的姿態解算算法

秦紅磊，陳萬通，金 天，叢 麗

(北京航空航天大學電子信息工程學院，100191北京，qhlmmm@sina．com)

一種基于GPS單差模型的姿態解算算法

秦紅磊，陳萬通，金 天，叢 麗

(北京航空航天大學電子信息工程學院，100191北京，qhlmmm@sina．com)

針對采用雙差模型無法克服雙差觀測量的強相關性這一問題，以觀測量不相關的單差模型為基礎，提出一種可遞歸處理的定姿算法．該方法用碼觀測量輔助載波相位觀測值，增加了信息量，通過遞歸估計模糊度浮點解及其方差協方差矩陣，從而可以利用LAMBDA算法估計整周模糊度．該算法計算量較小，且具有數值穩定性．實驗結果表明，該算法的初始化時間僅占雙差模型所用時間的1/50左右．該算法可以有效縮短初始化階段，能夠高效地用于動態姿態解算．

全球定位系統;姿態解算;整周模糊度;短基線

全球定位系統(GPS)具有全球性、全天候和連續的精密三維定位能力．應用GPS載波相位測量技術來確定載體姿態，因具有精度高、長期穩定的準確性、低成本和低能量耗損等優點而日益成為導航領域研究的熱點［1］．但由于載波是一種周期性的正弦信號，進行相位測量時存在未知的整周模糊度問題，必須先求得初始時刻的整周模糊度，才能獲得高精度的基線坐標［2］．采用雙差載波相位觀測方程能夠有效地減小電離層和對流層誤差、軌道誤差、衛星和接收機時鐘誤差，是姿態測量中通常使用的基本模型［3］．靜態情況下，基線坐標在任何時刻保持不變，可以利用最小二乘求解浮點解和LAMBDA算法估計初始整周模糊度;動態情況下，基線坐標在不同時刻都對應著不同的未知量，經典最小二乘法和LAMBDA算法不能直接應用，而采用“直接收斂法”由于雙差觀測量的強相關性，通常會有十幾分鐘的初始化時間，不利于實際姿態測量應用．

本文避開雙差模型，采用觀測量不相關的單差模型，遞歸求解模糊度浮點解及其方差協方差矩陣，使得可以應用LAMBDA算法進行整周模糊度估計，并在有約束條件的情況下，通過擴大候選值空間后再利用約束條件“識別”模糊度的方法，將初始化時間縮減至十幾秒，并給出了詳細的算法步驟和相應的實驗對比．

1 基本觀測方程

1.1 碼觀測方程

接收機s對衛星i的碼觀測方程為

1.2 載波相位觀測方程

接收機s對衛星i的載波相位觀測方程為

式中:φis為載波相位測量的小數部分;Nis為未知的整周模糊度;λ為L1載波的波長;vis為載波相位測量噪聲．

1.3 碼單差方程

對于短基線而言，采用單差方法可以消去電離層和對流層誤差、衛星軌道誤差和衛星時鐘誤差［4］．若在接收機A和接收機B之間針對同一顆衛星i進行碼觀測量的單差，則碼單差方程為

1.4 載波相位單差方程

碼觀測量的缺點是噪聲較大，因此在高精度測量中主要是利用載波相位觀測量，其觀測精度相對于碼觀測量要高2 ～3個數量級［5－6］．對于以A、B兩個天線為端點的短基線，其對衛星i的單差載波相位觀測方程表述為

根據式(5)，則整理方程為

1.5 載波相位雙差方程

將2個單差方程相減，可以消去接收機鐘差項．假設衛星1仰角最高，作為參考星，則由式(5)可知參考星的單差方程為

式(5)與式(7)相減可以得到衛星i對參考星的雙差載波相位觀測方程為

式中:▽ΔφiAB為A、B兩個天線到衛星i的雙差載波相位觀測值的小數部分;▽ΔNiAB為未知的雙差整周模糊度．

2 雙差觀測模型求解載體姿態

2.1 雙差觀測模型

對于m顆衛星而言，則存在m－1個載波相位雙差方程，寫成矩陣形式為

一般認為噪聲服從均值為零的正態分布，對于雙差載波相位觀測值向量y而言，其各個分量存在一定的相關性，y的方差協方差矩陣為

式中:σ2φ為接收機載波相位測量噪聲方差，一般認為 σφ=0.025周［7］;Em－1為 m － 1 階矩陣，其所有元素均為1;Im－1為m－1階單位矩陣．可見雙差觀測值y各個分量之間存在強相關性．

式(9)是虧秩的，單歷元的情況下無法得到唯一的解，必須將多個歷元的觀測方程進行聯立，通過最小二乘估計未知參數．但最小二乘估計的基線和模糊度參數均為實數解，而模糊度是整數，解決此問題通常是先求解模糊度的浮點解，然后根據其方差協方差矩陣估計模糊度固定解，再利用固定解重新修正基線坐標．

2.2 靜態觀測模型

對于靜態基線，b是保持不變的，對于連續的q個歷元有

則其數學模型和統計模型可以看成

即得到基線向量坐標和模糊度浮點解及其相應的方差協方差矩陣為

然后通過目標函數最小來求得模糊度固定解為

LAMBDA算法是求解式(15)的有效算法．

2.3 動態觀測模型

對于動態情況，每個歷元k都對應一個不同的基線坐標bk，則其批處理模型為

則式(16)中矩陣階數會隨著k的增加而增加，特別是求解各歷元的浮點解時，矩陣求逆計算量會非常大，無法滿足實時處理的要求．一種解法是將方程組進行解耦變換［9］，即

移項有Bi(yi－vi)=Biz，但是觀測噪聲是未知的，若忽略噪聲，可得到模糊度的近似浮點解為

這種做法雖然可以避免大矩陣求逆求得模糊度的浮點解，也易遞歸實現，但是需要多個歷元來等待模糊度浮點解的收斂，即直接收斂法．其收斂的速度曲線基本呈指數衰減的趨勢，在開始的2～3 min內，向量NDD可迅速收斂到距真實值2～3個波長λ的位置，以后的收斂速度就變得十分緩慢．要到400～600 s后，NDD才逼近到距真實值NDD＜0.5λ的位置，從而獲得正確的NDD．這主要是由于短時間內雙差觀測量的強相關性造成的，由于式(19)的特殊結構無法準確求得浮點解的方差協方差矩陣，所以無法直接使用LAMBDA算法進行估計．

3 單差觀測模型求解姿態模型

3.1 求解單差整周模糊度

假設存在m顆衛星，則根據式(6)和式(4)將分別存在m個載波相位單差方程與m個碼單差方程，寫成矩陣形式分別為

式中E為接收機到衛星的設計矩陣．為了將噪聲統計特性歸一化［10］，令 σr≡ σφ/σρ，其中:σφ，σρ分別為載波和碼噪聲的標準差;σr為其比值，一般取為0. 001，將其乘在式(21)的兩端，則有

則式(23)等價于

其中a的各個分量均為整數．將k個歷元的上述方程進行聯立，其矩陣形式為

通過Householder變換進行QR分解可以得到正交陣Tk，式(28)兩邊同時乘以TTk，使得:

LAMBDA算法對初始模糊度進行Z變換，變換結果為

在變換域求解到z，在通過反變換求得原始的雙差整周模糊度為

3.2 模糊度檢驗

由于各種誤差的作用，并不能保證所求出的模糊度就是正確的．因此需要對模糊度候選值進行檢驗，以得到正確的固定解．工程中最為常用的一種是比率檢驗，即式中Z矩陣及其轉置的作用是對搜索空間進行不同方向的“擠壓”，使搜索空間由拉長的高維橢球近似于球，以提高搜索效率，經過變換，式(34)的目標函數等價為

式中:as為次優模糊度固定解;α為常數，其如何選取該關鍵值仍是一個開放的問題［11］，目前多通過經驗選取．為了確保模糊度驗證的可靠性，本文采用的方法為如果連續n個歷元的基線坐標都通過了基線長度檢驗且初始整周模糊度的固定解都一致，就可以認為該模糊度求解正確，即看是否滿足

否則，將采用更多個歷元繼續遞歸求解，直到可以確定初始整周模糊度．

3.3 求解基線坐標和鐘差

求得整周模糊度a，代入式(23)即有

3.4 算法遞歸實現

在動態載體應用中，數據處理常采用遞歸的方式，為了便于實現模糊度浮點解的遞歸估計，利用第k歷元估計的結果，則第k+1歷元有

通過 Householder變換實現 QR分解，即式(41)兩邊同時乘以正交陣TTk+1，即

用k+1個歷元的數據得到的初始整周模糊度的浮點估計值為，那么

式(45)與式(31)形式統一，即在各歷元都可以采用LAMBDA算法對浮點解進行整周模糊度估計．另外，遞歸算法使用QR分解來完成最小二乘計算，具有數值穩定性，而若使用遞歸最小二乘法，則會出現隨著時間的增長，不可避免的截斷誤差或舍入誤差導致結果偏離真值．

4 附加約束條件下的算法改進

4.1 歷元數與模糊度正確性的關系

4.2 擴大候選值空間

減少歷元數將導致模糊度浮點解的精度降低，雖然低精度的浮點解經變換后可能落到正確z的鄰近整數點的歸屬域中，但其一定在z的周圍．假設其落入z'的歸屬域，那么以z'為中心一定半徑的覆蓋范圍內的全部歸屬域所對應的整數點會包含z，如果不包含，則擴大覆蓋半徑到一定程度，必定包含z．如果能夠對通過擴大候選值空間得到的所有候選點進行符合某些準則的判別和檢驗，從而篩選出正確的整周模糊度，那么則對k的精度也放低．幸運的是，LAMBDA算法能夠按照“最佳程度”依次給出多個候選值，這里用“正確候選點出現的位置”表示正確的候選值在多個候選值中的次序．例如若第3個候選值是真實的模糊度，可記作“正確候選點出現的位置為3”．

4.3 根據約束條件篩選候選值

基線長度或俯仰角等約束信息可以輔助篩選正確的整周模糊度，令LAMBDA算法給出的候選值個數為N，具體篩選步驟為:對于某個候選值a'，代入式(40)，求得對應的基線坐標記為b'k，則其基線長度為

式中:δbl為預先設定的閾值;bl為基線的真實長度．在姿態測量應用中，δbl的選取是基線長度約束的關鍵．因為天線相位中心的變化和噪聲影響，不能將δbl的值設置的過小，同時過大的δbl又無法有效地識別a是否為正確的整周模糊度．該值一般為基線真實長度的1%．

如果由a'得到的b'k滿足式(47)，則可利用b'k求解的姿態角繼續進行約束檢驗，即對于地面載體，滿足|θ|≤10°;對于有慣性器件輔助的載體，則應滿足|θ－|≤δθ并且|ψ－|≤δψ，其中:和分別為慣性器件給出的俯仰角和航向角觀測值;δθ和δψ分別為設定的誤差變化范圍．如果通過該檢驗，則可認為a'在很大概率上為正確的模糊度．如果連續幾個歷元計算初始整周模糊度相同，則基本可以確定其為正確解．

如果由a'得到的b'k沒有通過基線長度檢驗或者沒有通過姿態角檢驗，則需要試驗下一個模糊度候選值，直到篩選出正確的整周模糊度;如果所有模糊度候選值均不滿足，則需擴大候選值個數中的N值，或者需要更多歷元繼續按照上述方法檢測，直至找出滿足約束條件的模糊度候選解及其解算得到的姿態角．

5 實際數據測試

為了比較新算法的優勢，進行了單基線姿態測量的試驗．數據采集地點為北京航空航天大學新主樓樓頂，即盡可能避免信號被遮擋，OEM板型號為NavAtel公司的Super Star II，其載波相位輸出頻率為1 Hz，天線型號為GPS－701GG，其內含扼流圈，具有多徑抑制功能．基線長度為2 m，保持水平放置，持續跟蹤的可見星數目為7顆．為了得到初始時刻的雙差整周模糊度真值，以作為參考基準，首先采用靜態觀測模型算法進行求解，通過1 200個歷元的觀測數據 求 得 NDD= ［ 1， － 9， 7，－ 8，－ 8，7］T，b =［0.868 3，1.805 8，0.017 0］T，ψ =64.320 3°，θ=0.484 9°．

采用動態觀測模型的直接收斂法求解，利用式(19)得到的模糊度浮點解的收斂情況如圖 1，2所示，其中:“?”為該時刻浮點解四舍五入后沒有映射到正確的整周模糊度上;“·”為映射成功．可以看出直接收斂法確定模糊度分量分別至少需要 452， 465， 543， 335， 478，517 s才會收斂到正確整周模糊度的0.5λ范圍內，也即需要近550 s才能保證所有的模糊度分量浮點解四舍五入后收斂到正確的整周模糊度上，進而求解出姿態．

圖1 模糊度第 1， 2，3分量收斂情況

圖2 模糊度第 4， 5，6分量收斂情況

采用基于單差模型的算法進行解算，首先假設無基線長度和姿態角約束等先驗信息．實驗結果如圖3所示．結果顯示，僅用31 s的觀測數據，就可以得到正確的姿態解，解算過程中求得初始時刻的單差整周模糊度為

以第1顆星為參考星，可以得到對應的雙差整周模糊度N(1)DD，與真實值一致．

圖3 姿態角解算結果

對比圖3和圖4可知，基于單差模型的算法求解姿態的時間要遠小于直接收斂法，其原因是單差模型的觀測量各個分量之間不相關，同時利用了多個歷元的數據遞歸求得初始時刻整周模糊度的浮點解及其方差協方差矩陣，進而使用LAMBDA算法準確地估計出單差整周模糊度，從而縮短了初始化時間．

圖4 基于約束條件的算法結果

由于基線長度事先可以測量，對于地面載體，俯仰角也存在約束，則可采用附加約束條件下的算法改進的算法求解．為了避免計算量過大，候選值個數N選定為500．實驗結果如圖4所示，此時初始化時間僅為10 s，小于直接收斂法所用時間的1/50．對比圖3可知，比不用約束條件時的初始化時間還要少20 s，從而說明基于約束條件的算法更加有效．圖5顯示了前100個歷元識別到的正確模糊度的位置，在第31 s之前，模糊度正確解是通過約束條件篩選出來的，出現位置雖然較遠，但可以看出隨著時間的推移其逐漸“向前”靠攏，直至第31 s以后，第1個候選值均為正確解．

圖5 利用基線長度篩選候選點

6 結論

1)采用觀測值無關的單差模型，克服了雙差模型由于觀測值存在相關性而導致的模糊度浮點解收斂速度太慢的缺點．

2)可利用LAMBDA算法提高模糊度估計的準確性．

3)該算法可遞歸實現，數值計算穩定，計算量較小，能夠滿足實時應用．

［1］JOHN R．Current issues in the use of the global positioning system aboard satellites［J］．Acta Astronautica， 2000，47(2/9):377－387．

［2］ XIA Kewen，ZHANG Xinying，GAO Jinyong，et al．Study on GPS attitude determination technology based on QPSO algorithm［C］//Proceedings of the 7th World Congress on Intelligent Control and Automation．Washington DC:IEEE，2008:1869－1873．

［3］ CHANG Xiaowen，PAIGE C C，YIN Lan．Code and carrier phase based short baseline GPS positioning:Computational aspects ［J］．GPS Solutions， 2004，7(4):230－240．

［4］HOFMANN-WELLENHOF B．GNSS-Global Navigation Satellite Systems GPS，GLONASS，Galileo，and More［M］．New York:Springer Wien，2008:227－234．

［5］CHANG Xiaowen，PAIGE C C．An orthogonal transformation algorithm for GPS Positioning［J］．SIAM Journal on Scientific Computing， 2002，24(5):1710－1732．

［6］CHANG Xiaowen，PAIGE C C．An algorithm for combined code and carrier phase based GPS positioning［J］．BIT Numerical Mathematics， 2003，43(5):915 －927．

［7］ MISRA P，ENGE P．Global Positioning System:Signals，Measurements，and Performance［M］．2nd ed．Lincoln MA:Ganga-Jamuna Press，2006:238 －239．

［8］VERHAGEN S．The GNSS Integer Ambiguities:Estimation and Validation［D］．Delft:Delft University of Technology，2003:34－35．

［9］俞文伯，高國江，趙剡，等．單頻GPS動態相對定位的模糊度逼近搜索解法［J］．北京航空航天大學學報， 2002，28(2):242－244．

［10］CHANG Xiaowen，GUO Ying．Huber＇s M-estimation in relative GPS positioning:Computational aspects［J］．Journal of Geodesy， 2005，79(6):351－362．

［11］VERHAGEN S．Integer ambiguity validation:An open problem［J］．GPS Solution， 2004，8(1):36－43．

An attitude determination algorithm based on GPS single-differenced mode

QIN Hong-lei，CHEN Wan-tong，JIN Tian，CONG Li

(School of Electronic and Information Engineering，Beihang University，100191 Beijing，China，qhlmmm@sina．com)

It is difficult to overcome the high correlation of double differential observation using double-differenced model and the initialization procedure costs a long time with makes it difficult to apply in practice．To deal this problem，a recursive algorithm based on the single-differenced model is presented in which the observations are not correlated．The new method increases observation information with the aid of code for original carrier phase observation．At the same time the float ambiguity solutions and its variance-covariance matrix can be obtained recursively and it makes the LAMBDA method be used to estimate the integer ambiguities．The method is also less computation with numerical stability．The experiment results show that the time of initialization is reduced to about 1/50 of that using the double-differenced method，so this method can be effectively applied in practical attitude determination for shortening initialization procedure．

GPS;attitude determination;integer ambiguity;short baseline

V249.3

A

0367－6234(2011)09－0105－07

2010－01－08．

國家高技術研究發展計劃資助項目(2009AA12Z313)．

秦紅磊(1975—)，男，博士，副教授．

(編輯 張 紅)