地基GNSS-IR風速反演原理及方法初探

王笑蕾，何秀鳳，陳 殊，張 勤，宋敏峰

1. 河海大學地球科學與工程學院，江蘇 南京 211100； 2. 長安大學地質工程與測繪學院，陜西 西安 710054

隨著全球導航衛星系統(Global Navigation Satellite Systems,GNSS)的不斷發展，過去被認為是誤差源的多路徑效應，已經被證實可以用來監測反射面物理參數，并逐步發展為GNSS反射遙感(GNSS reflectometry,GNSS-R)技術。其中的地基分支被稱作GNSS干涉遙感(GNSS-interferometry reflectometry,GNSS-IR)技術[1]。目前，GNSS-IR已經被證實可以對水位[2-4]、雪深[5-9]、植被指數[10]、土壤濕度[11-13]、地表凍融[14]等相關參數進行反演。為了擴展地基GNSS技術的監測對象和應用范圍，本文提出了一種基于GNSS-IR的風速反演方法。

利用反射信號進行風速探測在星載GNSS-R中已經實現，并有相關測風衛星播發[15]。星載GNSS-R測風的思路為：不同的風速會導致不同的海面粗糙度，進而導致反射信號的不同。依據反射信號的特征，就可以反推風速。這樣的思路對于地基GNSS-IR同樣適用，即利用反射信號特征推演粗糙度及風場信息。目前，關于GNSS-IR海面粗糙度反演及海面風場反演的研究較少。文獻[16]利用信噪比(signal-to-noise ratio,SNR)的衰減參數對海面有效波高進行反演。文獻[17]利用信噪比的多徑截止高度角參數，提出了一種GNSS-IR風向反演方法。然而，該方法對于站點要求較高，需要站點四面環海，接收到全方位角海域反射信號，才能夠確定風向。相關研究表明：SNR衰減因子及多徑截止高度角可以反映海面粗糙度及風場信息，而且多徑截止高度角參數更為準確[16-17]。

因此，本文將利用多徑截止高度角參數進行風速反演。而獲取該參數的方法目前主要基于曲線擬合方法和經驗閾值，兩種方法對于SNR的質量要求較高[17]。由于小波分析可以探測SNR序列的細節信息，已經有相關研究將小波分析應用到SNR的細節信息探測中。文獻[18—19]利用小波分析提取了SNR序列的瞬時頻率，并利用瞬時頻率來反演瞬時潮位。文獻[18]中提到，利用小波分析方法可以探測有效瞬時頻率消失的高度角，即多路徑振蕩消失的截止高度角。因此，本文將利用小波分析方法探測截止高度角，并研究截止高度角與風速之間的相關關系。

1 原 理

直射信號和反射信號會發生干涉，從而在SNR序列中表現出干涉振蕩。在只考慮一次反射的情況下，低高度角、去除直射信號分量的殘余SNR序列會呈現較為明顯的干涉振蕩，可以用式(1)表示[20]

(1)

式中，h為反射面到天線相位中心的高度，常被稱為有效高度(reflector height,RH)，根據式(1)，h與振蕩頻率f之間的關系為h=(λf)/2；λ為信號波長；k=2π/λ；d為衰減因子；e為高度角；Amp為振幅。然而，在考慮到反射面并不光滑，存在粗糙度時，信號會存在散射現象。此時的殘余SNR序列可以寫為

(2)

式中，下標i表示第i個散射分量的相關參數；M為散射分量的總數；δSNRs可以寫為相干分量(共性部分)δSNRcoh和非相干分量(非共性部分)δSNRincoh；相干分量表現為反射主導的分量，非相干量表現為噪聲n，即式(2)可寫為

δSNRscatter=δSNRcoh+δSNRincoh=Amp·

(3)

式中，當相干分量δSNRcoh更顯著時，SNR序列表現有干涉振蕩；當非相干量(噪聲)δSNRincoh更顯著時，SNR序列表現無干涉振蕩。

1.1 散射模型

(4)

式中

ρ0=sin(vxX)sin(vyY)

vx=k(sinθd-sinθrcosφr)x

vy=k(sinθrsinφr)y

(5)

(6)

當eecutoff，非相干分量更顯著，SNR振蕩消失。根據式(6)，σh越大時，ecutoff越?。沪襤越小時，ecutoff越大。不同的σh會造成不同的ecutoff。從而，根據ecutoff，可以估計σh。

1.2 波浪模型

不同的風速會造成不同的海面粗糙度。海浪譜可以用來描述一定區域的海浪結構，它反映了各分量波的分布情況。若只考慮能量S相對于頻率ω的分布，稱一維譜或波浪頻譜。根據Pierson-Moskowitz (P-M)波浪譜模型，波浪頻譜為[23]

(7)

式中，g是重力加速度；U19.5是距海面19.5 m處的風速；Υ表示海浪隨風速變化的響應程度因子。在成熟開闊海域條件下，Υ=1；在近岸海域，海浪處于破碎狀態，可以根據海面情況調整Υ，0<Υ<1。利用P-M波浪譜仿真Υ=1和Υ=0.5時，不同風速下的σh；并根據該σh和式(6)計算對應的ecutoff，得到風速與ecutoff的對應關系。圖1(a)給出了在理想開闊成熟海面(Υ=1)和理想散射條件假設下的風速與ecutoff之間的數據關系；圖1(b)給出了在理想近岸海面(Υ=0.5)和理想散射條件假設下的風速與ecutoff之間的數據關系?？梢钥闯?，風速與ecutoff有單調遞減的數學關系。

2 方 法

2.1 小波分析提取截止高度角

小波分析是一種通過對時間和頻率進行局域變換，實現信號多尺度細化分析的方法，可以對序列的細節信號進行放大，起到“放大器”的作用[18-19]。設ψ(t)是基本小波函數，對于輸入SNR序列δS(sin(e))，其小波變換為

(8)

式中，a和b分別是尺度縮放因子和平移因子。要獲得小波分析結果，需要計算每個移位參數b的卷積，并對a重復此過程。MATLAB Wavelet Toolbox中的函數——centfrq和scal2frq，給出了a、b參數與頻率的映射關系。得到頻率后，根據h=(λf)/2(式(1))，將頻率轉換為RH。圖2給出了HKQT站2018年DOY 258，PRN 1衛星的SNR序列、對應的小波分析譜圖及根據譜圖能量最大(峰值)選取的RH值，細節詳見文獻[18—19]。由于SNR噪聲總是表現為高頻噪聲，所以有效RH值和噪聲非常容易區分[18-19]。例如，圖2中，當RH大于4時，認為是能反映海面變化的有效RH值，對應的SNR序列振蕩明顯，相干能量顯著；當RH小于4時，認為是噪聲信息，對應的SNR序列無明顯振蕩，非相干能量顯著。即RH小于4的起始高度角，可以判定為多徑振蕩消失的截止高度角。文獻[11]中，判定截止高度角的方法是基于最小二乘曲線擬合和特定閾值的；該方法對數據質量要求高且閾值需要根據經驗確定，不同站點的閾值可能并不相同，方法較為復雜。而本文的小波分析判定方法，簡單易行且不需要基于經驗參數。

圖1 風速與ecutoff在不同相關長度下的對應關系Fig.1 Relationship between ecutoff and wind speed with different T

圖2 HKQT站2018年DOY 258，PRN 1衛星的SNR序列、小波譜圖及峰值RH序列Fig.2 SNR sequence, wavelet spectrum and peak RH sequence of PRN 1 satellite on DOY 258 of 2018 for HKQT station

2.2 截止高度角基準統一

風速引起海面粗糙度變化，粗糙度引起截止高度角變化。理論上，風速越大，粗糙度越大，截止高度角越??；相反地，風速越小，粗糙度越小，截止高度角越大。這里用來描述截止高度角的越小或越大是一種定性描述，需要轉換到定量描述。因此，本文構建了一個截止高度角變化量δecutoff，用以定量描述截止高度角的大??；δecutoff=ecutoff,wind_speed=1-ecutoff(其中，ecutoff,wind_speed=1為基準)，它描述的是以平靜海面(風速為1 m/s)時的截止高度角為基準，ecutoff相對于基準的變化。圖3給出了依據圖1獲得的δecutoff與風速之間的關系。圖3(a)給出了在理想開闊成熟海面(Υ=1)和理想散射條件假設下的風速與δecutoff之間的數學關系；圖3(b)給出了在理想近岸海面(Υ=0.5)和理想散射條件假設下的風速與δecutoff之間的數學關系。圖3中，δecutoff與風速呈現指數函數的數學關系。

圖3 風速與δecutoff在不同相關長度下的對應關系Fig.3 Relationship between δecutoff and wind speed with different T

在仿真試驗中，反射面粗糙度變化是影響截止高度角的唯一因素。然而，在實際觀測中，截止高度角與許多因素有關，其中包括粗糙度、反射率因子等反射面特性，也包括信號解調方式、信號頻率、信號帶寬等信號特性[24]。為了提高反演精度，需要保證粗糙度是影響截止高度角變化的唯一因素。由于GPS衛星的周期近似為1 d(11 h 58 min)，因此，同一衛星每天同一軌跡被接收機捕獲的同一信號，理論上具有相同的反射區及相同的解調方式、信號頻率、信號帶寬；為了實現唯一變量影響，需要計算同一衛星同一信號同一軌跡的截止高度角的變化。同時，由于風速為1 m/s時的截止高度角在實際中并不容易獲得，因此，選取同一衛星同一信號同一軌跡的截止高度角的均值作為基準，進而計算δecutoff。將第i個GPS衛星、第j個GPS信號，在第n天、第k軌跡獲得的截止高度角記為ecutoffi,j,n,k，則定義該截止高度角對應第i個GPS衛星、第j個GPS信號、第k軌跡的δecutoffi,j,n,k為

(9)

圖4 GNSS-IR風速反演Fig.4 GNSS-IR wind speed retrieval

3 算 例

3.1 站點及數據

站點HKQT(114.21°E，22.29°N)位于香港鰂魚涌，隸屬香港衛星定位參考站網(https:∥www.geodetic.gov.hk/)，站點環境如圖5(a)所示。接收機類型為TRIMBLE NETR5。該站可以接收多模多頻GNSS信號的1 Hz采樣的觀測數據?？紤]到反射區重返的問題，僅使用GPS數據進行風速反演。該站接收到的GPS信號有L1C/A、L2P、L2C及L5C；由于GPS L2P為加密信號，在信號解調過程中，容易引起一些GNSS-IR反演偏差[24-25]，故不在本算例中使用。根據文獻[26]提供的基于Google Earth的反射區繪制工具，繪制反射區(圖5(b))，根據水域范圍，確定有效海域方位角[-60°,105°]；反射區相關原理參見文獻[2,5,20]?？紤]到平靜水面條件下，截止高度角為25°～35°，因此，用于小波分析的數據弧段，選擇高度角區域[5°,40°]。距離該站約1 km處，有一氣象站點，編號為45 007，提供每小時一次的風速、降水等氣象參數。同時，距離HKQT站點2 m處有一驗潮站Quarry Bay可提供實測的潮位數據。為了驗證GNSS-IR技術能否探測強風風速，算例時間選擇了“天鴿”臺風登港時間(2018年9月16日)前后日期2018年DOY 249—DOY 270，及“山竹”臺風登港時間(2017年8月23日)前后日期2017年DOY 224—DOY 243。

3.2 風速反演結果

圖6給出了2018年DOY 247—DOY 270間(臺風“山竹”前后時間)的截止高度角變化參數與實測風速序列。圖6(c)中，GPS L5C信號的截止高度角變化δecutoff與實測風速對應關系良好，相關系數達到0.85，高度相關。而圖6(a)中，GPS L1C/A信號的截止高度角變化δecutoff與實測風速對應關系很差，相關系數僅為0.14，微弱相關；圖6(b)中，GPS L2C信號的截止高度角變化δecutoff與實測風速對應關系一般，相關系數為0.40，低度相關。許多研究表明[24,27]，GPS L5C是GPS中信噪比質量最好的信號。而在該算例中，GPS L5C是δecutoff與風速相關性最好的信號。臺風“山竹”在2018年9月7日于太平洋西北方位生成，2018年9月16日(DOY 258—DOY 259)臺風向我國廣東沿海逼近，2018年9月16日凌晨，香港處于臺風風圈以內，一直持續到2018年9月16日中午。圖6(c)中，在該時間內，實測風速劇增至32 m/s，截止高度角變化量劇增至18°；包括在風速攀升階段，截止高度角也能很好地刻畫風速在低風速至高風速下的變化情況。值得提出的是，在小于5 m/s的情況下，記錄的風場信息，風向一直不穩定；這種不持續風向的風場條件會導致GNSS-IR風速探測的不準確。

圖5 HKQT站Fig.5 HKQT station

圖7給出了2017年DOY 224—DOY 243(臺風“天鴿”前后時間)的截止高度角變化參數與實測風速序列。圖7(c)中，GPS L5C信號的截止高度角變化δecutoff與實測風速對應關系良好，相關系數達到0.70，高度相關。而圖7(a)中，GPS L1C/A信號的截止高度角變化δecutoff與實測風速對應關系很差，相關系數僅為0.05，不相關；圖7(b)中，GPS L2C信號的截止高度角變化δecutoff與實測風速對應關系一般，相關系數為0.42，低度相關。在該算例中，GPS L5C是δecutoff與風速相關性最好的信號。臺風“天鴿”在2017年8月20日于太平洋西北方位生成，受其影響，香港局部區域于2017年8月23日(DOY 234—DOY 235)遭受10級大風。圖7(c)中，在該時間內，實測風速劇增至23 m/s，截止高度角變化量劇增至16°。除去臺風期間，香港2019年8月26日(DOY 238)持續有3～4級東風，在該時間內，實測風速增至15 m/s，截止高度角變化量劇增至10°。同時，在風速攀升階段，截止高度角也能很好地刻畫風速在低風速至高風速下的變化情況。在小于5 m/s的情況下，記錄的風場信息，風向一直不穩定；這種不持續風向的風場條件會導致GNSS-IR風速探測的不準確。

圖6 HKQT站2018年DOY 250—DOY 270 δecutoff序列與實測風速序列Fig.6 δecutoff vs. measured wind speed during DOY 250 and DOY 270 of year 2018 for HKQT station

圖7 HKQT站2017年DOY 225—DOY 243 δecutoff序列與實測風速序列Fig.7 δecutoff vs. measured wind speed during DOY 225 and DOY 243 of year 2017 for HKQT station

為了更好地比較截止高度角變化δecutoff與風速的關系，圖8給出了GPS L5Cδecutoff與風速之間的對應關系。由于GPS L1C/A與GPS L2C信號δecutoff與風速之間的相關性較低，故未在圖8中分析。圖8表明，δecutoff與風速之間存在一種近似指數函數的變化；圖8中的風速s與δecutoff的擬合指數函數(擬合函數的RMSE為3.09 m/s)為

s=f(δecutoff)=3.70×e0.12×δecutoff

(10)

同時，2018年和2017年不同時間段內的函數關系相同。這說明，同一個站點，δecutoff與風速之間的數學對應關系是一定的。該數學關系與圖3(b)中T=4時根據仿真獲得的數學關系近似。因此，在對海面情況進行調查后，可以根據散射模型和波浪譜模型獲得δecutoff與風速之間的數學關系，再通過實測的δecutoff，反推獲得相關的風速值。

圖8 HKQT站δecutoff與風速Fig.8 δecutoff VS. measured wind speed at HKQT station

4 結論與討論

本文提出了一種基于小波分析的GNSS-IR風速反演方法，通過算例證實了δecutoff與風速存在特定的指數函數關系，證實了利用GNSS-IR技術進行風速探測的可能性，擴充了GNSS-IR的監測對象，拓寬了GNSS-IR的應用范圍。地基GNSS-IR風速反演和星載GNSS-R風速反演的方法不同，但是核心相同——風速會造成不同的海面粗糙度，使反射信號特性不同；從而可以利用反射信號的特性反推粗糙度，進而反推風速。為了保證風速反演的準確性，需要盡量保證風速是引起粗糙度的唯一因素，這就對站點的海面環境提出了要求。對于海面變化劇烈(如潮位振幅較大或者有船行波影響)的站點，本身海面的粗糙度就較大，理論上是無法進行風速反演的；對于某些粗糙度的變化主要為風導致(一般海面變化較小)的站點，GNSS-IR風速反演技術理論上可用。

值得注意的是，利用截止高度角參數、基于波浪譜模型和散射模型的仿真在文獻[17](研究粗糙度及有效波高反演)中均已提及。本文的貢獻主要在：

(1) 使用了更好的截止高度角參數提取方法——小波分析，來反演風速。

(2) 提出了截止高度角基準統一方法。文獻[17]認為不同衛星不同軌跡相同信號的截止高度角的基準相同；而經過算例分析，本文發現不同衛星不同軌跡相同信號的截止高度角的基準不同。文獻[17]認為基準相同，可能是由于文獻[17]中的試驗區域環境極好(站點很高、海域開闊、四面環海)，不同軌跡的反射區環境大致相同導致。而本文算例中的站點環境更為常見，所用的截止高度角基準統一方法適應性更強。

(3) 證實了GNSS-IR技術可以監測風速。

本文只是GNSS-IR探測風速的初步探索，想要推進該技術的實際應用，還需要對下列問題進行研究。

(1)δecutoff與風速之間的數學關系。根據理論，式(10)所得到的數學關系，僅適用于附近海域(或相似特性的海域)及GPS L5C信號，對于其他特性的海域及其他信號并不適用。數學關系的準確與否決定了能否根據δecutoff獲得準確的風速參數。該數學關系與海面特性(主要指海面粗糙度隨風速的變化情況)及信號特性(包括信號調制和解調方式、信號頻率等)有關。數學關系的建立可以從兩方面入手。①使用更好的波浪譜模型，得到更好的仿真數學關系。目前效果比較好的波浪譜模型大都基于成熟廣闊的海域。本文所用到的沿岸波浪譜是最簡單的波浪譜，并不準確。實際上，沿岸波浪受海洋深度、海洋特性、沿岸海岸線、風速風向等影響，相關波浪譜的建立十分復雜。波浪譜模型的準確建立，對于建立δecutoff與風速之間的數學關系非常重要。②根據大量實測數據，建立數學關系。由于利用仿真建立的數學關系是基于理想情況下，可能與實際情況不符。在大量算例的支撐下，可以根據大量實測δecutoff和實測風速數據，利用擬合方法或神經網絡方法，建立二者之間的數學關系；從而根據任一站點的實測δecutoff推算風速。

(2) 風速反演閾值的研究。風速反演必然存在上限，即，對于大于某一特定上限的風速，GNSS-IR方法無法探測。在上限風速情況下，引起的海面粗糙度已經使截止高度角達到最低(本文為5°)。上限風速與信號特性、信號解調方式、海面特性等有關。對于本算例用到的站點情況，上限風速在25 m/s～30 m/s。

(3) 信號研究。本文的算例表明，GPS L5C可用，而GPS L1C/A、GPS L2C信號并不可用。但是在某些環境較好的情況下(站點較高、海域開闊)，GPS L1C/A、GPS L2C信號也可能能夠反映風速變化。需要相關研究，以獲得不同信號在不同環境下的截止高度角參數與風速之間的數學關系。

另外，本文提出的相關理論和方法除了進行風速反演外，還隱藏著另一個重要的潛在應用——推算波浪特性。波浪特性對于海岸工程建設十分重要，波浪資料的獲得目前主要依賴于浮標；而浮標站點數量有限，監測范圍小。根據本文的理論及試驗，在有實測風速的情況下，根據SNR的δecutoff值，可以繪制δecutoff與風速圖(如圖8所示)；而根據不同的風速波浪譜，可以得到不同的仿真δecutoff與風速圖(如圖1及圖3所示)；將實際δecutoff與風速圖與仿真δecutoff與風速圖進行對比，根據符合程度，可以推算該地區在風速影響下的波浪特性。例如本文根據實測圖8及仿真圖3，得到HKQT站點附近海域的海浪表現近似為T=4，Υ=0.5時P-M波浪譜模型(式(7))。今后，本團隊將繼續研究GNSS-IR波浪特性反演的相關理論和方法，并挖掘GNSS-IR其他潛在應用。