多模多頻GNSS-IR水位反演中的頻間偏差分析及改正

王笑蕾，何秀鳳，宋敏峰，陳 殊，牛紫瑾

河海大學地球科學與工程學院，江蘇 南京 211100

準確的水位監(jiān)測對于水資源調控、水災監(jiān)控及氣候氣象研究十分重要。傳統(tǒng)水位計監(jiān)測技術存在造價較高或需人工維護的缺點，而且需要附加技術才能實現基準統(tǒng)一[1]。同時，衛(wèi)星測高技術無法獲得近岸水面高度且時間分辨率較低[2]。隨著全球導航衛(wèi)星系統(tǒng)(GNSS)的不斷發(fā)展與完善，一種GNSS干涉遙感(GNSS-interferometry reflectomety,GNSS-IR)技術被發(fā)現可以用來進行水位(包括潮位)[3-7]、雪深[8-11]、土壤濕度[12-14]及地表凍融[15]等環(huán)境參數的監(jiān)測。其中，GNSS-IR水位監(jiān)測技術只基于沿岸的測量型接收機，根據其信噪比(signal-to-noise ratio,SNR)或其他觀測數據中的干涉特性，便可完成水位監(jiān)測[3-4]。利用新興的GNSS-IR技術，可以實現全自動、小成本、高時間分辨率、長期連續(xù)的水位監(jiān)測，并且監(jiān)測結果自動固定在穩(wěn)定框架下[5]。

文獻[3—4]提出GNSS-IR技術，并由此發(fā)展出一套經典的GNSS-IR水位反演理論。目前，地基GNSS-IR技術主要有兩類誤差源：反射面高度變化和大氣折射效應。文獻[4]發(fā)現海面起伏會引起反演結果中的高度變化誤差，并提出了相應的改正算法。文獻[16]利用5個GNSS站進行了潮位反演，并且利用經驗潮波系數來改進誤差改正算法。文獻[17]基于潮波分析模型，提出了一種不需要經驗潮波系數的海面高度擬合方法，進一步完善了高度變化誤差的改正算法。大氣折射效應則會引起兩類誤差：一類是由于大氣折射引起的信號彎曲帶來的高度角彎曲誤差；另一類是由于大氣折射引起的對流層延遲誤差。文獻[18]針對大氣折射引起的高度角彎曲誤差，利用折射改正公式來改正高度角偏差。文獻[19]利用VMF1映射函數模型和GPT2w對流層延遲模型來改正對流層延遲誤差。

2020年6月23日，北斗三號全球導航衛(wèi)星系統(tǒng)星座建設完成。同時，GPS和GLONASS正在推進現代化，Galileo在不斷發(fā)展，各個區(qū)域系統(tǒng)也在持續(xù)改進和建設。多系統(tǒng)發(fā)展使得GNSS-IR技術擁有了豐富的信號源，而四系統(tǒng)的多模多頻信號已經均被證明可以監(jiān)測水位。文獻[3]證明了GPS信號的水位可監(jiān)測性。文獻[10]證實了GLONASS系統(tǒng)的水位可監(jiān)測性。文獻[21]證明了北斗數據的水位可監(jiān)測性。文獻[22]利用GPS、GLONASS和北斗信號進行了水位監(jiān)測試驗。文獻[23]利用GPS、GLONASS、Galileo和北斗系統(tǒng)進行了水位反演，分析了四系統(tǒng)的監(jiān)測性能，并提出了一套多模多頻融合反演算法。多模多頻信號的引入，提升了時間采樣，也為提高精度提供了可能。但是多模多頻數據也帶來了更多的誤差，需要對這些誤差進行正確認識及合理處理，以獲得最優(yōu)反演性能。

為了更好地挖掘GNSS-IR水位反演誤差，本文研究4個GNSS連續(xù)跟蹤站(HKQT、SW50、SW51和SW52)的多模多頻反演結果，發(fā)現除了上述已經發(fā)現的3類誤差，還存在一種明顯的頻間偏差。而目前對GNSS-IR頻間偏差的研究極少，尚未達成將其歸為GNSS-IR誤差源的共識。為了進一步挖掘頻間偏差的特性、加深對該誤差的認識，本文分析此類頻間偏差的特性，并提出相關的誤差改正方法。

1 GNSS-IR原理及系統(tǒng)誤差

1.1 GNSS-IR原理

當僅存在一次多路徑反射條件下，直射信號與反射信號之間的路徑差為D=2h(sine)[3-5]；其中，h為反射面與天線中心之間的垂直距離，稱為有效高度(reflector height,RH)；e為高度角。因此，二者之間的相位差Δφ為[3-5]

Δφ=2πD/λ=4πh(sine)/λ

(1)

式中，λ為信號波長。相位差中隱藏了一個頻率信息f

(2)

式中，ωφ為角頻率。如果忽略輸入數據對應時間中h的變化，即在反射面在一定時間內靜止的假設下，式(2)通過求導可寫為

(3)

1.2 GNSS-IR系統(tǒng)誤差

根據已有研究，本文匯總了已有共識的3類GNSS-IR系統(tǒng)誤差，并給出相關的誤差改正或削弱方法。

1.2.1 高度變化誤差

(4)

(5)

1.2.2 高度角彎曲誤差

大氣折射效應會引起信號彎曲，從而導致計算高度角與實際高度角之間的偏差，而這一高度角偏差又會進一步導致反演結果的偏差。文獻[18]發(fā)現該高度角彎曲誤差，并提出利用大氣折射改正公式(6)改正高度角偏差δe的方法

(6)

式中，T為溫度，單位為℃；P為氣壓，單位為mb。也可利用其他的大氣折射公式進行高度角偏差改正[27]。利用改正后的高度角進行反演計算，即可削弱此誤差。該誤差與高度角有關，高度角越小，誤差越大；該誤差會引起反演結果的尺度誤差及均值誤差。對于5°～20°高度角區(qū)間，反演誤差量級在亞厘米至毫米級，可視精度要求情況進行改正。

1.2.3 對流層延遲誤差

大氣折射效應會引起信號延遲，由于反射信號比直射信號傳播路徑更長，因而受到的延遲影響更大，從而造成反演結果的誤差。文獻[19,24—26]利用GPT2w模型和VMF1映射模型來修正該誤差。先計算反射信號和直射信號間的相對延遲值τT

(7)

(8)

式中，ΔhT為延遲誤差改正量。該誤差與RH大小及高度角有關，該誤差會引起反演結果的尺度誤差及均值誤差。一般情況下(有效高度10 m以內，最低高度角選擇5°)，影響量級在厘米級至毫米級；可視精度要求情況進行改正。對于某些RH較大的站點(例如安置于燈塔上部、高壩頂部、高樓頂部等的站點)或者極低高度角的情況，該誤差必須改正。

2 頻間偏差分析

2.1 站點介紹

本文選取了4個GNSS連續(xù)跟蹤站。其中3個站點(SW50、SW51和SW52)來自山東雙王城水庫的GNSS大壩監(jiān)測系統(tǒng)，另一個站點(HKQT)來自香港衛(wèi)星定位參考站網(satellite positioning reference station network，SatRef)。

雙王城水庫位于山東濰坊壽光市，水庫大壩軸線總長為9.636 km，壩高約12.5 m，最大庫容量達到6150萬m3，是南水北調東線膠東干線工程的重要調蓄水庫。為了監(jiān)測水庫大壩的穩(wěn)定性，在大壩周圍建立了3個監(jiān)測站(SW50、SW51和SW52)和1個基準站(SW43)，如圖1所示。該GNSS網配備CHC N72接收機及CHCC220GR天線，記錄每15 s一次的GPS(L1和L2)和北斗(B1、B2和B3)觀測值。相關信號對應的頻率、波長及SNR類型見表1。監(jiān)測站SW50、SW51和SW52均可接收到來自水庫水面的反射信號。本文選取5°～15°高度角區(qū)間SNR弧段進行反演。同時， 根據圖2中的反射區(qū)分布，選定對應的水域方位角為：SW50的方位角區(qū)域為60°～165°，SW51的方位角區(qū)域為130°～230°，SW52的方位角區(qū)域為270°～350°。算例選取2017-11-22—2018-01-10期間的數據進行分析。在SW50附近有一個水位監(jiān)測站，提供每天一次的水庫水位測量值。

圖1 雙王城水庫GNSS監(jiān)測網Fig.1 GNSS net of Shuangwangcheng Reservoir

圖2 雙王城站點信號反射區(qū)Fig.2 FFZs of sites of Shuangwangcheng Reservoir

HKQT站位于香港鲗魚涌，屬于香港衛(wèi)星定位參考站網，站點環(huán)境如圖3所示。HKQT站配有Trimble NetR9接收機及Trimble 59800.0型天線。本文使用其接收到的四系統(tǒng)——GPS(L1C/A、L2P、L2C和L5)、GLONASS(G1C、G1P、G2C和G2P)、Galileo(E1、E5、E7和E8)和北斗(B1和B2)5 s采樣的觀測數據進行反演。相關信號對應的頻率、波長及SNR類型見表1。算例選取5°～15°高度角區(qū)間SNR弧段進行反演。同時，根據圖4中的反射區(qū)分布，選定對應的海域方位角為-60°～105°。算例選取2017年DOY 224—DOY 244期間的數據進行分析。距離HKQT站點2m處有一隸屬Intergovern-mental Oceanographic Commission(IOC,http:∥www.ioc-sealevelmonitoring.org/)的驗潮站Quarry Bay可提供實測的潮位數據。

圖3 HKQT站點環(huán)境Fig.3 HKQT site surroundings

2.2 雙王城站點反演序列

圖4 HKQT站點信號反射區(qū)Fig.4 FFZs of HKQT

表1 使用信號對應的頻率、波長、碼及SNR類型

圖5中，L1、B1、B2及B3的反演與水位測量對應極好；而L2的反演序列表現出雙序列現象，這是由于L2信號在LSP譜圖中表現出的雙波峰現象導致[18]。由表2可知，L1、B1、B2及B3的反演精度很好，RMSE在1.74～11.63 cm，且大部分集中在3～6 cm。這種高精度結果主要是由于水庫表面平靜、反射條件好、SNR多徑特性明顯導致。表2中，對于3個監(jiān)測站，L2信號的常數比L1信號大0.19 m，B1的常數與L1相似，B2的常數比L1信號大0.20～0.21 m，B3的常數比L1信號大0.13～0.15 m；這意味著，L2信號的RH比L1信號整體小0.19 m，B1的RH與L1相似，B2的RH比L1信號小0.20～0.21 m，B3的RH比L1信號小0.13～0.15 m。比較各個信號間的偏差值及表1中各個信號的波長值，可以發(fā)現偏差大小與波長呈顯著線性關系(CORR=98.82%)。因此，偏差可表示為a×δλ，其中a為常數。根據擬合，雙王城站點3個站反演結果頻間偏差的a=3.6。這種偏差量級在分米級，明顯不是相位中心偏移或是隨機誤差導致；而且其呈現了明確的系統(tǒng)特性(與波長相關)。因此，可以確定有一類與波長相關的系統(tǒng)誤差確實存在。而波長λ等于光速與信號頻率的比值，即該偏差也與頻率呈現系統(tǒng)特性，故其為一種頻間偏差。這種頻間偏差可能是由電磁偏差、菲涅耳反射系數的組合相互作用及天線方向圖和兩個圓極化之間的功率偏移引起的[28]，這些硬件造成的誤差，理論上與頻率線性相關。

圖5 雙王城水庫站點反演序列 Fig.5 Retrievals of sites of Shuangwangcheng Reservoir

表2 雙王城站點RH反演結果統(tǒng)計表

2.3 香港鲗魚涌站點反演序列

經過與2.2節(jié)所述相同的處理后，獲得了HKQT站各個信號的RH值。根據1.1節(jié)所述原理，可以通過站點坐標及基準轉換參數獲得特定常數，然后通過常數減去RH就可以得到潮位反演值。IOC(Intergovernmental Oceanographic Commission,http:∥www.ioc-sealevel monitoring.org/)提供了HKQT站的三維坐標結果以及相應框架轉換參數，根據相關數據計算獲得用于估計潮位的常數值為7.797 m。然而，根據2.2節(jié)的分析，頻間偏差是存在的。因此，通過RH與潮位觀測值之間的對應關系，計算獲得對應的常數。在計算對應常數的過程中，發(fā)現不同的信號對應的常數不同，且差距較大。本文統(tǒng)計了不同信號對應的常數、反演個數、反演值與實測值間的RMSE和CORR。同時為了更好地表現信號間的頻間偏差，統(tǒng)計了各信號常數與L1信號常數之間的差距，相關結果記錄見表3。各個信號的反演序列及反演序列與實測序列的偏差如圖6所示。

表3 香港鲗魚涌站點RH反演結果統(tǒng)計表

圖6中，GPS、GLONASS、Galileo和北斗各個信號的反演序列與潮位實測序列對應良好。根據表3，HKQT站各個信號的反演精度在低分米級，RMSE在15～23 cm。表3中，波長與L1近似的信號，其對應的常數與L1信號的常數相同或近似；而與L2波長近似的信號，其對應的常數較L1的常數低0.12 m。這意味著，與L1波長近似的信號的RH與L1信號的RH近似，L2信號及與其波長近似的信號的RH比L1信號整體小0.12 m。比較表3中各個信號間的偏差值及表1中各個信號的波長值，可以發(fā)現：偏差大小與波長呈顯著線性關系(CORR=96.97%)。偏差可表示為a×δλ，其中a為常數。根據擬合，HKQT站點對應的a=2.3。

圖6 HKQT站點反演序列及偏差Fig.6 Retrievals and bias of HKQT

3 頻間偏差改正

3.1 誤差改正方法

(9)

(10)

(11)

式(10)或式(11)可簡寫為

(12)

3.2 改正結果

利用3.1節(jié)所述方法，對第2節(jié)獲得的各個信號的RH結果進行處理，獲得各個窗口的hi輸出結果，并將其轉換至水位基準，獲得對應的水位序列及其不確定度。值得注意的是，由于雙王城站點的GPS L2P信號有雙波峰現象，因此相關RH反演值不參與改正方程組的建立。本文相關站點誤差改正后的多模多頻融合水位反演結果如圖7所示。

圖7 各站誤差改正后的多模多頻融合水位反演序列Fig.7 Multi-GNSS combined retrievals with error correction of sites

經過誤差改正處理后，各個測站的水位反演結果與實測水位結果的對應關系較好(CORR>99%)。各個測站的誤差改正后得出多模多頻融合水位反演序列的RMSE，SW50為1.56 cm，SW51為0.67 cm，SW52為0.99 cm，HKQT為6.98 cm。與第2節(jié)所述的各個信號的反演結果的RMSE相比，改正后的融合反演值的精度提高了30%～80%。同時，本文計算了未顧及頻間偏差， 只顧及其他誤差的改正結果，得出相關的多模多頻融合水位反演序列的RMSE，SW50為3.15 cm，SW51為5.90 cm，SW52為2.24 cm，HKQT為18.97 cm。與顧及頻間偏差的改正結果相比，未顧及頻間偏差的反演結果的RMSE高1.5～12 cm。具體而言，SW50和SW52站點改正頻間偏差后提高了約1.5 cm，SW51站提高了約5 cm，而HKQT站提高了約12 cm。這是由于SW51站的北斗反演值個數顯著大于SW50和SW52站，因此，在反演結果改正中，頻間偏差對SW51的誤差影響更大；而HKQT站，由于其多系統(tǒng)多頻信號均可用，因此在反演結果改正中，頻間偏差對融合反演結果的誤差影響遠大于雙王城各站點。圖8顯示了HKQT站未顧及頻間偏差的改正結果，可以看出，未改正頻間偏差時，融合反演結果的精度變差、粗差增多。

圖8 HKQT站未顧及頻間偏差的多模多頻融合水位反演結果及其不確定度Fig.8 Multi-GNSS combined retrievals without consideration of inter-frequency bias and corresponding uncer-tainties at HKQT

4 結論與討論

本文分析了頻間偏差的相關特性，發(fā)現了頻間偏差與波長之間存在線性關系這一規(guī)律；同時，本文顧及高度變化誤差、對流層延遲誤差及頻間偏差，提出了相關的誤差改正方法。本文結果表明：顧及頻間偏差的改正結果比未顧及頻間偏差的反演結果的RMSE高1.5～12 cm；同時，改正后的融合反演值的精度較未改正的獨立信號的反演值提高了30%～80%；精度的提高得益于多模多頻信號提供的大量冗余數據及對各類誤差(包括系統(tǒng)誤差、粗差和隨機誤差)的正確處理。值得注意的是，改正方法中的相關方程組是經典平差函數模型的表達方式，可借鑒測量數據處理方法進行模型擴充和解算方法擴充；在有新的誤差被發(fā)現后，也可以根據誤差特性進行擴展。第3.1節(jié)所述誤差改正模型具有很好的擴展性。

除了系統(tǒng)誤差，本文還表現出了一類L2P信號獨有的反演結果特性。雙王城站點中，L2的反演序列表現出雙序列現象，這是由于L2信號在LSP譜圖中表現出的雙波峰現象導致[18]。而HKQT站則并未表現出L2P雙波峰現象導致的反演結果很差的情況。對于此類表現在L2P上的獨特情況——有時出現、有時不出現“雙波峰”的情況，其根本原因還沒有定論，猜測與P碼調制、L1干擾或對L2的半無碼跟蹤有關。這類獨特現象還需要更深一步的研究。

此外，本文只研究了GNSS-IR水位反演的情況，后續(xù)可以對雪深、冰厚等參數進行GNSS-IR反演，以探討在其他參數反演中是否存在頻間偏差。同時，對于形成頻間偏差的具體原因還未被發(fā)現；而該誤差成因對于更好地認識及改正頻間偏差非常重要。第3.1節(jié)中，利用a×δλ來表示頻間偏差大小，而表2和表3中，某些信號并不完全符合a×δλ的數值，有約1～3 cm的差距。文獻[28]指出，頻間偏差的產生可能與電磁偏差、菲涅耳反射系數的組合相互作用及功率偏移等接收機硬件有關，而這些硬件造成的誤差理論上與頻率線性相關。因此筆者認為，對于特定站點，a取決于接收機硬件，它對于不同信號是唯一的，相關差距應該是隨機誤差導致的。但是，該結論僅基于現有研究推斷而來，a與信號特征是否有關依然待研究。

目前，GNSS-IR技術仍處在快速發(fā)展階段，還有許多的誤差特性還未理清，例如地形引起的誤差[30]，不同高度角SNR質量引起的誤差[31]等，也有許多的誤差現象還沒有被發(fā)現成因。需要后續(xù)研究者對此進行深入研究，以實現更好的GNSS-IR監(jiān)測性能，推進GNSS-IR技術的實際應用進程。

致謝:感謝香港衛(wèi)星定位參考站網提供的GNSS數據，感謝IOC(Intergovernmental Oceano-graphic Commission, http:∥www.ioc-sealevel monitoring.org/)提供的實測潮位數據。