基于优化滤波的水平非均匀蒸发波导反演算法*

董道广，芮国胜，康　健，田文飚

（海军航空工程学院信号与信息处理山东省重点实验室，山东　烟台　264000）

基于优化滤波的水平非均匀蒸发波导反演算法*

董道广，芮国胜，康健，田文飚

（海军航空工程学院信号与信息处理山东省重点实验室，山东烟台264000）

针对水平非均匀蒸发波导反演中大气折射率剖面水平分布规律不确知和大范围海域气象参数观测难的问题，提出一种基于优化滤波的水平非均匀蒸发波导折射率剖面参数反演算法，该算法采用粒子群优化估计折射率剖面参数初始状态和描述剖面参数水平分布规律的状态转移函数，结合GPS海面散射信号波导传播的DMFT步进计算正演模型和实测GPS散射信号功率建立系统方程和观测方程，并采用粒子滤波迭代估计水平步进上的非均匀折射率剖面参数。该算法的精估计环节有效修正了粒子群估计误差在后续反演过程中导致的误差积累，保证了反演精度，并使得算法具备距离分段条件下逐段反演的能力，为远距离海域蒸发波导折射率剖面参数的反演提供了可行的方法。

蒸发波导，粒子群优化，粒子滤波，反演

0　引言

蒸发波导是出现概率最高的海上大气波导，过去人们对其反演问题的研究更多停留在均匀蒸发波导，事实上，蒸发波导因海上温度、湿度、气压及风速等气象因子在大尺度区域分布上的差异性而呈现出空间分布的非均匀性。

Karimian［1-3］采用现代滤波技术对时变大气波导进行了卓有成效的反演。由于滤波反演算法的收敛性依赖于系统方程和初始状态，上述困难严重影响到了反演的精度。为了解决上述两方面问题，本文提出一种基于优化滤波的GPS海面散射信号反演水平非均匀蒸发波导算法。本文算法解决的是空间上的反演问题。该算法将前人提出的离散混合傅立叶步进计算方法（DMFT）作为GPS海面散射信号蒸发波导传播的正演工具，通过粒子群优化估计滤波反演算法的初始状态和系统方程，依据估计值建立系统方程和观测方程，并将正演GPS海面散射信号波导传播的DMFT步进计算和粒子滤波的迭代状态估计同步处理。对于远距离、大范围的反演，该算法采取分段优化的方法来修正估计误差引起的反演误差积累。

1　反演建模

1.1折射率剖面参数化建模

折射率剖面参数既是反演输入量，也是反演输出量，故必须对折射率剖面作参数化建模。经典的蒸发波导剖面模型有MGB［4-5］模型和P-J模型等，本文采用P-J模型，该模型适用于热中性气象条件，是 Paulus［6］在 Jeske［7］的研究基础上基于Monin-Obukhov相似理论提出的，其数学解析式为

其中，z为高度，δ为波导高度，z0为空气动力学粗糙度因子，M0为海面修正折射率。图1给出了该剖面模型的示意，根据该模型可确定两个参数，分别为波导高度δ和波导强度ΔM。

图1　P-J模型参数化建模示意

在海上复杂的气象参数环境下，实际折射率剖面不如模型剖面平滑，所以实现对实际蒸发波导折射率剖面的高精度拟合要求模型具备灵活的参数调整能力，以使模型能够遍历和构成的二维可行域。本文以P-J模型为依据，对折射率剖面作二维参数化建模，现导出其数学解析式，给定参数δ和ΔM，设z=δ以下高度范围内的修正大气折射率为含辅助参变量的函数，如式（2）

z=δ处的修正大气折射率为满足，

解式（3）可得

z≥δ范围内的大气修正折射率仍如式式（1）所示，为使式（1）和式（2）在z=δ处计算的数值相等，该范围内修正折射率为

其中

式（2）～式（6）即为反演所用P-J二参数剖面模型。

1.2观测方程及系统方程建模

观测方程描述折射率剖面参数与正演GPS海面散射信号的函数关系，系统方程描述折射率剖面参数在水平方向上按距离步进的更新。为方便分析计，将折射率剖面参数记为，根据抛物方程的DMFT步进算法［8-9］，正演GPS海面散射信号步进传播损耗的自变量即当前步进处的垂直修正折射率剖面，据此将该步进的观测方程简记为式（7）

式中，L代表传播损耗，下标n为当前步进的距离索引，上标j为参数δ和ΔM在当前步进剖面可行域内众多参数中的序号，v是加性白噪声，方差为R。根据式（7）即可获得GPS海面散射信号的正演功率序列

其中，下标sim表示正演，其后的数字标号表示不同步进的距离索引。GPS海面直射波本质上是球面波，由于卫星到海面距离较大，可将直射波近似为平面波［10-11］，GPS散射功率接收机接收到的是一个功率-延迟序列，该序列是由海面散射信号在水平方向上不同距离处的接收延迟引起的，可记为

其中，下标obs表示实测，其后的数字标号意义同式（9）。反演的本质即通过式（8）和式（9）的拟合，求取拟合度最高的折射率剖面参数δopt和ΔMopt。

其中，F是当前距离步进到下一距离步进的剖面参数状态转移函数，是为了方便分析而增加的一个加性白噪声。

2　反演算法

2.1初始状态及系统方程估计

初始状态即，是GPS海面散射正演起算点的折射率剖面参数。系统方程的建立依赖于状态转移函数F，假设F有确定的解析形式且海上折射率剖面在水平方向上的变化是连续而平滑的，则可将作级数展开，并取二阶近似，如式（11）

综上，初始状态及系统方程的估计即是对参数δ1，ΔM1，α1，α2，β1和β2的估计，考虑到国内已有粒子群优化反演蒸发波导较成功的应用［12］，本文采用粒子群优化给出初始状态的估计流程：

1）结合式（7）～式（9），建立各距离步进上度量GPS散射信号正演和实测功率拟合度的目标函数

5）设fmin为停止迭代更新的阈值，反复执行上述4步，直到达到迭代次数上限或满足fgbest小于fmin，停止迭代，获得初始状态估计值；

α1，α2，β1和β2估计流程类似于上述流程，所不同的是其目标函数是经多步DMFT步进计算获得的正演功率序列构成的，如式（15）

2.2水平非均匀折射率剖面参数反演

水平非均匀折射率剖面参数的反演由粗估计和精估计两部分组成，粗估计基于2.1节估计结果建立观测方程和系统方程，采用粒子滤波完成一段距离范围内的参数反演；精估计对粗估计最后一个距离步进上的δN，Δ，α1，α2，β1和β2的粒子群优化，通过线性插值修正2.1节的估计误差在后续粒子滤波过程中造成的误差积累，并平滑反演结果。现给出粗估计的步骤：

1）根据式（7）～式（8）建立观测方程，根据式（10）～式（11）建立系统方程；

2）设初始状态的概率密度函数为均匀分布类型，依此分布随机产生S个状态粒子，，…，，结合系统方程对每个粒子作时间更新，获得下一步进上的先验状态粒子，，…；

3）由观测方程计算每个先验状态粒子在下一步进上的正演功率，i∈［1，S］，并结合该步进的实测功率Pobs，2计算每个先验状态粒子的似然概率密度，如式（16）～式（17）

4）基于归一化的似然概率密度对先验状态粒子重采样，获得后验状态粒子，对后验状态粒子取均值即得下一步进上的状态粒子估计；

5）反复执行上述4步，直至获得0～R1距离范围内的剖面参数，…。

假设距离对应的步进索引为N，现在粗估计R1步骤的基础上给出精估计步骤：

1）参照2.1节粒子群优化步骤重新估计所在步进上的剖面参数状态矢量，得；

2）设距离r对应的步进索引为n，对r～R1距离范围内各步进上的剖面参数作线性插值，获得，…，，如式（18）

对于较大距离范围内的反演，可将距离区间分为数段，依次执行上述反演步骤。

3　算法性能分析及仿真实验

针对波导折射率剖面水平分布规律不确知和大范围海域气象参数观测难的问题，基于优化滤波的反演算法将折射率剖面参数的水平分布以现代滤波理论的系统方程形式表达，对其状态转移函数进行二次有理多项式近似，采用粒子群优化估计GPS海面散射区域的初始折射率剖面参数和状态转移函数。由于水平方向上的折射率剖面参数是连续而平滑的，二阶近似的误差仅为状态转移函数对距离步进二阶导的高阶无穷小，可以忽略不计。综上分析，状态转移函数的优化精度是影响算法性能关键因素。合理设置粒子群规模是提高优化估计精度的有效方法，以α1和β1的优化为例，设4个系数分别为α1=6.0×10-3，α2=2.0×10-5，β1=5.0×10-3，β2=2.0× 10-5，将粒子群规模分别取30、40、50及60，逐个进行优化操作，获得优化结果如表1所示。分析表1结果，可见种群规模的扩大可以提高优化结果的精度。

表1　状态转移函数参数粒子群优化性能表

即便如此，估计误差仍是不可避免的，随着距离步进上的迭代估计，反演误差仍会逐渐累积，本文算法通过2.2节的精估计修正反演结果，消除误差积累。如2.2节所述，远距离的反演问题可将距离区间分段处理，分段的基本原则是区间分段不宜过大，线性差值区间不妨取为分段区间的后半段，考虑到不同分段策略对反演精度的影响不是本文重点，不予详细讨论。若δ为真实波导高度，δ为反演波导高度，则反演相对误差可定义为

波导强度的反演相对误差定义方法相同。以50 km范围内的反演为例，本文将其分为0 km～16 km、16 km～32 km、32 km～40 km、40 km～48 km，状态转移函数的系数设置同前，得反演结果的相对误差分布图2，从图中可以看出，4个优化点处的反演相对误差均为极小值点，由于状态转移函数优化误差的影响，每个距离区间的前半段内的反演相对误差随距离增加而变大，由于采取了线性插值，区间后半段内的反演相对误差被有效降低，整个50 km范围内的反演相对误差不超过6%。

图2　蒸发波导状态参数反演相对误差

基于上述分析，本文通过仿真实验检验算法的反演性能，为方便研究，仅以α1不同取值条件下的反演为例，设粒子群规模为60，α1=6.0×10-3，α2、β1和β2取值不变，获得反演结果如下页图3所示；设α1=4.5×10-3，α2、β1和β2取值不变，获得反演结果如下页图4所示，其中δ为波导高度，ΔM为波导强度，由图中可见，反演结果的相对误差均不超过5%。

分析图4仿真结果，精估计对粒子群优化估计误差造成的反演误差积累起到了重要的修正作用，通过线性插值，平滑了各距离区间内折射率剖面参数的反演结果。相对于本文算法，没有误差修正的算法在反演初始阶段尚能提供较高的反演精度，随着距离的增加，前期粒子群优化的估计误差使得粒子滤波的反演结果出现明显的误差积累，直至完全偏离真实波导折射率剖面参数。

4　结论

在波导剖面参数水平方向分布连续平滑的条件下将状态转移函数以级数展开的二阶近似表示，将折射率剖面初始状态观测问题转化为参数估计问题，采用粒子群优化估计初始状态和状态转移函数，该算法的不足之处是时效性有待提高，下一步可致力于降低反演运算量，提高实时性等方面的研究工作。

图3　α1=6.0×10-3条件下的反演结果

图4　α1=4.5×10-3条件下的反演结果

［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］田斌,孔大伟,周沫.蒸发波导迭代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］张青洪，廖成，盛楠.抛物方程法的非均匀网格技术研究［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］庞佳玮，杜晓燕，张水莲.粒子群优化算法反演无线电波导传输损耗分布［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-），男，山东济南人，硕士。研究方向：蒸发波导反演、雷达探测技术。