GNSS信噪比數據反演大壩水位的方法與應用

何騫，竇邵華

(廣州市城市規劃勘測設計研究院，廣東 廣州 510060)

1 引 言

隨著全球衛星導航系統(GNSS)的發展和研究的不斷深入，學者們發現GNSS接收機不僅接收來自衛星的直射信號，同時也會接收經天線周圍環境反射后的衛星信號，這兩種信號之間會發生干涉效應，從而影響了衛星導航定位的精度，即多路徑效應[1]。1993年，歐空局科學家Martin-Neira提出了全球定位系統(GPS)反射測量技術，開啟了GNSS反射測量這一全新研究領域[2]。Bilich等人對GPS信噪比(SNR)觀測值進行了分離，并研究了反射信號與反射環境之間的關系[3]。Larson等人提出并發展了GPS反射測量技術，其主要利用GPS信噪比數據對雪深、海平面、土壤濕度以及植被覆蓋等地表環境進行監測[4～7]。針對基于SNR數據進行GNSS反射測量這一技術，國內學者也積極研究，并且在潮位變化監測、大壩水位監測、不同GNSS信號的反演精度以及反演算法等方面取得了諸多成果[8～11]。

為了讓大家了解學習GNSS反射測量技術，本文從理論原理和實現方法兩個方面進行了詳細介紹，并使用深圳茜坑水庫的監測數據對該項技術進行了實驗驗證。

2 理論原理

2.1 菲涅爾反射區

介紹GNSS反射測量原理之前，先介紹一下菲涅爾反射區。因為了解了菲涅爾反射區，可以讓我們在進行水位反演之前先對數據進行合理的預處理。

根據惠更斯-菲涅爾原理，菲涅爾反射區是在收發天線之間接收信號反射的區域，其中接收點信號最強的區域稱為第一菲涅爾反射區，如圖1所示。對于波長為λ的電波，由電波的直線路徑與折線路徑的行程差為λ/2的折點形成的區域即第一菲涅爾反射區。

圖1 菲涅爾反射區示意圖

為了更好地將菲涅爾反射區應用于GNSS反射測量，國外學者Larson對第一菲涅爾反射區的表達式進行了修正[12]，可以表示為：

(1)

Cy=0

(2)

(3)

(4)

式(1)～式(4)中：(Cx，Cy)為反射區中心點C的坐標，a為反射區橢圓的長半軸，b為反射區橢圓的短半軸，e為衛星高度角。

2.2 衛星反射點軌跡計算

結合菲涅爾反射區我們可以確定用于GNSS反射測量的衛星的高度角范圍。衛星反射點軌跡則能更加直觀地反映可以用于GNSS反射測量的衛星數據分布。

通過衛星、接收機以及反射點三者之間的簡單幾何關系，我們可以將衛星反射點軌跡表示為：

x=hcote·sinθ

(5)

y=hcote·cosθ

(6)

式(5)和式(6)中：x為反射點軌跡在X軸的分量，y為反射點軌跡在Y軸的分量，θ為衛星方位角。

2.3 GNSS反射測量原理

GNSS反射測量技術利用了來自衛星的直射信號和來自地表的反射信號之間的延遲，根據衛星、接收機天線和信號反射點之間的幾何關系來反演地表參數。利用單天線GNSS信噪比數據進行反射測量的幾何關系示意圖如圖2所示。

（三）區域認知。地理的各個知識點是具有區域性的，每個地區當中的自然環境、條件各不相同，相對應的，其人文事物也就有所不同，因此，學生應該具有區域認知能力，在學習到新的知識的時候，應該將知識與相應的區域相聯系起來，加強對于區域內地理相關知識的理解。

圖2 單天線GNSS反射測量幾何關系示意圖

圖2中，e為衛星高度角，即信號入射角，h為接收機天線相位中心到反射面的垂直距離，δ為GNSS衛星直射信號和反射信號之間的程差。

由于程差的存在，GNSS衛星直射信號和反射信號之間會存在一定的相位延遲ψ，可以表示為：

(7)

式(7)中，λ為載波波長。從式(7)中可以看出，直射信號和反射信號之間的相位延遲與衛星高度角相關，因為衛星高度角是隨時間變化的，所以直射信號和反射信號之間的相位延遲也是隨時間變化而變化的。

對于傳統測量型GNSS接收機，其接收到的復合信號可以表示為：

(8)

式(8)中，AC為復合信號的振幅，Ad為直射信號的振幅，Ar為反射信號的振幅。

由于衛星高度角較低時，接收機天線接收到的反射信號更多，但此時在多路徑效應和天線增益模式的共同影響下，直射信號的變化趨勢決定了復合信號整體的變化趨勢，因此為了能夠有效地提取到可以用于反演地表參數的反射信號，需要去除接收到的信噪比數據的趨勢項，本文采用的方法為通過低階多項式去除趨勢項。結合式(7)，去除趨勢項之后的信噪比殘差序列可以表示為：

(9)

(10)

由于式(10)中的頻率f包含了接收機天線相位中心至反射面的垂直距離h，因此對去除趨勢項之后的信噪比殘差序列進行頻譜分析，得到振幅最大值所對應的頻率值即可計算出接收機天線相位中心至反射面的垂直距離。因為式(10)中t為與衛星高度角相關的正弦函數值，為非等間隔采樣，因此需要采樣Lomb-Scargle頻譜分析方法進行頻譜分析處理。

2.4 Lomb-Scargle頻譜分析

傅立葉變換作為常用的時頻變換分析方法，要求數據連續且分布均勻，而Lomb-Scargle頻譜分析不要求數據連續均勻分布，并且可以有效地從非等間隔采樣的數據中提取較弱的周期信號，還可以減弱由于序列的不連續性導致的虛假結果。因此Lomb-Scargle頻譜分析方法很適合用于提取去除趨勢項之后的信噪比殘差序列的最大振幅值。

對于Lomb-Scargle頻譜分析方法，離散非均勻變化的時域數據X(tj)，j=1，2，3，…，N，其功率譜可定義為關于頻率f的函數如下[13]：

(11)

式(11)中，Px(f)為頻率為f的周期信號的功率，X(tj)為離散的實驗數據，tj為X(tj)所對應的時間，N為實驗數據的統計量，τ為時間平移的不變量。

3 實現方法

基于上述理論原理，GNSS信噪比反射測量數據處理平臺的主要流程設計如圖3所示。詳細流程設計如圖4所示。Matlab作為一門面向科學與工程計算的高級語言，編程效率高，使用方便，并且具有完備的圖形處理能力，因此本數據處理平臺采用Matlab語言實現。

圖3 主要流程

圖4 詳細流程

GNSS信噪比反射測量數據處理平臺主要包含數據預處理、反演分析以及精度評估三部分。

在數據預處理部分，重點是對GNSS觀測數據中的信噪比數據以及對應的衛星方位角和高度角進行提取。其中，信噪比數據可以直接從衛星的觀測數據中進行提取，衛星的高度角和方位角可以利用衛星的精密星歷和測站的位置計算得到。計算菲涅爾反射區時對于不同的測站我們需要設置合適的天線到反射面的高度才可以準確地計算出相應的菲涅爾反射區范圍。

在反演分析部分，首先是根據預處理部分的結果人工判斷選擇合適的衛星高度角和方位角范圍內的數據，其次，根據不同的情況，我們需要設置合適的多項式階數去去除信噪比數中的趨勢項，本文實驗采用的是3階多項式，具體設置可以根據反演結果的精度進行變動。為了避免使用過短的弧段進行反演分析，可以設置當滿足衛星高度角條件和方位角條件且連續觀測的信噪比數據大于30個時才用于反演分析，反之直接剔除。反演分析時用到的Lomb-Scargle頻譜分析方法，本文直接調用了MATLAB自帶的lomb函數進行分析。當計算得到反演結果之后，我們還需要通過設置最小振幅、最大振幅與平均振幅的比值等排除掉反演結果不理想的數據，具體參數值不同測站可以根據不同設置值情況的反演結果進行設定，無固定值設置。

在精度評估部分，需要對反演結果值的精度以及整體的趨勢進行評估，即單個反演結果與實際水位觀測結果的差值以及連續時間的反演結果與整體水位觀測結果的吻合性。這一部分的評估我們主要采用了均方根誤差(RMS)以及相關性分析。

4 實驗與結果分析

4.1 實驗數據

本文選取了深圳茜坑水庫變形監測系統XK03號監測站采集的2018年第180天～240天共61天的GPS數據。XK03號測站的位置示意圖如圖5所示。用于驗證大壩水位反演精度的數據為每天早上8點由人工進行采集的水位數據。

圖5 XK03測站位置示意圖

4.2 實驗結果分析

為了選取可以用于反演分析的最佳衛星高度角區間，首先模擬當天線相位中心到水面的垂直距離為 6 m時，XK03測站周邊的菲涅爾反射區情況如圖6所示。

圖6 XK03測站菲涅爾反射區

從圖6可以看出，當衛星高度角在5°～15°范圍內變化時，接收機天線可以接收到周圍 120 m范圍內的反射信號。結合圖5，本次實驗數據采用衛星方位角在220°～350°之間，衛星高度角在5°～15°之間的GPS觀測數據。

以衛星高度角和衛星方位角范圍確定限制區域，可以繪制衛星反射點軌跡圖如圖7所示，圖7反映了可以用于反演分析的GPS衛星編號。

圖7 GPS衛星反射點軌跡圖

以G21衛星為例，其觀測值中信噪比數據隨時間的變化情況如圖8所示，可以看出，信噪比數據整體趨勢呈拋物線形式，因此可以通過二次多項式擬合的方法去除趨勢項。去除趨勢項之后的信噪比殘差數據如圖9所示?？梢钥闯鲂旁氡葰埐钚蛄信c衛星高度角之間存在正弦函數關系。通過Lomb-Scargle頻譜分析方法我們可以得到圖10所示的結果。圖10中頻譜振幅最大值所對應的高度即為我們所要反演的天線相位中心至水面的垂直距離。

圖8 G21原始信噪比數據隨時間變化圖

圖9 G21信噪比殘差序列隨衛星高度角變化圖

圖10 Lomb-Scargle頻譜分析結果

將反演水位結果與人工實測水位結果進行對比可以得到圖11和圖12。圖13給出了反演水位結果與人工實測水位結果的相關系數圖。

圖11 反演水位變化量與實測水位變量

圖12 反演水位變化量與實測水位變量差值

圖13 相關性分析

由于GNSS反演得到的水位數據很多，并且大壩水位一天之中變化相對較為緩慢，為了便于進行統計分析，我們將一天之中的反演水位高度求平均之后作為當天的反演水位高度。結果表明，反演水位變化量和實測水位變量的相關系數達0.99，兩者差值的RMS值為 0.11 m。

5 結 語

GNSS反射測量技術作為全球衛星導航系統的一個全新研究領域，本文對該項技術的理論原理及實現方法進行了詳細介紹，并以深圳茜坑水庫的變形監測數據為基礎進行了實驗驗證。實驗結果表明，本文的方法可以很好地實現GNSS反射測量，成果豐富，驗證了GNSS在大壩水位反演方面的可行性，并且反演的大壩水位變化量與人工實測水位變化量的相關系數達到0.99，兩者差值的RMS為0.11 m。