干堅(jiān)定，周 適，段太生，郭 平，晏 勇，劉立正

（1.中鐵二局第五工程有限公司，四川 成都 610091；2.中鐵二局集團(tuán)有限公司，四川 成都 610031）

北斗衛(wèi)星導(dǎo)航系統(tǒng)（BDS）最后一顆全球組網(wǎng)衛(wèi)星（GEO衛(wèi)星）已于2020年6月23日9：43在西昌衛(wèi)星發(fā)射中心成功發(fā)射，這是由我國自主研發(fā)的衛(wèi)星導(dǎo)航系統(tǒng)。對(duì)于全球?qū)Ш叫l(wèi)星系統(tǒng)（GNSS）而言，GPS雖然仍是當(dāng)前技術(shù)最先進(jìn)和最成熟的導(dǎo)航定位系統(tǒng)，但BDS通過近3年來高密度的衛(wèi)星發(fā)射，目前已具備和GPS共同作用、兼容并存的實(shí)力。從衛(wèi)星數(shù)量和衛(wèi)星類型上看，北斗的衛(wèi)星數(shù)量和亞太地區(qū)的衛(wèi)星分布比GPS更具優(yōu)勢，國內(nèi)大量專業(yè)人士對(duì)BDS進(jìn)行了相關(guān)深入研究。BDS即將組網(wǎng)完成，極大地增強(qiáng)了國人的自信心和民族自豪感，也激發(fā)了眾多相關(guān)從業(yè)人員學(xué)習(xí)和研究BDS的興趣。

1 相鄰歷元相位平滑偽距的數(shù)學(xué)模型

楊元喜等[1]提出，對(duì)北斗所有衛(wèi)星進(jìn)行統(tǒng)計(jì)，北斗B1I和B2I偽距野外測量的精度約33cm，B1和B2載波相位野外測量的精度約2mm。由于載波相位觀測的精度遠(yuǎn)高于偽距觀測[2]，因此采用偽距值進(jìn)行單站定位計(jì)算時(shí)，有必要對(duì)偽距值進(jìn)行平滑處理。偽距平滑的成熟算法包括載波相位平滑偽距算法和基于卡爾曼濾波的偽距平滑算法[3]。文章采用簡易數(shù)學(xué)公式計(jì)算，BDS共設(shè)為3個(gè)頻率，由于外業(yè)使用的接收機(jī)只采集B1、B2兩個(gè)頻率的觀測數(shù)據(jù)，因此文章中的數(shù)學(xué)模型是針對(duì)北斗雙頻（B1、B2）展開的。對(duì)BDS不同頻率分別列出偽距平滑公式，由于相鄰歷元時(shí)間較短，因此偽距中包含的噪聲（電離層延遲、對(duì)流層延遲、多路徑效應(yīng)）可以忽略不計(jì)，整周模糊度N也可消去，具體公式如下[4]：

同理，可得B2頻率下的各參數(shù)含義。利用相位平滑后的偽距進(jìn)行后續(xù)單站定位計(jì)算，若監(jiān)測到載波相位發(fā)生整周跳變，則相關(guān)歷元不采用相位平滑的偽距值，而采用原始偽距觀測值。

2 BDS單站單歷元定位的數(shù)學(xué)模型

在t歷元時(shí)刻，接收機(jī)共接收到n顆北斗衛(wèi)星，通過北斗廣播星歷文件可計(jì)算每顆北斗衛(wèi)星的空間坐標(biāo)，選取高度角最高的衛(wèi)星作為參考衛(wèi)星r，其余衛(wèi)星記作衛(wèi)星s，將偽距觀測方程線性化后，通過星間單差（其他衛(wèi)星和參考衛(wèi)星作差），可列出B1頻率的偽距觀測方程，此時(shí)接收機(jī)鐘差可消去，并忽略多路徑效應(yīng)誤差，具體公式如下[5]：

式中：x0、y0、z0為測站點(diǎn)初始坐標(biāo)，初始坐標(biāo)可設(shè)為(0,0,0)；(xr,yr,zr)為參考衛(wèi)星r坐標(biāo)；(xs,ys,zs)為其他衛(wèi)星s坐標(biāo)。由于t歷元時(shí)刻共接收n顆衛(wèi)星，由式（2）可列出n-1個(gè)方程。

同理可列出B2頻率的星間單差偽距觀測方程，具體公式如下：

采用不同頻率觀測值的線性無電離層組合，由于北斗B1頻率值為1561.098MHz，B2頻率值為1207.14MHz，可寫出北斗B1和B2頻率的波長比值，具體公式如下：

因此式（2）乘以7632，式（4）乘以-5902，相加后移項(xiàng)可得到如下公式：

式（6）已消去大部分電離層延遲，方程左邊為常數(shù)項(xiàng)，其中對(duì)流層延遲，可采用Saastamoinen模型進(jìn)行改正計(jì)算，對(duì)流層改正模型如下[6]：

式中：z為衛(wèi)星天頂距，rad；p為大氣壓，hPa；T為開氏溫度，K；e為水氣分壓，hPa。

式中：a0、a1、a2分別為衛(wèi)星鐘的鐘差、鐘速、鐘漂。

式（6）可簡寫為如下公式[8]：

1個(gè)歷元接收到n顆北斗衛(wèi)星，可列出n-1個(gè)誤差方程，式（9）中，常數(shù)項(xiàng)矩陣為L(n-1)×1，系數(shù)陣為B(n-1)×3，改正數(shù)矩陣為V(n-1)×1，未知數(shù)矩陣為x3×1。構(gòu)建誤差方程式后，根據(jù)最小二乘原理，可計(jì)算未知數(shù)x，具體公式如下[9]：

由于測站初值設(shè)為(0,0,0)，通過式（10）計(jì)算得出的x還不是最終結(jié)果，需進(jìn)行迭代運(yùn)算，因不知道迭代次數(shù)，C#程序中可采用while循環(huán)語句，循環(huán)結(jié)束的判別標(biāo)志可令未知數(shù)x＜1cm，根據(jù)實(shí)測數(shù)據(jù)驗(yàn)證，一般情況下迭代次數(shù)不超過5次即可滿足條件。

式（10）的權(quán)陣P可通過衛(wèi)星高度角不同進(jìn)行定權(quán)，高度角越大權(quán)值越高。式（6）的觀測值為B1和B2頻率下的偽距觀測值，B1、B2兩個(gè)頻率的偽距值是獨(dú)立的。觀測值組合中誤差可設(shè)為σ，其取值與衛(wèi)星高度角有關(guān)。觀測值方差陣D可用如下公式表示：

權(quán)陣P可通過方差陣D求逆運(yùn)算得到，代入式（10）可計(jì)算未知參數(shù)x。

3 BDS和GPS組合單站單歷元定位的數(shù)學(xué)模型

在BDS單衛(wèi)星系統(tǒng)下單歷元單站定位的數(shù)學(xué)模型中采用星間單差、雙頻觀測值無電離層組合的模式進(jìn)行計(jì)算。若采用BDS和GPS組合衛(wèi)星系統(tǒng)，則星間單差計(jì)算方式有兩種，一種是采用GPS和BDS各系統(tǒng)選擇自己的參考衛(wèi)星分開組合星間單差；另一種是僅選擇一個(gè)參考衛(wèi)星組混合單差[10]?？紤]到組混合單差需要引入系統(tǒng)間偏差，算法較為復(fù)雜，故采用第一種計(jì)算方式，即BDS和GPS分別各選擇系統(tǒng)內(nèi)的一顆衛(wèi)星作為參考衛(wèi)星，按照上述公式進(jìn)行計(jì)算。由于GPS雙頻觀測值（不采用GPS第3個(gè)民用頻率L5的數(shù)據(jù)進(jìn)行計(jì)算，因L5觀測值目前部分GPS衛(wèi)星還沒有，僅采用L1、L2兩個(gè)頻率）L1頻率值為1575.42MHz，L2頻率值為1227.60MHz，易得GPS兩個(gè)頻率的波長關(guān)系，77λ1=60λ2。將BDS單站定位計(jì)算公式中式（6）、式（11）、式（12）中的系數(shù)763替換成77，590替換成60，即可得到GPS單站定位的公式。

現(xiàn)將兩者組合，設(shè)1個(gè)歷元共接收n顆北斗衛(wèi)星、m顆GPS衛(wèi)星的觀測數(shù)據(jù)，可按照上述公式寫出組合衛(wèi)星系統(tǒng)的誤差方程式，具體公式如下：

除未知參數(shù)x不變外，其余矩陣的行數(shù)有所擴(kuò)充，北斗衛(wèi)星觀測值可列出n-1個(gè)方程，GPS衛(wèi)星觀測值可列出m-1個(gè)方程，組合形成n+m-2個(gè)方程。

此時(shí)方差陣D可用分塊矩陣表示，具體公式如下：

式中：DBDS為BDS觀測值的方差陣，采用式（12）計(jì)算；DGPS為GPS觀測值的方差陣，將式（12）的北斗雙頻組合系數(shù)換成GPS的組合系數(shù)即可；x為未知參數(shù)，按式（10）計(jì)算。

4 BDS單站多歷元定位的數(shù)學(xué)模型

文章采用單站定位的數(shù)學(xué)模型未知參數(shù)為3個(gè)，即測站點(diǎn)的未知數(shù)(dx,dy,dz)；而傳統(tǒng)單站定位的數(shù)學(xué)模型未知參數(shù)一般有4個(gè)，包括接收機(jī)鐘差。文章通過星間單差（其他衛(wèi)星和參考衛(wèi)星觀測方程作差）消去了接收機(jī)鐘差，不將其作為未知參數(shù)，由于是靜態(tài)觀測，測站點(diǎn)坐標(biāo)不會(huì)改變，而接收機(jī)鐘差可能會(huì)發(fā)生變化，故這里消去了接收機(jī)鐘差，再利用多歷元觀測數(shù)據(jù)計(jì)算單站定位的坐標(biāo)時(shí)，便于采用法方程疊加的原理進(jìn)行計(jì)算，具體公式如下：

式中：Xk-1為前k-1個(gè)歷元計(jì)算的測站點(diǎn)坐標(biāo)值；Xk為k個(gè)歷元計(jì)算的測站點(diǎn)坐標(biāo)值；Bi、Pi、Li分別為第i個(gè)歷元所形成的系數(shù)陣、權(quán)陣和常數(shù)項(xiàng)矩陣，不同歷元觀測的北斗衛(wèi)星數(shù)量可能不同，因此各歷元的系數(shù)陣、權(quán)陣、常數(shù)項(xiàng)矩陣的行數(shù)可能不同。

式（15）是單獨(dú)BDS多歷元單站定位的數(shù)學(xué)模型，通過單歷元GPS和BDS組合衛(wèi)星系統(tǒng)的單站定位計(jì)算模型，易知GPS和BDS組合衛(wèi)星系統(tǒng)下多歷元單站定位的計(jì)算公式，亦采用法方程累加的原理，構(gòu)建的各矩陣由于觀測方程的增加而需相應(yīng)擴(kuò)充。

當(dāng)計(jì)算出測站點(diǎn)坐標(biāo)時(shí)，精度指標(biāo)也可同時(shí)計(jì)算得出。即可通過相應(yīng)公式計(jì)算單位權(quán)中誤差σ0，求出測站在N、E、U方向的點(diǎn)位精度MN、ME、MU，平面精度因子HDOP值，高程精度因子VDOP值及空間精度因子PDOP值等精度指標(biāo)。點(diǎn)位精度MN、ME、MU表征計(jì)算得到的測站點(diǎn)各方向上的精度，DOP值表征衛(wèi)星分布的合理性，DOP值與測站和觀測衛(wèi)星構(gòu)成的單位多面體的體積成反比，DOP值越小，則相應(yīng)的解算精度越高[8]。

5 實(shí)測數(shù)據(jù)的坐標(biāo)及精度分析

利用2019年5月在山東菏澤地區(qū)某控制點(diǎn)觀測約70min的靜態(tài)數(shù)據(jù)，采樣間隔10s，高度截止角設(shè)為15°，共采集357個(gè)歷元數(shù)據(jù)。對(duì)每個(gè)歷元進(jìn)行單站定位計(jì)算，只采用BDS觀測值，計(jì)算得到的測站空間坐標(biāo)X、Y、Z坐標(biāo)如圖1～圖3所示。

由圖1～圖3可知，單獨(dú)利用BDS偽距觀測值進(jìn)行單站定位計(jì)算，X坐標(biāo)變化幅度為7.01m，（X最大值為-2291357.8834m，最小值為-2291364.9011m），Y坐標(biāo)變化幅度為10.21m，（Y最大值為4670024.9080m，最小值為4670014.6962m），Z坐標(biāo)變化幅度為4.52m（Z最大值為3678370.6212m，最小值為3678366.1014m）。其中Y坐標(biāo)變化幅度最大，這是因?yàn)橹袊蟛糠值貐^(qū)Y坐標(biāo)（空間直角坐標(biāo)）主要影響高程方向，衛(wèi)星定位高程方向的精度一般比平面精度低。各歷元測站X坐標(biāo)值大部分落在-2291360～-2291362m，各歷元測站Y坐標(biāo)值大部分落在4670017～4670021m，各歷元測站Z坐標(biāo)值大部分落在3678367～3678369m。結(jié)合圖1～圖3分析，測站定位坐標(biāo)在第150～180個(gè)歷元時(shí)間段內(nèi)突變較大，其余時(shí)刻測站坐標(biāo)各分量變化不大。參考RINEX數(shù)據(jù)O文件給出的測站近似坐標(biāo)，與BDS各歷元計(jì)算的單站定位坐標(biāo)平均值進(jìn)行較差比較，如表1所示。

由表1可知，O文件提供的近似坐標(biāo)與利用BDS觀測值數(shù)據(jù)計(jì)算各歷元單站定位坐標(biāo)進(jìn)行比較，坐標(biāo)差別不大，各分量較差均在1m以下。

現(xiàn)將單獨(dú)利用BDS和單獨(dú)利用GPS計(jì)算測站各歷元單站定位坐標(biāo)進(jìn)行較差比較，如圖4～圖6所示。

圖1 BDS偽距觀測值計(jì)算每個(gè)歷元測站點(diǎn)X坐標(biāo)（單位：m）

圖2 BDS偽距觀測值計(jì)算每個(gè)歷元測站點(diǎn)Y坐標(biāo)（單位：m）

圖3 BDS偽距觀測值計(jì)算每個(gè)歷元測站點(diǎn)Z坐標(biāo)（單位：m）

表1 RINEX數(shù)據(jù)的O文件測站近似坐標(biāo) 單位：m

由圖4～圖6可知，單獨(dú)采用BDS系統(tǒng)和單獨(dú)采用GPS系統(tǒng)進(jìn)行單站定位，各歷元計(jì)算測站點(diǎn)的坐標(biāo)較差差別不大，大部分歷元X坐標(biāo)較差在2m以下，Y坐標(biāo)較差在4m以下，Z坐標(biāo)較差在3m以下。將BDS和GPS各歷元單站定位坐標(biāo)取平均值進(jìn)行比較，如表2所示。

由表2可知，在分別利用BDS和GPS系統(tǒng)計(jì)算各歷元單站定位坐標(biāo)平均值時(shí)，分量較差均在2m以下。在BDS系統(tǒng)下單站定位在N、E、U方向的點(diǎn)位精度MN、ME、MU如圖7～圖9所示。

由圖7～圖9可知，BDS系統(tǒng)下單站定位精度為米級(jí)，大部分歷元的E方向的ME和N方向的MN精度值在1.5m以下。該實(shí)例數(shù)據(jù)的E方向點(diǎn)位精度最高，U方向點(diǎn)位精度最低，在第170～180個(gè)歷元計(jì)算U方向的精度最低，MU出現(xiàn)峰值。從圖1～圖3各歷元單站定位計(jì)算的坐標(biāo)各分量X、Y、Z值也可以看出，在第170～180個(gè)歷元區(qū)段內(nèi)，計(jì)算坐標(biāo)值有明顯差異，這與圖9歷元U方向上的點(diǎn)位精度MU變化特征是相符的。

由于單站定位的精度在米級(jí)，文章也只采用了廣播星歷計(jì)算衛(wèi)星坐標(biāo)，再按照數(shù)學(xué)模型進(jìn)行單站定位計(jì)算。從BDS和GPS計(jì)算的精度指標(biāo)結(jié)果分析，單獨(dú)采用BDS和GPS系統(tǒng)進(jìn)行單站定位精度相當(dāng)，均能達(dá)到米級(jí)[11-13]，若采用BDS和GPS混合衛(wèi)星系統(tǒng)進(jìn)行單站定位計(jì)算，由于廣播星歷計(jì)算的衛(wèi)星坐標(biāo)精度約為2m[13]，單站定位精度也只能達(dá)到米級(jí)，但在個(gè)別時(shí)段衛(wèi)星數(shù)量較少時(shí)進(jìn)行單站定位解算，此時(shí)混合衛(wèi)星系統(tǒng)下單站定位的解算精度能得到更好的保證。若需要提高定位精度，達(dá)到厘米級(jí)的定位精度[14]，衛(wèi)星星歷需采用精密星歷，觀測值需采用載波相位觀測值[15]，并準(zhǔn)確計(jì)算整周模糊度，各種誤差分別采用各自的誤差改正模型進(jìn)行改正計(jì)算。

圖4 BDS和GPS計(jì)算各歷元測站點(diǎn)X坐標(biāo)較差（單位：m）

圖5 BDS和GPS計(jì)算各歷元測站點(diǎn)Y坐標(biāo)較差（單位：m）

圖6 BDS和GPS計(jì)算各歷元測站點(diǎn)Z坐標(biāo)較差（單位：m）

表2 BDS和GPS各歷元單站定位坐標(biāo)的平均值較差比較表 單位：m

圖7 BDS計(jì)算各歷元單站定位的N方向的點(diǎn)位精度MN（單位：m）

圖8 BDS計(jì)算各歷元單站定位的E方向的點(diǎn)位精度ME（單位：m）

圖9 BDS計(jì)算各歷元單站定位的U方向的點(diǎn)位精度MU（單位：m）

6 結(jié)論

文章對(duì)BDS利用偽距觀測值進(jìn)行單站定位的數(shù)學(xué)模型進(jìn)行了介紹，并通過C#編程實(shí)現(xiàn)計(jì)算過程。通過實(shí)例數(shù)據(jù)分析，得出以下結(jié)論：（1）BDS和GPS單獨(dú)采用一種衛(wèi)星系統(tǒng)進(jìn)行單站定位的精度相當(dāng)，精度都能達(dá)到米級(jí)。若采用BDS、GPS混合衛(wèi)星系統(tǒng)進(jìn)行單站定位計(jì)算，由于廣播星歷計(jì)算衛(wèi)星坐標(biāo)的精度約為2m，因此BDS、GPS混合衛(wèi)星系統(tǒng)下單站定位的坐標(biāo)精度仍為米級(jí)。（2）采用文章的數(shù)學(xué)模型僅使用BDS進(jìn)行單站定位計(jì)算，將計(jì)算結(jié)果與RINEX的O文件中提供的近似坐標(biāo)進(jìn)行比較，文章實(shí)例數(shù)據(jù)計(jì)算的坐標(biāo)較差在1m以下。BDS單站定位計(jì)算的點(diǎn)位精度，平面的點(diǎn)位精度MN、ME高于高程方向的點(diǎn)位精度MU。（3）采用多歷元靜態(tài)觀測數(shù)據(jù)進(jìn)行單站定位計(jì)算時(shí)，由于數(shù)學(xué)模型消去了接收機(jī)鐘差未知參數(shù)，僅保留測站的3個(gè)未知參數(shù)，在各歷元測站點(diǎn)進(jìn)行靜態(tài)觀測并保持不動(dòng)，便于采用法方程累加的原理計(jì)算多歷元單站定位的測站坐標(biāo)。