一种高频天地波雷达海洋表面流反演方法与流程

本发明属于雷达技术领域，特别涉及一种高频天地波雷达海洋表面流反演方法，用于解决天地波超视距雷达面临的电离层状态未知以至于难以实现精准定位的问题。

背景技术：

高频天地波超视距雷达是一种适用于大范围海态监测的遥感新体制雷达设备，因其探测距离远、信噪比高等特点，近年来已经被广泛应用于目标探测、海洋动力学参数反演等领域。

天地波超视距雷达工作于天波反射/地波绕射的新传播模式，由内陆放置的发射机发射高频无线电波并斜入射到电离层，经过电离层反射到达地(海)面，再经地(海)面绕射传播到地波接收站。地波接收站可以放置在远离发射站的海岸或舰船上，相比传统的高频地波雷达具有很大的灵活性。然而，由于电离层具有分层、非平稳、时变等特性，导致无线电波的反射点位置难以确定，从而给探测与反演工作带来了一定的困难。

目前国内已有关于高频天地波雷达提取海流方面的研究，而这些研究都是建立于平面无倾斜电离层反射模型的基础之上，假设无线电波的波束照射在同一反射高度的平面电离层，通过解析几何来确定海面散射元的位置。这种简易模型忽略了由于电离层倾斜和高度变化引起的探测误差，会使散射元定位不准确，严重影响探测精度。这将使得高频天地波雷达的应用受到了极大的限制，即只有在电离层非常平稳无倾斜的状态下该简易模型才是有效的，而在实际情况中难以实现。高频段(3～30mhz)的无线波通常会在e/f层(200～300km)的高度发生反射，电离层的小角度倾斜和反射高度的微小变化都将引起无线电波在地(海)投射距离的巨大变化。因此，基于平面无倾斜电离层反射模型的天地波海态探测研究结果目前仍存在较大误差。

技术实现要素：

本发明针对背景技术存在的问题，提出了一种高频天地波雷达海洋表面流反演新方法，利用傅里叶级数展开法重新提取天地波模式下的真实一阶bragg频率，并利用最小二乘约束求解出天地波超视距雷达在发生海面散射时难以确定的俯仰角和双基地角，从而实现精准定位。

为达到上述目的，本发明采用如下技术方案：

一种高频天地波雷达海洋表面流反演方法，包括如下步骤，

步骤1：海洋回波数据预处理。

步骤2：将一阶bragg散射频率fb和由径向流产生的多普勒偏移量δf分别用傅里叶级数展开为关于群距离r和波达方向θ的函数。

步骤3：根据海洋回波的散射传播机理和多普勒效应建立了线性方程组，利用最小二乘拟合求解超定方程，提取一阶bragg频率。

步骤4：双站遍历搜索公共散射元，求解双基地角和俯仰角。

步骤5：联立双站的径向流结果以及投影关系，合成天地波模式下矢量流结果。

上述技术方案的特点在于：

步骤1中对实验数据作前期预处理，包括距离匹配滤波、通道校准、谱区划分及doa(directionofarrival)估计等关键技术，得到群距离r、波达方向θ和海洋回波的谱点信息f⁺/f^-。

步骤2中将一阶bragg散射频率fb和由径向流产生的多普勒偏移量δf看作在海洋表面上关于r和θ的分布函数，由若干不同频率的正弦波叠加组合而成。根据二维傅里叶级数展开公式可以得到：

其中，p1,p2是δf傅里叶展开的阶数，q1,q2是fb傅里叶展开的阶数,ε1和ε2分别为δf和fb傅里叶级数展开的残余量，可忽略；

其中，k＝1,2,3,4

x1,mn＝amn,x2,mn＝bmn,x3,mn＝cmn,x4,mn＝dmn；

y1,mn＝αmn,y2,mn＝βmn,y3,mn＝γmn,y4,mn＝ηmn；

a1,mn＝cos(mω1r)cos(nω2θ),a2,mn＝cos(mω1r)sin(nω2θ),。

a3,mn＝sin(mω1r)cos(nω2θ),a4,mn＝sin(mω1r)sin(nω2θ)

步骤3中根据海洋回波的散射传播机理和多普勒效应建立线性方程组。在实际的海洋探测过程中，海洋回波谱分布在零频两侧，且满足如下关系式：

即

其中，f⁺和f^-分别是正负一阶海洋回波的谱点，我们可以通过高频超视距雷达获得。

将步骤2中fb和δf的傅里叶级数表达式带入上式方程组，可得到如下超定方程组：

其中，a和β分别是δf和fb的傅里叶级数展开系数矩阵,下标的正负号分别对应正谱点和负谱点，

为保证线性方程组有解，方程组的个数不能少于未知变量的个数。然而，在长期的雷达观测试验中发现，探测结果会存在一定的随机误差，因此我们需要使用尽可能多的谱点将随机误差平均减弱。此时，问题就变成了超定方程组求解，这里采用最小二乘的方法求解：

minα||λα-f||²

其中，

得到

[xy]^t＝(λ^tλ)^-1λ^tf

为获得方程组的最优解，傅立叶展开阶数的选择也至关重要。高阶最小二乘解经常会带来数值计算的不稳定性，阶数太小又会导致收敛误差过大。在实际的探测试验中，δf比fb变化更剧烈，所以我们考虑适当减小fb的展开阶数，增大δf的展开阶数。综合以上几点，加上仿真测试，我们得到fb和δf的展开阶数经验值为：

p1＝3，p2＝7，q1＝q2＝3

于是，我们可以提取真实的一阶bragg频率fb为：

fbinversion＝[b+yb-y]^t

步骤4中假设照射在同一散射元上的电磁波发生各个方向的散射，两个雷达接收站分别接收到该散射元散射的回波，两个回波在空中的传播路径认为是相同的，沿海面散射的路径分别指向雷达接收站，如附图2所示。利用两站雷达数据群距离的差值进行匹配。假设对应散射元s，雷达接收站a的群距离为g1，雷达接收站b的群距离为g2。散射元s到雷达接收站a的实际距离为r1，到雷达接收站b的实际距离为r2。在假设空中传播路径相同的情况下，群距离与散射元实际距离应满足下式：

r1-r2＝g1-g2

由此可以确定两个天波雷达接收站的公共散射元s。傅里叶级数展开法可以确定两个雷达接收站的一阶bragg频率fb1和fb2。假设照射公共散射元的是同一列波束，则双站在公共散射元的地方具有相同的俯仰角γ，且双基地角满足：

β2＝β1+θ3

其中θ3为散射元到两个雷达接收站连线的夹角，如附图2。

天地波模式下的一阶bragg频率表达式如下：

其中，f0为雷达工作频率，c为电磁波传播速度，γ为掠射角，β为双基地角。

基于天地波模式下的一阶bragg频率表达式，可以联立求解俯仰角γ和双基地角β1和β2。进而通过解析几何求解散射元距离雷达接收站的实际距离r1或r2，实现坐标配准。

步骤5中首先计算出单站bragg散射方向的矢量流分量，立体几何关系如附图3。其中为天波照射海面的入射波，为反射波，则bragg散射方向的海流分量流vb可用下式计算：

其中vr为单站雷达接收到的径向流速，vb为vb在入射波与反射波角平分线上的分量，α为vb与vb的夹角，通过解析立体几何，可得到

双站接收时，将两个分量进行合成即可得到矢量流，如附图4所示。其中标号1和2分别代表两个接收站的标识，通过解如下方程组，即可求得矢量流的大小和方向。

其中κ1和κ2分别为两个接收站bragg散射方向的分量流vb1和vb2与矢量流方向的夹角，v是矢量流，κ是分量流vb1和vb2的夹角。

与现有技术相比，本发明的优势在于：

本发明提出的高频天地波雷达海洋表面流反演方法，通过傅里叶级数展开模型和双站搜索公共散射元的方法，实现精准定位，提高了天地波模式下的海流探测精度。该方法摆脱了平面无倾斜电离层的假设，求解出受电离层影响后的真实一阶bragg频率，进而求解海面上各散射元的俯仰角和双基地角。与传统方法相比，解决了电离层倾斜、反射点高度不同引起的距离跳变、定位不准确的问题，从而大大提升高频天地波雷达在海洋动力学探测方面的实际应用价值。

附图说明

图1是高频天地波雷达系统工作示意图；

图2是双站接收天地波模式的信道传播示意图；

图3是电磁波在海面发生bragg散射时海流分量关系示意图；

图4是矢量流合成示意图；

图5是东山站天波海洋回波距离多普勒谱图；

图6是天地波雷达与地波雷达矢量流结果对比图；图6(a)是2017年4月4日10:04海流对比图；图6(b)是2017年4月4日10:24海流对比图。

具体实施方式

下面结合实例对本发明作进一步的详细描述，此处所描述的实施示例仅用于说明和解释本发明，并不用于限定本发明。

一种高频天地波雷达海洋表面流反演方法，该实现方法主要包括雷达回波信号预处理和算法的具体实现，具体包括如下步骤，

步骤1：海洋回波数据预处理。武汉大学于2017年4月4进行了高频天地波雷达组网实验。天波发射站位于武汉，地波接收站分别位于福建东山和赤湖附近的一艘船只上。

对两个接收站的回波数据分别作处理。包括回波数据的匹配滤波、通道校准、谱区划分、doa(directionofarrival)估计等，从而得到群距离r、波达方向θ和海洋回波的谱点信息f⁺/f^-。附图5为东山站接收到的2017年4月4日10:04的天地波海洋回波距离多普勒谱。船载回波数据还包括晃动补偿与由船体金属引起的天线方向图校正的前期处理。

步骤2：将一阶bragg散射频率fb和由径向流产生的多普勒偏移量δf分别用傅里叶级数展开为关于群距离r和波达方向θ的函数。根据二维傅里叶级数展开公式可以得到：

其中，p1,p2是δf傅里叶展开的阶数，q1,q2是fb傅里叶展开的阶数,ε1和ε2分别为δf和fb傅里叶级数展开的残余量，可忽略；

其中，k＝1,2,3,4

x1,mn＝amn,x2,mn＝bmn,x3,mn＝cmn,x4,mn＝dmn；

y1,mn＝αmn,y2,mn＝βmn,y3,mn＝γmn,y4,mn＝ηmn

a1,mn＝cos(mω1r)cos(nω2θ),a2,mn＝cos(mω1r)sin(nω2θ),。

a3,mn＝sin(mω1r)cos(nω2θ),a4,mn＝sin(mω1r)sin(nω2θ)

将步骤1中得到的回波数据的群距离r、频点f+/f-和方位信息θ带入以上展开式。

步骤3：根据海洋回波的散射传播机理和多普勒效应建立线性方程组，最小二乘拟合求解超定方程，提取一阶bragg频率。在实际的海洋探测过程中，海洋回波谱分布在零频两侧，且满足如下关系式：

其中，f⁺和f^-分别是正负一阶海洋回波的谱点，我们可以通过高频超视距雷达获得。

将步骤2中fb和δf的傅里叶级数表达式带入步骤上式方程组，可得到如下超定方程组：

其中，a和β分别是δf和fb的傅里叶级数展开系数矩阵,下标的正负号分别对应正谱点和负谱点，

利用最小二乘求解超定方程，得到：

[xy]^t＝(λ^tλ)^-1λ^tf

其中，

为获得方程组的最优解，傅立叶展开阶数的选择也至关重要。高阶最小二乘解经常会带来数值计算的不稳定性，阶数太小又会导致收敛误差过大。在实际的探测试验中，δf比fb变化更剧烈，所以我们考虑适当减小fb的展开阶数，增大δf的展开阶数。综合以上几点，加上仿真测试，我们得到fb和δf的展开阶数经验值为：

p1＝3，p2＝7，q1＝q2＝3

于是，我们可以提取真实的一阶bragg频率fb为：

fbinversion＝[b+yb-y]^t

步骤4：双站遍历搜索公共散射元，求解双基地角和俯仰角。利用两站雷达数据群距离的差值进行匹配。假设对应散射元s，雷达接收站a的群距离为g1，雷达接收站b的群距离为g2。散射元s到雷达接收站a的实际距离为r1，到雷达接收站b的实际距离为r2。在假设空中传播路径相同的情况下，群距离与散射元实际距离应满足下式：

r1-r2＝g1-g2

由于天波群距离的分辨率δd＝10km，要达到上式完全相等会出现散射元堆叠与模糊的现象。于是我们修正约束条件为：

其中δr＝r1-r2,δg＝g1-g2。

通过遍历谱点，初步确定两个天波雷达接收站的公共散射元s。傅里叶级数展开法可以确定两个雷达接收站的一阶bragg频率fb1和fb2。假设照射公共散射元的是同一列波束，则双站在公共散射元的地方具有相同的俯仰角γ，且双基地角满足：

β2＝β1+θ3

其中θ3为散射元到两个雷达接收站连线的夹角，如附图2。

天地波模式下的一阶bragg频率表达式如下：

其中，f0为雷达工作频率，c为电磁波传播速度，γ为掠射角，β为双基地角。

联立fb1和fb2，求解俯仰角γ、双基地角β1和β2，以及实际散射元距离r。利用直达波反射高度h和r计算无倾斜时俯仰角γ0，约束条件返回修正为γ0-3＜γ＜γ0+3。进而通过解析几何求解散射元距离雷达接收站的实际距离r1或r2，实现坐标配准。

步骤5：联立双站的径向流结果以及投影关系，合成天地波模式下的矢量流。首先由径向流计算出单站bragg散射方向的矢量流分量vb，由下式计算：

其中vr为单站雷达接收到的径向流速，vb为vb在入射波与反射波角平分线上的分量，α为vb与vb的夹角，通过解析立体几何，可得到

双站接收时，将两个分量进行合成即可得到矢量流。其中标号1和2分别代表两个接收站的标识，通过解如下方程组，即可求得矢量流的大小和方向。

其中κ1和κ2分别为两个接收站bragg散射方向的分量流vb1和vb2与矢量流方向的夹角，v是矢量流，κ是分量流vb1和vb2的夹角。

附图6为2017年4月4日10:04和10:24天波矢量流与地波矢量流在空间海平面上的对比结果。其公共覆盖区域的流速大小均方根误差分别为16.4cm/s和12.8cm/s，流向均方根误差分别为15.7°和38.8°。

以上实验结果表明，本发明所提的一种高频天地波雷达海洋表面流反演方法对高频天地波雷达海流探测是适用的，并提高探测精度，具备与地波雷达相当的探测水平。