吳迪,張波,李博聞
(北京航空航天大學 電子信息工程學院,北京100083)
自20世紀90年代以來,利用全球導航衛(wèi)星系統(tǒng)(Global Navigation Satellite System,GNSS)反射信號進行遙感探測的技術成為導航領域的研究熱點之一[1]。該技術廣泛應用于海面測風、海面測高、海冰探測、海鹽探測、海面溢油、土壤濕度等諸多領域,具有重要的研究意義和廣闊的應用前景[2-5]。
隨著遙感技術的發(fā)展,關于GNSS反射信號模擬仿真的研究引起了越來越多的關注。2000年,Zavorotny和Voronovich[6]基于雙基雷達方程和基爾霍夫幾何光學(Kirchhoff Approximation-Geometric Optics,KA-GO)近似法,提出GNSS散射信號的時延-多普勒相關功率模型,即經典的Z-V模型,該模型描述了全球定位系統(tǒng)(Global Positioning System,GPS)的散射信號功率隨幾何及環(huán)境參數的變化。Schiavulli等[7]為改進基于該模型的仿真技術,提出了一種基于面元的散射模型,生成隨機三維海面,利用基爾霍夫物理光學(Kirchhoff Approximation-Physical Optics,KA-PO)近似法計算GNSS散射信號相關功率和二維時延-多普勒相關功率分布圖像(Delay-Doppler Maps,DDM),研究GNSS反射信號的極化特性。伊朗學者Ghavidel等基于Z-V模型模擬反射信號,利用數值模擬計算電磁偏離,分析了其與海浪譜、入射/散射角的關系[8],同時引入降雨、涌浪和表面洋流對電磁偏離的影響[9]。為了更好地描述GNSS散射信號相關功率的統(tǒng)計特性,Park等[10]為歐洲PAU/PARIS端對端性能模擬器(PAU/PARIS End-to-end Performance Simulator,P2EPS)系統(tǒng)增加了大氣誤差和斑點噪聲的影響,并考慮了地形地貌等地球物理參數對散射功率的影響,完善了其探測陸地參數的性能[11-12]。李博聞等[13]提出反向建模思路,將反射信號等效為n條具有不同時延、多普勒以及幅度的信號疊加。
現有仿真研究大多基于Z-V模型,采用風驅海浪譜分析風驅海面對GNSS反射信號的影響[14-15],而未考慮涌浪、降雨等復雜環(huán)境的影響。其中,降雨可以改變大氣邊界層的結構,撞擊海面產生毛細重力波環(huán),導致表面粗糙度的改變,這種改變將導致電磁場偏離增大[16]。在測高應用中,降雨可以引入分米級的電磁場偏離誤差[9]。涌浪的存在也會使海面粗糙度發(fā)生變化,造成后向散射截面增大,且在L波段以及入射角較小的情況下產生的影響更大[17]。這些復雜環(huán)境影響導致的海面粗糙度變化在反射信號風速反演等應用中都會引入誤差,對真實的復雜環(huán)境下的GNSS反射信號模擬及遙感技術的研究至關重要。
本文針對涌浪、降雨2種海面影響因子,提出在復雜環(huán)境下的GNSS海面反射信號的建模方法。首先,對Elfouhaliy海浪譜、涌浪譜、降雨譜進行仿真,從海浪譜的角度分別論述了涌浪、降雨對海面統(tǒng)計特征及GNSS海面反射信號的影響。其次,利用該影響對風驅海浪譜進行修正,建立復雜環(huán)境下的GNSS海面反射信號模型,并基于該模型進行仿真得到DDM。然后,將仿真的DDM 與英國技術演示衛(wèi)星(United Kingdom TechDemoSat-1,UK TDS-1)實測數據的處理結果進行比對驗證。最后,對所建立的反射信號模型、處理的結果和不足之處進行了總結。
海面受海浪的影響,上下起伏,可視為隨機過程,一般采用海浪譜描述其統(tǒng)計特征。其中Elfouhaily海浪譜[18]考慮了風區(qū)對能量的影響,廣泛應用于海面風場反演。其可表示為[18]
式中:ME(k)表示海浪譜各向同性部分;fE(k,φ)為對應的方位函數;k為波數;φ為波向。
在不同風速下,仿真得到的Elfouhaily海浪譜如圖1所示。圖中:Kc為截止波數,其將波譜分為大小尺度粗糙度2個部分,小于Kc的部分為大尺度粗糙度,大于Kc的部分為小尺度粗糙度。Brown[19]提出大尺度粗糙度的截止波數可以表示為
圖1 不同風速下的Elfouhaily海浪譜Fig.1 Elfouhaily wave spectrum at different wind speeds
式中:θi為入射角;λs為信號波長。
對于GPS L1信號,截止波長約為0.6m,可計算此時Kc約為10 rad/m。由圖1可知,隨著風速增大,風浪的總能量增大,同時譜峰向低頻、大尺度粗糙度方向推移,使得前向散射系數減小,漫反射增強。
涌浪不同于風浪,是由其他風場傳播而來的海浪,持續(xù)運動且不再受當地風影響。據海洋觀測數據顯示,僅風浪或僅涌浪出現的情況極少,大多數情況以混合浪形式出現[20]。因此,考慮涌浪的影響才能模擬更真實的GNSS反射信號。涌浪可被建模為一個窄帶高斯過程,即[21]
式中:〈h2〉為涌浪的高度變化;σkx和σky為譜的標準偏差;kx、ky分別為波數k在x、y方向的分量;kxm、kym分別為涌浪在x、y方向的譜峰波數。為了突出涌浪的影響,一般取〈h2〉為4 m2,波長Λm為400 m[9],且滿足:
對該涌浪譜進行仿真,結果如圖2所示。從圖2中可以看出,涌浪驅動的海浪譜位于截至波數Kc左側,主要形成對GNSS海面反射信號影響較大的大尺度粗糙海面。且在低風速時,涌浪的影響更為顯著,該情況下進行風速反演,涌浪的存在將帶來巨大的反演誤差,從而降低風速反演的精度。
圖2 風驅、涌浪驅動的Elfouhaily海浪譜Fig.2 Wave spectrum driven by Elfouhaily wind and swell
降雨時,飛濺的雨滴進入海水產生渦旋使表面粗糙度發(fā)生變化,變化大小與降雨率、雨滴大小分布有關。Bliven等[22]通過實驗模擬降雨擾動水面的場景,并對散射計數據進行了分析,實驗結果證明可將降雨譜建模為對數高斯環(huán)形波譜。其表示為
式中:Cg(k)為群速度;Sfp為冪律譜模型;g為重力加速度;τw為水面張力;ρ為水的密度;Δf為帶寬;fp為峰值頻率。
利用該降雨譜進行仿真,結果如圖3所示。降雨驅動的海浪譜主要對小尺度粗糙海面產生影響,對GNSS反射信號的影響較小,且對風速不敏感。僅在低風速強降雨情況下,降雨才會對GNSS反射信號產生一定影響。
圖3 風驅、降雨驅動的海浪譜Fig.3 Wave spectrum driven by wind and rain
接收機端接收到的GNSS反射信號可以看作各個散射單元散射信號的加權和[23]。DDM 作為GNSS反射信號的基礎觀測量,可用基于雙基雷達方程的GNSS反射信號時延-多普勒二維相關功率模型來表示[24]:
式中:Pt為發(fā)射機功率;Gt為發(fā)射天線增益;λ為GNSS信號波長;Tc為相干積分時間;Gr_x,y為散射單元Sx,y方向上的接收天線增益;σp,q_x,y為p極化入射波對應q極化反射波時散射單元Sx,y的散射系數;Rt_x,y、Rr_x,y分別為GNSS衛(wèi)星和接收天線到散射單元Sx,y的距離;Λ和S分別為偽隨機碼自相關函數和多普勒濾波函數;τ和f分別為鏡面反射點處的時延、多普勒;fx,y為散射單元Sx,y處的反射信號分量的多普勒;N(τ,f)為噪聲項。
當多普勒頻率固定時,利用式(6)可推導GNSS反射信號一維時延相關功率波形(Delay Waveform,DW)的表達式:
基于風驅、涌浪驅及降雨驅海浪譜,建立受涌浪、降雨擾動的雙基電磁散射模型,從而計算GNSS反射信號的時延-多普勒二維相關功率,完成GNSS反射信號的建模。其研究框圖如圖4所示。具體建模步驟如下:
1)建立散射參考坐標系(Scattering Reference Frame,SRF)[25],根據GNSS衛(wèi)星和接收機狀態(tài)(位置、速度)計算鏡面反射點。
2)在10000平方公里的正方形區(qū)域內,將表面劃分為1 km×1 km的網格,基于Elfouhaily海浪譜仿真每個網格所對應的二維隨機海面,如圖5所示(風速為5m/s,風向φw=45°,色標為每個像素點的海面高度,單位為m)。該尺寸足以囊括最高風速條件下與涌浪相關的長波[9]。然后進行空間域到時延-多普勒域的映射,計算每個網格的時延、多普勒等參數。
圖4 復雜環(huán)境下GNSS反射信號建模研究框圖Fig.4 Block diagram of GNSS reflected signal modeling in complex environment
圖5 基于Elfouhaily海浪譜的二維隨機海面Fig.5 Two-dimensional random sea surface based on Elfouhaily wave spectrum
3)結合風驅、涌浪驅、降雨驅海浪譜,基于雙尺度模型建立受涌浪、降雨擾動的電磁散射模型,完成散射系數的計算。
4)加入噪聲,利用式(6)計算DDM。
5)進行非相干累加。
相比基于風驅海浪譜的建模,為得到適合于涌浪、降雨這些復雜環(huán)境因素影響下的反射信號模型,其中關鍵的步驟在于受擾動的雙基電磁散射模型的建立。
當GNSS信號入射到海洋表面時,發(fā)生散射。由于海表粗糙度較大,反射信號的相干分量減少,可被認為是漫反射。散射強度的大小通常用散射系數描述。目前,常用的電磁散射模型有基于KAGO模型、微擾動法(Small Perturbation Method,SPM)模型、雙尺度(Two Scale Model,TSM)模型、小斜率近似(Small Slope Approximation,SSA)模型[1]。
根據1.2節(jié)及1.3節(jié)的分析,涌浪與風速類似,主要影響海浪譜的低波數部分,即主要對大尺度粗糙海面產生影響。而降雨則主要影響小尺度粗糙海面,對GNSS信號的前向散射影響較小。因此可以采用雙尺度模型,其中大尺度即風速和涌浪的散射系數用KA-GO近似法計算,小尺度即降雨的散射系數用SPM法計算。綜上所述,引入涌浪、降雨影響后的散射系數為
式 中:σKA-GO,wind、σKA-GO,swell、σSPM,rain分別為風速 和涌浪影響的大尺度粗糙度海面的散射系數及降雨影響的小尺度粗糙度海面的散射系數。其中,σKA-GO,wind、σKA-GO,swell用KA-GO近似法計算:
其中:積分上限Kc為截止波數,可由式(2)計算得到。
GNSS直射信號經海面反射后,菲涅爾反射系數為
式中:下標r、l、v、h分別表示右旋圓極化、左旋圓極化、垂直線極化、水平線極化;ε為海面復介電常數;θ為衛(wèi)星高度角。
σSPM,rain用SPM法計算,
式中:θs為散射角;Srain(kx,ky)為降雨驅海浪譜。
利用上述雙尺度模型計算受涌浪、降雨影響下的散射系數的流程如圖6所示。圖中坐標系將由SRF轉換為地心地固坐標系(Earth-Centered,Earth-Fixed,ECEF)。
GNSS反射信號主要被2種噪聲干擾:熱噪聲和斑點噪聲[24]。加噪后的二維時延-多普勒相關功率可被建模為
式中:Nsk(τ,f)為乘性斑點噪聲;Nt(τ,f)為加性熱噪聲,服從(0,σ2t)的高斯分布;σ2t為噪聲功率[26]。
圖6 散射系數的計算流程圖Fig.6 Flowchart of scattering coefficient calculation
根據概率理論,2個標準高斯分布的平方和服從自由度為2的中心卡方分布。因此,熱噪聲可建模為
式中:Nt_I、Nt_Q分別為同相、正交支路的熱噪聲;χ2(2)為自由度為2的中心卡方分布。
斑點噪聲是乘性噪聲,一般假設其服從均值為1,方差為σ2sk的指數分布。
根據式(15)對二維時延-多普勒相關功率進行計算,并對得到的DDM 進行非相干累加,以降低熱噪聲和斑點噪聲,從而改善波形。
利用已知的空間幾何、散射場景、接收場景等參數,根據式(6)計算GNSS反射信號的二維時延-多普勒相關功率。
設置仿真場景為星載條件下,GNSS衛(wèi)星的位置和速度分別為[-11 178 791.991 294,-13 160 191.204 988,20 341 528.127 540]m、[2523.258023,-361.592839,1163.748104]m/s,低軌衛(wèi)星的位置和速度分別為[-4069896.7033860330,-3583236.963735 0840,4 527 639.271 758 164 0]m、[-4 738.074 234 206 3,-1 796.252 568 996 4,-5 654.995 201 365 7]m/s,海面風速為5m/s,波長為400m,接收天線增益為12 dB,天線3 dB波束寬度為20°,當降雨率超過200mm/d時為特大暴雨[27],因此為突出降雨的影響,將降雨率設為100mm/h,同樣為突出涌浪的影響,將涌浪高度方差設為4m2[9]。
該系統(tǒng)主要采用JSP語言編寫,開發(fā)環(huán)境為MyEclips,服務器采用Apache Tomcat,數據庫采用MSSQL Server2005。
仿真得到復雜環(huán)境下的DDM、DW 分別如圖7和圖8所示。從圖7中可以看出,仿真得到的DDM整體均呈現典型的“馬蹄形”分布,但在不同的影響因素下,DDM的峰值功率發(fā)生顯著變化,如在涌浪影響下,峰值功率由圖7(a)中的1.2×10-27下降至圖7(b)中的8×10-28。比較圖7(a)、(b)及圖8可知,加入涌浪影響后,海面粗糙度發(fā)生較大變化,海面GNSS信號的漫反射增強,而鏡面反射的GNSS信號減弱,導致DDM峰值功率明顯減小,DW 拖尾延長,幅度明顯下降。比較圖7(a)、(c)及圖8可以發(fā)現,相比于涌浪,降雨情況下DDM 和DW 的形狀、峰值功率均無明顯變化,雖然降雨也會改變海面粗糙度,但其對GNSS信號的前向散射功率影響非常小,與理論預期相符,這也是可以利用GNSS反射信號透過臺風墻對內部風速進行探測的原因。
為了驗證建立的復雜環(huán)境影響下的GNSS反射信號模型的準確性,利用UK TDS-1于2014-09-01 20:24:00采集的真實海面GNSS反射信號數據進行處理。該組數據的海況如下:海面風速為5m/s,涌浪高度方差為4m2,降雨率為0 mm/h。其他數據參數及仿真參數如表1所示。
按照上述參數對GNSS反射信號進行仿真,計算引入涌浪后的DDM、DW,并與UK TDS-1數據處理得到的DDM、DW 進行對比,結果分別如圖9和圖10所示。其中,GNSS衛(wèi)星的位置和速度分別為[-21 842 415.002 320 271,-15180714.434890306,-970098.63007445715]m、[258.9350681325609,-174.18019867283397,-3 155.347 262 329 083 2]m/s,低軌衛(wèi)星的位置和速度分別為[-6 806 318.464 608 931,-1262592.5946520383,-1 103923.093 804 532 4]m、[847.567 995 738 401 81,1 803.738 038 472 3278,-7 368.059 680 401 838]m/s。
圖9顯示,引入涌浪影響后的模型仿真得到的DDM與UK TDS-1數據處理的DDM波形的形狀及峰值功率之間有很好的一致性。從圖10可以看出,在延遲為(-1~15)C/A碼片范圍內,引入涌浪影響因子后的DW 與UK TDS-1數據處理的DW 波形更為接近,二者的相關系數達到0.92。而在相同場景下,未修正模型的DW與UK TDS-1數據處理結果的相關系數為0.88,明顯低于修正后的模型。因此,引入涌浪影響因子后的GNSS反射信號模型不僅在物理層面上更接近真實情況,并且模型的精度有一定提高,優(yōu)于未修正模型。所以引入涌浪影響后的GNSS海洋反射信號模型能夠更真實的模擬GNSS反射信號。證明該復雜環(huán)境下的GNSS反射信號模型具有一定的可行性,可以有效的用于GNSS海面反射信號的仿真。
圖7 仿真的復雜環(huán)境下的GNSS反射信號二維時延-多普勒相關功率分布圖像Fig.7 Simulated two-dimensional DDM of GNSS reflected signal in complex environment
圖8 仿真的復雜環(huán)境下的GNSS反射信號一維時延相關功率波形Fig.8 Simulated one-dimensional DW of GNSS reflected signal in complex environment
表1 UK TDS-1真實數據參數Table 1 UK TDS-1 real data parameters
圖9 模型仿真及UK TDS-1數據處理所得二維時延-多普勒相關功率分布圖像Fig.9 Two-dimensional DDM obtained from model simulation and UK TDS-1 data processing
圖10 模型仿真及UK TDS-1數據處理所得一維時延相關功率波形Fig.10 One-dimensional DW obtained from model simulation and UK TDS-1 data processing
本文在風驅海面下GNSS反射信號模型的基礎之上,從海浪譜角度分析涌浪、降雨對GNSS反射信號的影響,提出利用雙尺度模型計算引入涌浪、降雨影響因子后散射系數的方法,建立涌浪、降雨影響下的GNSS海面反射信號模型。并根據該建模流程進行仿真,將得到的DDM 和DW 與UK TDS-1實測數據進行比對,結果表明:
1)提出的涌浪、降雨影響下的GNSS反射信號模型簡單,算法清晰,易于程序實現,并加入噪聲影響,提高仿真精度。
2)涌浪的存在使得海面粗糙度發(fā)生變化,形成對GNSS反射信號影響較大的大尺度粗糙海面。而降雨對反射信號的影響較小,與理論分析結果相符。因此在風速反演中加入涌浪的影響將會有效提高反演精度。
3)引入涌浪影響因子后仿真得到的DDM、DW 與實測數據結果相關系數達到0.92,優(yōu)于未修正模型與實測數據的對比結果。驗證了該建模方法的可行性和有效性,為模擬真實復雜環(huán)境下的反射信號和GNSS反射信號星載探測應用的研究奠定基礎。
此外,為使本文提出的模型能用于模擬更為真實的反射信號,后續(xù)考慮引入泡沫、洋流等對GNSS反射信號的影響以及相應軟件模擬器的實現。