GNSS-R潮位監(jiān)測抗差估計

王澤明，李浩軍，孫亞峰

1.同濟(jì)大學(xué)測繪與地理信息學(xué)院，上海 200092；2.上海市城市建設(shè)設(shè)計研究總院(集團(tuán))有限公司，上海 200125

臨海國家的生產(chǎn)生活受海平面升降影響巨大，因此對海平面進(jìn)行實時監(jiān)測向來是眾多研究領(lǐng)域的重點[1]。監(jiān)測海平面高度變化的傳統(tǒng)方法主要包括驗潮站實測與激光測高技術(shù)。作為觀測精度最高的驗潮站數(shù)據(jù)，仍然存在著地理分布不均和數(shù)據(jù)不公開等問題。具有良好測高能力的激光測高技術(shù)在近岸地區(qū)的水位監(jiān)測方面存在著精度低和衛(wèi)星重訪問周期長等問題[2]。全球?qū)Ш叫l(wèi)星系統(tǒng)反射技術(shù)(global navigation satellite system reflectometry，GNSS-R)是一種基于GNSS反射信號的全新遙感技術(shù)，它通過分析反射信號的波形、振幅、相位、頻率等信息，提取出隱含在其中的反射面的物理特性[3]。該技術(shù)的優(yōu)勢在于反射點具有密集的空間分布，可以覆蓋測站周圍整個區(qū)域，且接收機(jī)天線距反射面高度越高，覆蓋范圍越大，因而不再局限于諸如浮標(biāo)、驗潮站等單一測量點的測量模式[4]。與傳統(tǒng)遙感測量技術(shù)相比，GNSS-R技術(shù)不需要專門的雷達(dá)發(fā)射機(jī)，并且具有全天候、覆蓋范圍廣、時間分辨率高、成本低等優(yōu)點[5]，目前已應(yīng)用于地表土壤濕度、植被變化、海面風(fēng)速、雪深和潮位等的反演[6-11]。

目前，GNSS-R技術(shù)在水面高度反演領(lǐng)域已經(jīng)取得眾多成果。文獻(xiàn)[12]在Crater湖上開展的基于1 Hz采樣率的試驗獲得了2 cm精度的水面高度。文獻(xiàn)[13]設(shè)計了一款基于GNSS反射信號的潮汐計量儀，在Onsala空間天文臺展開的試驗獲得了4 cm精度的水面高度。文獻(xiàn)[14]基于機(jī)載GNSS-R技術(shù)，對已有反演模型進(jìn)行優(yōu)化，獲得了5 cm高程精度和5 km空間分辨率。文獻(xiàn)[15]對信噪比序列和頻譜分析閾值進(jìn)行約束，基于最小二乘方法，結(jié)合多個GNSS星座，對時間跨度為3個月的數(shù)據(jù)進(jìn)行潮位反演，試驗所得相關(guān)系數(shù)高達(dá)0.97，同時指出GNSS-R技術(shù)在反演海面高度方面有很強的實用價值。文獻(xiàn)[16]提出一種圓周回歸估計算法和一種接收機(jī)體系結(jié)構(gòu)，試驗所得海面高度反演精度為厘米級。

針對驗潮站參考面至天線相位中心間垂直距離(下文簡稱“測站高”)的求解，以往絕大多數(shù)研究采用對應(yīng)歷元水面高度與潮位取平均的方法，本文提出了基于抗差估計的優(yōu)化方法，采用4個IGS GNSS連續(xù)運行跟蹤站HNLC、SC02、TDAM、TPW2的觀測數(shù)據(jù)進(jìn)行數(shù)據(jù)處理。結(jié)果表明，本文方法對水面高度粗差反演值的剔除、測站高計算精度的提升，以及對最終反演潮位結(jié)果都具有重要的意義。同時，驗證了基于先驗潮位數(shù)據(jù)求得的測站高用于后續(xù)無驗潮站或驗潮站數(shù)據(jù)中斷情況下的GNSS-R潮位反演的可行性。

1 精準(zhǔn)測站高計算原理

1.1 SNR反演水面高度原理

信噪比(signal noise ratio，SNR)是GNSS接收機(jī)觀測數(shù)據(jù)中信號與噪聲的比例，由衛(wèi)星直射、反射信號疊加干涉形成?？捎墒?1)計算得出[17]

(1)

式中，Ad表示直射信號的振幅；Am表示由多路徑效應(yīng)引起的局部周期性變化[18]；Ψ表示直反射信號間的相位差。因反射信號對GNSS定位精度造成嚴(yán)重影響，通常被接收機(jī)天線的扼流圈和天線增益所抑制[19]，且反射信號在被反射時會發(fā)生能量衰減，故直射信號與反射信號的振幅滿足Ad?Am的關(guān)系[18]，因此式(1)可近似為

(2)

式(2)表明，直射信號構(gòu)成了SNR的主體，表征著整體變化趨勢；而反射信號對SNR的貢獻(xiàn)體現(xiàn)在小幅度的周期性震蕩擾動，且該擾動僅在低衛(wèi)星高度角時較為顯著[20]。

相位差的產(chǎn)生源于直反射信號之間存在路徑延遲，二者間的關(guān)系可表示為[21-22]

(3)

式中，λ表示衛(wèi)星信號的載波波長；δ表示直反射信號的路徑延遲。如圖1所示，路徑延遲δ可表示為

圖1 GNSS-R反演模型

δ=2hsinθ

(4)

式中，h表示自反射面起算的接收機(jī)天線相位中心的高度；θ為衛(wèi)星高度角。

根據(jù)式(3)、式(4)，對相位差Ψ求時間的導(dǎo)數(shù)，可獲得多路徑效應(yīng)振蕩頻率為[15]

(5)

(6)

(7)

由于LSP分析提取的是水面高度，故在與驗潮站實測數(shù)據(jù)對比分析前，需要進(jìn)行如下轉(zhuǎn)換[19,26]

htide(t)=[h(t0)+Htide(t0)]-h(t)

(8)

式中，htide(t)為t時刻反演所得潮位高度；h(t0)為某一參考時刻t0反演所得水面高度；Htide(t0)為t0時刻驗潮站實測潮位；h(t)為t時刻反演所得水面高度?？梢钥闯?，h(t0)+Htide(t0)表示的正是基于驗潮站參考基準(zhǔn)起算的測站高。對于該高度的求取，多數(shù)研究將試驗期間所有反演結(jié)果取均值代入h(t0)，并將與反演結(jié)果對應(yīng)時刻的實際潮位高度取均值代入Htide(t0)，計算得測站高。由于GNSS-R技術(shù)反演水面高度利用的是水面反射的衛(wèi)星信號，而測站周圍的建筑、植被、行船等地物同樣會將部分衛(wèi)星信號反射傳入接收機(jī)；此外，大氣折射效應(yīng)會使衛(wèi)星信號發(fā)生彎曲，從而導(dǎo)致計算高度角與實際高度角不同，引起潮位反演結(jié)果的偏差[27]；海面粗糙度大時，高頻海浪會導(dǎo)致絕大部分反射信號來自波峰，低頻海浪會導(dǎo)致絕大部分反射信號來自波谷，導(dǎo)致獲取的水面高度并非是基于平均海平面起算的水面高度[15]；同時風(fēng)和雨等因素也會嚴(yán)重影響反演結(jié)果。因此將試驗期間所有反演所得水面高度取均值代入求得的h(t0)，必定是包含著受多個誤差源影響的水面高度均值，以此求出的測站高并非其最優(yōu)解。

1.2 抗差估計求解精準(zhǔn)測站高原理

抗差估計又稱為穩(wěn)健估計，是處理離群值的一種有效途徑，即使在分析模型與實際模型存在偏差時仍能得到穩(wěn)定可靠的結(jié)果，其利用等價權(quán)使傳統(tǒng)的最小二乘估計具備抗干擾能力[28]。權(quán)函數(shù)是抗差估計的核心問題，就抗差目的而言，所選權(quán)函數(shù)應(yīng)具備3個階段：當(dāng)觀測誤差很小時，采用原始權(quán)進(jìn)行最小二乘處理(正常段)；當(dāng)觀測誤差較大但并不顯著時，采用降權(quán)估計處理(可疑段)；當(dāng)觀測值顯著異常時，采用零權(quán)淘汰處理(淘汰段)。正常段可以準(zhǔn)確揭示觀測數(shù)據(jù)的內(nèi)在分布模式，利用極大似然估計即可獲得其最優(yōu)估值；可疑段雖不能準(zhǔn)確揭示數(shù)據(jù)的內(nèi)在分布模式，但若能合理地利用其中內(nèi)含的基本分布特征，對提高估值精度仍大有裨益；淘汰段提供的是有害信息，影響結(jié)果的正常估計，應(yīng)對其予以剔除。本文權(quán)函數(shù)選擇IGG Ⅲ方案，該方案較Huber、Tukey、Andrews權(quán)函數(shù)具備完整的3個階段，較Hampel權(quán)函數(shù)處理方便，實用性更強[29]。

對于受多種誤差源影響的GNSS-R水面高度反演值，采用抗差估計方法，將可疑段的反演值依其標(biāo)準(zhǔn)化殘差分別降權(quán)，對顯著異常反演值不承認(rèn)其信息，取零權(quán)淘汰，以此削弱反演粗差，確保最終求解的精準(zhǔn)測站高不受其影響。

對于任意時間段內(nèi)反演所得水面高度，存在如下關(guān)系式

(9)

(10)

l=hGNSS+htide-hstation_0

(11)

式中,hstation_0為測站高初始值，可由試驗期間所有反演所得水面高度與實際潮位求和后取平均獲得。

由抗差估計理論可得如下極值條件

(12)

式中,Pi為第i個觀測量權(quán)值；ρ為凸函數(shù)；νi為第i個殘差。對式(12)求極值，由等價權(quán)理論得

(13)

(14)

(15)

式中，σνi為νi的標(biāo)準(zhǔn)差。

本文將GNSS-R技術(shù)與抗差估計方法相結(jié)合，直接目的是求解出測站高的最優(yōu)估值，進(jìn)而提升GNSS-R技術(shù)的反演精度。采用抗差估計的方法削弱粗差對反演結(jié)果的影響，本質(zhì)上并非從消除誤差源的角度入手，而是從反演結(jié)果入手，篩選可靠的水面高度信息代入測站高度最優(yōu)值的求解，以期達(dá)到與消除誤差源的技術(shù)路線相近的優(yōu)化效果，而操作步驟卻更加簡單易行，圖2給出了本文的技術(shù)流程。

圖2 技術(shù)路線

2 數(shù)據(jù)處理與分析

本文使用的試驗數(shù)據(jù)來自2018年1至2月IGS GNSS連續(xù)運行跟蹤站HNLC、SC02、TDAM、TPW2。各測站的經(jīng)緯度、接收機(jī)類型、天線類型、整流罩類型、采樣間隔見表1，各測站的周圍環(huán)境、地理位置如圖3所示。為驗證GNSS-R潮位反演的精度，分別選取緊鄰上述4個測站的火奴魯魯驗潮站、星期五港驗潮站、加爾維斯頓驗潮站、阿斯托里亞驗潮站的潮位數(shù)據(jù)，4個驗潮站均以平均低低潮面為潮位基準(zhǔn)面，提供每6 min更新一次的潮位數(shù)據(jù)。

圖3 各測站信息

表1 HNLC、SC02、TDAM、TPW2測站信息

2.1 抗差估計求解精準(zhǔn)測站高

以2018年DOY 01 PRN27號衛(wèi)星L2波段的部分?jǐn)?shù)據(jù)為例進(jìn)行反演分析。剔除衛(wèi)星時間序列中的直射信號，提取衛(wèi)星高度角5°～15°、方位角205°～330°范圍內(nèi)的SNR殘差序列并進(jìn)行線性化處理，LSP分析結(jié)果中頻譜振幅峰值所對應(yīng)的頻率即水體反射頻率。圖4(a)、(b)分別給出了PRN27號衛(wèi)星的SNR殘差序列和LSP頻譜分析結(jié)果。

圖4 PRN27號衛(wèi)星提取水體反射頻率

本節(jié)采用傳統(tǒng)法(將試驗期間所有反演得到的水面高度的均值與對應(yīng)時刻實際潮位高度的均值相加求得測站高)和抗差估計法對HNLC、SC02、TDAM、TPW2這4個測站2018年1月的數(shù)據(jù)進(jìn)行反演。高度角、方位角、反演波段、抗差因子的選取見表2。圖5給出了4個測站采用傳統(tǒng)法和抗差估計法反演所得潮位與時空同步的驗潮站實測潮位的對比結(jié)果。圖6給出了4個測站采用上述兩種方法的相關(guān)性分析對比結(jié)果。表3給出了4個測站傳統(tǒng)法和抗差估計法反演結(jié)果的比較。通過對比可以發(fā)現(xiàn)，對于4個測站，采用兩種方法求得的測站高分別相差1.02、1.31、16.2、1.22 cm；RMSE方面，4個測站采用抗差估計方法后分別降低了26.7%、34.4%、84%、31.6%；相關(guān)系數(shù)方面，4個測站采用抗差估計方法后分別提高了5.6%、0.69%、41.5%、1.04%，反演潮位更加貼合實測潮位，潮位殘差更小。這是因為抗差估計方法有效削弱了水面高度反演粗差的權(quán)重，阻止其參與精準(zhǔn)測站高度的計算，以此獲得測站高度最優(yōu)解?？梢钥闯鯰DAM測站使用兩種方法反演結(jié)果差異明顯，由傳統(tǒng)法反演潮位與實測潮位序列圖可知，除反演潮位整體較實測潮位略微偏低以外，絕大多數(shù)反演結(jié)果與實測結(jié)果吻合較好，只有極少數(shù)時刻反演值嚴(yán)重偏離實測值。由于TDAM站位于加爾維斯頓港口，常有船只來往，因此推測上述極少數(shù)偏離值是由船體噪聲混入水體反射信號造成的。該部分偏離值為數(shù)不多但誤差極大，造成了16.2 cm的測站高差別。而抗差估計方法可以準(zhǔn)確識別并剔除該部分反演粗差，利用剩余正確反演值求解測站高度無偏估值，從而達(dá)到去噪目的。如圖5(e)、(f)所示，利用經(jīng)抗差估計去噪后求得的測站高進(jìn)行潮位反演，反演潮位整體不再比實測潮位偏低，契合度大大提升，可見經(jīng)抗差估計去噪后求得的測站高度更接近TDAM測站高度的真值。

表3 各測站1月數(shù)據(jù)反演結(jié)果

圖5 各測站傳統(tǒng)法與抗差估計法反演潮位時間序列與殘差

圖6 各測站傳統(tǒng)法與抗差估計法相關(guān)性分析

表2 各測站反演參數(shù)選擇

2.2 精準(zhǔn)測站高用于無驗潮站場合的潮位反演

基于抗差估計的GNSS-R潮位反演方法的實際應(yīng)用價值在于能夠求出測站高度的無偏估值，在忽略測站地基極其微小沉降的情況下，可將求得的測站高度視為定值，供此后若干年間繼續(xù)使用，不再需要驗潮站參與潮位高度的反演。因此本節(jié)將上一節(jié)采用傳統(tǒng)法和抗差估計法求解各測站2018年1月數(shù)據(jù)所得的測站高度作為已知值，分別直接代入2月數(shù)據(jù)中進(jìn)行潮位反演(下文分別稱為傳統(tǒng)-直接法和抗差估計-直接法)。反演潮位由代入的測站高度與反演所得水面高度直接獲得，無須驗潮站的參與，驗潮站僅用于提供真實潮位以檢驗兩種直接法反演潮位的精度。圖7、圖8分別為抗差估計-直接法反演結(jié)果與實測潮位的比較，以及相關(guān)性分析結(jié)果。由于上述兩種直接法的差別僅在于所使用測站高不同，因此二者反演結(jié)果除整體相差一個定值(測站高之差)外，圖像走勢等要素完全相同，因此不再給出傳統(tǒng)-直接法反演結(jié)果與實測潮位的比較及相關(guān)性分析結(jié)果)。表4給出了4個測站兩種直接法反演的潮位與時空同步的驗潮站實測潮位的比較結(jié)果。通過對比發(fā)現(xiàn)，采用抗差估計法求解1月數(shù)據(jù)所得的測站高度比采用傳統(tǒng)法的更適合2月數(shù)據(jù)的反演，4個測站的抗差估計-直接法的RMSE均小于傳統(tǒng)-直接法的RMSE。同樣的，對于TDAM測站，抗差估計-直接法較傳統(tǒng)-直接法有著明顯的優(yōu)勢，相關(guān)系數(shù)大幅提升，RMSE大幅減小，這也間接說明對于TDAM測站1月數(shù)據(jù)采用抗差估計求解的測站高度較傳統(tǒng)法更接近其真值。

圖7 各測站抗差估計-直接法潮位時間序列與殘差

圖8 各測站抗差估計-直接法相關(guān)性分析

表4 各測站2月數(shù)據(jù)反演結(jié)果

由以上試驗可知，采用抗差估計求解IGS GNSS連續(xù)運行跟蹤站數(shù)據(jù)所得的精準(zhǔn)測站高度，可以應(yīng)用于之后無驗潮站或驗潮站數(shù)據(jù)中斷情況下潮位的獲取。4個算例中，TDAM測站明顯體現(xiàn)出抗差估計方法的優(yōu)越性，這表明測站觀測數(shù)據(jù)受周圍地物噪聲影響越大，抗差估計的去噪優(yōu)勢越明顯，即便是用于之后無驗潮站參與情況下的潮位反演，抗差估計-直接法依舊顯著優(yōu)于傳統(tǒng)-直接法。此外，表4中除TDAM測站外的其余測站兩種直接法的反演精度較為接近，原因在于當(dāng)無驗潮站參與時，不易判斷反演所得水面高度是否包含粗差，因此均未對其進(jìn)行剔除；但同時，由于抗差估計-直接法所使用的測站高更接近其真值，因此在反演2月潮位數(shù)據(jù)時，其RMSE總是小于傳統(tǒng)-直接法。另一方面，2月潮位的反演精度普遍低于1月潮位的反演精度，其根本原因在于，抗差估計-直接法使用的測站高度是求解該測站1月數(shù)據(jù)所得的無偏估值，由于各測站1月的樣本數(shù)量有限，無法保證1月樣本所得的測站高無偏估值等于2月樣本的無偏估值，即不同的時間序列求解出的測站高無偏估值往往存在差異，因此要充分發(fā)揮抗差估計法反演潮位優(yōu)勢的前提是，必須具備充足的樣本量，即求解測站高度無偏估值時選用GNSS觀測數(shù)據(jù)的時間序列越長，所得結(jié)果越接近測站高度真值。相應(yīng)地，對于某一驗潮站，其擁有的歷史潮位數(shù)據(jù)越多，越能發(fā)揮抗差估計的優(yōu)勢以求解其附近IGS GNSS連續(xù)運行跟蹤站的測站高度的最優(yōu)解，也就越能發(fā)揮岸基IGS GNSS測站監(jiān)測潮位的能力。

2.3 抗差估計抵抗人工噪聲

為更好地體現(xiàn)抗差估計計算精準(zhǔn)測站高、提高GNSS-R反演精度的優(yōu)勢，本文還設(shè)計了一組用傳統(tǒng)法和抗差估計法兩種方法處理觀測數(shù)據(jù)中包含人工噪聲的試驗。仍以SC02測站2018年1月的觀測數(shù)據(jù)為例，向信噪比文件中加入人工噪聲，添加的方案是：每隔若干個歷元，在原SNR值的基礎(chǔ)上添加某一隨機(jī)數(shù)。表5給出了每隔20個歷元，在原SNR值中添加不同大小的噪聲后兩種方法的反演結(jié)果對比。由前文分析可知，可以預(yù)測不同反演模式下的精度應(yīng)滿足：無人工噪聲的抗差估計法>無人工噪聲的傳統(tǒng)法>有人工噪聲的抗差估計法>有人工噪聲的傳統(tǒng)法，試驗結(jié)果證實了這一猜想。

表5 傳統(tǒng)法和抗差估計法處理恒定歷元間隔的人工噪聲

3 結(jié) 語

結(jié)合試驗對比結(jié)果，本文有如下主要結(jié)論：

(1)針對HNLC、SC02、TDAM、TPW2測站2018年1、2月的試驗數(shù)據(jù)，基于岸基IGS GNSS站的GNSS-R技術(shù)與驗潮站實測數(shù)據(jù)在整體趨勢上吻合較好，精度為dm級。

(2)抗差估計法能夠更準(zhǔn)確地求解出基于驗潮站潮位起算面的測站高度，能夠更好地減小反演的均方根誤差，提高相關(guān)系數(shù)，對于受周圍地物噪聲影響嚴(yán)重的測站優(yōu)化效果尤為顯著。

(3)以HNLC、SC02、TDAM、TPW2測站為例進(jìn)行的試驗研究，證明了基于先驗潮位數(shù)據(jù)求得的測站高度，可以繼續(xù)用于后期無驗潮站或驗潮站數(shù)據(jù)中斷情況下潮位的反演。

在導(dǎo)航定位中被視為嚴(yán)重影響定位精度并被人們設(shè)法抑制的反射信號，在近地水環(huán)境等領(lǐng)域的遙感監(jiān)測中卻發(fā)揮著重要作用。岸基CORS網(wǎng)絡(luò)為開展GNSS-R技術(shù)研究提供了很好的數(shù)據(jù)支撐，可有效填補多地監(jiān)測數(shù)據(jù)的空白。本文利用4個IGS GNSS連續(xù)運行跟蹤站反演海面潮位數(shù)，采用抗差估計的方法，使GNSS更好地應(yīng)用于潮位監(jiān)測，并證明了岸基IGS GNSS站可連續(xù)對潮位變化進(jìn)行監(jiān)測。