基于RMSprop的高精度GNSS-R鏡點快速預(yù)測算法

盛志超，周 勃，秦 瑾，顧斯祺，孟婉婷

(上海航天電子技術(shù)研究所，上海 201109)

0 引言

全球?qū)Ш叫l(wèi)星系統(tǒng)反射信號(Global Navigation Satellite System Reflectometry，GNSS-R)技術(shù)作為外源助動式新興遙感技術(shù)，以其寬刈幅、高時空分辨率、低功耗、低成本和低重量等特點受到了國內(nèi)外遙感領(lǐng)域的廣泛關(guān)注[1]。GNSS-R技術(shù)能夠進(jìn)行中尺度海洋環(huán)境要素的探測，海面測高是其主要應(yīng)用方向之一。如何實現(xiàn)5cm的海面測高精度，一直是GNSS-R研究的重點和難點。

GNSS-R本質(zhì)上是雙基雷達(dá)，GNSS衛(wèi)星發(fā)射的GNSS信號經(jīng)過地球表面反射后，被GNSS-R接收機(jī)接收和處理[2]。但GNSS信號經(jīng)過大氣、電離層和反射面的衰減和散射之后，信號及其微弱，為了能夠更好的捕獲反射信號，GNSS-R接收機(jī)往往只接收來自以鏡面反射點為中心的閃耀區(qū)的信號[3]。鏡面反射點預(yù)測的精確與否對于高精度海面測高應(yīng)用尤為重要，一方面，鏡面反射點是整個測高反演過程中的坐標(biāo)基準(zhǔn)；另一方面，鏡面反射點的精度直接影響了反射信號的信噪比，進(jìn)而影響海面高度的反演精度。

目前，常用的鏡面反射點預(yù)測算法有:Gleason算法[4]、S.C.Wu算法[5]、線段二分法[6]。其中Gleason算法利用向量共線的思想，算法效果不理想，精度差，迭代次數(shù)多；S.C.Wu算法精度高，但搜索完成之后需要進(jìn)行復(fù)雜的校正；線段二分法模型簡單，迭代次數(shù)少，但沒有解決徑向和法向不統(tǒng)一問題[7]。

本文針對鏡面反射點預(yù)測問題，分析利用梯度下降法和AdaGrad算法預(yù)測鏡點存在的不足，提出了一種基于RMSprop的高精度鏡面反射點預(yù)測算法，該算法相比于現(xiàn)有的算法，通過平滑梯度平方，獲得更高的預(yù)測精度，適合對精度要求很高的GNSS-R海面測高應(yīng)用中。

1 鏡面反射點定義

GNSS-R技術(shù)中GNSS衛(wèi)星、地球表面和GNSS接收機(jī)的空間幾何位置關(guān)系如圖 1所示。在理想情況下，鏡面反射點滿足以下條件：① 鏡面反射點在地球表面；② 發(fā)射機(jī)-鏡面反射點-接收機(jī)是所有反射路徑中的最短路徑；③ 滿足斯涅耳定律——入射角等于反射角；④ 鏡面反射點與接收機(jī)、發(fā)射機(jī)和地心共面[8]。

圖1 GNSS-R空間幾何位置關(guān)系Fig.1 Spatial geometric relations of GNSS-R

圖2是GNSS空間幾何位置關(guān)系的二維視圖。

圖2 幾何關(guān)系二維視圖Fig.2 2D view of geometric relations

圖中，R表示接收機(jī)，T表示發(fā)射機(jī)，S表示鏡面反射點，O表示地心，M表示OC所在直線與RT所在直線的交點，U表示R關(guān)于OM的對稱點，C表示M關(guān)于RU的對稱點，αr、αt分別表示信號的入射角和反射角。

根據(jù)上面對鏡面反射點條件的描述，鏡面反射點預(yù)測問題就是一個帶非線性約束條件的非線性優(yōu)化問題，數(shù)學(xué)表達(dá)式可以歸納為：

(1)

式中，f(S)是鏡面反射點預(yù)測的目標(biāo)函數(shù)；R表示接收機(jī)；T表示發(fā)射機(jī)；S表示地球表面的鏡面反射點；Sx、Sy、Sz為S在地心坐標(biāo)系下對應(yīng)的三軸坐標(biāo)；Rx、Ry、Rz為R在地心坐標(biāo)系下對應(yīng)的三軸坐標(biāo)；Tx、Ty、Tz為T在地心坐標(biāo)系下對應(yīng)的三軸坐標(biāo)；Ox、Oy、Oz為地心坐標(biāo)系原點的三軸坐標(biāo)；αt、αr分別表示信號的入射角和反射角。

2 基于RMSprop的鏡面反射點預(yù)測算法原理

2.1 梯度下降法及其問題

梯度下降法就是利用梯度的方向是函數(shù)變化最快的方向這一準(zhǔn)則進(jìn)行最優(yōu)化求解的[9][10]。針對反射路徑長度f(S)，也可以用梯度下降法來實現(xiàn)對鏡面反射點的預(yù)測：

f(S)=|S-T|+|R-S|。

(2)

即信號從GNSS衛(wèi)星發(fā)射經(jīng)鏡面反射點反射后到達(dá)GNSS-R接收機(jī)所走過的路徑長度。

梯度下降法預(yù)測鏡點的具體求解流程如圖 3所示。首先求得一個粗略的鏡面反射點估計值S1，然后求在該點處反射路徑長度函數(shù)f(S)的梯度，記為f(S1)，隨后利用f(S1)更新鏡面反射點位置S2，更新方程為：

S2=S1-K×f(S1)，

(3)

圖3 梯度下降法示意Fig.3 Schematic diagram of gradient descent method

式中，K為學(xué)習(xí)率。

再計算S2處f(S)的梯度f(S2)，更新S3的位置……,對于第n次鏡面反射點預(yù)測，Sn更新表達(dá)式為：

Sn=Sn-1-K×f(Sn)。

(4)

梯度下降法會一直計算到前后2個預(yù)測值之間的距離小于某個閾值，則認(rèn)為到達(dá)鏡面反射點，算法結(jié)束。

在梯度下降迭代的過程中，學(xué)習(xí)率決定了每一步沿梯度負(fù)方向前進(jìn)的長度。梯度下降法采用了固定的學(xué)習(xí)率，并不會隨著預(yù)測的深入進(jìn)行自適應(yīng)調(diào)節(jié)。所以梯度下降法非常依賴學(xué)習(xí)率K的選擇，學(xué)習(xí)率太小，收斂慢，學(xué)習(xí)率太大，會導(dǎo)致算法在鏡面反射點附近發(fā)生震蕩，無法收斂[11]。

2.2 Adagrad算法及其問題

為了解決梯度下降法學(xué)習(xí)率無法調(diào)節(jié)的問題，引入了Adagrad算法[12]。Adagrad算法利用累積梯度平方和對學(xué)習(xí)率進(jìn)行約束，即將式(4)的鏡面反射點更新表達(dá)式改為：

(5)

隨著算法遞推的深入，梯度平方和越來越大，每次移動的步長越來越小，最后到達(dá)預(yù)測值。但是，在Adagrad算法執(zhí)行的中后期，分母上梯度平方和的累加將會越來越大，使得每次更新步長趨近于0，預(yù)測會提前結(jié)束。

2.3 基于RMSprop的鏡面反射點預(yù)測算法

為了解決Adagrad算法中學(xué)習(xí)率下降過快的問題，需要降低式(5)分母的平方和項，RMSprop通過將平方和變?yōu)橐苿蛹訖?quán)平均平方和來解決，算法鏡面反射點的更新表達(dá)式為：

(6)

式中，τn=τn-1×β+(1-β)×f(Sn)2；τn是加權(quán)平均梯度和；β是加權(quán)平均指數(shù)。

式(6)中，τn完成的是對梯度的平方的一次平滑處理。在更新鏡面反射點Sn+1時，先用梯度f(Sn)除以的值，相當(dāng)于對梯度做了一次歸一化。在鏡面反射點預(yù)測的過程中，豎直方向的梯度比水平方向上的梯度大很多，此時豎直方向的τn也很大，當(dāng)f(Sn)除以之后，歸一化梯度變小，豎直方向的步長變小了；而水平方向梯度小，相應(yīng)的τn也小，水平方向上歸一化的梯度就大了，步長也變大了。

基于RMSprop的鏡面反射點預(yù)測算法作為梯度下降法的優(yōu)化算法，一方面能夠解決Adagrad算法遞推過程學(xué)習(xí)率下降過快的問題，另一方面可以調(diào)整不同維度上的步長，加快收斂速度，最終實現(xiàn)對鏡面反射點的預(yù)測。

3 基于RMSprop的鏡面反射點預(yù)測算法流程

根據(jù)GNSS信號中攜帶的星歷數(shù)據(jù)，可以解算出在WGS-84坐標(biāo)系統(tǒng)中接收機(jī)RWGS-84和發(fā)射機(jī)TWGS-84的位置，作為已知條件利用如下步驟進(jìn)行鏡面反射點的預(yù)測。

基于RMSprop的高精度GNSS-R鏡面反射點預(yù)測算法流程如圖 4所示，其處理流程主要包括準(zhǔn)球映射、鏡面反射點粗估和鏡面反射點精估3個內(nèi)容。

圖4 算法總流程Fig.4 SOverall flow of algorithm

3.1 準(zhǔn)球映射

將WGS-84坐標(biāo)系統(tǒng)映射成單位圓，WGS-84坐標(biāo)系統(tǒng)下的接收機(jī)位置RWGS-84和發(fā)射機(jī)位置TWGS-84將分別被映射成R和T。

3.2 鏡面反射點的粗估

根據(jù)如圖 2所示的空間幾何位置關(guān)系，可以得到如下數(shù)學(xué)關(guān)系：

|RM|=|MU|=|UC|，

(7)

|RM|/|RT|=|RM|/(|RM|+|MT|)=|CU|/(|CU|+|MT|)=

|SR|/(|SR|+|ST|)≈|OR|/(|OR|+|OT|)。

(8)

M點的估計值為：

M=R+RM≈R+[|OR|/(|OR|+|OT|)]×RT。

(9)

S1是M的星下點，則：

S1=(OM/|OM|)×Re，

(10)

式中，Re為地球半徑，因為進(jìn)行了準(zhǔn)球映射，此時地球是單位圓，即Re=1。

3.3 RMSprop鏡面反射點精估

設(shè)置全局學(xué)習(xí)率ε，鏡面反射點的初始值S1，小常數(shù)δ(典型值10-8，用于被小數(shù)除時的數(shù)值穩(wěn)定)，加權(quán)平均系數(shù)β。

① 初始化梯度累積變量τ1=0；

f(Sn)=(Sn-T)/|T-Sn|+(Sn-R)/|R-Sn|，

(11)

③ 計算累積梯度平方和τn:

τn=τn-1×β+(1-β)×f(Sn)2，

(12)

(13)

④ 計算更新增量ΔS；

⑤ 更新沿梯度方向下降后的預(yù)測值S′:

S′=Sn+ΔS，

(14)

⑥ 取S′的星下點作為當(dāng)前鏡面反射點預(yù)測值Sn+1:

Sn+1=(OS′/|OS′|)×Re(S′)，

(15)

⑦ 計算當(dāng)前鏡面反射點預(yù)測值與之前的歐氏距離作為更新的幅度，如果大于預(yù)設(shè)門限則返回步驟①，反之則退出循環(huán)，以當(dāng)前鏡面反射點預(yù)測值Sn+1作為鏡面反射點最優(yōu)值進(jìn)行反射信號的捕獲跟蹤。

AMP=|Sn+1-Sn|。

(16)

4 仿真與分析

本文利用Matlab和STK仿真工具進(jìn)行如圖 5所示的聯(lián)合仿真。利用Matlab對STK進(jìn)行軌道參數(shù)配置后，STK完成衛(wèi)星可見性分析，將分析結(jié)果導(dǎo)入Matlab，在Matlab中執(zhí)行鏡面反射點計算并分析算法性能。本次仿真實驗中，GNSS星座采用GPS星座和北斗星座共58顆GNSS衛(wèi)星，LEO衛(wèi)星的軌道高度為500 km，接收入射角最大45°的反射信號，仿真了在軌運行24 h(2019年3月17日—18日)的數(shù)據(jù)，采樣時間間隔為1 s。利用GPS的PRN1衛(wèi)星與LEO所得到的9段數(shù)據(jù)共12 341個時間點進(jìn)行下面的仿真。

圖5 聯(lián)合仿真結(jié)構(gòu)框架Fig.5 Co-simulation framework

4.1 粗估的S1的精度

利用本文算法預(yù)測每個時刻的鏡面反射點，將最終預(yù)測結(jié)果作為鏡面反射點預(yù)測值，把粗估的鏡面反射點值與預(yù)測值進(jìn)行比較，分析粗估的精度。

鏡面反射點粗估性能如圖6所示。

(a) 鏡面反射點粗估誤差

(b) 各時刻GNSS衛(wèi)星與LEO衛(wèi)星距離圖6 鏡面反射點粗估性能Fig.6 Performance of specular point by coarse estimation

圖6(a)中，粗估值與預(yù)測值相距最近65 m，最遠(yuǎn)172 320 m，均值68 290 m。圖 6(b)是各個時刻GNSS衛(wèi)星與LEO衛(wèi)星的距離。圖6(a)中曲線的起伏與圖6(b)中曲線的起伏基本一致，說明粗估的精度與GNSS衛(wèi)星和LEO衛(wèi)星所處的幾何位置有關(guān)，二者距離越近，衛(wèi)星高度角就越高，粗估的精度就越高；二者距離越遠(yuǎn)，衛(wèi)星高度角就越低，粗估的精度就越低。

4.2 單個數(shù)據(jù)仿真性能比較

選取了上述12 341個時間點內(nèi)某一時間點GNSS衛(wèi)星和LEO衛(wèi)星的位置，分別利用Gleason算法、S.C.Wu算法、線段二分法和本文算法進(jìn)行鏡面反射點預(yù)測，分析其收斂性能和算法精度。收斂性能的評判依據(jù)是算法的迭代次數(shù)，算法精度的依據(jù)是GNSS信號從GNSS衛(wèi)星發(fā)射經(jīng)鏡面反射點反射到達(dá)LEO衛(wèi)星所走過的路徑長度。

4.2.1 算法精度

圖 7是4種算法路徑長度的比較。從圖中可以看出，線段二分法、S.C.Wu算法和Gleason算法得到的反射路徑長度非常接近。其中，S.C.Wu算法得到的反射路徑比線段二分法和Gleanson算法得到的反射路徑長度更小。而本文提出的算法得到的反射路徑長度最短。這說明本文提出的算法獲得的鏡面反射點與真實鏡面反射點距離最近，計算精度最高。

圖7 4種算法得到的路徑長度Fig.7 Path length obtained by four algorithms

4.2.2 收斂性能

將仿真中本文算法計算得到的鏡面反射點作為參考，對算法的收斂性能進(jìn)行分析。圖 8為4種鏡面反射點算法收斂性能的比較圖。通過圖 8分析可以得到，線段二分法在大約10次迭代后收斂至最小點，其曲線下降過程比較震蕩；S.C.Wu算法在大約10次迭代后，曲線達(dá)到最小值，隨后曲線迅速上升且無法收斂，該現(xiàn)象是由未經(jīng)過坐標(biāo)校正引起的；Gleason算法在收斂過程很慢，曲線較為平滑，經(jīng)過大約100次迭代后，曲線開始收斂，其迭代次數(shù)是最多的；本文提出的算法在經(jīng)過約10次迭代后達(dá)到預(yù)測點，迭代次數(shù)最少，運算速度最快。

圖8 4種算法的收斂性能Fig.8 Convergence performance of four algorithms

4.3 統(tǒng)計數(shù)據(jù)仿真性能比較

將12 341個時間點每個時刻GNSS衛(wèi)星和LEO衛(wèi)星的位置分別利用線段二分法和基于RMSprop的鏡面反射點預(yù)測算法進(jìn)行預(yù)測，將結(jié)果進(jìn)行統(tǒng)計意義上的收斂性能比較和精度比較。

4.3.1 收斂性能

圖 9是本文算法和線段二分法的效率統(tǒng)計圖。其中橫軸表示對單個數(shù)據(jù)點進(jìn)行計算，收斂到預(yù)測值所需的迭代次數(shù)?？v軸表示需要相應(yīng)迭代次數(shù)的數(shù)據(jù)點個數(shù)。通過比較可以得到，本文算法中所有數(shù)據(jù)點計算的迭代次數(shù)區(qū)間為9～15，而線段二分法中所有數(shù)據(jù)點計算的迭代次數(shù)區(qū)間為9～19，且數(shù)據(jù)點落在17～18區(qū)間的最多。利用均值和方差對本文提出的算法和線段二分法進(jìn)行統(tǒng)計意義分析，得到本文算法與線段二分法的收斂性能統(tǒng)計表，如表1所示。

表1 本文算法與線段二分法的收斂性能統(tǒng)計表Tab.1 Convergence performance list of the proposed algorithm and the dichotomy of line segment

(a) RMSprop算法的迭代次數(shù)統(tǒng)計

(b) 線段二分法的迭代次數(shù)統(tǒng)計圖9 2種算法的迭代次數(shù)直方圖Fig.9 Iteration times histogram of two algorithms

分析表 1可以得到，本文提出的算法運算效率略優(yōu)于線段二分法。

4.3.2 算法精度

圖10是本文算法和線段二分法的路徑長度比較圖。

(a) 本文算法和線段二分法路徑長度比較

(b) 2種算法路徑長度差值圖10 2種算法的路徑長度比較Fig.10 Path length comparison between two algorithms

圖10(a)中，2種算法的路徑長度曲線基本重合，其中RMSprop算法的路徑長度要小一點，為了進(jìn)行更直觀的分析，本文將線段二分法路徑長度與RMSprop算法路徑長度進(jìn)行差值，得到如圖10(b)所示的路徑長度差值圖，可以看出，RMSprop算法要比線段二分法算出來的路徑短18km左右。經(jīng)統(tǒng)計分析可以得到，在觀察數(shù)據(jù)點樣本中，任何時刻，本文算法的精度都更高。

綜上所述，無論是算法精度還是收斂性能，本文算法都明顯優(yōu)于傳統(tǒng)的Gleason算法、S.C.Wu算法和線段二分法。

5 結(jié)束語

本文通過研究鏡面反射點定義、梯度下降法及其優(yōu)化算法的相關(guān)知識，提出了一種基于RMSprop的鏡面反射點預(yù)測算法。作為機(jī)器學(xué)習(xí)中梯度下降法優(yōu)化算法的一種，該算法解決了梯度下降法的收斂問題，也不會和AdaGrad算法一樣提前收斂，專門用于對鏡面反射點精確度要求很高的GNSS-R海面測高應(yīng)用中。為了驗證算法的正確性，利用Matlab和STK進(jìn)行聯(lián)合仿真，通過與常用算法的比較，結(jié)果表明，本文提出的算法可以對鏡面反射點進(jìn)行正確預(yù)測，該算法不僅精度高，而且收斂性能好，運算效率高，適用于GNSS-R海面測高應(yīng)用。