顧及系統偏差的BDS鐵路邊坡監測方法

何義磊,張云龍,張冠軍,胡錦民,陳旭升

(1. 中國鐵路設計集團有限公司,天津 300308; 2. 天津市軌道交通導航定位及時空大數據技術重點實驗室,天津 300251)

當前我國鐵路事業正迅猛發展,截至2021年末,全國鐵路營業總里程已突破15萬km,鐵路工程基礎設施建設和運營維護的需求日漸增加。我國幅員遼闊,山地眾多,特別是西北黃土高原,土質松軟,在極端天氣或施工擾動情況下,滑坡和泥石流等突發地質災害嚴重威脅鐵路工程施工建設與運營,邊坡的穩定性對地面設施、人身安全及工程進度影響巨大[1]。受內外多種不確定性因素的影響,運動規律難以掌控,除了采取必要的防護措施外,還應進行實時監測,以掌握邊坡變化趨勢并及時預警,為鐵路工程建設與維護提供技術支撐和安全保障[2]。

GNSS具有自動化程度高、全天候觀測、測點間無需通視、可快速獲取高精度三維信息等特點[3],被廣泛應用于邊坡、橋梁、大壩等大型工程的變形監測[4],可滿足西北部山區惡劣的工程環境和鐵路遮擋環境下快速高精度的監測需求。由于GPS成熟較早,諸多學者已通過靜態相對定位技術、RTK技術、精密單點定位(precise point positioning,PPP)技術等多種數據處理方式,對GPS變形監測技術進行研究,并取得了良好效果[5-8]。靜態相對定位和RTK技術均需基準站提供差分信息或進行雙差。然而在土質松軟地區,基準站的穩定性將嚴重影響定位精度,選擇可獨立觀測且定位精度相當的PPP技術對獲取高精度位置信息具有重要意義。文獻[7]采用PPP技術對高層建筑的傾斜和沉降進行長期監測,通過穩定的局部參考框架和當地季節性地面變形模型相結合的方式得出,華北地區PPP平面和高程精度分別達到2～3 mm和6～9 mm(24 h連續觀測)。文獻[8]研究了復雜艱險地區的鐵路場景下,顧及相位偏差的GPS PPP定位性能,結果表明水平和高程方向定位精度分別優于10、15 cm(90 min連續觀測)。

北斗衛星導航系統(BeiDou Navigation Satellite System,BDS)于2020年7月31日全面建成,現階段處于北斗二號(BDS-2,15顆)/北斗三號(BDS-3,30顆)共存的局面,并向全球用戶提供高精度的導航、定位、授時服務,其衛星數量、信號質量、定位性能均顯著提升[9-11]。BDS的廣泛應用和GNSS接收機的國產化是必然趨勢,然而有研究表明多系統組合PPP中存在系統偏差(inter-system bias, ISB)。文獻[12]基于無電離層組合解算GPS+BDS ISB,發現ISB大小與接收機類型相關,大小為10～100 ns;文獻[13]針對7 d的MGEX(multi-GNSS experiment)跟蹤網觀測數據,基于GPS+BDS組合PPP模型解算GPS、BDS短期ISB,并構建二次多項式加周期函數的短期ISB預報模型,建議預報時長為1 d,利用預報ISB作為先驗信息進行組合PPP,發現在N、E、U方向上的定位精度和收斂效率有所提升。文獻[14]針對MGEX跟蹤網2014—2016年各7 d的數據,基于BDS+GPS組合PPP模型分別解算4個跟蹤站的ISB值,以分析ISB的長期穩定性,結果表明ISB序列的天變化較為穩定,而3年的ISB周平均值和標準差(standard deviation,STD)具有顯著差異。因此,為提升遮擋環境下的定位性能,需削弱ISB的影響。

本文以西北地區某鐵路隧道南側山體的邊坡北斗監測工程為例,詳細介紹利用BDS PPP技術進行邊坡變形監測的方法,并對BDS-2、BDS-3間ISB進行估計、分析和建模,將預報ISB作為先驗約束加入BDS-2+BDS-3組合PPP中,以檢驗PPP定位精度和收斂時間的提升水平,并對案例中的監測點進行穩定性分析。

1 顧及ISB的BDS PPP基本理論

圖1為2021年12月1日BDS在軌運行衛星星下點軌跡。BDS-3衛星在BDS-2衛星播發B1I和B3I信號的基礎上,新增觀測數據質量更佳的B1C和B2a信號[15],同時應用了新的衛星姿態模型、新型原子鐘等技術。因此本文將兩者視為兩個不同的衛星導航系統,重點分析BDS-3衛星B1C和B2a信號與BDS-2衛星組合PPP的性能水平。

圖1 2021年12月1日BDS在軌運行衛星星下點軌跡

1.1 顧及ISB的GPS+BDS-2+BDS-3組合PPP模型

接收機R與衛星S在原始偽距和載波相位觀測值的觀測方程為[16]

δion+δrel+δmult+εR(P)

(1)

δion+δrel+δmult+εR(Φ)

(2)

式中,S為GPS+BDS-2+BDS-3衛星;ρ為衛星和接收機天線相位中心間的幾何距離;c為真空中的光速;dtR為接收機鐘差;dtS為衛星鐘差;bR、bS、Br、Bs分別為接收機和衛星的硬件延遲;δion為電離層延遲誤差;δtrop為對流層延遲誤差;δrel為相對論效應誤差;δmult為多路徑效應誤差;N為載波相位整周模糊度;ε(P)和ε(Φ)分別為偽距和載波相位觀測噪聲等未模型化的參數[16]。

基于雙頻偽距和載波相位觀測值的無電離層組合,式(1)—式(2)可分別轉換為

(3)

(4)

利用精密衛星軌道和鐘差,式(3)—式(4)可分別轉換為

(5)

(6)

(7)

ISB參數(ISBG,B2和ISBG,B3)可與其他參數一同被估計,兩者相減可得到BDS-2和BDS-3間的系統偏差ISBB2,B3,即

ISBB2,B3=ISBG,B3-ISBG,B2

(8)

在ISB被忽略或已知時,可通過式(1)—式(7)實現GPS+BDS-2+BDS-3組合PPP。

1.2 基于N次多項式+ARIMA的ISB模型

針對ISBB2,B3序列中未含有周期項,提出一種基于N次多項式+自回歸移動平均模型(autoregressive integrated moving average model,ARIMA)的組合建模方法[16]。利用N次多項式對原始ISB＿O序列進行最小二乘擬合[16],即

ISB＿Ft=antn+an-1tn-1+…+a1t+a0

(9)

式中,ISB＿Ft表示第t個歷元ISBB2,B3的擬合值;an和a0為N次項和常量項的系數。通過計算ISB＿O與ISB＿F殘差序列的RMS,確定最佳擬合次數。

利用ARIMA對殘差序列進行建模,ARIMA(p,d,q)模型可表示為[17]

(10)

式中,ISB＿Rt為ISB＿Ot和ISB＿Ft之差;A、B和p、q分別為ARIMA模型的系數和階數;v為誤差項,{vi}～WN(0,δ2),WN為白噪聲。

ISB＿Pt+1預報值可表示為

ISB＿Pt+1=a8(t+1)8+…+a1(t+1)+a0+

(11)

式中,各項系數均可由式(9)、式(10)得到。為了減弱歷史數據的影響,采用滑動窗口的策略。

2 BDS-2、BDS-3間ISB特性分析及模型構建

2.1 數據獲取與處理

選取2021年12月1—31日(年積日335—365)西北地區某鐵路隧道南側山體高邊坡BJ01和BC01北斗監測站數據,接收機類型為SEPT MOSAIC-X5C,天線類型為HGGCYH8373 HGGS,可同時接收GPS、BDS-2、BDS-3信號,采樣間隔30 s。圖2為BJ01和BC01北斗監測站示意圖,其中BJ01建在穩定的開闊區域,BC01建在高邊坡頂端,周邊植被叢生,遮擋嚴重。該邊坡為高山峽谷地貌,地形陡峭,坡度一般為30°～50°。邊坡多為泥石流堆積體和發育滑坡,南側受336國道施工擾動,為避免滑坡威脅鐵路安全,需對其進行自動化實時監測。

圖2 BJ01和BC01北斗監測站示意

2.2 BDS-2+BDS-3組合PPP中ISB特性分析

對估計的BDS-2、BDS-3間系統偏差ISBB2,B3日穩定性進行分析,探究其短期變化規律。圖3為BJ01和BC01監測站年積日335—348為期14 d的原始ISBB2,B3序列。

圖3 BJ01和BC01監測站ISBB2,B3時序

由圖3可以看出,BJ01和BC01北斗監測站ISBB2,B3在0～700 ns之間波動,一天內變化較為平穩,呈多條連續且平順的曲線,但天與天之間存在跳變,呈階梯狀。分別計算兩個監測站單日的ISBB2,B3平均值和標準差,其中年積日344時兩站單日ISBB2,B3平均值均為最大,分別為660、658 ns,14 d的ISBB2,B3平均值分別470、467 ns,因此ISBB2,B3的存在已嚴重影響到PPP定位精度。而14 d的ISBB2,B3標準差的平均值均為10 ns,表明ISBB2,B3在一天內的變化較為穩定。兩站天與天之間的ISBB2,B3存在跳變,且表現出良好的一致性,該跳變是由估計衛星精密鐘差時更換參考星所致[16]。因此需要先消除跳變的影響,跳變修復的公式為

(12)

Δd=ISB(Ti)-ISB(Ti-1)

(13)

圖4為跳變修復后的ISBB2,B3序列,可以看出,原序列中存在的跳變被修復,但由某些未得到修復的小跳變的積累導致年積日343后BJ01和BC01站的ISB序列產生明顯差異。

圖4 BJ01和BC01監測站修復跳變后ISBB2,B3時序

2.3 基于N次多項式+ARIMA的BDS-2、BDS-3間ISB模型

利用快速傅里葉變換(fast Fourier transform, FFT)對修復后的ISBB2,B3序列進行頻譜分析。圖5為BJ01和BC01站ISBB2,B3序列的功率譜密度(power spectral density,PSD)。可以看出,只有一個明顯信號頻率為零,表明序列ISBB2,B3不含有周期項。

圖5 監測站修復跳變后ISBB2,B3的PSD

為得到最佳的多項式系數,利用最小二乘法估計BJ01和BC01站1～100次多項式擬合的系數。次數為8時,擬合系數最佳。采用八次多項式(octave polynomial,OP)對原始ISBB2,B3序列進行擬合。表1為年積日308—321,BJ01和BC01站的基于OP擬合的各次項和常數項系數。圖6為BJ01和BC01站的擬合結果。

表1 年積日335—348時BJ01和BC01北斗監測站的ISBB2,B3八次多項式擬合系數

圖6 監測站八次多項式擬合ISBB2,B3結果

采用ARIMA對ISBB2,B3原始序列ISB＿O和八次多項式擬合的序列ISB＿P之間的殘差序列進行建模。使用Eviews 9.0軟件構建ARIMA(p,d,q)模型,得到BJ01和BC01站的ARIMA(p,d,q)模型參數、殘差方差、AIC和SC信息準則值,見表2。圖7為年積日335—348時BJ01和BC01站的ISBB2,B3殘差序列擬合結果,ISB＿f表示擬合值。

表2 BJ01和BC01站的ARIMA模型參數和結果

圖7 ISBB2,B3殘差擬合結果

本文以平均絕對誤差(mean absolute error,MAE)和均方根誤差(root mean square error,RMSE)作為ISB建模與預測的精度指標。MAE表示為

(14)

式中,ISB＿Ot、ISB＿Ft、ISB＿ft分別為t時刻ISB的原始估計值、八次多項式、ARIMA模型擬合值。評估ISB的預報精度時用ISB＿Pt代替(ISB＿Ft+ISB＿ft)擬合值。

圖8為基于OP+ARIMA的ISBB2,B3擬合結果。可以看出,ISB＿F與ISB＿O非常吻合,殘差序列ISB＿R較為穩定,大多數情況下近似為0。BJ01和BC01站擬合殘差的MAE分別為0.17、0.19 ns,RMSE分別為0.50、0.54 ns,表明本文ISB建模方法具有較高的擬合精度。

圖8 北斗監測站ISBB2,B3擬合結果

為了驗證OP+ARIMA模型的預報精度,利用以年積日335—348兩周的原始ISB＿O序列為基礎構建的模型,預報年積日349的ISB＿P,并與同期估計的ISB＿O進行比較。圖9為年積日349時BJ01和BC01站ISBB2,B3的預報結果。可以看出,預報早期ISB＿O與ISB＿P具有較好的一致性,但隨著時間的推移,預報殘差ISB＿R變大,準確性有所降低。兩站ISB＿R的MAE分別為1.322、2.365 ns,RMSE分別為1.623、3.145 ns,表明本文ISB模型具有較高的精度。為了保證預報的ISB精度,建議以14 d為時長建模,預報時長為1 d,若進行長時間預報,滑動窗口大小為14 d。

圖9 年積日349北斗監測站ISBB2,B3預報結果

2.4 顧及ISB約束的BDS PPP定位性能分析

為了驗證OP+ARIMA模型ISB預報值的可用性,針對年積日349時BJ01和BC01北斗監測站作靜態PPP試驗,根據兩站的地理位置,可分別模擬開闊環境下和遮擋環境下的PPP驗證;同時將上述預報的ISB視為先驗約束,ISB預報殘差的RMSE作為先驗精度。按如下兩種方案進行數據處理,并比較其在N、E、U方向的收斂時間和定位精度。收斂時刻以定位誤差小于0.1 m的時刻,且其后20個時刻的誤差均滿足要求為準則[13]。

(1)GPS、BDS-2和BDS-3組合PPP時,不考慮之間的ISB,將兩者視為同一系統。

(2)GPS、BDS-2和BDS-3組合PPP時,將預報的BDS-2、BDS-3間系統偏差ISB＿P作為先驗約束。

圖10為BJ01北斗監測站0～3 h的PPP收斂時序。可以直觀地看出,N、E、U方向均在短時間內達到收斂要求,隨后曲線保持平穩,在附加ISB先驗約束后收斂時間縮短,N方向優于E和U方向。表3和圖11為BJ01和BC01北斗監測站在顧及ISB約束前后的定位精度和收斂時間的對比。可以看出,在增加預報ISB作為約束條件后,BJ01站在N、E、U方向上的定位精度分別提升了6.25%、4.23%、4.66%,收斂時間分別縮短了5.71%、7.14%、11.11%;BC01站在N、E、U方向上的定位精度分別提升了9.24%、16.02%、8.46%,收斂時間分別縮短了9.30%、15.59%、15.00%。

表3 BJ01和BC01站在兩種方案下的定位精度和收斂時間對比 (%)

圖10 年積日349 BJ01北斗監測站PPP收斂時序

圖11 年積日349 BJ01和BC01北斗監測站顧及ISB約束的PPP定位精度和收斂時間

結果表明,在觀測環境較好的情況下增加ISB的約束可提升定位精度和解算效率,而在衛星數較少或現場存在遮擋的情況下增加ISB的約束,可顯著提升定位精度和解算效率,在鐵路高邊坡自動化監測中具有良好的適用性。

2.5 監測點變形分析

利用本文顧及ISB約束的BDS-2+BDS-3組合PPP,對BJ01和BC01北斗監測站在年積日335—365間的點位進行估計,分析監測點的變形情況。以年積日328—334(7 d)的坐標平均值為參考值,圖12為BJ01和BC01兩監測站的坐標變化。由于靜態相對定位技術可消除衛星端、接收機端、大氣等誤差影響,而受多路徑影響較大,以此驗證顧及ISB約束的BDS-2+BDS-3組合PPP的定位精度。圖12為BJ01和BC01北斗監測站靜態相對定位結果及其與PPP的比較。

圖12 年積日335—365 BJ01和BC01北斗監測站坐標變化

由圖12可知,BJ01北斗監測站在年積日335—365一個月內基本穩定,無異常形變,在N、E、U方向的最大形變量分別為4.10、7.05、7.01 mm,RMSE分別為1.85、3.08、4.13 mm。對于BC01北斗監測站,在年積日352后E和U方向發生明顯形變,年積日365時N、E、U方向的形變量分別為0.97、-22.33、-19.17 mm,年積日335—365內RMSE分別為2.60、11.10、9.64 mm。由于BJ01站基本穩定,BC01和BJ01站間靜態相對定位結果即為BC01站的變形量。如圖13所示,靜態相對定位結果與BC01站PPP結果表現一致,在年積日352后E和U方向發生明顯變形,年積日365時N、E、U方向的變形量分別為0.30、-21.40、-17.70 mm。靜態相對定位結果與BC01和BJ01站間PPP坐標差在N、E、U方向的RMSE分別為1.26、2.11、2.22 mm。表明本文方法適用于遮擋環境下鐵路邊坡長期靜態變形監測,監測精度可達毫米級。

圖13 年積日335—365 BJ01和BC01北斗監測站相對定位

3 結 論

本文詳細介紹了顧及BDS-2、BDS-3間ISB的鐵路邊坡自動化監測方法;以某鐵路隧道邊坡工程監測為例,在PPP中估計了BDS-2、BDS-3間的ISB,并分析了ISB短期特性;提出了八次多項式+自回歸移動平均模型的ISB建模方法,驗證了模型精度和該方法對PPP定位性能的改進;分析了BJ01和BC01監測點為期一個月的穩定性,主要結論如下。

(1)提出了GPS+BDS-2+BDS-3組合PPP估計ISB的解算方法,分析了BDS-2和BDS-3間的系統偏差ISBB2,B3的短期穩定性,其14 d的平均值量級約為470 ns,但單日內變化較小,天與天之間存在數百納秒的跳變,在BDS-2(B1I,B3I)和BDS-3(B1C,B2a)組合PPP中,需考慮其影響。

(2)分析表明ISB序列中不含周期項,提出OP+ARIMA的ISB建模方法,以為期14 d的ISB序列進行建模,在預報時長為1 d時,該方法具有較高的擬合精度和較好的一致性,預報精度隨時間推移有所降低。

(3)將預報的ISB作為先驗約束,應用于BDS-2+BDS-3組合PPP中,可在復雜的觀測環境下提升解算效率和定位精度。

(4)對于BJ01北斗監測站在年積日335—365一個月內基本穩定,BC01北斗監測站在年積日352后發生明顯形變,與靜態相對定位結果相近,監測精度可達毫米級,能滿足西北山區鐵路邊坡長期靜態變形監測的需求。