劉睿國(guó)，劉立龍，吳 晗，薛張芳，吳昊艦，張 志

(1.桂林理工大學(xué) 測(cè)繪地理信息學(xué)院，廣西 桂林 541006；2.廣西空間信息與測(cè)繪重點(diǎn)實(shí)驗(yàn)室，廣西 桂林 541006)

0 引言

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

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

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

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

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

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

Ad?Ar,

(1)

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

(2)

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

(3)

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

2 SSA的原理與方法

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

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

假設(shè)混有噪聲的SNR序列Li的序列長(zhǎng)度為N，對(duì)其進(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)

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

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

(7)

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

3 實(shí)驗(yàn)分析

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

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

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

(8)

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

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

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

(a)原始SNR序列

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

圖5 SSA獲得的殘差序列

圖6 Lomb-Scargle 譜分析

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

圖7 2種方法反演海平面高度與實(shí)測(cè)海平面高度

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

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

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

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

4 結(jié)束語

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

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

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