CN110297247B - 基于低秩矩阵稀疏恢复的气象雷达风电场杂波抑制方法 - Google Patents

Abstract

本发明公开了一种基于低秩矩阵稀疏恢复的气象雷达风电场杂波抑制方法，该方法利用气象信号的空间相关性，首先将同时包含气象信号和风力涡轮机杂波的距离单元置零，并在该距离单元两侧对称取40个距离单元，然后逐脉冲将距离向量重构为满足零元素随机分布的低秩快拍矩阵，最后利用非精确增广拉格朗日乘子法(IALM)最小化核范数有效恢复气象数据。仿真实验结果表明，本发明能有效抑制风力涡轮机杂波(WTC)和噪声，提高气象信号输出信噪比，适合工程应用。

Description

基于低秩矩阵稀疏恢复的气象雷达风电场杂波抑制方法

技术领域

本发明涉及一种基于低秩矩阵稀疏恢复的气象雷达风电场杂波抑制方法，属于气象雷达杂波抑制领域。

背景技术

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

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

发明内容

本发明所要解决的技术问题是提供一种基于两维联合插值的气象雷达风电场杂波抑制方法，将稀疏优化理论引入气象雷达WTC抑制，研究基于矩阵补全的气象雷达小型风电场杂波抑制方法，其中矩阵补全理论(MC)作为稀疏恢复理论的两个重要研究方向之一，能够规避上述几种方法中造成气象信息缺失的问题，高精度矩阵补全被WTC干扰的气象信号，提高气象信息的预测精度。

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

本发明提供基于低秩矩阵稀疏恢复的气象雷达风电场杂波抑制方法，该方法的具体步骤如下：

步骤1，输入气象雷达回波信号，其中，第i个距离单元第m个脉冲下输入信号为：

x_i(m)＝s_i(m)+c_i(m)+w_i(m)+n_i(m)

式中，i＝1,...,L，L为距离单元数，m＝1,...,M，M为相干积累脉冲数，s_i(m)、c_i(m)、w_i(m)和n_i(m)分别为第i个距离单元第m个脉冲下的气象信号、地杂波信号、风力涡轮机杂波WTC信号和噪声信号；

步骤2，随机采样的低秩快拍矩阵构建，具体为：

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

由X_L×M构建随机采样的低秩快拍矩阵，其构建准则为：逐次将观测矩阵X_L×M第m个脉冲下的向量[x₁(m),x₂(m),...,x_L(m)]^T构建成快拍矩阵

其中，m₁和m₂分别表示快拍矩阵的行数和列数，m₁×m₂＝L，m₁＝m₂，快拍矩阵的第p行、第q列的元素

回波信号x_i(m)构建的低秩快拍矩阵X为：

由X定义气象信号s_i(m)第i个距离单元置零后构建的低秩快拍矩阵S为：

由X定义地杂波信号c_i(m)在第i个距离单元置零后构建的低秩快拍矩阵C为：

由X定义噪声信号n_i(m)在第i个距离单元置零后构建的低秩快拍矩阵N为：

WTC信号w_i(m)构建的低秩快拍矩阵W为零矩阵；

步骤3，通过矩阵补全模型恢复抑制WTC后的气象信号：

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

步骤4，利用非精确增广拉格朗日乘子法IALM求解矩阵补全模型，输出抑制WTC后的信号。

作为本发明的进一步技术方案，步骤4中拉格朗日函数为：

其中，Y＝Y₀+μ(X-S-N-W)为拉格朗日乘子矩阵；Y₀为拉格朗日乘子矩阵初值，取值为0；μ为惩罚因子，||·||_F表示F范数，||·||_F表示F范数，

tr(·)表示取矩阵的迹，

表示取复数的实部，<·,·>表示矩阵的内积。

作为本发明的进一步技术方案，利用非精确增广拉格朗日乘子法IALM求解的步骤为：

利用非精确增广拉格朗日乘子法IALM求解的步骤为：

1)令Y₀＝0、W₀＝0，N＝0，μ₀＞0，ρ＞1，k＝0，η＝10^-3，其中W₀＝0表示需要抑制的风力涡轮机杂波的初始值；

2)利用公式(U,Σ,V^H)＝svd(X-N_k-W_k+μ_k ^-1Y_k)和

更新S：

其中S_k+1和S_k分别表示气象信号S的第k+1和k次更新，W_k表示WTC信号W的第k次更新，N_k表示噪声N的第k次更新，Y_k表示拉格朗日乘子矩阵Y的第k次更新，μ_k表示惩罚因子μ的第k次更新；

3)更新W：

其中

表示Ω以外的指标集；

4)

5)更新Y_k：Y_k+1＝Y_k+μ_k(X-S_k+1-N_k+1-W_k+1)；

6)更新μ_k至μ_k+1：μ_k+1＝ρμ_k；

7)若下式不成立，则算法未收敛，令k←k+1,转步骤2，否则转步骤8：

||S_k-S_k-1||_F/||S_k||_F≤η；

8)结束循环，输出：

作为本发明的进一步技术方案，利用IALM逐次输出补全后恢复出的各脉冲下回波信号的快拍矩阵

依次提取出各矩阵中的第一列与最后一行，将第一列与最后一行构成M个L×1维的向量

然后将上述M个L×1维的向量组成一个M×L维的回波信号恢复矩阵

取

中的第i行向量作为稀疏恢复出的第i个距离单元的气象信号

有益效果

本发明利用IALM算法能够有效的恢复气象数据，抑制了风力涡轮机杂波，提高回波信号的信噪比，实用性强，运算量较小，具有很好的工程应用前景。

附图说明

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

图2是快拍矩阵特征值分布图；

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

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

图5是不同信噪比下信号恢复误差曲线图；

图6是IALM算法恢复前后气象信号速度值。

具体实施方式

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

本发明主要研究基于快拍矩阵的气象雷达风电场杂波矩阵补全抑制算法，图1是方法处理流程。其主要步骤如下：

步骤一、气象雷达回波信号的建模，具体为：

气象雷达接收信号主要包含气象信号、地杂波信号和噪声信号。假定第i个距离单元同时包含WTC信号，第m个脉冲下输入信号记为：

x_i(m)＝s_i(m)+c_i(m)+w_i(m)+n_i(m),m＝1,...,M (1)

其中，s_i(m)为气象信号、c_i(m)为地杂波信号、w_i(m)为WTC信号、n_i(m)为噪声信号，M为相干积累脉冲数，可取M＝64；

步骤二、随机采样的低秩快拍矩阵构建，具体为：

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

其中X_L×M的维数为L×M，L为距离单元数，可取L＝81，M为脉冲数；

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

逐次将观测矩阵X_L×M第m个脉冲下的向量[x₁(m),x₂(m),...,x_L(m)]^T构建成快拍矩阵

m₁,m₂满足m₁×m₂＝L，m₁＝m₂＝9，分别表示快拍矩阵的行与列数。快拍矩阵第p行、第q列的元素用k_p,q表示，即

满足

回波信号x_i(m)构建的低秩快拍矩阵X为：

由X定义气象信号s_i(m)第i个距离单元置零后构建的低秩快拍矩阵S为：

由X定义地杂波c_i(m)在第i个距离单元置零后构建的低秩快拍矩阵C为：

由X定义噪声n_i(m)在第i个距离单元置零后构建的低秩快拍矩阵N为：

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

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

步骤三、低秩矩阵稀疏恢复气象信号，具体为：

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

其中min(·)表示最小化处理，||·||_F表示F范数：

tr(·)表示取矩阵的迹，P_Ω表示投影到仅在指标集Ω非零的稀疏矩阵子空间上的映射，它使得矩阵在Ω中的元不变，Ω以外的元置零，用公式表示如下：

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

其中Y＝Y₀+μ(X-S-N-W)为拉格朗日乘子矩阵，Y₀为初值，取值为0，μ＞0表示惩罚因子，||·||_*表示核范数，为矩阵所有奇异值的和，

表示取复数的实部，＜X,Y＞＝tr(X^HY)表示矩阵的内积。

首先定义一个收缩运算符为：

IALM算法步骤如下所示：

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

1)令Y₀＝0、W₀＝0，N＝0，μ₀＞0，ρ＞1，k＝0，η＝10-³，其中W₀＝0表示需要抑制的风力涡轮机杂波的初始值，取值为0；

2)当未满足9)中收敛公式时，利用3)、4)解出

其中S_k表示气象信号S的第k次更新，W_k表示WTC信号W的第k次更新，N_k表示噪声N的第k次更新，Y_k表示拉格朗日乘子矩阵Y的第k次更新，μ_k表示惩罚因子μ的第k次更新；

3)(U,Σ,V^H)＝svd(X-N_k-W_k+μ_k ^-1Y_k)；

4)

5)更新W：

其中

表示Ω以外的指标集；

6)

7)更新Y_k：Y_k+1＝Y_k+μ_k(X-S_k+1-N_k+1-W_k+1)；

8)更新μ_k至μ_k+1：μ_k+1＝ρμ_k；

9)若下式不成立，则算法未收敛，令k←k+1,转步骤2，否则转步骤10：

||S_k-S_k-1||_F/||S_k||_F≤η；

10)结束循环，输出：

利用IALM迭代算法逐次输出补全后恢复出的各脉冲下回波信号的快拍矩阵

k为迭代次数，依次提取出各矩阵中的第一列与最后一行，将第一列与最后一行构成M个L×1维的向量

然后将M个L×1维的向量组成一个M×L维的回波信号恢复矩阵

取矩阵

中的第i行向量为：

其中

为稀疏恢复出的第i个距离单元的气象信号。

下面通过MATLAB对算法进行测试，验证矩阵补全算法恢复气象信号的有效性，雷达***仿真参数如表1所示；假设回波信号存在于第1到第100个距离单元，共有64个脉冲，采样矩阵第25个距离单元存在WTC，对其置零处理，依次将各脉冲下的向量构造为快拍矩阵。

表1快拍矩阵仿真参数

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

已知快拍矩阵为10行、10列的矩阵，将补全后的矩阵重新排列为L×1维矩阵，依次对剩余63个脉冲信号进行仿真，提取出的各脉冲下第25个向量元素构成了WTC抑制后的回波信号。为了清晰的了解本发明算法对WTC以及噪声抑制效果，图3给出了添加不同功率噪声情况下，气象信号的多普勒谱，由噪声和缺陷信号重构误差引起的基底噪声在矩阵补全抑制WTC之前波动较大。图3中仿真结果显示，在信噪比SNR＝30dB的情况下，矩阵补全前后气象信号中心频率约为360Hz，接近中心频率时，噪声干扰的影响很微弱，气象信号恢复精度更高。峰值旁瓣对高干扰噪声有一定的抑制作用，噪声功率降低约5～10dB。多普勒频率f_d＝300Hz、400Hz时噪声降低了15～20dB，抑制效果最佳。

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

根据图4的仿真数据可得基于快拍矩阵构建补全的气象信号幅值的均方根误差RMSE＝1.71e-2，定义方差Var来衡量抑制前后信号的离散程度

矩阵补全抑制前时域下气象信号幅值的方差Var1＝7.1107e-005，抑制后时域下气象信号幅值的方差Var2＝1.3303e-005，可以看出离散程度大大降低，与真实值的偏离程度减小，噪声干扰引起的波动降低，信噪比得到改善提高。仿真结果、数据分析表明，该算法大大提高了矩阵补全的性能，在抑制WTC的同时实现了气象信号的精确恢复。

为了分析了输入信号信噪比对本算法性能的影响，将不同信噪比下的信号恢复误差绘制出如图5所示曲线图，其中信噪比变化范围为0dB～30dB，从图中可以看出，对缺失数据进行矩阵补全的均方根误差随着输入信号的信噪比的增加而降低，输入信号的信噪比越高，进行奇异值分解时，噪声因素对奇异值分解的影响越小，矩阵补全恢复精度越高。

图6为IALM算法恢复前后气象信号速度的值，恢复后的气象信号速度为：

其中∠为取相位，

表示回波采样序列的一阶自相关参数，其中

为恢复后的气象信号的共轭转置。径向速度估计在低信噪比下波动较大，与真值存在较大偏差。但随着信噪比的增加，径向速度估计的误差逐渐减小，最终收敛到真值。

该发明首先利用WTC检测算法，剔除干扰数据，其次依据气象信号的空间相关性，构建出满足缺失元素随机分布的低秩补全矩阵，并通过IALM迭代算法解决约束秩最小化问题有效恢复气象数据。仿真实验结果表明，该算法能有效抑制WTC和噪声干扰，提高回波信号信噪比，具有很好的工程应用前景。

Claims (1)

1.基于低秩矩阵稀疏恢复的气象雷达风电场杂波抑制方法，其特征在于，该方法的具体步骤如下：

步骤1，输入气象雷达回波信号，其中，第i个距离单元第m个脉冲下输入信号为：

x_i(m)＝s_i(m)+c_i(m)+w_i(m)+n_i(m)

式中，i＝1,...,L，L为距离单元数，m＝1,...,M，M为相干积累脉冲数，s_i(m)、c_i(m)、w_i(m)和n_i(m)分别为第i个距离单元第m个脉冲下的气象信号、地杂波信号、风力涡轮机杂波WTC信号和噪声信号；

步骤2，随机采样的低秩快拍矩阵构建，具体为：

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

由X_L×M构建随机采样的低秩快拍矩阵，其构建准则为：逐次将观测矩阵X_L×M第m个脉冲下的向量[x₁(m),x₂(m),...,x_L(m)]^T构建成快拍矩阵

其中，m₁和m₂分别表示快拍矩阵的行数和列数，m₁×m₂＝L，m₁＝m₂，快拍矩阵的第p行、第q列的元素

回波信号x_i(m)构建的低秩快拍矩阵X为：

由X定义气象信号s_i(m)第i个距离单元置零后构建的低秩快拍矩阵S为：

由X定义地杂波信号c_i(m)在第i个距离单元置零后构建的低秩快拍矩阵C为：

由X定义噪声信号n_i(m)在第i个距离单元置零后构建的低秩快拍矩阵N为：

WTC信号w_i(m)构建的低秩快拍矩阵W为零矩阵；

步骤3，通过矩阵补全模型恢复抑制WTC后的气象信号：

s.t.P_Ω(X)＝P_Ω(S+N)

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

步骤4，利用非精确增广拉格朗日乘子法IALM求解矩阵补全模型，输出抑制WTC后的信号，具体为：

拉格朗日函数为：

其中，Y＝Y₀+μ(X-S-N-W)为拉格朗日乘子矩阵；Y₀为拉格朗日乘子矩阵初值，取值为0；μ为惩罚因子，||·||_F表示F范数，||·||_F表示F范数，

tr(·)表示取矩阵的迹，

表示取复数的实部，<·,·>表示矩阵的内积；

利用非精确增广拉格朗日乘子法IALM求解的步骤为：

1)令Y₀＝0、W₀＝0，N＝0，μ₀＞0，ρ＞1，k＝0，η＝10^-3，其中W₀＝0表示需要抑制的风力涡轮机杂波的初始值；

2)利用公式(U,∑,V^H)＝svd(X-N_k-W_k+μ_k ^-1Y_k)和

更新S：

其中S_k+1和S_k分别表示气象信号S的第k+1和k次更新，W_k表示WTC信号W的第k次更新，N_k表示噪声N的第k次更新，Y_k表示拉格朗日乘子矩阵Y的第k次更新，μ_k表示惩罚因子μ的第k次更新；

3)更新W：

其中

表示Ω以外的指标集；

4)

5)更新Y_k：Y_k+1＝Y_k+μ_k(X-S_k+1-N_k+1-W_k+1)；

6)更新μ_k至μ_k+1：μ_k+1＝ρμ_k；

7)若下式不成立，则算法未收敛，令k←k+1,转步骤2，否则转步骤8：

||S_k-S_k-1||_F/||S_k||_F≤η；

8)结束循环，输出：

利用IALM逐次输出补全后恢复出的各脉冲下回波信号的快拍矩阵

依次提取出各矩阵中的第一列与最后一行，将第一列与最后一行构成M个L×1维的向量

然后将上述M个L×1维的向量组成一个M×L维的回波信号恢复矩阵

取

中的第i行向量作为稀疏恢复出的第i个距离单元的气象信号

Images