[發明專利]一種基于球諧展開的GNSS單點定位方法

有效

• 申請號: 202110283670.9
• 申請日: 2021-03-17
• 公開（公告）號: CN113093242B
• 公開（公告）日: 2022-03-11
• 發明（設計）人: 郭金運;郭恒洋;楊洲銘;邢云鵬;劉新;孔巧麗
• 申請（專利權）人: 山東科技大學
• 主分類號: G01S19/42
• 分類號: G01S19/42
• 代理公司: 青島錦佳專利代理事務所(普通合伙) 37283
• 代理人: 朱玉建

1.一種基于球諧展開的GNSS標準單點定位方法，其特征在于，包括如下步驟：

I.1.建立利用偽距觀測值進行的標準單點定位觀測方程，如公式(1)所示：

式中，i表示第i觀測歷元，j表示衛星編號；ρ表示偽距觀測值，表示衛星j在第i觀測歷元的偽距觀測值；表示接收機的位置與衛星j在第i觀測歷元的位置之間的幾何距離；c表示光速；t_r表示測站接收機鐘差，(t_r)_i表示第i觀測歷元測站的接收機鐘差；t_s表示衛星鐘差，表示衛星j在第i觀測歷元的鐘差；表示衛星j在第i觀測歷元的信號傳播路徑上的對流層延遲誤差；表示衛星j在第i觀測歷元的信號傳播路徑上的電離層延遲誤差；ε_ρ表示偽距觀測數據殘差；

定義衛星j在第i觀測歷元觀測瞬間的空間三維坐標為測站的近似坐標為(X₀,Y₀,Z₀)，則衛星到測站近似位置的幾何距離表示為：

I.2.基于球諧展開表示與測站和衛星之間的高度角和方位角有關的誤差項，誤差項包括對流層延遲誤差以及電離層延遲誤差；基于球諧展開的標準單點定位觀測方程表示為：

式中，n為球諧展開的階數，m為球諧展開的次數，Nmax為球諧展開的最大階數；表示n階m次的締合勒讓德多項式；C_nm和S_nm分別表示n階m次的球諧展開的系數，C_nm和S_nm為球諧展開的待求參數；分別表示測站與衛星j在第i觀測歷元之間的高度角和方位角；為了方便表示球諧展開，令：

則公式(2)簡化為：

由簡化后的標準單點定位觀測方程(3)，得到標準單點定位誤差方程(4)，公式如下：

其中，v_ρ表示偽距觀測數據的改正值；

表示衛星j在第i觀測歷元的偽距觀測數據的改正值；

I.3.基于公式(4)得到標準單點定位誤差方程的線性化表達式，并進行簡化處理，然后對簡化后的標準單點定位誤差方程的線性化表達式進行滑動解算；

I.3.1.對標準單點定位誤差方程(4)在測站的近似坐標(X₀,Y₀,Z₀)處進行泰勒級數展開并保留一階項，則標準單點定位誤差方程的線性化表達式如公式(5)所示；

公式(5)中，令：

其中，表示測站近似坐標與衛星j在第i觀測歷元坐標計算出的待求參數dX的系數，待求參數dX是測站近似坐標X₀的改正數；表示測站近似坐標與衛星j在第i觀測歷元坐標計算出的待求參數dY的系數，待求參數dY是測站近似坐標Y₀的改正數；表示測站近似坐標與衛星j在第i觀測歷元坐標計算出的待求參數dZ的系數，待求參數dZ是測站近似坐標Z₀的改正數；表示由測站近似坐標與衛星j在第i觀測歷元坐標計算出的球諧展開的待求參數C_nm的系數；表示由測站近似坐標與衛星j在第i觀測歷元坐標計算出的球諧展開的待求參數S_nm的系數；d(t_r)_i表示在第i觀測歷元測站接收機鐘差的待求參數；

則簡化后的標準單點定位誤差方程的線性化表達式，如公式(6)所示；

公式(6)中，i＝1,2,…,epoch，epoch表示最大觀測歷元數；表示由衛星j在第i觀測歷元偽距觀測值、衛星j在第i觀測歷元的鐘差以及測站近似坐標與衛星j在第i觀測歷元之間的幾何距離計算出的常數項，其中，

設定滑動窗口選擇i個歷元，每個歷元最多能夠觀測到j顆衛星，則在這個滑動窗口下所有可視衛星的偽距觀測值組成imax個方程，其中，imax＝i×j；

定義標準單點定位誤差方程的系數矩陣、偽距觀測值的改正數向量矩陣、常數項矩陣以及未知參數矩陣分別為矩陣B、V、L、X，則它們的表達式分別如下：

X＝[dX dY dZ (dt_r)₁ … (dt_r)_i C₀₀ … C_NmaxNmax S_NmaxNmax]^T；

則標準單點定位誤差方程的線性化表達式(6)用一個矩陣形式表達，如公式(7)所示：

V＝BX-L (7)

未知參數X的最小二乘估值為：X＝(B^TPB)^-1B^TPL (8)

公式(8)中，P為單位權矩陣；

I.3.2.判斷公式(7)中系數矩陣B是否為病態，經過判斷，若系數矩陣B是病態，則執行步驟I.3.3；否則，系數矩陣B是良態，執行步驟I.3.4；

I.3.3.首先利用截斷奇異值分解的方法，為最小二乘譜修正迭代提供待求解的未知參數X的初值，然后基于最小二乘譜修正迭代計算未知參數X的值；

設系數矩陣B∈R^n×m，R^n×m表示n行m列的實數陣；則系數矩陣B的奇異值分解為：

B＝USV^T (9)

公式(9)中，U∈n×n，V∈m×m，U、V均為正交矩陣，S∈n×m為一對角矩陣；

系數矩陣B∈R^n×m的截斷奇異值矩陣B_k定義為：

公式(10)中，矩陣S中最小的(r-k)個非奇異值用零代替，即被截斷，且k≤r；r表示系數矩陣B的秩，k表示矩陣S中保留的奇異值數量；

u_i表示矩陣U對應的向量，v_i表示矩陣V對應的向量，σ_i表示矩陣S中保留的奇異值；

計算矩陣S中奇異值的均值，將該均值作為截斷奇異值的閾值，其中，矩陣S中奇異值大于截斷奇異值的閾值的保留，小于截斷奇異值的閾值的進行零化處理；

則公式(7)的截斷奇異值解為：

公式(11)中，

根據最小二乘原理V^TPV＝min，公式(8)寫成以下形式：

(B^TPB)X＝(B^TPL) (12)

將公式(12)的等式左端和右端同時加未知參數KX，并化簡得到公式(13)：

(B^TPB+KI)X＝B^TPL+KX (13)

由公式(13)整理得，最小二乘譜修正迭代公式為：

X＝(B^TPB+KI)^-1(B^TPL+KX) (14)

公式(14)中，K為任意實數；

基于最小二乘譜修正迭代計算未知參數X的解，轉到步驟I.3.5；

I.3.4.利用最小二乘法求解未知參數X的解，轉到步驟I.3.5；

I.3.5.根據計算得到的未知參數X的解，得到未知參數X中包含的參數值，即測站的位置參數(dX,dY,dZ)、接收機鐘差dt_r以及球諧展開的系數C_nm和S_nm；

將未知參數X中包含的位置參數(dX,dY,dZ)，分別加到測站的近似坐標(X₀,Y₀,Z₀)上，即標準單點定位求出的地心地固坐標系下的測站坐標(X,Y,Z)；

I.3.6.將測站坐標由地心地固坐標系轉換到站心坐標系，實現GNSS標準單點定位。