CN111624556A - 基于形态成分分析的气象雷达wtc抑制方法 - Google Patents

Abstract

本发明公开了一种基于形态成分分析的气象雷达WTC抑制方法。该方法首先对污染距离单元进行地杂波抑制预处理，其次在傅里叶变换域和短时傅里叶变换域对预处理输出信号分别进行稀疏表示，并采用分裂增广拉格朗日收缩算法进行不同形态成分的稀疏分离，进而抑制WTC并恢复该距离单元上的气象信号。仿真实验结果表明，本发明能有效抑制地杂波和WTC，并高精度恢复气象信号，适合工程应用。

Description

基于形态成分分析的气象雷达WTC抑制方法

技术领域

本发明属于气象雷达杂波抑制领域，具体涉及一种风力涡轮机杂波的抑制方法。

背景技术

风力发电在满足人们对电能日益增长的需求的同时，风电场数量的增多、涡轮机尺寸和转速的增加，使得涡轮机叶片的高速旋转引起的运动杂波给各类雷达的目标检测带来新挑战，尤其是气象雷达。和常见的杂波不同，WTC(wind turbine clutter,风力涡轮机杂波)的多普勒频谱严重展宽甚至产生重叠，气象目标常常淹没在WTC中，严重影响了气象信息的预测精度。为解决这一问题，美国、德国和西班牙等多国科学家先后提出三种WTC抑制方法，分别为：基于插值恢复的WTC抑制方法、基于距离-多普勒谱回归的WTC抑制方法以及基于信号分离的WTC抑制方法。

基于插值恢复的WTC抑制方法依据气象回波在相邻距离单元中的空间连续分布特性，采用多重二次插值算法，对任意维度矩阵进行插值恢复。插值恢复气象信号算法虽然能够有效地增强气象回波信号在各距离单元之间的空间连续性，但恢复的气象回波也不可避免地存在信息与分辨率的损失。

为提高上述插值恢复算法的精度，基于距离-多普勒谱回归(RDR)的WTC抑制方法逐渐进入人们的视野。RDR算法的核心内容是将被WTC污染的距离单元内的数据采用多项式拟合的方法估计其气象回波的平均多普勒频移和多普勒带宽，在此基础之上，通过二次拟合相邻距离单元气象回波多普勒频谱窗的分布函数，提高对被WTC污染的距离单元内气象目标的多普勒频谱拟合精度，进而提高气象信息的预测精度。

虽然上述两种算法均能有效地抑制WTC，但是他们的适用范围仅仅是被WTC污染的距离单元较少的小型风电场。2014年F.Uysal,I.Selesnick等学者针对气象雷达聚束工作模式，创造性地提出采用稀疏优化技术来分离多普勒频率时变的WTC和多普勒频率时不变的气象回波。该算法不再局限于从被WTC污染的距离单元及相邻距离单元入手抑制WTC，而是根据WTC和气象目标分别在时-频域和多普勒域具有不同的形态，设计对应的稀疏字典矩阵，并进行稀疏表示，采用形态成分分析方法(MCA)实现两信号的盲分离。虽然F.Uysal的研究为中、大型风电场气象雷达WTC的抑制提供了思路，但是如何精确选取求解各成分信号稀疏系数过程中各参数，优化算法，以达到快速收敛、减少运算量，实现分量信号无损分离的目的，该方法并未明确指出。

发明内容

本发明所要解决的技术问题是：针对WTC的多普勒频谱严重展宽甚至产生重叠，气象目标常常淹没在WTC中，严重影响了气象信息的预测精度的问题，本发明提出一种基于MCA(morphological component analysis，形态成分分析)的风力涡轮机杂波抑制方法，利用MATLAB仿真平台，在设计滤波器滤除地杂波的基础上，通过分析成分信号特征，在傅里叶变换域和短时傅里叶变换域对其分别进行稀疏表示，并采用分裂增广拉格朗日收缩算法进行稀疏分离，进而抑制WTC，高精度恢复气象信号。

本发明为解决上述技术问题采用以下技术方案：

本发明提出一种基于形态成分分析的气象雷达风力涡轮机杂波抑制方法，包括如下步骤：

步骤一、进行地杂波抑制预处理，得到含有气象信号、风力涡轮机杂波WTC和噪声的回波信号；

步骤二、进行气象信号和WTC稀疏域字典矩阵构造，得到回波形态成分信号；

步骤三、对回波形态成分信号进行稀疏分离，得到分离的气象信号和WTC；

步骤四、进行气象信号恢复，得到气象信号。

进一步的，本发明所提出的基于形态成分分析的气象雷达风力涡轮机杂波抑制方法，所述步骤一、进行地杂波抑制预处理，得到含有气象信号、风力涡轮机杂波WTC和噪声的回波信号，具体为：

滤波后得到输出信号y的表达式为

y＝x_{1_1}+x_{2_1}+n₁

其中x_{1_1}为滤波后的气象目标信号，x_{2_1}为滤波后的WTC信号，n₁为滤波后的高斯噪声。

进一步的，本发明所提出的基于形态成分分析的气象雷达风力涡轮机杂波抑制方法，所述步骤二、进行气象信号和WTC稀疏域字典矩阵构造，得到回波形态成分信号，具体为：

滤波后的雷达回波信号y，其形态成分信号可表示为

x_{1_1}＝F₁a₁,x_{2_1}＝F₂a₂

其中，F₁,F₂分别表示傅里叶变换阵和短时傅里叶变换阵，a₁,a₂分别为气象信号和WTC各自在F₁,F₂上的稀疏系数，从而有

y＝F₁a₁+F₂a₂。

进一步的，本发明所提出的基于形态成分分析的气象雷达风力涡轮机杂波抑制方法，步骤三、对回波形态成分信号进行稀疏分离，得到分离的气象信号和WTC，具体为：

求解稀疏系数a₁和a₂：采用l₁-范数作为稀疏诱导正则器估计稀疏系数a₁和a₂，

稀疏系数a₁和a₂的BP求解公式为

such that y＝F₁a₁+F₂a₂

其中，λ₁,λ₂为正则化参数，

为最优求解系数；

假定F₁和F₂均满足Parseval框架，即

F₁F₁ ^H＝I_K,F₂F₂ ^H＝I_K

其中，F_i ^H是F_i的右逆矩阵，I_K为阶数与采样脉冲数K相等的单位阵；

采用分裂增广拉格朗日收缩算法进行不同形态成分的稀疏分离，首先应用变量分离技术，引入满足u_i-a_i＝0的辅助变量u₁,u₂；其次，利用增广拉格朗日法ALM和交替方向乘子法ADMM交替最小化u_i,a_i，即

Initialize:μ＞0,c_i,i＝1,2

Repeat

c_i＝c_i-(u_i-a_i),i＝1,2

Until convergence.

其中，μ为步长，c_i为拉格朗日乘子，初始值为零向量；

将初始值a₁＝F₁ ^Hy，a₂＝F₂ ^Hy带入上述SALSA求解得到最优稀疏系数

进而得到不同形态成分的雷达回波信号y^*为

进一步的，本发明所提出的基于形态成分分析的气象雷达风力涡轮机杂波抑制方法，步骤四、基于MCA恢复的气象信号

为：

与现有技术相比，本发明具有以下技术效果：

1、基于形态成分分析的气象雷达WTC抑制这一方法，能够对中、大型风电场气象雷达WTC进行抑制，提高气象雷达的工作性能。

2、仿真实验结果说明，基于形态成分分析的气象雷达WTC抑制这一方法，在抑制地杂波的前提下，通过SALSA稀疏分离雷达回波形态成分，能有效抑制WTC，恢复高精度的气象信号。

附图说明

图1为本发明信号处理流程图。

图2为风力涡轮机与雷达的位置关系图。

图3为高通滤波器频率响应图。

图4为滤波前雷达回波信号(气象信号+WTC+地杂波+噪声)功率谱。

图5为滤波前雷达回波信号(气象信号+WTC+地杂波+噪声)时-频图。

图6为滤波后雷达回波信号(气象信号+WTC+噪声)时-频图。

图7为气象信号功率谱。

图8为WTC时-频图。

图9为分离的气象信号功率谱图。

图10为分离的WTC时-频图。

具体实施方式

下面结合附图对本发明的技术方案进行详细说明：

本发明主要研究基于形态成分分析的气象雷达WTC抑制方法，图1是方法的处理流程。假定第l个距离单元同时接收到气象信号、WTC信号和地杂波信号，雷达回波中的噪声为高斯噪声，定义回波信号x∈C^K，其表达式为x＝x₁+x₂+x₃+n，其中K为采样脉冲数，x₁为气象目标信号，x₂为WTC信号，x₃为地杂波信号，n为高斯噪声，所述气象雷达风力涡轮机杂波抑制方法，其主要步骤如下：

步骤一、雷达回波信号仿真，具体为：

在本发明中，雷达回波信号包括：气象信号、风力涡轮机信号和地杂波信号，这三种回波信号的仿真如下：

1.1气象信号：气象目标通常情况下可看作是以恒定或微变速度运动的散射单元。在某一时刻t采样的气象信号可以表示为：

其中，A_t为气象目标的幅度，f_c为雷达频率，R_t为目标与雷达之间的距离，v_t为目标运动速度，θ_t为目标与雷达视线夹角，c为光速。

1.2WTC信号：风力涡轮机主要由机舱、桅杆和叶片组成。机舱和桅杆静止不动，对分析WTC特性基本没有影响，而旋转的风轮机叶片产生的回波则是构成WTC的主要组成。本发明主要构建由旋转叶片产生的杂波模型。常见的风轮机叶片都具有长而宽的特点，因此，可用散射点叠加模型对叶片进行分析。图2为风力涡轮机与雷达的位置关系图。

在这里，设雷达与风机轴心距离为R_w，雷达视线的俯仰角为β，为方便分析，设置方位角为零。风机叶片个数为M，长度为L，并以速度v_w匀速旋转，叶片起始角为零。在某一时刻t接收到的理想风轮机叶片回波信号可表示为：

其中，幅度A_w和相位历程函数Ω_k(t)分别为

1.3地杂波信号：常见地杂波多为地面静止物体的回波信号，其多普勒频移f_clu＝0Hz，则某一时刻t采样的地杂波信号可表示为

x₃＝exp(-j2πf_clut)。

步骤二、地杂波抑制预处理，具体为：

为减少雷达回波信号中地杂波对形态成分分离算法的影响，设计一高通滤波器以滤除地杂波。为灵活控制滤波器各参数以保证较好的滤波效果，本发明采用窗函数设计法设计该滤波器。以凯塞窗为例，设计滤波器通带截止频率60Hz，阻带截止频率30Hz，阻带最小衰减为60dB，通带波纹为±0.00868dB。滤波后得到输出信号y的表达式为

y＝x_{1_1}+x_{2_1}+n₁

其中x_{1_1}为滤波后的气象目标信号，x_{2_1}为滤波后的WTC信号，n₁为滤波后的高斯噪声。

步骤三、气象信号和WTC稀疏域字典矩阵构造，具体为：

气象目标通常情况下可看作是以恒定或微变速度运动的散射单元。由运动目标多普勒频移特性可知，在某一距离单元上，单一气象目标的最大多普勒频移保持不变，即气象目标在多普勒域具有稀疏性。WTC中主要包含以恒定转速旋转的叶片回波。由旋转物体的多普勒特性可知，叶片回波的多普勒频移随时间的变化而变化，即WTC在时-频域具有稀疏性。MCA的适用条件要求，待处理信号由几个不同形态的成分线性混合，且每个形态成分分别在一个特定的变换域能够稀疏表示，而在其他的变换域不具有稀疏性。滤波后的雷达回波信号y满足该条件，其形态成分信号可表示为

x_{1_1}＝F₁a₁,x_{2_1}＝F₂a₂

其中，F₁,F₂分别表示傅里叶变换阵和短时傅里叶变换阵，a₁,a₂为气象信号和WTC各自在F₁,F₂上的稀疏系数。从而有

y＝F₁a₁+F₂a₂

步骤四、回波信号形态成分稀疏分离

由步骤三可知，MCA算法的核心不再是求解x_{1_1}和x_{2_1}，而是求解稀疏系数a₁和a₂。由于稀疏重构算法要求字典矩阵均是过完备矩阵，则上述公式的求解实际是一个病态问题，常采用l₁-范数作为稀疏诱导正则器估计稀疏系数a₁和a₂。

稀疏系数a₁和a₂的BP求解公式为

such that y＝F₁a₁+F₂a₂

其中，λ₁,λ₂为正则化参数，

为最优求解系数。

在这里，假定F₁和F₂均满足Parseval框架，即

F₁F₁ ^H＝I_K,F₂F₂ ^H＝I_K

其中，F_i ^H,i＝1,2是F_i的右逆矩阵，I_K为阶数与采样脉冲数K相等的单位阵。

SALSA首先应用变量分离技术，引入满足u_i-a_i＝0,i＝1,2的辅助变量u₁,u₂。其次，利用ALM(augmented Lagrangian method，增广拉格朗日法)和ADMM(alternatingdirection method of multipliers,交替方向乘子法)交替最小化u_i,a_i，即

Initialize:μ＞0,c_i,i＝1,2

Repeat

c_i＝c_i-(u_i-a_i),i＝1,2

Until convergence.

其中，步长μ可取0.5，它的数值影响算法的迭代速度，但不影响算法的最终迭代结果。c_i,i＝1,2为拉格朗日乘子，初始值为零向量，一般可取λ₁＝0.8，λ₂＝0.2。

将初始值a₁＝F₁ ^Hy,a₂＝F₂ ^Hy带入上述SALSA求解得到最优稀疏系数

进而得到不同形态成分的雷达回波信号y^*为

步骤五、气象信号恢复

基于MCA恢复的气象信号

为

下面通过计算机仿真验证本方法的有效性。雷达系统仿真参数如表1所示，气象信号仿真参数如表2所示，风力涡轮机系统仿真参数如表3所示。在本仿真实验中，首先将雷达回波信号高通滤波，滤除静止的地杂波。滤波的前提是尽量将通过频率段信号的幅度影响降到最低。

图3为本发明设计的高通滤波器频率响应图。从图中不难看出，通带范围内高通滤波器的功率基本稳定在0dB，这说明通过该滤波器后信号的幅度基本不会发生变化。

图4、图5分别为滤波前雷达回波信号的功率谱和时-频图。从图4中并不能够看出信号中是否包含气象信号，而图5则可以清晰地看到，气象信号是实际存在的，只是淹没在WTC中。

高通滤波后的雷达回波信号时-频图如图6所示。从图中可以看出，地杂波已被滤除。

由于MCA要求待处理信号由几个不同形态成分线性混合，且每个形态成分分别在某一特定变换域可稀疏表示，而在其他变换域不具有稀疏性，这就需要对回波信号中的各成分进行分析，找到各自的稀疏域。

表1雷达系统仿真参数

频率f<sub>c</sub>(GHz)	3
		脉冲重复频率f<sub>r</sub>(Hz)	3000
采样脉冲数K	15000

表2气象信号仿真参数

信号幅度A<sub>t</sub>(V)	0.003
		目标与雷达的距离R<sub>t</sub>(m)	600
运动速度v<sub>t</sub>(m/s)	10
		目标与雷达视线夹角θ<sub>t</sub>(°)	15

表3风力涡轮机仿真参数

叶片长度L(m)	60
		叶片个数M	3
叶片转速v<sub>w</sub>(r/min)	12
		俯仰角β(°)	45
雷达与风机轴心距离为R<sub>w</sub>(m)	707

图7、图8分别为气象信号的功率谱和WTC的时-频图，由图可以得知，多普勒频移时不变的气象信号在傅里叶变换域具有稀疏性，而多普勒频移时变的WTC在时-频域具有稀疏性。因此本发明设定满足MCA的两个变换分为傅里叶变换和短时傅里叶变换。

求解形态成分信号在各稀疏域的稀疏系数是信号分离的主要任务。常用l₁-范数作为稀疏诱导正则器求解稀疏系数，这就是常见的BP问题。由于最小化稀疏系数的l₁-范数是一个凸优化的问题，无法通过具体的线性方程组求得解析解，只能通过迭代算法得到最优解，本发明采用全局收敛且精度较高的SALSA解决该问题。

图9、图10为采用SALSA进行信号分离后的气象信号的功率谱和WTC的时-频图。从图9中可以看出，功率谱中只包含有气象信号，说明WTC已被有效抑制。而图10中WTC的时-频图也表明，WTC得到了良好的分离。

对比图7和图9，在信噪比SNR＝18dB时，原气象信号的功率为32.94dB，滤波处理后恢复的气象信号功率为32.92dB，信号损失0.02dB。仿真结果说明，预处理后基于形态成分分析的WTC抑制方法能够有效地抑制杂波，恢复高精度的气象信号。

本发明首先将雷达回波信号高通滤波，滤除地杂波，其次针对MCA适用条件，在傅里叶变换域和短时傅里叶变化域稀疏表示气象信号和WTC，最后运用SALSA稀疏分离雷达回波形态成分，有效抑制WTC，高精度恢复气象信号，适合工程应用。

以上所述仅是本发明的部分实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的公开保护范围。

Claims (5)

1.基于形态成分分析的气象雷达风力涡轮机杂波抑制方法，其特征在于，包括如下步骤：

步骤一、进行地杂波抑制预处理，得到含有气象信号、风力涡轮机杂波WTC和噪声的回波信号；

步骤二、进行气象信号和WTC稀疏域字典矩阵构造，得到回波形态成分信号；

步骤三、对回波形态成分信号进行稀疏分离，得到分离的气象信号和WTC；

步骤四、进行气象信号恢复，得到气象信号。

2.根据权利要求1所述的基于形态成分分析的气象雷达风力涡轮机杂波抑制方法，其特征在于，

所述步骤一、进行地杂波抑制预处理，得到含有气象信号、风力涡轮机杂波WTC和噪声的回波信号，具体为：

滤波后得到输出信号y的表达式为

y＝x_{1_1}+x_{2_1}+n₁

其中x_{1_1}为滤波后的气象目标信号，x_{2_1}为滤波后的WTC信号，n₁为滤波后的高斯噪声。

3.根据权利要求2所述的基于形态成分分析的气象雷达风力涡轮机杂波抑制方法，其特征在于，

步骤二、进行气象信号和WTC稀疏域字典矩阵构造，得到回波形态成分信号，具体为：

滤波后的雷达回波信号y，其形态成分信号可表示为

x_{1_1}＝F₁a₁,x_{2_1}＝F₂a₂

其中，F₁,F₂分别表示傅里叶变换阵和短时傅里叶变换阵，a₁,a₂分别为气象信号和WTC各自在F₁,F₂上的稀疏系数，从而有

y＝F₁a₁+F₂a₂。

4.根据权利要求3所述的基于形态成分分析的气象雷达风力涡轮机杂波抑制方法，其特征在于，

步骤三、对回波形态成分信号进行稀疏分离，得到分离的气象信号和WTC，具体为：

求解稀疏系数a₁和a₂：采用l₁-范数作为稀疏诱导正则器估计稀疏系数a₁和a₂，

稀疏系数a₁和a₂的BP求解公式为

such that y＝F₁a₁+F₂a₂

其中，λ₁,λ₂为正则化参数，

为最优求解系数；

假定F₁和F₂均满足Parseval框架，即

F₁F₁ ^H＝I_K,F₂F₂ ^H＝I_K

其中，F_i ^H是F_i的右逆矩阵，I_K为阶数与采样脉冲数K相等的单位阵；

采用分裂增广拉格朗日收缩算法SALSA进行不同形态成分的稀疏分离，首先应用变量分离技术，引入满足u_i-a_i＝0的辅助变量u₁,u₂；其次，利用增广拉格朗日法ALM和交替方向乘子法ADMM交替最小化u_i,a_i，即

Initialize:μ＞0,c_i,i＝1,2

Repeat

c_i＝c_i-(u_i-a_i),i＝1,2

Until convergence.

其中，μ为步长，c_i为拉格朗日乘子，初始值为零向量；

将初始值a₁＝F₁ ^Hy，a₂＝F₂ ^Hy带入上述SALSA求解得到最优稀疏系数

进而得到不同形态成分的雷达回波信号y^*为

5.根据权利要求4所述的基于形态成分分析的气象雷达风力涡轮机杂波抑制方法，其特征在于，

步骤四、进行气象信号恢复，得到气象信号

为：

Images