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

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

0 引言

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

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

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

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

1 鏡面反射點(diǎn)定義

GNSS-R技術(shù)中GNSS衛(wèi)星、地球表面和GNSS接收機(jī)的空間幾何位置關(guān)系如圖 1所示。在理想情況下，鏡面反射點(diǎn)滿足以下條件：① 鏡面反射點(diǎn)在地球表面；② 發(fā)射機(jī)-鏡面反射點(diǎn)-接收機(jī)是所有反射路徑中的最短路徑；③ 滿足斯涅耳定律——入射角等于反射角；④ 鏡面反射點(diǎn)與接收機(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表示鏡面反射點(diǎn)，O表示地心，M表示OC所在直線與RT所在直線的交點(diǎn)，U表示R關(guān)于OM的對(duì)稱點(diǎn)，C表示M關(guān)于RU的對(duì)稱點(diǎn)，αr、αt分別表示信號(hào)的入射角和反射角。

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

(1)

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

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

2.1 梯度下降法及其問題

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

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

(2)

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

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

S2=S1-K×f(S1)，

(3)

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

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

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

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

(4)

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

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

2.2 Adagrad算法及其問題

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

(5)

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

2.3 基于RMSprop的鏡面反射點(diǎn)預(yù)測(cè)算法

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

(6)

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

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

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

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

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

基于RMSprop的高精度GNSS-R鏡面反射點(diǎn)預(yù)測(cè)算法流程如圖 4所示，其處理流程主要包括準(zhǔn)球映射、鏡面反射點(diǎn)粗估和鏡面反射點(diǎn)精估3個(gè)內(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 鏡面反射點(diǎn)的粗估

根據(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點(diǎn)的估計(jì)值為：

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

(9)

S1是M的星下點(diǎn)，則：

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

(10)

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

3.3 RMSprop鏡面反射點(diǎn)精估

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

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

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

(11)

③ 計(jì)算累積梯度平方和τn:

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

(12)

(13)

④ 計(jì)算更新增量ΔS；

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

S′=Sn+ΔS，

(14)

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

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

(15)

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

AMP=|Sn+1-Sn|。

(16)

4 仿真與分析

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

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

4.1 粗估的S1的精度

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

鏡面反射點(diǎn)粗估性能如圖6所示。

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

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

圖6(a)中，粗估值與預(yù)測(cè)值相距最近65 m，最遠(yuǎn)172 320 m，均值68 290 m。圖 6(b)是各個(gè)時(shí)刻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 單個(gè)數(shù)據(jù)仿真性能比較

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

4.2.1 算法精度

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

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

4.2.2 收斂性能

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

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

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

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

4.3.1 收斂性能

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

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

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

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

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

4.3.2 算法精度

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

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

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

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

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

5 結(jié)束語

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