GNSS信噪比數(shù)據(jù)反演大壩水位的方法與應(yīng)用

何騫，竇邵華

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

1 引 言

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

為了讓大家了解學(xué)習(xí)GNSS反射測量技術(shù)，本文從理論原理和實現(xiàn)方法兩個方面進(jìn)行了詳細(xì)介紹，并使用深圳茜坑水庫的監(jiān)測數(shù)據(jù)對該項技術(shù)進(jìn)行了實驗驗證。

2 理論原理

2.1 菲涅爾反射區(qū)

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

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

圖1 菲涅爾反射區(qū)示意圖

為了更好地將菲涅爾反射區(qū)應(yīng)用于GNSS反射測量，國外學(xué)者Larson對第一菲涅爾反射區(qū)的表達(dá)式進(jìn)行了修正[12]，可以表示為：

(1)

Cy=0

(2)

(3)

(4)

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

2.2 衛(wèi)星反射點軌跡計算

結(jié)合菲涅爾反射區(qū)我們可以確定用于GNSS反射測量的衛(wèi)星的高度角范圍。衛(wèi)星反射點軌跡則能更加直觀地反映可以用于GNSS反射測量的衛(wèi)星數(shù)據(jù)分布。

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

x=hcote·sinθ

(5)

y=hcote·cosθ

(6)

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

2.3 GNSS反射測量原理

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

（三）區(qū)域認(rèn)知。地理的各個知識點是具有區(qū)域性的，每個地區(qū)當(dāng)中的自然環(huán)境、條件各不相同，相對應(yīng)的，其人文事物也就有所不同，因此，學(xué)生應(yīng)該具有區(qū)域認(rèn)知能力，在學(xué)習(xí)到新的知識的時候，應(yīng)該將知識與相應(yīng)的區(qū)域相聯(lián)系起來，加強對于區(qū)域內(nèi)地理相關(guān)知識的理解。

圖2 單天線GNSS反射測量幾何關(guān)系示意圖

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

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

(7)

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

對于傳統(tǒng)測量型GNSS接收機，其接收到的復(fù)合信號可以表示為：

(8)

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

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

(9)

(10)

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

2.4 Lomb-Scargle頻譜分析

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

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

(11)

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

3 實現(xiàn)方法

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

圖3 主要流程

圖4 詳細(xì)流程

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

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

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

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

4 實驗與結(jié)果分析

4.1 實驗數(shù)據(jù)

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

圖5 XK03測站位置示意圖

4.2 實驗結(jié)果分析

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

圖6 XK03測站菲涅爾反射區(qū)

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

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

圖7 GPS衛(wèi)星反射點軌跡圖

以G21衛(wèi)星為例，其觀測值中信噪比數(shù)據(jù)隨時間的變化情況如圖8所示，可以看出，信噪比數(shù)據(jù)整體趨勢呈拋物線形式，因此可以通過二次多項式擬合的方法去除趨勢項。去除趨勢項之后的信噪比殘差數(shù)據(jù)如圖9所示。可以看出信噪比殘差序列與衛(wèi)星高度角之間存在正弦函數(shù)關(guān)系。通過Lomb-Scargle頻譜分析方法我們可以得到圖10所示的結(jié)果。圖10中頻譜振幅最大值所對應(yīng)的高度即為我們所要反演的天線相位中心至水面的垂直距離。

圖8 G21原始信噪比數(shù)據(jù)隨時間變化圖

圖9 G21信噪比殘差序列隨衛(wèi)星高度角變化圖

圖10 Lomb-Scargle頻譜分析結(jié)果

將反演水位結(jié)果與人工實測水位結(jié)果進(jìn)行對比可以得到圖11和圖12。圖13給出了反演水位結(jié)果與人工實測水位結(jié)果的相關(guān)系數(shù)圖。

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

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

圖13 相關(guān)性分析

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

5 結(jié) 語

GNSS反射測量技術(shù)作為全球衛(wèi)星導(dǎo)航系統(tǒng)的一個全新研究領(lǐng)域，本文對該項技術(shù)的理論原理及實現(xiàn)方法進(jìn)行了詳細(xì)介紹，并以深圳茜坑水庫的變形監(jiān)測數(shù)據(jù)為基礎(chǔ)進(jìn)行了實驗驗證。實驗結(jié)果表明，本文的方法可以很好地實現(xiàn)GNSS反射測量，成果豐富，驗證了GNSS在大壩水位反演方面的可行性，并且反演的大壩水位變化量與人工實測水位變化量的相關(guān)系數(shù)達(dá)到0.99，兩者差值的RMS為0.11 m。