胡媛，鐘李程，陳行楊，顧旺旺，劉衛(wèi)

( 1.上海海洋大學(xué)工程學(xué)院上海， 201306；2.上海海事大學(xué)商船學(xué)院上海， 200135)

0 引言

近年來(lái)，全球氣候變暖，兩極冰川融化導(dǎo)致海平面的上升.隨著海平面的上升，海拔低的沿海地帶將被淹沒(méi)，數(shù)百萬(wàn)居民流離失所.此外，風(fēng)暴、極端海浪和颶風(fēng)等自然災(zāi)害頻繁發(fā)生，因此，監(jiān)控海平面的變化對(duì)于沿海人民的生活有著重大的意義.傳統(tǒng)監(jiān)測(cè)海平面的方法是通過(guò)岸邊驗(yàn)潮儀進(jìn)行監(jiān)測(cè)，但是傳統(tǒng)驗(yàn)潮儀測(cè)量的是相對(duì)海面高度，其測(cè)量結(jié)果是絕對(duì)海面和地殼垂直方向的形變共同組成的相對(duì)海面高度.測(cè)量絕對(duì)海面高度則需要通過(guò)全球衛(wèi)星導(dǎo)航系統(tǒng)(GNSS)抵消地殼垂直位移的影響，并且傳統(tǒng)驗(yàn)潮儀在風(fēng)浪的作用下，測(cè)量?jī)x器可能會(huì)出現(xiàn)零點(diǎn)漂移現(xiàn)象[1].

自Hall 等[2]提出全球衛(wèi)星導(dǎo)航反射信號(hào)(GNSS-R)的概念以來(lái)，GNSS-R 被廣泛應(yīng)用于遙感方面.如土壤濕度檢測(cè)[3-8]、降雪深度檢測(cè)[9-13]、海冰檢測(cè)[14-17]、海面溢油檢測(cè)[18-21]，同樣GNSS-R 也被運(yùn)用于監(jiān)測(cè)海平面.1993年，Martin-Neira[22]首次提出被動(dòng)反射和干涉測(cè)量系統(tǒng)，即PARIS(passive reflectometry and interferometry system)的概念，PARIS的主要理論思想是利用從海面反射過(guò)來(lái)的GPS信號(hào)測(cè)量海面高度，并且同時(shí)測(cè)量多個(gè)反射點(diǎn)，利用多個(gè)衛(wèi)星信號(hào)來(lái)達(dá)到很高的測(cè)量精度.

自PARIS概念提出，許多研究人員在GNSS-R海面測(cè)高方面做了大量的研究.目前為止，利用GPS信號(hào)測(cè)量海面高度可以大致分為兩類(lèi)，相位差分析[23]和信噪比(SNR)分析[24].相位差分析是使用兩根天線(xiàn)來(lái)確定直射信號(hào)和來(lái)自海面的反射信號(hào)之間的相位差，包括碼相位差法和載波相位差法[25]，從而得出兩個(gè)信號(hào)的路徑差來(lái)確定海面高度.其中直射信號(hào)為右圓極化(RHCP)信號(hào)，反射信號(hào)的極化特性通常情況下從RHCP轉(zhuǎn)變?yōu)樽髨A極化(LHCP)，這種改變由反射面的介電特性和發(fā)射機(jī)-反射面-接收機(jī)的幾何結(jié)構(gòu)決定[26].該技術(shù)需要運(yùn)用雙天線(xiàn)的接收機(jī)，這種接收機(jī)結(jié)構(gòu)復(fù)雜，成本大并且需要專(zhuān)業(yè)人員安裝.SNR 分析方法利用單根天線(xiàn)接收直射信號(hào)和反射信號(hào)所產(chǎn)生的干涉模式技術(shù)(IPT)[25]，利用所獲得的SNR 與衛(wèi)星仰角的函數(shù)關(guān)系來(lái)確定海面高度.運(yùn)用SNR 數(shù)據(jù)分析只需要單根天線(xiàn)來(lái)測(cè)量海面高度，并且對(duì)風(fēng)浪條件具有更好的魯棒性.但是目前該方法的測(cè)量精度還不如利用相位差測(cè)量方法的精度高，所以許多研究者正在研究運(yùn)用SNR 數(shù)據(jù)提高測(cè)量精度的算法.本文圍繞SNR 數(shù)據(jù)測(cè)量海面高度中運(yùn)用的算法進(jìn)行展開(kāi)，介紹在利用SNR 數(shù)據(jù)反演海面高度所運(yùn)用的算法，為以后運(yùn)用SNR 數(shù)據(jù)進(jìn)行海面測(cè)高提供一些參考.

1 基于SNR 測(cè)高的基本原理

1.1 GNSS-R 測(cè)高原理

利用GNSS-R 進(jìn)行海面測(cè)高的基本模型如圖1所示，接收機(jī)接收來(lái)自衛(wèi)星的直射信號(hào)和經(jīng)海平面反射的反射信號(hào).因?yàn)榻邮諜C(jī)與衛(wèi)星之間的距離足夠遠(yuǎn)，所以可以認(rèn)為直射信號(hào)和反射信號(hào)為平行關(guān)系.由圖1所示的幾何關(guān)系，可以得到

圖1 GNSS-R 海面測(cè)高幾何關(guān)系圖

式中：Δd為直射信號(hào)與反射信號(hào)的路程差；h為接收機(jī)到海面的垂直高度；θ為衛(wèi)星仰角.

1.2 SNR 信號(hào)的數(shù)據(jù)特點(diǎn)

因?yàn)槔肧NR 測(cè)量海面高度只用了一根天線(xiàn)，所以接收機(jī)的信號(hào)為反射信號(hào)和直射信號(hào)發(fā)生干涉而形成.根據(jù)文獻(xiàn)[27]可知，SNR 與接收機(jī)之間存在以下關(guān)系：

由式(5)中可知，基于SNR 數(shù)據(jù)反演海面高度，只需要求解出去除趨勢(shì)項(xiàng)后的SNR 數(shù)據(jù)的頻率.所以尋求合適的算法去除SNR 信號(hào)的趨勢(shì)項(xiàng)和提取去除趨勢(shì)項(xiàng)后的SNR 信號(hào)的頻率是利用SNR 信號(hào)反演海面高度的關(guān)鍵步驟.去除趨勢(shì)項(xiàng)和提取頻率對(duì)反演的海面高度精度起著至關(guān)重要的作用.

2 基于SNR 信號(hào)的測(cè)高方法

2.1 最小二乘法在反演海面高度的應(yīng)用

2.1.1運(yùn)用最小二乘法去除趨勢(shì)項(xiàng)

最小二乘法是一種被廣泛應(yīng)用于消除趨勢(shì)項(xiàng)的方法[28].如式(2)所示，SNR 數(shù)據(jù)具有一個(gè)二階多項(xiàng)式的趨勢(shì)項(xiàng)，SNR 信號(hào)的采樣數(shù)據(jù)為{xk}(k=1,2,3,···,n).設(shè)趨勢(shì)項(xiàng)多項(xiàng)式函數(shù)

從而得到趨勢(shì)項(xiàng)多項(xiàng)式的系數(shù)，得出趨勢(shì)項(xiàng)擬合曲線(xiàn)，再用原始的SNR 數(shù)據(jù)減去得到的趨勢(shì)項(xiàng)，達(dá)到去除趨勢(shì)項(xiàng)的目的，結(jié)果如圖2所示.

圖2 SNR 信號(hào)和擬合的趨勢(shì)項(xiàng)

2.1.2非線(xiàn)性最小二乘法擬合SNR 信號(hào)

通過(guò)非線(xiàn)性最小二乘法擬合得到的頻率即為去勢(shì)后信噪比數(shù)據(jù)的頻率，根據(jù)式(5)便可以反演出海面高度.

李惟等使用法國(guó)圖盧茲大學(xué)CESBIO實(shí)驗(yàn)室的實(shí)測(cè)數(shù)據(jù)對(duì)上述方法進(jìn)行驗(yàn)證，同時(shí)與LSP方法進(jìn)行對(duì)比.具體結(jié)果如表1所示.

表1 LSP 法和擬合法反演結(jié)果m

2.2 LSP在反演海面高度的應(yīng)用

SNR 信號(hào)在時(shí)間上是均勻采樣的，但是實(shí)際接收機(jī)接收的SNR 信號(hào)因變量是衛(wèi)星仰角的正弦值，這導(dǎo)致SNR 信號(hào)是非均勻分布的，并且接收信號(hào)時(shí)存在噪聲干擾.在運(yùn)用傅里葉變換時(shí)，這些噪聲和不均勻分布會(huì)對(duì)頻譜分析產(chǎn)生影響，造成頻譜圖與SNR 信號(hào)實(shí)際結(jié)果偏差很大.Larson K M 等[24]提出運(yùn)用Lomb標(biāo)準(zhǔn)化周期圖(LSP)方法對(duì)SNR 進(jìn)行分析，獲取該信號(hào)的頻率.

LSP是由Lomb[31]于1976年提出，并在1986年由Scargle[32]進(jìn)行修改完善，LSP有sin 和cos兩組基，并且在時(shí)域方面做了歸一化的處理.因此，LSP方法能夠有效地提取時(shí)序中的弱周期信號(hào)，并且能夠減少非均勻采樣產(chǎn)生的虛假信號(hào)所帶來(lái)的影響.

運(yùn)用LSP方法獲取去除趨勢(shì)項(xiàng)后SNR 的頻率后，根據(jù)式(5)反演得出海面高度，如圖3所示.

圖3 LSP分析結(jié)果

2.3 小波算法在海面測(cè)高中的應(yīng)用

小波分析與傅里葉變換和LSP 相比，傅里葉變換和LSP是全局的變換，得到的是一維的頻域信息，而小波分析是一個(gè)時(shí)間和頻率的局域變換.因此小波分析能夠有效地從信號(hào)中提取信息，解決傅里葉變換所不能解決的諸多問(wèn)題[33].Wang 等[34]使用小波分析對(duì)SNR 信號(hào)進(jìn)行去噪處理，并使用不同長(zhǎng)度的SNR信號(hào)進(jìn)行驗(yàn)證.結(jié)果表明：在較短的SNR 信號(hào)中，使用小波分析去噪精度與未使用小波去噪相比提升了20%.Chen 等[35]和蘇曉容等[36]提出了運(yùn)用小波算法去除噪聲.王杰等[37]提出了運(yùn)用小波算法提取去除趨勢(shì)后的SNR 信號(hào)的頻率來(lái)反演海面高度的方法.

因?yàn)樵趯?shí)際采樣中，SNR 信號(hào)是一個(gè)離散非連續(xù)的信號(hào)，并且使用連續(xù)小波變換計(jì)算量太大，所以在處理實(shí)際的SNR 信號(hào)時(shí)使用的是離散小波變換(DWT).將伸縮系數(shù)a和平移因子b進(jìn)行離散化處理，具體如下：

2.3.1小波去噪原理

1989年，Mallat[38]提出了一種快速DWT 算法，使用小波濾波器分解和重構(gòu)信號(hào)，具體如下：

式中：t是離散的時(shí)間序列號(hào)；f(t)表示原始信號(hào)；j是分解信號(hào)的層數(shù)；H和G分別是小波的低頻和高頻濾波器.

根據(jù)上述描述，將原始信號(hào)通過(guò)高通和低通濾波器，將信號(hào)分解成低頻部分和高頻部分，然后再將分解的低頻部分再分解成低頻部分和高頻部分，從而實(shí)現(xiàn)信號(hào)的多層分解.小波分解的原理圖如圖4所示.

將分解出來(lái)的高頻系數(shù)進(jìn)行相應(yīng)小波系數(shù)處理去除噪聲部分，然后再對(duì)分解的信號(hào)進(jìn)行重構(gòu)處理，這樣就可以得到去除噪聲后的SNR 信號(hào).

2.3.2小波算法提取頻率

運(yùn)用小波變換提取去趨勢(shì)后SNR 信號(hào)頻率，是將趨勢(shì)后的SNR 信號(hào)與經(jīng)過(guò)拉伸或壓縮后的小波進(jìn)行比較，通過(guò)在不同尺度m>0和不同位置n比較原信號(hào)與處理后的小波，可以得到一個(gè)關(guān)于兩個(gè)變量的函數(shù).對(duì)于去除趨勢(shì)后的SNR 信號(hào)，小波變換可以表示為

針對(duì)上述方法，王杰使用美國(guó)SC02站數(shù)據(jù)進(jìn)行驗(yàn)證.運(yùn)用小波分析反演海面高度，標(biāo)準(zhǔn)差略高于LSP法，但是其可以極大地提高反演海面高度的時(shí)間分辨率.具體結(jié)果如表2所示.

表2 LSP法、加窗LSP法和小波法反演結(jié)果

2.4 B樣條函數(shù)和卡爾曼濾波在海面測(cè)高的應(yīng)用

2.4.1 B樣條函數(shù)在海面測(cè)高的應(yīng)用

2016年Strandberg[29]提出使用B樣條函數(shù)對(duì)海面高度進(jìn)行建模，該方法考慮到監(jiān)測(cè)的海平面是連續(xù)變化的，可以將海面高度變化假設(shè)為一個(gè)平滑函數(shù).將海面高度變化作為分段線(xiàn)性模型進(jìn)行建?？赡苁亲顬楹?jiǎn)單的方法，但是分段線(xiàn)性模型在計(jì)算一階導(dǎo)數(shù)時(shí)會(huì)導(dǎo)致節(jié)點(diǎn)處的不連續(xù)，不能對(duì)海面高度提供足夠的可變性用來(lái)考慮次日和長(zhǎng)期的海平面變化.B樣條函數(shù)的基本函數(shù)定義為

2.4.2卡爾曼濾波在海面測(cè)高的應(yīng)用

之前通過(guò)SNR 數(shù)據(jù)反演海面高度的方法是提取接收機(jī)接收的一段時(shí)間內(nèi)的SNR 信號(hào)數(shù)據(jù)，提取出含有高度信息的頻率，然后反演海面高度，該方法只能得出該時(shí)間內(nèi)的平均海面高度，并不能提供實(shí)時(shí)的海面高度數(shù)據(jù).卡爾曼濾波[39]是當(dāng)前應(yīng)用最廣的一種動(dòng)態(tài)數(shù)據(jù)處理方法，但原始的卡爾曼濾波是建立在線(xiàn)性系統(tǒng)的基礎(chǔ)上的，對(duì)非線(xiàn)性系統(tǒng)效果并不理想[40].在2000年，Merwe等[41]提出無(wú)損卡爾曼濾波處理(UKF)更新步驟中的非線(xiàn)性.在2019年，Strandberg等[42]提出了運(yùn)用無(wú)損卡爾曼濾波實(shí)現(xiàn)實(shí)時(shí)監(jiān)測(cè)海面高度.該方法采用了一個(gè)預(yù)測(cè)方程，可以預(yù)測(cè)狀態(tài)向量和協(xié)方差矩陣的變換，即：

采用二次B樣條函數(shù)來(lái)描述海面高度隨時(shí)間的變化，我們需要估計(jì)一組縮放系數(shù)αi就可以在任何時(shí)間預(yù)測(cè)海面高度，即

B樣條函數(shù)定義如式(19)～(20)所示，如式(34)所示，想要知道某個(gè)時(shí)間點(diǎn)所對(duì)應(yīng)的高度，只需要知道該時(shí)間點(diǎn)附件的4個(gè)節(jié)點(diǎn).隨著卡爾曼濾波隨著時(shí)間的推移，4個(gè)節(jié)點(diǎn)B樣條系數(shù)和相應(yīng)的協(xié)方差會(huì)不斷的更新，直到引入tj+1處的新節(jié)點(diǎn).然后將最后一個(gè)節(jié)點(diǎn)從狀態(tài)向量中刪除，將前3個(gè)節(jié)點(diǎn)往后移，為新的節(jié)點(diǎn)釋放空間.這樣卡爾曼濾波便可以連續(xù)且實(shí)時(shí)的估計(jì)海面的高度.

針對(duì)上述方法，Strandberg 使用GTGU 站數(shù)據(jù)進(jìn)行驗(yàn)證.卡爾曼濾波法誤差為1.98 cm，并且隨著時(shí)間的推移，誤差縮小為1.48 cm.同樣的實(shí)驗(yàn)數(shù)據(jù)，使用LSP方法的誤差為8.84 cm.

3 結(jié)束語(yǔ)

GNSS-R 遙感作為一個(gè)新興的領(lǐng)域，在各個(gè)領(lǐng)域中有著廣泛的應(yīng)用.本文主要介紹了GNSS-R 在海面測(cè)高的研究進(jìn)展，著重介紹了使用SNR 數(shù)據(jù)反演海面高度的研究進(jìn)展.主要利用了最小二乘法，LSP法、小波算法、卡爾曼濾波、B樣條函數(shù)等算法在處理SNR 數(shù)據(jù)進(jìn)行海面測(cè)高的應(yīng)用.目前，基于SNR數(shù)據(jù)測(cè)量的方法比之前的基于相位差的方法所需要的設(shè)備更為簡(jiǎn)單且安裝更加簡(jiǎn)便，更具有應(yīng)用潛力.當(dāng)前，人工神經(jīng)網(wǎng)絡(luò)發(fā)展迅速，在諸多領(lǐng)域得到廣泛的應(yīng)用.神經(jīng)網(wǎng)絡(luò)特別適合變量之間關(guān)系復(fù)雜、耦合度高難以直接數(shù)學(xué)描述的應(yīng)用場(chǎng)景，在GNSS-R 遙感技術(shù)的其他領(lǐng)域，如海冰遙感、土壤濕度檢測(cè)、海面溢油、海面風(fēng)場(chǎng)反演等領(lǐng)域都有成功應(yīng)用.但目前為止，神經(jīng)網(wǎng)絡(luò)并未在海面測(cè)高方面有所應(yīng)用，對(duì)于提取特征頻率的方面還可以利用神經(jīng)網(wǎng)絡(luò)的優(yōu)勢(shì)來(lái)提高反演海面高度的精度.