CN110174651A - 基于低秩Hankel矩阵补全的气象雷达风电场杂波抑制方法 - Google Patents

Abstract

本发明公开了一种基于低秩汉克尔（Hankel）矩阵补全的气象雷达风电场杂波抑制方法，利用气象信号空间相关性，首先将同时包含气象信号和风力涡轮机杂波的距离单元置零，并在该距离单元两侧对称取32个距离单元，然后逐脉冲将距离向量重构为满足零元素随机分布的低秩Hankel矩阵，最后通过非精确增广拉格朗日乘子法（IALM）解决约束秩最小化问题有效恢复气象数据。仿真实验结果表明，该发明在抑制风力涡轮机杂波（WTC）的同时也抑制了噪声干扰，提高了气象信号的信噪比，具有很好的工程应用前景。

Description

基于低秩Hankel矩阵补全的气象雷达风电场杂波抑制方法

技术领域

本发明涉及气象雷达信号处理领域，尤其涉及一种基于低秩Hankel矩阵补全的气象雷达风电场杂波抑制方法。

背景技术

为应对能源的潜在危机和全球变暖等生态环境的恶化，世界各国都在积极开发可再生清洁能源，风力发电作为可再生清洁能源的一种重要形式受到了全世界的高度关注。近年来，全球范围内的风电场规模和数量正呈指数式增长，风轮机叶片的转速与长度不断增加，但研究表明，风电场风力涡轮机由于叶片高速旋转引起的运动杂波对雷达、通信导航等电子设备会产生严重影响，给各类雷达目标检测带来了新的挑战，而现有的杂波抑制技术均无法有效滤除风力涡轮机杂波(WTC)，对气象信息的预测精度产生了严重的影响，因此风力涡轮机杂波已成为当今气象雷达杂波抑制的核心问题。

现有杂波抑制技术如时域滤波方法、频域滤波方法、基于功率谱特征的滤波方法由于风力涡轮机的高速旋转造成的频谱展宽，使得气象信息损失严重，导致无法有效的抑制WTC，极大地影响了气象信息的预测精度。欧美科学家在详细分析了气象雷达不同工作模式下风力涡轮机杂波与气象回波的时、频域分布特性之后，提出的基于多重二次插值恢复、距离－多普勒谱回归、递归稀疏重构等风力涡轮机杂波抑制算法，也受风电场规模、风机转速、气象雷达工作模式等实际条件限制，上述算法均无法同时兼顾风力涡轮机杂波抑制与气象信息无损恢复。此外传统WTC抑制方法分别只单独处理每个距离单元的数据，而不需要利用其他距离单元的信息。且不能有效地抑制噪声信号。

发明内容

本发明所要解决的技术问题是针对背景技术中所涉及到的缺陷，提供一种基于低秩Hankel矩阵补全的气象雷达风电场杂波抑制方法。

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

基于低秩Hankel矩阵补全的气象雷达风电场杂波抑制方法，包括如下步骤：

步骤1)，建立气象雷达回波信号模型：

气象雷达接收信号包含气象信号、地杂波信号和噪声信号，令第i个距离单元包含WTC信号、其他距离单元均不包含WTC信号，第m个脉冲下输入信号记为x_i(m)＝s_i(m)+c_i(m)+w_i(m)+n_i(m),m＝1,..M.,，其中，s_i(m)为气象信号、c_i(m)为地杂波信号、w_i(m)为WTC信号、n_i(m)为噪声，M为脉冲数，s_i(m)与n_i(m)为目标信号，t_i(m)＝s_i(m)+n_i(m)；

步骤2)，随机采样构建低秩MC矩阵：

已知观测矩阵X_L×M，在第i个距离单元两侧各取32距离单元，并将第i个距离单元中的回波信号[x_i(1),x_i(2),...,x_i(M)]置零，得：

由X_L×M构建出低秩随机采样Hankel矩阵：逐次将观测矩阵X_L×M第m个脉冲下的向量[x₁(m),x₂(m),...,x_L(m)]^T构建成行为m₁、列为m₂的Hankel低秩矩阵，其中L为距离单元数，m₁、m₂满足m₁+m₂-1＝L；令Hankel矩阵的第p行、第q列的元素为h_p,q，则h_p,q＝h_p-1,q+1，且满足

则回波信号x_i(m)构建的Hankel矩阵X_H为：

由X_H定义气象信号s_i(m)第i个距离单元置零后构建的Hankel矩阵S_H为：

由X_H定义地杂波c_i(m)在第i个距离单元置零后构建的Hankel矩阵C_H为：

由X_H定义噪声n_i(m)在第i个距离单元置零后构建的Hankel矩阵N_H为：

已知WTC只存在于第i个距离单元，即观测矩阵中WTC信号为W_L×M：

在对W_L×M第i个距离单元中的信号置零后，WTC信号w_i(m)构建的Hankel矩阵W_H为零矩阵，即W_H＝0；

步骤3)，矩阵补全抑制WTC、恢复气象信号：

根据低秩补全理论，矩阵补全模型依据部分矩阵元素，通过约束秩最小化问题实现对未知元素的补全，将WTC剔除后，通过下面的矩阵补全模型恢复气象数据抑制WTC：

||P_Ω(N_H)||_F≤δ

其中，min(·)表示最小化处理，||·||_F表示F范数：P_Ω表示投影到仅在指标集Ω非零的稀疏矩阵子空间上的映射，它使得矩阵在Ω中的元不变，Ω以外的元置零，用公式表示为

利用非精确增广拉格朗日乘子法(IALM)解矩阵补全模型的最优化问题，将经典矩阵补全算法扩展到复数域，对应得拉格朗日函数L(S_H,W_H,N_H,Y,μ)表示为：

其中，Y＝Y₀+μ(X_H-S_H-N_H-W_H)为拉格朗日乘子矩阵，初值为0，μ＞0表示惩罚因子，||·||_*表示核范数，为矩阵所有奇异值的和，表示取复数的实部，＜X,Y＞＝tr(X^HY)表示矩阵的内积；

非精确增广拉格朗日乘子法的具体步骤如下：

输入：X_ij观测样本，(i,j)∈Ω，矩阵X_H∈R_m×n；

1.初始化：Y₀＝0,(W_H)₀＝0,N_H＝0,μ₀＝0,ρ>1,k＝0,η＝10^-3；

2.P_Ω(A)未收敛，求解

3.更新WTC矩阵：

4.更新噪声矩阵；

5.更新Y；Y_k+1＝Y_k+μ_k(X_H-(S_H)_k+1-(N_H)_k+1-(W_H)_k+1)；

6.更新μ、k；μ_k+1＝ρμ_k，k←k+1，若P_Ω(A)收敛，结束迭代；

输出：

利用非精确增广拉格朗日乘子法逐次输出补全后恢复出的各脉冲下目标信号的Hankel矩阵k为迭代次数，依次提取出各Hankel矩阵中的第一列与最后一行，将其构成M个L×1维的向量然后将其组成一个M×L维矩阵

取矩阵中的第i行向量为 即为基于Hankel矩阵补全恢复出的第i个距离单元的气象信号。

本发明采用以上技术方案与现有技术相比，具有以下技术效果：

本发明将MC理论引入气象雷达WTC抑制，研究基于矩阵补全(MC)的气象雷达小型风电场杂波抑制方法，能够规避现有技术的几种方法中造成气象信息缺失的问题，且通过高精度矩阵补全被WTC干扰的气象信号，提高了气象信息的预测精度。

附图说明

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

图2为Hankel矩阵特征值分布图；

图3无噪干扰气象信号、噪声干扰下矩阵补全前后气象信号功率谱；

图4无噪干扰气象信号、噪声干扰下矩阵补全前后气象信号幅值；

图5为IALM算法恢复前后气象信号速度值；

图6为不同信噪比下信号恢复误差曲线图。

具体实施方式

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

本发明公开了一种基于低秩Hankel矩阵补全的气象雷达风电场杂波抑制方法，如图1所示，其步骤如下：

步骤1)，建立气象雷达回波信号模型：气象雷达接收信号包含气象信号、地杂波信号和噪声信号，令第i个距离单元包含WTC信号、其他距离单元均不包含WTC信号，第m个脉冲下输入信号记为x_i(m)＝s_i(m)+c_i(m)+w_i(m)+n_i(m),m＝1,...,M，其中，s_i(m)为气象信号、c_i(m)为地杂波信号、w_i(m)为WTC信号、n_i(m)为噪声，M为脉冲数，s_i(m)与n_i(m)为目标信号，t_i(m)＝s_i(m)+n_i(m)。

步骤2)，随机采样构建低秩MC矩阵：

已知观测矩阵X_L×M，在第i个距离单元两侧各取32距离单元，并将第i个距离单元中的回波信号[x_i(1),x_i(2),...,x_i(M)]置零，得：

由X_L×M构建出低秩随机采样Hankel矩阵，其构建准则为：

逐次将观测矩阵X_L×M第m个脉冲下的向量[x₁(m),x₂(m),...,x_L(m)]^T构建成行为m₁、列为m₂的Hankel低秩矩阵，其中L为距离单元数，m₁、m₂满足m₁+m₂-1＝L。令Hankel矩阵的第p行，第q列的元素为h_p,q，有h_p,q＝h_p-1,q+1，且满足：

则回波信号x_i(m)构建的Hankel矩阵X_H为：

由X_H定义气象信号s_i(m)第i个距离单元置零后构建的Hankel矩阵S_H为：

由X_H定义地杂波c_i(m)在第i个距离单元置零后构建的Hankel矩阵C_H为：

由X_H定义噪声n_i(m)在第i个距离单元置零后构建的Hankel矩阵N_H为：

已知WTC只存在于第i个距离单元，即观测矩阵中WTC信号为W_L×M：

在对W_L×M第i个距离单元中的信号置零后，WTC信号w_i(m)构建的Hankel矩阵W_H为零矩阵，即W_H＝0。

步骤3)，矩阵补全抑制WTC、恢复气象信号：

根据低秩补全理论，矩阵补全模型可依据部分矩阵元素，通过约束秩最小化问题实现对未知元素的补全，将WTC剔除后，可以通过下面的矩阵补全模型恢复气象数据抑制WTC：

||P_Ω(N_H)||_F≤δ

其中，min(·)表示最小化处理，||·||_F表示F范数：P_Ω表示投影到仅在指标集Ω非零的稀疏矩阵子空间上的映射，它使得矩阵在Ω中的元不变，Ω以外的元置零，用公式表示如下：

利用非精确增广拉格朗日乘子法(IALM)解矩阵补全模型的最优化问题，经典的矩阵补全算法只适用于实矩阵，本文所构建矩阵为复矩阵，将经典矩阵补全算法扩展到复数域，对应得拉格朗日函数L(S_H,W_H,N_H,Y,μ)可以表示为：

其中Y＝Y₀+μ(X_H-S_H-N_H-W_H)为拉格朗日乘子矩阵，初值为0，μ＞0表示惩罚因子，||·||_*表示核范数，为矩阵所有奇异值的和，表示取复数的实部，＜X,Y＞＝tr(X^HY)表示矩阵的内积。

IALM算法步骤如下所示：

输入：X_ij观测样本，(i,j)∈Ω，矩阵X_H∈R_m×n

1.初始化：Y₀＝0,(W_H)₀＝0,N_H＝0,μ₀＝0,ρ>1,k＝0,η＝10^-3

2.P_Ω(A)未收敛，求解

3.更新WTC矩阵：

4.更新噪声矩阵：

5.更新Y：Y_k+1＝Y_k+μ_k(X_H-(S_H)_k+1-(N_H)_k+1-(W_H)_k+1)

6.更新μ、k：μ_k+1＝ρμ_k，k←k+1，若P_Ω(A)收敛，结束迭代

输出：

利用IALM迭代算法逐次输出补全后恢复出的各脉冲下目标信号的Hankel矩阵k为迭代次数，依次提取出各Hankel矩阵中的第一列与最后一行，将其构成M个L×1维的向量然后将其组成一个M×L维矩阵

取矩阵中的第i行向量为： 即为基于Hankel矩阵补全恢复出的第i个距离单元的气象信号。

下面通过MATLAB对算法进行测试，验证矩阵补全算法恢复气象信号的有效性，在本仿真实验中，设第1到第65个距离单元存在回波信号，共有64个脉冲，取采样矩阵的第10列，第33个距离单元存在WTC，对其置零处理，将其构造为Hankel矩阵形式，将补全后的矩阵重新排列L×1维矩阵，依次对剩余63个脉冲信号进行仿真，提取出的各脉冲下第33个向量元素构成了WTC抑制后的回波信号。仿真参数如下表所示：

表1Hankel矩阵仿真参数

脉冲重复频率	1000Hz
		载频	5.5GHz
雷达高度	1000m
		信噪比	35dB
脉冲数	64
		距离单元数	65
矩阵行数	33
		矩阵列数	33

如图2所示为Hankel矩阵特征值分布图，一个相干积分时间(CIT)中秩r远小于脉冲数M，即气象信号第m个脉冲下的列向量s_l(m)构建成的Hankel矩阵S满足低秩性要求。在接下来的仿真测试中，为了清晰的了解本发明算法对WTC以及噪声抑制效果，图3给出了添加不同功率噪声情况下，气象信号的多普勒谱，由噪声和缺陷信号重构误差引起的基底噪声在矩阵补全抑制WTC之前波动较大。图3中仿真结果显示，在信噪比SNR＝35dB的情况下，矩阵补全前后气象信号中心频率约为360Hz，接近中心频率时，噪声干扰较弱，气象信号恢复精度较高。峰值旁瓣对高干扰噪声有较明显的抑制作用，噪声功率降低约15dB，旁瓣处噪声的波动起伏降低，提高了信噪比。

气象信号幅值如图4所示，可以清楚的显示出矩阵补全前后的信号幅值的变化情况。通过矩阵补全恢复的幅值与原始值相比百分误差为2.5416e-002％。为了定量分析MC算法的恢复性能，定义如下均方根误差作为性能指标：

根据图4的仿真数据可得基于Hankel矩阵构建补全的气象信号幅值的均方根误差RMSE＝7.4e-3，可以看出补全后的气象信号幅值与真实值的偏离程度减小，噪声干扰引起的波动降低，信噪比得到改善提高，仿真结果表明，该算法大大提高了矩阵补全的性能，在抑制WTC的同时实现了气象信号的精确恢复。

图5为IALM算法恢复前后气象信号速度的值，恢复后的气象信号速度为：其中∠为取相位，表示回波采样序列的一阶自相关参数，其中为恢复后的气象信号的共轭转置，从图6可以看出，在低信噪比情况下，MC恢复出的气象信号速度方差大于高信噪比的情况，其速度估计波动较大，信噪比提高时，其速度估计波动较小，速度估计值最为精确。图6为不同信噪比下信号恢复误差曲线图，分析了输入信号信噪比对本算法性能的影响，其中信噪比变化范围为0dB～30dB，由图可知，对缺失数据进行矩阵补全的均方根误差随着输入信号的信噪比的增加而降低，输入信号的信噪比越高，进行奇异值分解时，噪声因素对奇异值分解的影响越小，矩阵补全恢复的精度越高。

本发明首先该发明利用气象信号空间相关性，首先检测出回波信号中同时含有气象信号和风电场杂波的距离单元，将其置零，然后依次将单个脉冲下的向量重构为满足零元素随机分布的低秩Hankel矩阵，最后通过IALM迭代算法解决约束秩最小化问题有效恢复气象数据。该算法实用性强，具有很好的工程应用前景。

本技术领域技术人员可以理解的是，除非另外定义，这里使用的所有术语(包括技术术语和科学术语)具有与本发明所属领域中的普通技术人员的一般理解相同的意义。还应该理解的是，诸如通用字典中定义的那些术语应该被理解为具有与现有技术的上下文中的意义一致的意义，并且除非像这里一样定义，不会用理想化或过于正式的含义来解释。

以上所述的具体实施方式，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施方式而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (1)

1.基于低秩Hankel矩阵补全的气象雷达风电场杂波抑制方法，其特征在于，包括如下步骤：

步骤1)，建立气象雷达回波信号模型：

气象雷达接收信号包含气象信号、地杂波信号和噪声信号，令第i个距离单元包含WTC信号、其他距离单元均不包含WTC信号，第m个脉冲下输入信号记为x_i(m)＝s_i(m)+c_i(m)+w_i(m)+n_i(m),m＝1,..M., ，其中，s_i(m)为气象信号、c_i(m)为地杂波信号、w_i(m)为WTC信号、n_i(m)为噪声，M为脉冲数，s_i(m)与n_i(m)为目标信号，t_i(m)＝s_i(m)+n_i(m)；

步骤2)，随机采样构建低秩MC矩阵：

已知观测矩阵X_L×M，在第i个距离单元两侧各取32距离单元，并将第i个距离单元中的回波信号[x_i(1),x_i(2),...,x_i(M)]置零，得：

由X_L×M构建出低秩随机采样Hankel矩阵：逐次将观测矩阵X_L×M第m个脉冲下的向量[x₁(m),x₂(m),...,x_L(m)]^T构建成行为m₁、列为m₂的Hankel低秩矩阵，其中L为距离单元数，m₁、m₂满足m₁+m₂-1＝L；令Hankel矩阵的第p行、第q列的元素为h_p,q，则h_p,q＝h_p-1,q+1，且满足

则回波信号x_i(m)构建的Hankel矩阵X_H为：

由X_H定义气象信号s_i(m)第i个距离单元置零后构建的Hankel矩阵S_H为：

由X_H定义地杂波c_i(m)在第i个距离单元置零后构建的Hankel矩阵C_H为：

由X_H定义噪声n_i(m)在第i个距离单元置零后构建的Hankel矩阵N_H为：

已知WTC只存在于第i个距离单元，即观测矩阵中WTC信号为W_L×M：

在对W_L×M第i个距离单元中的信号置零后，WTC信号w_i(m)构建的Hankel矩阵W_H为零矩阵，即W_H＝0；

步骤3)，矩阵补全抑制WTC、恢复气象信号：

根据低秩补全理论，矩阵补全模型依据部分矩阵元素，通过约束秩最小化问题实现对未知元素的补全，将WTC剔除后，通过下面的矩阵补全模型恢复气象数据抑制WTC：

||P_Ω(N_H)||_F≤δ

其中，min(·)表示最小化处理，||·||_F表示F范数：P_Ω表示投影到仅在指标集Ω非零的稀疏矩阵子空间上的映射，它使得矩阵在Ω中的元不变，Ω以外的元置零，用公式表示为

利用非精确增广拉格朗日乘子法(IALM)解矩阵补全模型的最优化问题，将经典矩阵补全算法扩展到复数域，对应得拉格朗日函数L(S_H,W_H,N_H,Y,μ)表示为：

其中，Y＝Y₀+μ(X_H-S_H-N_H-W_H)为拉格朗日乘子矩阵，初值为0，μ＞0表示惩罚因子，||·||_*表示核范数，为矩阵所有奇异值的和，表示取复数的实部，＜X,Y＞＝tr(X^HY)表示矩阵的内积；

非精确增广拉格朗日乘子法的具体步骤如下：

输入：X_ij观测样本，(i,j)∈Ω，矩阵X_H∈R_m×n；

1.初始化：Y₀＝0,(W_H)₀＝0,N_H＝0,μ₀＝0,ρ>1,k＝0,η＝10^-3；

2.P_Ω(A)未收敛，求解

3.更新WTC矩阵：

4.更新噪声矩阵；

5.更新Y；Y_k+1＝Y_k+μ_k(X_H-(S_H)_k+1-(N_H)_k+1-(W_H)_k+1)；

6.更新μ、k；μ_k+1＝ρμ_k，k←k+1，若P_Ω(A)收敛，结束迭代；

输出：

利用非精确增广拉格朗日乘子法逐次输出补全后恢复出的各脉冲下目标信号的Hankel矩阵k为迭代次数，依次提取出各Hankel矩阵中的第一列与最后一行，将其构成M个L×1维的向量然后将其组成一个M×L维矩阵

取矩阵中的第i行向量为 即为基于Hankel矩阵补全恢复出的第i个距离单元的气象信号。