基于優(yōu)化濾波的水平非均勻蒸發(fā)波導反演算法*

董道廣，芮國勝，康　健，田文飚

（海軍航空工程學院信號與信息處理山東省重點實驗室，山東　煙臺　264000）

基于優(yōu)化濾波的水平非均勻蒸發(fā)波導反演算法*

董道廣，芮國勝，康健，田文飚

（海軍航空工程學院信號與信息處理山東省重點實驗室，山東煙臺264000）

針對水平非均勻蒸發(fā)波導反演中大氣折射率剖面水平分布規(guī)律不確知和大范圍海域氣象參數(shù)觀測難的問題，提出一種基于優(yōu)化濾波的水平非均勻蒸發(fā)波導折射率剖面參數(shù)反演算法，該算法采用粒子群優(yōu)化估計折射率剖面參數(shù)初始狀態(tài)和描述剖面參數(shù)水平分布規(guī)律的狀態(tài)轉(zhuǎn)移函數(shù)，結(jié)合GPS海面散射信號波導傳播的DMFT步進計算正演模型和實測GPS散射信號功率建立系統(tǒng)方程和觀測方程，并采用粒子濾波迭代估計水平步進上的非均勻折射率剖面參數(shù)。該算法的精估計環(huán)節(jié)有效修正了粒子群估計誤差在后續(xù)反演過程中導致的誤差積累，保證了反演精度，并使得算法具備距離分段條件下逐段反演的能力，為遠距離海域蒸發(fā)波導折射率剖面參數(shù)的反演提供了可行的方法。

蒸發(fā)波導，粒子群優(yōu)化，粒子濾波，反演

0　引言

蒸發(fā)波導是出現(xiàn)概率最高的海上大氣波導，過去人們對其反演問題的研究更多停留在均勻蒸發(fā)波導，事實上，蒸發(fā)波導因海上溫度、濕度、氣壓及風速等氣象因子在大尺度區(qū)域分布上的差異性而呈現(xiàn)出空間分布的非均勻性。

Karimian［1-3］采用現(xiàn)代濾波技術對時變大氣波導進行了卓有成效的反演。由于濾波反演算法的收斂性依賴于系統(tǒng)方程和初始狀態(tài)，上述困難嚴重影響到了反演的精度。為了解決上述兩方面問題，本文提出一種基于優(yōu)化濾波的GPS海面散射信號反演水平非均勻蒸發(fā)波導算法。本文算法解決的是空間上的反演問題。該算法將前人提出的離散混合傅立葉步進計算方法（DMFT）作為GPS海面散射信號蒸發(fā)波導傳播的正演工具，通過粒子群優(yōu)化估計濾波反演算法的初始狀態(tài)和系統(tǒng)方程，依據(jù)估計值建立系統(tǒng)方程和觀測方程，并將正演GPS海面散射信號波導傳播的DMFT步進計算和粒子濾波的迭代狀態(tài)估計同步處理。對于遠距離、大范圍的反演，該算法采取分段優(yōu)化的方法來修正估計誤差引起的反演誤差積累。

1　反演建模

1.1折射率剖面參數(shù)化建模

折射率剖面參數(shù)既是反演輸入量，也是反演輸出量，故必須對折射率剖面作參數(shù)化建模。經(jīng)典的蒸發(fā)波導剖面模型有MGB［4-5］模型和P-J模型等，本文采用P-J模型，該模型適用于熱中性氣象條件，是 Paulus［6］在 Jeske［7］的研究基礎上基于Monin-Obukhov相似理論提出的，其數(shù)學解析式為

其中，z為高度，δ為波導高度，z0為空氣動力學粗糙度因子，M0為海面修正折射率。圖1給出了該剖面模型的示意，根據(jù)該模型可確定兩個參數(shù)，分別為波導高度δ和波導強度ΔM。

圖1　P-J模型參數(shù)化建模示意

在海上復雜的氣象參數(shù)環(huán)境下，實際折射率剖面不如模型剖面平滑，所以實現(xiàn)對實際蒸發(fā)波導折射率剖面的高精度擬合要求模型具備靈活的參數(shù)調(diào)整能力，以使模型能夠遍歷和構(gòu)成的二維可行域。本文以P-J模型為依據(jù)，對折射率剖面作二維參數(shù)化建模，現(xiàn)導出其數(shù)學解析式，給定參數(shù)δ和ΔM，設z=δ以下高度范圍內(nèi)的修正大氣折射率為含輔助參變量的函數(shù)，如式（2）

z=δ處的修正大氣折射率為滿足，

解式（3）可得

z≥δ范圍內(nèi)的大氣修正折射率仍如式式（1）所示，為使式（1）和式（2）在z=δ處計算的數(shù)值相等，該范圍內(nèi)修正折射率為

其中

式（2）～式（6）即為反演所用P-J二參數(shù)剖面模型。

1.2觀測方程及系統(tǒng)方程建模

觀測方程描述折射率剖面參數(shù)與正演GPS海面散射信號的函數(shù)關系，系統(tǒng)方程描述折射率剖面參數(shù)在水平方向上按距離步進的更新。為方便分析計，將折射率剖面參數(shù)記為，根據(jù)拋物方程的DMFT步進算法［8-9］，正演GPS海面散射信號步進傳播損耗的自變量即當前步進處的垂直修正折射率剖面，據(jù)此將該步進的觀測方程簡記為式（7）

式中，L代表傳播損耗，下標n為當前步進的距離索引，上標j為參數(shù)δ和ΔM在當前步進剖面可行域內(nèi)眾多參數(shù)中的序號，v是加性白噪聲，方差為R。根據(jù)式（7）即可獲得GPS海面散射信號的正演功率序列

其中，下標sim表示正演，其后的數(shù)字標號表示不同步進的距離索引。GPS海面直射波本質(zhì)上是球面波，由于衛(wèi)星到海面距離較大，可將直射波近似為平面波［10-11］，GPS散射功率接收機接收到的是一個功率-延遲序列，該序列是由海面散射信號在水平方向上不同距離處的接收延遲引起的，可記為

其中，下標obs表示實測，其后的數(shù)字標號意義同式（9）。反演的本質(zhì)即通過式（8）和式（9）的擬合，求取擬合度最高的折射率剖面參數(shù)δopt和ΔMopt。

其中，F(xiàn)是當前距離步進到下一距離步進的剖面參數(shù)狀態(tài)轉(zhuǎn)移函數(shù)，是為了方便分析而增加的一個加性白噪聲。

2　反演算法

2.1初始狀態(tài)及系統(tǒng)方程估計

初始狀態(tài)即，是GPS海面散射正演起算點的折射率剖面參數(shù)。系統(tǒng)方程的建立依賴于狀態(tài)轉(zhuǎn)移函數(shù)F，假設F有確定的解析形式且海上折射率剖面在水平方向上的變化是連續(xù)而平滑的，則可將作級數(shù)展開，并取二階近似，如式（11）

綜上，初始狀態(tài)及系統(tǒng)方程的估計即是對參數(shù)δ1，ΔM1，α1，α2，β1和β2的估計，考慮到國內(nèi)已有粒子群優(yōu)化反演蒸發(fā)波導較成功的應用［12］，本文采用粒子群優(yōu)化給出初始狀態(tài)的估計流程：

1）結(jié)合式（7）～式（9），建立各距離步進上度量GPS散射信號正演和實測功率擬合度的目標函數(shù)

5）設fmin為停止迭代更新的閾值，反復執(zhí)行上述4步，直到達到迭代次數(shù)上限或滿足fgbest小于fmin，停止迭代，獲得初始狀態(tài)估計值；

α1，α2，β1和β2估計流程類似于上述流程，所不同的是其目標函數(shù)是經(jīng)多步DMFT步進計算獲得的正演功率序列構(gòu)成的，如式（15）

2.2水平非均勻折射率剖面參數(shù)反演

水平非均勻折射率剖面參數(shù)的反演由粗估計和精估計兩部分組成，粗估計基于2.1節(jié)估計結(jié)果建立觀測方程和系統(tǒng)方程，采用粒子濾波完成一段距離范圍內(nèi)的參數(shù)反演；精估計對粗估計最后一個距離步進上的δN，Δ，α1，α2，β1和β2的粒子群優(yōu)化，通過線性插值修正2.1節(jié)的估計誤差在后續(xù)粒子濾波過程中造成的誤差積累，并平滑反演結(jié)果?，F(xiàn)給出粗估計的步驟：

1）根據(jù)式（7）～式（8）建立觀測方程，根據(jù)式（10）～式（11）建立系統(tǒng)方程；

2）設初始狀態(tài)的概率密度函數(shù)為均勻分布類型，依此分布隨機產(chǎn)生S個狀態(tài)粒子，，…，，結(jié)合系統(tǒng)方程對每個粒子作時間更新，獲得下一步進上的先驗狀態(tài)粒子，，…；

3）由觀測方程計算每個先驗狀態(tài)粒子在下一步進上的正演功率，i∈［1，S］，并結(jié)合該步進的實測功率Pobs，2計算每個先驗狀態(tài)粒子的似然概率密度，如式（16）～式（17）

4）基于歸一化的似然概率密度對先驗狀態(tài)粒子重采樣，獲得后驗狀態(tài)粒子，對后驗狀態(tài)粒子取均值即得下一步進上的狀態(tài)粒子估計；

5）反復執(zhí)行上述4步，直至獲得0～R1距離范圍內(nèi)的剖面參數(shù)，…。

假設距離對應的步進索引為N，現(xiàn)在粗估計R1步驟的基礎上給出精估計步驟：

1）參照2.1節(jié)粒子群優(yōu)化步驟重新估計所在步進上的剖面參數(shù)狀態(tài)矢量，得；

2）設距離r對應的步進索引為n，對r～R1距離范圍內(nèi)各步進上的剖面參數(shù)作線性插值，獲得，…，，如式（18）

對于較大距離范圍內(nèi)的反演，可將距離區(qū)間分為數(shù)段，依次執(zhí)行上述反演步驟。

3　算法性能分析及仿真實驗

針對波導折射率剖面水平分布規(guī)律不確知和大范圍海域氣象參數(shù)觀測難的問題，基于優(yōu)化濾波的反演算法將折射率剖面參數(shù)的水平分布以現(xiàn)代濾波理論的系統(tǒng)方程形式表達，對其狀態(tài)轉(zhuǎn)移函數(shù)進行二次有理多項式近似，采用粒子群優(yōu)化估計GPS海面散射區(qū)域的初始折射率剖面參數(shù)和狀態(tài)轉(zhuǎn)移函數(shù)。由于水平方向上的折射率剖面參數(shù)是連續(xù)而平滑的，二階近似的誤差僅為狀態(tài)轉(zhuǎn)移函數(shù)對距離步進二階導的高階無窮小，可以忽略不計。綜上分析，狀態(tài)轉(zhuǎn)移函數(shù)的優(yōu)化精度是影響算法性能關鍵因素。合理設置粒子群規(guī)模是提高優(yōu)化估計精度的有效方法，以α1和β1的優(yōu)化為例，設4個系數(shù)分別為α1=6.0×10-3，α2=2.0×10-5，β1=5.0×10-3，β2=2.0× 10-5，將粒子群規(guī)模分別取30、40、50及60，逐個進行優(yōu)化操作，獲得優(yōu)化結(jié)果如表1所示。分析表1結(jié)果，可見種群規(guī)模的擴大可以提高優(yōu)化結(jié)果的精度。

表1　狀態(tài)轉(zhuǎn)移函數(shù)參數(shù)粒子群優(yōu)化性能表

即便如此，估計誤差仍是不可避免的，隨著距離步進上的迭代估計，反演誤差仍會逐漸累積，本文算法通過2.2節(jié)的精估計修正反演結(jié)果，消除誤差積累。如2.2節(jié)所述，遠距離的反演問題可將距離區(qū)間分段處理，分段的基本原則是區(qū)間分段不宜過大，線性差值區(qū)間不妨取為分段區(qū)間的后半段，考慮到不同分段策略對反演精度的影響不是本文重點，不予詳細討論。若δ為真實波導高度，δ為反演波導高度，則反演相對誤差可定義為

波導強度的反演相對誤差定義方法相同。以50 km范圍內(nèi)的反演為例，本文將其分為0 km～16 km、16 km～32 km、32 km～40 km、40 km～48 km，狀態(tài)轉(zhuǎn)移函數(shù)的系數(shù)設置同前，得反演結(jié)果的相對誤差分布圖2，從圖中可以看出，4個優(yōu)化點處的反演相對誤差均為極小值點，由于狀態(tài)轉(zhuǎn)移函數(shù)優(yōu)化誤差的影響，每個距離區(qū)間的前半段內(nèi)的反演相對誤差隨距離增加而變大，由于采取了線性插值，區(qū)間后半段內(nèi)的反演相對誤差被有效降低，整個50 km范圍內(nèi)的反演相對誤差不超過6%。

圖2　蒸發(fā)波導狀態(tài)參數(shù)反演相對誤差

基于上述分析，本文通過仿真實驗檢驗算法的反演性能，為方便研究，僅以α1不同取值條件下的反演為例，設粒子群規(guī)模為60，α1=6.0×10-3，α2、β1和β2取值不變，獲得反演結(jié)果如下頁圖3所示；設α1=4.5×10-3，α2、β1和β2取值不變，獲得反演結(jié)果如下頁圖4所示，其中δ為波導高度，ΔM為波導強度，由圖中可見，反演結(jié)果的相對誤差均不超過5%。

分析圖4仿真結(jié)果，精估計對粒子群優(yōu)化估計誤差造成的反演誤差積累起到了重要的修正作用，通過線性插值，平滑了各距離區(qū)間內(nèi)折射率剖面參數(shù)的反演結(jié)果。相對于本文算法，沒有誤差修正的算法在反演初始階段尚能提供較高的反演精度，隨著距離的增加，前期粒子群優(yōu)化的估計誤差使得粒子濾波的反演結(jié)果出現(xiàn)明顯的誤差積累，直至完全偏離真實波導折射率剖面參數(shù)。

4　結(jié)論

在波導剖面參數(shù)水平方向分布連續(xù)平滑的條件下將狀態(tài)轉(zhuǎn)移函數(shù)以級數(shù)展開的二階近似表示，將折射率剖面初始狀態(tài)觀測問題轉(zhuǎn)化為參數(shù)估計問題，采用粒子群優(yōu)化估計初始狀態(tài)和狀態(tài)轉(zhuǎn)移函數(shù)，該算法的不足之處是時效性有待提高，下一步可致力于降低反演運算量，提高實時性等方面的研究工作。

圖3　α1=6.0×10-3條件下的反演結(jié)果

圖4　α1=4.5×10-3條件下的反演結(jié)果

［1］KARIMIAN A，YARDIM C,HODGKISS W S.Estimation of radio refractivity using a multiple angle clutter model［J］.Radio science，2012，47，RS0M07.

［2］KARIMIAN A，YARDIM C，GERSTOFT P.Multiple grazing angle sea clutter modeling［J］.IEEE Trans on Antennas and Propagation，2012，60（9）：4408-4417.

［3］KARIMIAN A，YARDIM C，GERSTOFT P.Refractivity estimation from sea clutter：An invited review［J］.Radio science，2011，46，RS6013.

［4］MUSSON G L，GAUTHIER S，BRUTH E.A simple method determine evaporation duct height in the sea surface boundary layer［J］.Radio science，1992，27（5）：635-644.

［5］田斌,孔大偉,周沫.蒸發(fā)波導迭代MGB模型適用性研究［J］.電波科學學報，2013，28（3）：590-594.

［6］PAULUS R A.Specification for environmental measurements to assess radar sensors［C］//Technical Document 1685，DTIC Document，Naval Ocean Systems Center.San Diego，California，1989.

［7］JESKE H.State and limits of prediction methods of radar wave propagation conditions over sea［C］//Proc.of Modern Topics in Microwave Propagation and Air-Sea Interaction, 1973：130-148.

［8］LEVY M.Parabolic equation methods for electromagnetic wave propagation［M］.London：IEEE Press，2000.

［9］張青洪，廖成，盛楠.拋物方程法的非均勻網(wǎng)格技術研究［J］.電波科學學報，2013，28（4）：635-640.

［10］BALVEID G C，WALTER F.Analysis of GPS signal propagation in tropospheric ducts using numerical methods［C］// Proc.of 11th URSI Comission F Open Symposition on Radio Wave Propagation and Remote Sensing，Rio de Janeiro,，RJ，Brazil,2007.

［11］BALVEDI G C，WALTER F.GPS signal propagation in tropospheric ducts［C］//Proc.of The XXIXth URSI General Assembly，Chicago，Illinois，USA，2008：9-16.

［12］龐佳瑋，杜曉燕，張水蓮.粒子群優(yōu)化算法反演無線電波導傳輸損耗分布［J］.電波科學學報，2013，28（1）：14-20.

A Horizontally Nonuniform Evaporation Duct Inversion Algorithm Based on Optimal Filtering

DONG Dao-guang，RUI Guo-sheng，KANG Jian，TIAN Wen-biao
（Signal and Information Processing Provincial Key Laboratory in Shandong，Naval Aeronautical and Astronautical University，Yantai 264000，China）

In the application of horizontally nonuniform evaporation duct，the meteorological parameters observation of remote waters and horizontal distribution regularity of atmospheric refractive index profile is hard to get.To solve the problem，this article puts forward a horizontally nonuniform evaporation duct refractive index profile parameters inversion algorithm based on optimal filtering.This algorithm estimates the initial state of refractive index profile and the state transfer function that describes the horizontal distribution regularity of profile parameters with PSO，establishes the system equation and observation equation based on the DMFT step calculation forward model of GPS sea surface scattering signal propagation in duct，then makes iterative estimation of nonuniform refractive index profile parameters on the horizontal step with particle filtering.The accurate estimation of this algorithm effectively modifies the error accumulation caused by the PSO estimation error，which ensures the inversion precision and enable the algorithm inverse in segments when the distance is divided.This algorithm provides a feasible method to inverse the evaporation duct refractive index profile parameters of remote waters.

evaporation duct，PSO，particle filtering，inversion

TN951

A

1002-0640（2016）07-0077-05

2015-06-05

2015-07-09
*

國家自然科學基金資助項目（41476089）

董道廣（1990-），男，山東濟南人，碩士。研究方向：蒸發(fā)波導反演、雷達探測技術。