奇異譜分析在GNSS-MR海平面高度反演中的應(yīng)用

劉睿國，劉立龍，吳 晗，薛張芳，吳昊艦，張 志

(1.桂林理工大學(xué) 測繪地理信息學(xué)院，廣西 桂林 541006；2.廣西空間信息與測繪重點實驗室，廣西 桂林 541006)

0 引言

隨著全球氣候變暖，南極和北極冰川融化速度加快，海平面逐年上升。為了更好地保護(hù)人類活動，特別是沿海地區(qū)的環(huán)境管理，監(jiān)測海平面的變化至關(guān)重要。GNSS多路徑反射測量(GNSS Multipath Reflectometry,GNSS-MR)作為海平面高度反演的一種新方法，尤其對于陸地邊緣海平面高度的變化監(jiān)測更為有利，該方法利用測量型接收機(jī)來接收GNSS直射信號和反射信號。直射信號和反射信號的相位不同引起多路徑效應(yīng)，相位的變化隨著衛(wèi)星的位置變化而變化，從而形成干涉圖樣，干涉圖樣以振幅的形式記錄在信噪比(Signal to Noise Ratio,SNR)中，用于進(jìn)行參數(shù)分析，具有全天候、實時自動化和成本低等優(yōu)點，因此在海平面高度反演[1-3]、有效波高反演[4]、積雪深度反演[5]、土壤濕度反演[6-7]以及形變監(jiān)測[8]等方面得到廣泛應(yīng)用。Larson等[9]利用GNSS接收機(jī)獲取的GPS信號數(shù)據(jù)進(jìn)行潮位反演，實驗表明GPS-MR的反演結(jié)果與驗潮站數(shù)據(jù)具有較好的一致性。吳繼忠等[10]利用GPS觀測數(shù)據(jù)對靜止的湖面高度進(jìn)行了反演，反演結(jié)果與實測結(jié)果比較吻合，GPS-MR反演的靜止湖面高度標(biāo)準(zhǔn)偏差達(dá)到3 cm。張雙成等[11]利用GPS觀測數(shù)據(jù)進(jìn)行海平面高度反演，結(jié)果表明GPS-MR反演結(jié)果與驗潮站實測值相差約10 cm。劉立龍等[1]對反演所需數(shù)據(jù)的選擇方法進(jìn)行了研究，依據(jù)觀測站觀測環(huán)境、有效反射區(qū)域反算的衛(wèi)星高度角可以對反演所需的數(shù)據(jù)進(jìn)行選擇，結(jié)果表明海平面高度反演結(jié)果與驗潮儀實測數(shù)據(jù)具有較強(qiáng)的相關(guān)性。陳發(fā)德等[3]利用GPS和BDS的SNR數(shù)據(jù)對海平面高度進(jìn)行了反演，并與驗潮站實測數(shù)據(jù)進(jìn)行對比分析，結(jié)果表明GPS，BDS聯(lián)合反演結(jié)果與驗潮站實測高度的均方根誤差(Root Mean Square Error,RMSE)達(dá)到0.286 m。但上述對于海平面高度反演的研究均采用二次多項式擬合的方法獲取SNR殘差序列，無法擬合出更為準(zhǔn)確的SNR信號趨勢。為了獲取更高質(zhì)量的反射信號，蘇曉容等[12]利用小波分析的方法研究GNSS-R技術(shù)監(jiān)測潮位高度變化，并對復(fù)雜環(huán)境下采集的信號進(jìn)行小波去噪，有效地提高了反演精度。王瑞芳[13]利用經(jīng)驗?zāi)B(tài)分解的方法對SNR殘差序列進(jìn)行分解，剔除其中的噪聲信號，提取出較為純凈的海面反射信號，在一定程度上提高了反演精度。但是以上方法在處理SNR數(shù)據(jù)時存在一定的局限性，如經(jīng)驗?zāi)B(tài)分解的模態(tài)混疊問題，導(dǎo)致獲取的SNR信號質(zhì)量相對不高，包含噪聲和干擾信息，影響GNSS-MR海平面測高的精度。

奇異譜分析(Singular Spectrum Analysis,SSA)是一種對時間序列進(jìn)行主成分分析的方法，SSA分解的空間結(jié)構(gòu)與時間尺度密切相關(guān)，在細(xì)節(jié)分量分析方面具有優(yōu)勢，它不受正弦波假定的約束，且無需先驗信息即可從含有噪聲的數(shù)據(jù)中提取出大量的可靠信息，并對其進(jìn)行重構(gòu)[14-15]。然而，鮮有研究SSA在海平面高度反演的應(yīng)用。本文通過SSA法以布雷斯特港的BRST測站為例進(jìn)行海平面高度反演，通過對比來驗證本文方法的可行性與有效性。

1 GNSS-MR監(jiān)測海平面高度變化原理

衛(wèi)星發(fā)射的直射信號和經(jīng)海平面反射的反射信號同時被GNSS接收機(jī)接收，因為直射信號和反射信號二者的相位不同，所以會引起多路徑效應(yīng)，多路徑效應(yīng)主要記錄在SNR中，用來評估觀測信號的質(zhì)量。利用GNSS-MR技術(shù)反演海平面高度的反演原理如圖1所示，H為GNSS天線相位中心與驗潮站高程基準(zhǔn)(平均低低潮面)間的垂直距離，hGNSS-MR為GNSS-MR反演海平面高度，θ為衛(wèi)星高度角。

圖1 GNSS-MR海平面高度反演原理

直射信號的振幅Ad與反射信號的振幅Ar存在如下關(guān)系：

Ad?Ar,

(1)

記合成信號的振幅為Ac，cosφ為直射信號與反射信號夾角的余弦值，其關(guān)系為：

(2)

因為直射信號的振幅Ad遠(yuǎn)遠(yuǎn)大于反射信號的振幅Ar，所以合成信號總體的變化趨勢由直射信號決定。為了獲取低高度角下因多路徑影響所生成的SNR殘差序列，通常采用二次多項式擬合的方法來消除趨勢項Ad。低高度角下因多路徑效應(yīng)所形成的SNR殘差序列可表示為：

(3)

式中，λ為載波波長；θ為衛(wèi)星高度角；h為垂直反射距離。對SNR殘差序列進(jìn)行Lomb-Scargle頻譜分析得到頻率f,由h=fλ/2可得天線相位中心到海平面的距離h，再把h轉(zhuǎn)換為和驗潮站相同基準(zhǔn)下的海平面高度，從而實現(xiàn)GNSS-MR對海平面高度的反演。

2 SSA的原理與方法

SSA是一種基于主成分分析的時間序列處理方法[16]，主要步驟為：構(gòu)建軌跡矩陣、奇異值分解、分組和重構(gòu)。具體過程如下[17-19]：

(1)構(gòu)建軌跡矩陣

假設(shè)混有噪聲的SNR序列Li的序列長度為N，對其進(jìn)行滯后排列處理，并選擇嵌入維度M，M為整數(shù)并且取值為1≤M≤N/2[20]，得到軌跡矩陣XM×K，其中K=N-M+1，軌跡矩陣X為：

(4)

(2)奇異值分解

X=X1+X2+…+Xd。

(5)

(3)分組

(4)重構(gòu)

(6)

利用特征值變化率選取有用信號和噪聲，將貢獻(xiàn)率最大的特征分量作為其趨勢項；特征值變化率突變拐點之后的有限個特征分量作為數(shù)據(jù)處理所需的殘差序列；剩下的特征分量即為噪聲。特征值變化率為:

r=(λk-λk+2)/2,k=1,2,…,M-2，

(7)

一般認(rèn)為r≤0.001時符合實際情況，此時p=k-1。

3 實驗分析

本文選用BRST觀測站(http:∥www.igs.org)的觀測數(shù)據(jù)進(jìn)行實驗分析，該測站位于法國西海岸的布雷斯特港(48°22′49.78″N，4°29′47.75″W)，該站點屬于法國國家地理研究所。BRST測站周圍環(huán)境如圖2所示，GNSS接收機(jī)為天寶NETR9型，天線型號為TRM57971.00，數(shù)據(jù)采樣間隔為30 s。位于BRST測站以北500 m的布雷斯特港驗潮站提供實測數(shù)據(jù)，可從REFMAR官網(wǎng)下載時間分辨率為1 min的驗潮數(shù)據(jù)。

圖2 BRST測站周圍環(huán)境

本實驗利用BRST觀測站2021年6月1日—15日連續(xù)15 d的GPS觀測數(shù)據(jù)，以均誤差(ME)、RMSE和相關(guān)系數(shù)(R)三項精度指標(biāo)來對比分析SSA法獲得的有效殘差序列和二次多項式擬合獲得的有效殘差序列反演海平面高度的精度：

(8)

式中，n為反演數(shù)據(jù)總數(shù)；X1(i)，X2(i)分別為驗潮站實測數(shù)據(jù)和反演海平面高度。

本文選取2021年6月15日GPS衛(wèi)星08號L2波段、衛(wèi)星高度角在5°～20°范圍內(nèi)、方位角在130°～270°范圍內(nèi)的原始SNR序列。奇異值貢獻(xiàn)率如圖3所示，當(dāng)特征值個數(shù)為2時，奇異值貢獻(xiàn)率變化速率突然降低，出現(xiàn)了明顯的拐點，根據(jù)特征值變化率以及大量的實驗，最終選擇第2～17個特征分量重構(gòu)有效殘差序列。經(jīng)SSA法分解和重構(gòu)的結(jié)果如圖4所示，橫軸表示高度角正弦值，縱軸表示振幅。其中IMF1為選取的貢獻(xiàn)率最大的特征分量所構(gòu)成的趨勢項；IMF2為通過特征值變化率選取的第2～17個特征分量所構(gòu)成的有效殘差序列；IMF3為選取的第18～25個特征分量所構(gòu)成的噪聲。

圖3 奇異值貢獻(xiàn)率

(a)原始SNR序列

將IMF2作為有效殘差序列，并對其進(jìn)行L-S譜分析得到頻率f。圖5是當(dāng)天用SSA得到的殘差序列。圖6是通過SSA得到的SNR殘差序列的Lomb-Scargle譜分析圖，橫軸表示天線相位中心到海平面的距離，縱軸表示SNR反射信號頻譜振幅，其中,振幅的最高峰值所對應(yīng)頻率為149.530 4 Hz，和實際的天線相位中心到海平面的垂直反射距離18.26 m相符合。

圖5 SSA獲得的殘差序列

圖6 Lomb-Scargle 譜分析

本實驗統(tǒng)計了GPS衛(wèi)星在2021年6月1日—15日的觀測數(shù)據(jù)，圖7為SSA法獲得的有效殘差序列和二次多項式擬合獲得的有效殘差序列反演的海平面高度結(jié)果對比圖。

圖7 2種方法反演海平面高度與實測海平面高度

驗潮站的實測海平面高度的變化有很明顯的日周期性，主要是因為受潮汐影響。從圖7中可以看出，SSA法獲得的有效殘差序列和二次多項式擬合獲得的有效殘差序列二者的反演結(jié)果和驗潮站的實測結(jié)果基本一致，驗證了這2種方法反演海平面高度變化的可行性，但仍有一些反演的海平面高度和驗潮站的高度存在較大偏差，這主要是由波浪在某一時間段發(fā)生急劇變化所導(dǎo)致，同時這也是利用GNSS-MR技術(shù)反演海平面高度的結(jié)果與驗潮站監(jiān)測結(jié)果存在偏差的主要原因。為了進(jìn)一步驗證SSA法的有效性，實驗選取在連續(xù)暴雨天氣和外界干擾較嚴(yán)重的環(huán)境下，GPS衛(wèi)星在2021年6月17日—21日的觀測數(shù)據(jù)進(jìn)行反演，反演結(jié)果如圖8所示。

圖8 連續(xù)暴雨天氣下，2種方法反演海平面高度與實測海平面高度

表1為SSA法獲得的有效殘差序列和二次多項式擬合獲得的有效殘差序列反演的海平面高度結(jié)果與驗潮站實測結(jié)果的精度統(tǒng)計。結(jié)果表明，利用SSA法很好地對原始SNR序列進(jìn)行了去噪處理，進(jìn)而得到與驗潮站監(jiān)測值更為接近的反演海平面高度。結(jié)合圖7、圖8和表1可知，通過SSA法獲得的有效殘差序列反演的海平面高度的RMSE比二次多項式擬合降低了6.5%，其相關(guān)系數(shù)為0.97，比二次多項式擬合提高了0.01；而在連續(xù)暴雨天氣和外界干擾較嚴(yán)重的環(huán)境下，通過SSA法獲得的有效殘差序列反演的海平面高度的RMSE比二次多項式擬合降低了14%，其相關(guān)系數(shù)為0.98，表明SSA法獲得的有效殘差序列用于海平面高度反演具有良好的精度。

表1 二次多項式擬合和SSA反演的海平面高度與驗潮站海平面高度的精度統(tǒng)計

4 結(jié)束語

本文利用BRST測站的數(shù)據(jù)進(jìn)行GNSS-MR海平面高度反演，通過對原始SNR序列進(jìn)行SSA，將有用的信號與噪聲分離，提取出高質(zhì)量的殘差序列，進(jìn)而利用GNSS-MR技術(shù)反演海平面高度。結(jié)果表明：

① SSA得到的有效殘差序列能夠準(zhǔn)確地反映原始SNR信號的變化趨勢，有效地對原始SNR信號進(jìn)行了去噪，從而得到精度更高的海平面高度。

② SSA方法能夠更加精確地提取SNR直反信號的干涉信號，反演的海平面高度與驗潮站高度對比，相關(guān)系數(shù)為0.97，RMSE降低了6.5%，而在連續(xù)暴雨天氣和外界干擾較嚴(yán)重的環(huán)境下，RMSE降低了14%，相關(guān)系數(shù)為0.98，驗證了SSA方法應(yīng)用于海平面測高的有效性。