[發明專利]一種用于實時監測雙頻多星座星基增強系統性能的方法

有效

• 申請號: 202010445180.X
• 申請日: 2020-05-23
• 公開（公告）號: CN111596315B
• 公開（公告）日: 2022-07-22
• 發明（設計）人: 邵搏;丁群;原彬;耿永超;熊帥;吳顯兵
• 申請（專利權）人: 中國電子科技集團公司第二十研究所
• 主分類號: G01S19/07
• 分類號: G01S19/07;G01S19/11;G01S19/08;G01S19/23
• 代理公司: 西北工業大學專利中心 61204
• 代理人: 金鳳

1.一種用于實時監測雙頻多星座星基增強系統性能的方法，其特征在于包括下述步驟：

步驟一：數據預處理；

機場附近的監測站采集所監測到全球衛星導航系統(Global Navigation SatelliteSystem,GNSS)衛星的觀測數據、GNSS導航電文和DFMC SBAS增強電文，監測站i觀測到衛星j的雙頻觀測數據如下：

其中，和分別為L1和L5頻點上的偽距觀測量；和分別為L1和L5頻點上的載波相位觀測量；為監測站i和衛星j間的幾何距離；為對流層延遲；b_i為監測站接收機時鐘與GNSS系統時之間的偏差；B^j為衛星時鐘與GNSS系統時之間的偏差；為電離層延遲，對偽距觀測量的影響是滯后，對載波相位觀測量的影響是超前；f₁＝1575.42MHz為載波L1的頻率，f₅＝1176.45MHz為載波L5的頻率；和為偽距觀測量上的觀測噪聲；N₁和N₅為整周模糊度，λ₁＝C/f₁和λ₅＝C/f₅分別為載波L1和L5的波長，C為光速；和為載波相位觀測量上的觀測噪聲；

利用監測站i觀測到衛星j的雙頻偽距觀測量和雙頻載波相位觀測量進行數據預處理，i＝1，2，…，M，具體步驟如下：

步驟1.1：周跳探測；

周跳探測利用前5個采樣時刻t-1、t-2、t-3、t-4和t-5的載波觀測量外推出當前時刻的觀測量，并將此觀測量與當前時刻接收機的載波相位觀測量比較，如果超出門限，則認為出現周跳；監測站i對衛星j的周跳監測有下式：

式(5)和式(6)中，a₀、a₁、a₂為擬合系數，[a₀，a₁，a₂]^T＝(F^TF)^-1F^TX_L1-L5，F取固定值；為t時刻觀測量組合，和分別為t時刻L1和L5頻點上的載波相位觀測量；為由多項式擬合得到的t時刻觀測量組合的擬合值；T_L1-L5＝0.055為探測門限；

如果L1和L5頻點同時出現相同的周跳時，無法檢測出周跳，則只利用L5相位觀測量進行單頻點周跳探測，利用公式(7)、(8)再檢測一次：

其中，b₀、b₁、b₂為擬合系數，[b₀，b₁，b₂]^T＝(F^TF)^-1F^TX_L5；T_L5＝0.35為探測門限；

步驟1.2：雙頻載波平滑

載波相位觀測量通過周跳探測后，則認為沒有周跳出現，利用載波相位觀測量對偽距觀測量進行平滑，首先對載波觀測量進行如下變化：

式(9)、(10)中

用平滑偽距觀測量中的噪聲；

其中，L_k表示L1或L5頻點，為相應頻點的偽距觀測量，為相應頻點平滑后的偽距觀測量，τ＝100為平滑時間；

步驟1.3：消除電離層延遲；

利用L1和L5頻點平滑后的偽距觀測量消除電離層延遲，消除電離層延遲后的偽距觀測量為：

步驟二：導航電文解算；

GNSS衛星導航電文中播發的軌道參數為：星歷參考時間t_oe，衛星軌道長半軸a_s的平方根，軌道偏心率e_s，t_oe時刻的軌道傾角i₀，周內時等于0時的軌道升交點赤經Ω0_，軌道近地角距ω，t_oe時刻的平近點角M₀，平均運動角速度校正值Δn，軌道傾角變化率i′，軌道升交點赤經變化率升交點角距余弦調和校正振幅C_uc，升交點角距正弦調和校正振幅C_us，軌道半徑余弦調和校正振幅C_rc，軌道半徑正弦調和校正振幅C_rs，軌道傾角余弦調和校正振幅C_ic，軌道傾角正弦調和校正振幅C_is；利用導航電文播發的軌道參數及衛星位置解算算法，得到衛星星歷位置利用衛星星歷位置和監測站基準位置[x_R，i，y_R，i，z_R，i]計算星歷距離

利用導航電文中播發的參考時間t_oe、參考時刻的衛星時鐘偏差a_f0、衛星時鐘漂移速度a_f1和衛星時鐘漂移速度的變化率a_f2計算t時刻的衛星時鐘偏差

步驟三：增強電文解算

步驟3.1：差分改正數解算；

星歷改正數解算從DFMC SBAS電文32讀取星歷位置改正信息，參考坐標系為地心地固坐標系，計算公式如下：

其中，[δx^j δy^j δz^j]^T為衛星j在時刻t的星歷改正數；為參考時刻t_D衛星j的星歷改正數，由電文32播發；為參考時刻t_D衛星j的星歷改正數變化率，由電文32播發；

鐘差改正數解算從DFMC SBAS電文32讀取鐘差改正信息，計算公式如下：

其中，為衛星j在時刻t的鐘差改正數，單位：秒；為參考時刻t_D衛星j的鐘差改正數，由電文32播發；為參考時刻t_D衛星j的鐘差改正數變化率，由電文32播發；C為光速；

步驟3.2：完好性參數解算；

首先利用電文32播發的第j顆衛星的協方差矩陣信息比例因子sc和矩陣元素其中，x,y為1,2,3,4且x≤y”，其中為矩陣E^j中的元素，計算衛星j的計算公式如下：

R^j＝SF^j·E^j (15)

其中，I^j為衛星j到監測站的4維方向矢量，前三維是單位方向矢量，第四維是1；C_covariance由電文37播發；

完好性參數解算利用DFMC SBAS電文32和37播發的信息進行計算，計算公式如下：

其中，利用電文34、35、36播發的DFREI以及電文37播發的DFREI映射表得出；為電文32第一比特播發時間；R_CORR和I_CORR由電文37播發，由電文32播發，當t-t_CORR≤I_CORR時，當t-t_CORR＞I_CORR時，(R_CORR)_sv＝R_CORR；為向下取整；

步驟四：定位解算；

步驟4.1：對流層延遲估計

使用對流層模型進行修正，對流層延遲估計計算如下：

其中，d_hyd與d_wet分別表示對流層的干分量和濕分量，為仰角，b＝acos[cos(φ^j-φ_i)×cos(δ^j-δ_i)]，φ^j和δ^j分別為衛星j所在位置的緯度和經度，φ_i和δ_i分別為監測站i所在位置的緯度和經度；

d_hyd與d_wet由監測站高度信息及五個氣象參數的估值計算：

其中，g＝9.80665m/s²，g_m＝9.784m/s²，H為監測站海拔，單位為米，k₁＝77.604K/mbar，k₂＝382000K²/mbar，R_d＝287.054J/kg/K；

氣象參數氣壓P(mbar)、溫度T(K)、水汽壓e(mbar)、溫度變化率β(k/m)、水汽變化率λ由監測站的氣象傳感器提供；插值公式如下：

其中，北緯時D_min＝28，南緯時D_min＝211，φ_i為監測站所在緯度和D為年積日，ξ₀(φ_i)和Δξ(φ_i)分別表示緯度φ處的氣象參數平均值和季節變化值；

計算P、T、e、β、λ時，將ξ分別替換為P、T、e、β、λ，按表2由下式插值得到：

表2對流層延遲的氣象參數表

如果φ_i≤15或φ_i＝30或φ_i＝45或φ_i＝60或φ_i≥75，直接利用ξ₀(φ_i)和Δξ(φ_i)在表2中對應的數值通過式(25)計算，其他情況下，以φ_i＝40為例，對應的φ_k＝30，φ_k+1＝45，利用ξ₀(φ_k+1)、ξ₀(φ_k)、Δξ(φ_k+1)和Δξ(φ_k)在表2中對應的數值通過式(26)和(27)計算ξ₀(φ_i)和Δξ(φ_i)；

步驟4.2：監測站位置解算；

衛星j星歷位置經星歷改正數[δx^j δy^j δz^j]改正后的位置如下：

衛星j時鐘偏差經鐘差改正數改正后的時鐘偏差如下：

其中，C為光速；

從平滑后的偽距觀測量中消除對流層延遲和衛星時鐘偏差得到偽距如下：

其中，[x_i y_i z_i]為監測站位置，t_i為監測接收機時鐘偏差；

該偽距方程為非線性方程，要用泰勒級數展開，并取一階量，將偽距方程轉化為線性方程；

其中，

為監測站位置估計值，[Δx_i Δy_i Δz_i]為監測站位置與估計值之間的差值，為用戶時鐘偏差估計值，Δt_i為監測接收機時鐘偏差與估計值之間的差值；

將式(31)進行變換得到：

其中，

式(32)對應的矩陣形式為：

Z＝HX (33)

其中，

X＝[Δx_i Δy_i Δz_i -C·Δt_i]^T，N為監測站觀測到的衛星數量；

利用最小二乘法得到：

X＝(H^TH)^-1H^TZ (34)

則監測站位置和時鐘偏差為：

取為[x_i y_i z_i]，通過多次迭代當時，得到監測站位置和時鐘偏差；

步驟4.3：監測站定位誤差解算

利用定位解算得到的監測站位置[x_i y_i z_i]，并結合測繪標定得到的真實位置[x_R，i，y_R，i，z_R，i]，得到地心地固(Earth-Centered Earth-Fixed,ECEF)坐標系下的監測站定位誤差為：

[Δx_R，i Δy_R，i Δz_R，i]＝[x_R，i，y_R，i，z_R，i]-[x_i y_i z_i] (37)

從ECEF坐標系到東北天(East North Up,ENU)坐標系的轉換矩陣為：

其中，φ_i和λ_i分別是監測站所處位置的地理緯度和經度；

得到ENU坐標系下的定位誤差如下：

[ΔE_i ΔN_i ΔU_i]＝P_i·[Δx_R，i Δy_R，i Δz_R，i]^T (39)

基于式(39)得到水平定位誤差(Horizontal Position Error,HPE)和垂直定位誤差(Vertical Position Error,VPE)如下：

VPE＝ΔU_i (41)

步驟五：保護級解算

首先計算監測站與可觀測衛星間的觀測矩陣G，該矩陣的第j行如下所示：

其中，為監測站與衛星j間的仰角；為監測站與衛星j間的方位角；

監測站與第j個可觀測的衛星間觀測偽距的噪聲方差為：

其中，為衛星j對流層投影系數；

監測站與衛星間觀測偽距的協方差矩陣為W，其對角線元素其余元素為0，通過G和W得到：

針對航路、終端區、非精密進近(Non-Precision Approach,NPA)飛行階段，水平保護級(Horizontal Protection Level,HPL)由式(45)計算：

HPL＝K_H，NPA·d_major (45)

其中，

針對一類垂直引導進近(APproach with Vertical guidance I,APV-I)、二類垂直引導進近(APV-II)和CAT-I飛行階段，HPL和垂直保護級(Vertical Protection Level,VPL)由下式計算：

HPL＝K_H，PA·d_major (46)

VPL＝K_Vd_U (47)

其中，K_H，PA＝6.0，K_V＝5.33；

步驟六：服務性能評估；

對于針對航路、終端區、NPA飛行階段，如果HPE≤HPL，表明DFMC SBAS服務正常；如果HPE＞HPL，表明DFMC SBAS服務不能用于導航；

對于APV-I、APV-II和CAT-I飛行階段，如果HPE≤HPL且VPE≤VPL，表明DFMC SBAS服務正常；如果HPE＞HPL或者VPE＞VPL，則DFMC SBAS服務不能夠用于引導飛機進行精密進近，當DFMC SBAS服務不可用時，機場管制員將DFMC SBAS服務不可用的信息告知機場附近準備進近著陸的飛機，飛機需要采用其他導航手段進行進近。