1、1 / 6參數估計理論在 GPS 解算中的應用錢鵬鵬何秀鳳(河海大學地球科學與工程學院江蘇南京 210098)1 概述GPS 在變形監測、線路勘測、大氣參數測量、車船及高速飛行器導航定位等領域有著十分 廣泛的應用。然而從 GPS 接收機獲取的數據不能被直接應用，必須經過適當的數據處理。GPS 基線方程的解算是數據處理的重要環節。以往幾乎所有的基線解算軟件都是將GPS 基線方程約束在線性空間解算，這不僅產生較大的模型誤差，而且導致函數特征的改變。然而 不論是傳統的幾何大地測量學，還是現代空間大地測量學、攝影與遙感測童學都面臨著復雜的非線性數據處理問題；一方面要求側量數據質量日益精化，另一方面新的

2、非線性模型空間 數據處理理論尚未形成。為克服這一矛盾，現有的測量數據處理方法正朝著以下兩個方向努 力：矚慫潤厲釤瘞睞櫪廡賴賃軔朧。(I)尋找非線性函數的先驗真值，在真值附近按泰勒級數展開，盡量減小因線性化而產生的模型誤差；計算過程中盡量避免函數求導，特別是高階求導，力求保證函數的既有特性，同時減少 計算工作量。本文試圖探索一種非線性函數的差分迭代方法。方法的最大優點是：在迭代求解時，無需計算 f(X)各分量f,X)的偏導數值，同時使函數求值次數顯著減少。以期為擺脫目前測量數據處理理論所面臨的困境做出些許貢獻。聞創溝燴鐺險愛氌譴凈禍測樅。2分線性差分迭代解算原理設有非線性方程：V=f(X)-L(

3、 2-1)式中，X 為參數向量，L 為觀測值向量。將(1)式按泰勒公式展開為1 2f(Xo) f(X)xo(X -Xo) -f (X)x(X -X。)2Rn(X)-L(2-2)由于先驗值Xo充分接近真值X，因此將(2-2)式截斷為這里，B(X) = f (X)X%。由(2-3)式構造 Gauss-Newton 迭代格式Xk 1= Xk-(B(xk)TB(xk)(B(xk)TP( f (xk) - L), k二0,1,2,.n(2-4)式中，P 為觀測權矩陣，B(xk)實際是一個 n 階 Jacobi 方陣。因每步迭代都要計算出n2個 偏導數值，這一方面計算工作量大， 特別當 f(X)的各分量f

4、/x)復雜時求導困難。利用xk附近的差分值代替xk處一階導數，這樣可以避免求導，同時極大地減少計算工作量。殘騖樓諍錈瀨濟溆塹籟婭驟東。由差分與倒數的關系知壯=hnfn( ),(xkXj(2-5)V = f(Xo) B(X)(X -Xo) -L(2-3)2 / 6取 n=1 時即為一階導數與一階差分的關系:：fk=hf()，(Xk,Xk i)(2-6)這里h=Xk-Xkj，為步長，恒值；經驗表明在計算中取中心差分計算導數更加準確即k二f(Xkh/2) - f(Xk-h/2)/2(2-7)因此可將(2-4)式化為Xk 1=xk-( fk/h)T：.fk/hJ(：fk/h)TP(f(xk-L)k =

5、0,1,2,.n(2-8)式中，P 為觀測值得權陣，上式表明，在迭代求解時無需計算f(X)的各分量fj(x)的偏導數值，這是一個顯著的優點，但每次迭代都要計算輔助差分點的函數值，如果能合理的選擇輔助差分點，可使函數求值次數顯著減少。釅錒極額閉鎮檜豬訣錐顧葒鈀。3輔助差分點的選取方法設有初始近似值向量Xo出發已求的第 k 次近似值Xk和相應的函數值f (XJ，為求得第K+1 次近似值Xk 1，采取下面闡述的輔助點選取方法；彈貿攝爾霽斃攬磚鹵廡詒爾膚。已知Xk, f (Xk)，設在已知點Xk附近有 2n 個點Xk1,Xk2,.,Xkn及yk1,yk2,對這 2n 個點的要求是：點Xk(i =1,2

6、，,n)可以相互重合，亦可以不是；yx互不相同且式兀。現取輔助點不重合的情形：謀蕎摶篋飆鐸懟類蔣薔點鉍雜。rXk,i =Xk工4Z-；-( 3-1)iyk,i =Xk_(i)+hk_(_i4)x：ek_(i4)其中：hkg4)=0(i =1,2,.,n;k一n)為給定數。此時已知函數值為f (Xk,J =f(Xzi4)- (3-2)if(yk ,i ) =f(Xk工4) *hk _(i 4) Xek_(i)為了使算法盡可能利用已知信息，對輔助點中的eki4)(i =1,2,.,n)向量作如下兩種定義1.令 k=qn,q=1,2,.時，取 Q 衛書二eqn書=編單書(e 為單位向量，i=1,2,

7、3,)則有-hk _n十000 1人Xk =0hk-n七0000003 / 6n _(i書,0, 1 ,0,0,此時：Xk形成如下的對角矩陣：廈礴懇蹣駢時盡繼價騷巹癩龔。相應的m n階矩陣：fk為:000hk一(3-3)eTk -(i 4)二0,4 / 6也fk=f(Xkm丄嚴ejf(XkJ.,f(Xk+hkXeQf(xQ(3-4)式中hk車)7 故5)，存在。于是可以唯一的 確定矩陣Ak-：fk(.x)-12令 q*k(q+1)n,(q=1,2,取ek= eqn寸=ej( j : n),其中e：二0,.0,1,0.,0,此時-Xk可表達為:hqn十hqn 2(3-5)h(q)n (j 1).

8、hqn這時相應的m n階矩陣.fk為:i-fki-fqnj f(xqn1hqn 1eqn 1)_f (xqn1),.,f(Xqn亠jhqn jeqn亠j) _f(xqnk),.,f(X(q)n (j 1)h(qd) n (j 1)e(qd)n -(j 1)f (x(q J)n(j1),.,f (xqnhqneqn_ f (xqn)（3-6）1_1式中（1jn）。因hk護，故（必2）存在，可以唯一的確定矩陣A=3k（AXk）.將式（11）、12）或（13）、（ 14）式代入（8）式可解算出X心。實際計算時可取hj=o（|f（Xj）|），特別時取：hk_n 1=hk_n 2 =. = hk=o(

9、f(xQ )。上述迭代解法，除了點Xk和X處的函數值需計算外， 其他均為已知信息。算法計算工作量 少，且簡單易行。4 非線性 GPS 基線方程的解算模型目前 GPS 定位的方法有測距碼和載波相位定位兩種，測距碼偽距定位的基本方程為：Pk=（xj-Xk）2（yj- yk）2（zj-Zk）2Ctk一ck-pOn卩爲vk（4-1） 式中，j表示衛星編號，k 表示測站點號，p 為偽距，（xj,yj,zj）為衛星坐標，（xk,yk,zk）為測站坐標，（4-1）式右邊的后 5 項依次為接收機鐘差、衛星鐘差、電離層誤差、大氣層誤 差、隨機觀測誤差（不考慮時間延遲影響），由于測距碼波長長并受 p 碼的制約，因

10、而定位精度低，一般只作單點定位測量，測地形GPS 接收機大都采用載波相位定位，載波相位測量 t1 時刻的偽距誤差方程為：煢楨廣鰳鯡選塊網羈淚鍍齊鈞。vk=-用1-1/-込1-k+c盯MNkj+也吧 +卯爲 + 旳kron( 4-2 )cc基線較短時（4-2）式中-討/C項為零。由上式構成 ti 時刻站間單差基線求差模式中， J5 / 6為單差虛擬觀測值；j 表示衛星編號；i、k 表示測站編號；f 為頻率；鵝婭盡損鶴慘歷蘢鴛賴縈詰聾。八ik（k-彳）kfFk - Mk * Vjk（4-3）cC 為光速；Fk表示接收機 i、k 的相對鐘差；N 表示相位整周差。lij的表達式如下：likj-f2?k

11、j冷一彳計一f仁Tk一譏Ti（ 4-4）cc站間一次差消除了衛星鐘差、軌道誤差、電離層與對流程延遲誤差。上述非差法（4-2）與求差法（4-3）本質是一致的，只是后者求解參數相對少且精度相對較高，因此載波相位測 量多用求差法解算。本節將重點討論單差、雙差模型的非線性解法，三差模型的解算原理與 此無異。分析（4-2）式：籟叢媽羥為贍債蟶練淨櫧撻曉。1:ion trop分別為電離層延遲和對流層延遲，經相應的模型改正或技術處理后作差，其影響可以忽略；2.N為整周模糊差，待定；3. ik為相位差，實際測量值；4.、：Tk為接收機站間鐘差，當 GPS 接收機采用外接原子頻標時，測站鐘差用多項式描述如下：r

12、k（ti） =ay（ti-t。）a2（ti-t。）2（4-5）5.lj為i、k兩站上 GPS 信號傳播距離變率，微小項，可忽略；考慮上述因素，取測站i 為已知點，k 為未知點。其 WGS 84 坐標分別為（人，,乙）,（xk,yk,zk）用 i 點坐標表示 k點的坐標預頌圣鉉儐歲齦訝驊糴買闥齙。Xk=Xi+ AXki,yk= yi+yki,Zk=z +Zki（4-6）將（4-6）代入偽距T的空間距離表達式，然后將（4-2）式在（備，,乙）處按多元泰勒級數 展開并取至一階項，顧忌到 j-i,j-k 信號路徑的相似性得到站間單差（ SD）基線誤差模型滲釤 嗆儼勻諤鱉調硯錦鋇絨鈔。Vjf （Xj-X

13、fJ2（yj-y2（xj-Zi）2- .（xj二Xk）2（yj匚yk）2（Zj匚Zk）2c（k）j-（f、Tk）j-（MJ（4-7）（4-7）式中，衛星 j 的坐標已知（從其星歷中獲得）如圖1。顯然如果在各歷元ti時刻兩站均同步觀測 n 顆衛星，且觀測時間內共觀測 m 個歷元，則可組成n m個方程：其中 3 個未 知坐標參數，n 個整周模糊度，m 個站間鐘差未知參數，當（3+n+m ） r,則可在 LSE 或 LAS 準則下求解。贓熱俁閫歲匱閶鄴鎵騷鯛漢鼉。n1 n,r r.1n.1構成二范和目標函數得解X -（ATA）JATl構成一范的目標函數f（X）=E Vikjp,依（2-8）式構造參數

14、解的迭代方程式中，X =C x y z, N1,Nm4）為未知參數向量，h 為均差值，X0中基線分向量取C/A 碼偽距定位所得的點位坐標，X0中整周模糊度向量取初始化值。迭代終止條件為：當Xk4t-Xk|v &（ E 為給定的微小正數）時迭代終止。壇搏鄉囂懺蔞鍥鈴氈淚躋馱釣。5 算例用 Trimble 4000SE 在 A、B 站采樣，多時段，90m/時段，采樣率 15s，截止高度角為15。觀測值為 C/A 碼偽距和L1載波相位。SAB= 1549.254m.觀測基本信息：蠟變黲癟報倀鉉錨鈰贅籜葦繯觀測起始時間：01:00:45:000 ( GPS 時為：522045.00) 觀測終止

15、時間：02:17:15:000 ( GPS 時為：526635.00)vk化耐 計農 V 化諾NikP_DDikP(4-8)Xk+=Xkf(Xk)hTfh3fhQTf(Xk)-Lhh(4-9)7 / 6表 1 基線 AB 觀測歷元統計衛星編號58917212330292515起始歷兀0000000030102終止歷兀10218470102102102102102102表 2 測站偽距定位結果(WGS-84 坐標系)測站xyzA-2191828.68745182411.77052993199.3386B-2190397.45585182831.28822993595.5923表 3 整周模糊度初始化值衛星號9172125302925模糊度27450242954796-703439498287462258319329848288表 4 雙差基線向量解(WGS-84 坐標系)分量dxdydzDis(m)坐標差-1453.2628-378.3641-380.76341549.23526結束語1迭代收斂速度依賴于收斂值得準確度，取CIA 碼 偽距 定位 所得 的 點位 坐標 及 而且 可靠 的 ,且加快收斂速度；2由差分該代替求導值，完全避免了對非線性方程的迭代求導，幾同時極大地減少了計算工 作量；3以向量空間參數距離最小為迭代終止條件；