CN112433215B - 基于先验知识加权的气象雷达风力涡轮机杂波抑制方法 - Google Patents

Abstract

本发明公开了一种基于先验知识加权的气象雷达风力涡轮机杂波抑制方法，首先将被风力涡轮机杂波(WTC)污染的距离单元置零，并逐脉冲将距离向量重构为随机采样的低秩托普利兹矩阵，利用气象信号观测矩阵的先验知识，通过对气象信号观测矩阵核范数中的每个奇异值赋予不同的权值，避免了用同一常数对所有的奇异值进行阈值化，更精确的获取原气象信号的托普利兹矩阵在低秩上的最优逼近矩阵，从而输出气象信号。仿真实验结果表明，本发明有效提高了气象信号的恢复精度，更好的抑制风力涡轮机杂波和噪声，适合工程应用。

Description

基于先验知识加权的气象雷达风力涡轮机杂波抑制方法

技术领域

本发明涉及气象雷达杂波抑制领域，尤其涉及一种基于先验知识加权的气象雷达风力 涡轮机杂波抑制方法。

背景技术

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

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

发明内容

本发明所要解决的技术问题是提供一种基于先验知识加权的气象雷达风力涡轮机杂波 抑制方法，将矩阵补全理论引入气象雷达WTC抑制，研究基于矩阵补全的气象雷达小型风 电场杂波抑制方法。矩阵补全方法能够规避上述几种方法中造成气象信息缺失的问题，高 精度的补全被WTC干扰的气象信号。本发明针对气象信号观测矩阵的先验知识，即较大奇 异值代表矩阵主要成分，在每次迭代运算中通过对气象信号观测矩阵核范数中的每个奇异 值赋予不同的权值，避免了用同一常数对所有的奇异值进行阈值化，更精确的获取原气象 信号的托普利兹矩阵在低秩上的最优逼近矩阵，抑制了风力涡轮机杂波及噪声信号，进一 步提高了补全精度。

为了实现以上目的，本发明采取的一种技术方案是：

本发明提供一种基于先验知识加权的气象雷达风力涡轮机杂波抑制方法，包括如下步 骤：

步骤一、气象雷达回波信号的建模，具体为：输入气象雷达回波信号，假定第i个距离 单元同时包含WTC信号，第i个距离单元第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为相干积累脉冲数；

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

在第i个距离单元两侧各取

个距离单元，其中L为距离单元数，并将第i个距离单 元中的回波信号[x_i(1)，x_i(2)，...，x_i(M)]置零，可得观测矩阵X_L×M：

由X_L×M构建出低秩随机采样Toeplitz矩阵T，其构建准则为：逐次将观测矩阵X_L×M第 m个脉冲下的向量[x₁(m)，x₂(m)，...，x_L(m)]^T构建成行为m₁、列为m₂的低秩托普利兹矩阵，其中， m₁及m₂满足m₁+m₂-1＝L，令托普利兹矩阵的第p行，第q列的元素为t_p，q，有t_p，q＝t_p+1，q-1， 且满足：

则回波信号x_i(m)在第i个距离单元置零后构建的低秩托普利兹矩阵T为：

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

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

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

WTC信号w_i(m)在第个i距离单元置零后构建的低秩托普利兹矩阵W为零矩阵；

步骤三、通过加权奇异值矩阵补全模型抑制WTC信号：

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

||S||_*，ω是气象信号Toeplitz矩阵S加权核范数：

σ_i(S)表示矩阵S 的第i个奇异值，r表示矩阵S的奇异值个数，ω＝[ω₁，ω₂，...，ω_r]为非负加权矢量满足：

其中，ε是为了避免出现奇异值为0的情况而添加的正极小值，v是调整参数满足

步骤四：利用加权非精确增广拉格朗日乘子法WIALM求解矩阵补全模型，逐脉冲输出抑 制WTC后矩阵补全后的气象信号

进一步地，步骤四中拉格朗日函数为：

其中，Y＝Y₀+μ(T-S-N-C-W)为拉格朗日乘子矩阵，Y₀为拉格朗日乘子矩阵初值， 取值为0；μ＞0表示惩罚因子，λ为正则化参数，设置为

||·||_*，ω表示加权核范 数，||·||₁表示矩阵的l₁范数，

表示取复数的实部，||·||_F表示F范数，

tr(·)表示取矩阵的迹，<X，Y>＝tr(X^HY)表示矩阵的内积。

进一步地，利用加权非精确增广拉格朗日乘子法WIALM求解的步骤为：

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

2)更新

首先利用公式(U，∑，V^H)＝svd(T-C_k-N_k -W_k+μ_k ^{- 1}Y_k)更新加权矢量ω，利用收缩运算符S_W(Σ)＝max(∑_i，i-ω_i，0)i＝1，2，...，r求解 S_k+1＝US_W(∑)V^H，其中r为矩阵(T-C_k-N_k-W_k+μ_k ^-1Y_k)奇异值数量；其中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)更新WTC矩阵W：

4)更新噪声矩阵N：

5)更新噪声矩阵C：

6)更新Y_k：Y_k+1＝Y_k+μ_k(T-S_k+1-N_k+1-C_k+1-W_k+1)；

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

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

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

9)结束循环，输出：

进一步地，利用WIALM迭代算法逐次输出补全后恢复出的各脉冲下回波信号的Toeplitz 矩阵为：

提取矩阵

的第一列与第一行，将其构成M个L×1维的列向量

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

取矩阵

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

本发明的有益效果在于：利用气象信号观测矩阵的先验知识，即较大奇异值代表矩阵 主要成分，通过对核范数中的每个奇异值赋予不同的权值，即用小值收缩大的奇异值，而 用大值收缩并滤除小的奇异值(若收缩后小于零，则置零)，以保护数据中的主要部分，忽 略不重要的或是产生噪声的部分，避免了用同一常数对所有的奇异值进行阈值化，更精确 的获取原矩阵在低秩上的最优逼近矩阵，抑制了风力涡轮机杂波及噪声信号，进一步提高 了气象信号矩阵补全的精度，运算量低，实用性强，具有很好的工程应用前景。

附图说明

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

图2为噪声干扰下IALM算法及本发明的WIALM算法恢复的气象信号功率谱对比图；

图3为IALM算法及本发明的WIALM算法恢复的气象信号幅值对比图；

图4为不同信噪比下IALM算法及本发明的WIALM算法恢复误差曲线对比图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地 描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本 发明中的实施例，本领域技术人员在没有作出创造性劳动前提下所获得的所有其他实施例， 都属于本发明保护的范围。

本发明主要研究一种基于先验知识加权的气象雷达风力涡轮机杂波抑制方法，图1是 信号处理流程，其主要步骤如下：

步骤一、气象雷达回波信号的建模，具体为：输入气象雷达回波信号，假定第i个距离 单元同时包含WTC信号，第i个距离单元第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；

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

在第i个距离单元两侧各取

个距离单元，其中L为距离单元数，可取L＝81，即在 第41个距离单元两边各取40个距离单元，并将第i个距离单元中的回波信号[x_i(1)，x_i(2)，...， x_i(M)]置零，可得观测矩阵X_L×M：

由X_L×M构建出低秩随机采样Toeplitz矩阵T，其构建准则为：逐次将观测矩阵X_L×M第 m个脉冲下的向量[x₁(m)，x₂(m)，...，x_L(m)]^T构建成行为m₁、列为m₂的低秩托普利兹矩阵，其中， m₁及m₂满足m₁+m₂-1＝L取m₁＝42，m₂＝40，令托普利兹矩阵的第p行，第q列的元素为t_p，q， 有t_p，q＝t_p+1，q-1，且满足：

则回波信号x_i(m)在第i个距离单元置零后构建的低秩托普利兹矩阵T为：

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

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

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

由于WTC信号w_i(m)仅仅存在第i个距离单元，即观测矩阵中WTC信号W_L×M：

因此在第i个距离单元置零后构建的低秩托普利兹矩阵W_T为零矩阵；

步骤三、通过加权奇异值矩阵补全模型抑制WTC信号：

根据低秩补全理论，对于不完全的数据矩阵，可依据部分矩阵元素，通过约束秩最小 化问题实现对未知元素的补全。矩阵补全模型可以表示为：

公式(2)中，T是观测矩阵，S是气象信号矩阵，N是噪声矩阵，C是地杂波矩阵， W是WTC矩阵；min表示取最小化，s.t.表示约束，rank(·)表示取矩阵的秩，P_Ω为投影算 子，表示投影到仅在指标集Ω非零的稀疏矩阵子空间上的映射，它使得矩阵在Ω中的元不 变，Ω以外的元置零，用公式表示如下：

在公式(2)中，rank(·)是一个非凸函数，所示的优化问题为NP-Hard问题，求解起来比较困难。因此可以将约束秩最小化问题松弛为约束核范数最小问题，将NP-Hard问题转换为凸优化问题，此时矩阵补全模型变为：

公式(4)中，||·||_*表示对矩阵·取核范数，即

其中r表示矩阵的奇异值 个数，σ_i(·)表示矩阵·的第i个奇异值。

传统的核范数最小化，对矩阵全部奇异值分配相等且固定的软阈值，以追求目标函数 的凸性，却忽略了低秩矩阵具有较大奇异值代表了矩阵的主要成分的先验知识。用同一常 数阈值核范数的奇异值，会导致大奇异值信息的丢失，从而使得恢复数据得到较低的峰值 信噪比，极大地限制了它处理许多实际问题(如去噪)的能力和灵活性。由于低秩矩阵奇异 值的前10％甚至1％的和通常占所有奇异值总和的99％以上，为了避免用同一常数对所有的奇 异值进行阈值化，在实践中引入了加权的思想，将公式(4)重构为：

在公式(5)中，min(·)表示最小化处理，s.t.表示约束，P_Ω表示投影到仅在指标集Ω非 零的稀疏矩阵子空间上的映射，|S||_*，ω是气象信号矩阵S加权核范数：

σ_i(S)表示矩阵S的第i个奇异值，r表示矩阵S的奇异值个数，ω＝[ω₁，ω₂，...，ω_r]为非负加 权矢量满足：

公式(6)中，ε是为了避免出现奇异值为0的情况而添加的正极小值，v是调整参数满 足

步骤四：利用加权非精确增广拉格朗日乘子法WIALM求解矩阵补全模型，逐脉冲输出抑 制WTC后矩阵补全后的气象信号

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

其中，Y＝Y₀+μ(T-S-N-C-W)为拉格朗日乘子矩阵，Y₀为拉格朗日乘子矩阵初值， 取值为0；μ＞0表示惩罚因子，λ为正则化参数，设置为

||·||_*，ω表示加权核范 数，||·||₁表示矩阵的l₁范数，

表示取复数的实部，||·||_F表示F范数，

tr(·)表示取矩阵的迹，<X，Y>＝tr(X^HY)表示矩阵的内积。

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

输入：T_j观测样本，(i，j)∈Ω，矩阵T∈R_m×n

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

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

其中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)(U，∑，V^H)＝svd(T-C_k-N_k-W_k+μ_k ^-1Y_k)求解ω＝[ω₁，ω₂，...，ω_r]，其中

4)S_k+1＝US_W(∑)V^H，其中收缩运算符S_W(Σ)＝max(∑_i，i-ω_i，0)i＝1，2，...，r，r为矩阵 (T-C_k-N_k-W_k+μ_k ^-1Y_k)奇异值数量；

5)更新WTC矩阵W：

其中

表示Ω以外的指标集；

6)更新噪声矩阵N：

7)更新噪声矩阵C：

8)更新Y_k：Y_k+1＝Y_k+μ_k(T-S_k+1-N_k+1-C_k+1-W_k+1)；

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

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

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

11)结束循环，输出：

及

在构建托普利兹矩阵之前，已经对气象信号c_i已经滤除，因此气象信号托普利兹矩阵C 及WTC信号托普利兹矩阵W在WIMLA算法迭代过程中一直部为零，因此经过WIMLA算法迭代结束后输出信号仅为

及

利用WIALM迭代算法逐次输出补全后恢复出的各脉冲下回波信号的Toeplitz矩阵为：

提取矩阵

的第一列与第一行，将其构成M个L×1维的列向量

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

取矩阵

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

下面通过MATLAB对算法进行测试，验证均矩阵补全算法恢复气象信号的有效性，雷达 系统仿真参数如表1所示；假设第1到第81个距离单元存在回波信号，共有64个脉冲，第41 个距离单元存在WTC，对其置零处理，将其构造为托普利兹矩阵形式，将补全后的矩阵重新 排列L×1维矩阵，依次对64个脉冲信号进行仿真，提取出的各脉冲下第41个向量元素构成 了WTC抑制后的回波信号。仿真参数如下表所示：

表1 托普利兹矩阵仿真参数

需要说明：本实例中对比的截断矩阵补全算法采用了IALM算法，具体过程如下：

步骤一：输入气象雷达回波信号，假定第i个距离单元同时包含WTC信号，第i个距离 单元第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)及n_i(m)分别为气象信号、地杂波信号、风力涡轮机杂波WTC信号和噪声信号，M为脉冲数，s_i(m)与n_i(m)为目标信号，记为z_i(m)＝s_i(m)+n_i(m)；

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

在第i个距离单元两侧各取

个距离单元，其中L为距离单元数，并将第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₂满足m₁+m₂-1＝L。令托普利兹矩阵的第p行，第q列的元素为t_p，q，有 t_p，q＝t_p+1，q-1，且满足：

则回波信号x_i(m)在第i个距离单元置零后构建的低秩托普利兹矩阵X_T为：

气象信号s_i(m)第i个距离单元置零后构建的低秩托普利兹矩阵S_T为：

地杂波信号c_i(m)在第i个距离单元置零后构建的低秩托普利兹矩阵C_T为：

噪声信号n_i(m)在第i个距离单元置零后构建的低秩托普利兹矩阵N_T为：

WTC信号w_i(m)在第i个距离单元置零后构建的低秩托普利兹矩阵W_T为零矩阵；

步骤三、通过截断矩阵补全模型抑制WTC信号：

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

截断式核范数||S_T||_o为气象信号托普利兹矩阵S_T降序排列的第1个到第0个奇异值相 加之和：

0对应的位置为气象信号S_L×M奇异值快速衰落的截断位置，气 象信号S_L×M为：

步骤四：利用非精确增广拉格朗日乘子法IALM求解矩阵补全模型，输出抑制WTC后目 标信号各脉冲下的托普利兹矩阵

拉格朗日函数为：

其中，Y_T＝Y_T0+μ(X_T-S_T-N_T-C_T-W_T)为拉格朗日乘子矩阵，Y_T0为拉格朗日乘子矩 阵初值，取值为0；μ为惩罚因子，||·||_F表示F范数，

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

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

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

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

2)更新

首先利用公式(U，∑，V^H)＝svd(X_T-C_Tk -N_Tk-W_Tk+μ_k ^-1Y_Tk)更新

从而求解

其中S_T(k+1)及S_Tk表示气象信号S_T的第k+1次及第k次更新，W_Tk表示WTC信号W_T的第k次更新，N_Tk表示噪声N_T的第k次更新，C_Tk表 示地杂波C_T的第k次更新，Y_Tk表示拉格朗日乘子矩阵Y_T的第k次更新，μ_k表示惩罚因子μ的 第k次更新；

3)更新W_T：

其中

表示Ω以外的指标集；

4)更新N_T：

5)更新C_T：

6)更新Y_T：Y_T(k+1)＝Y_Tk+μ_k(X_T-S_T(k+1)-N_T(k+1)-W_T(k+1)-C_T(k+1))；

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

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

||S_Tk-S_T(k-1)||_F/||S_Tk||_F≤η；

9)结束循环，输出：

及

步骤五：逐脉冲均值处理利用IALM算法恢复的目标信号的托普利兹矩阵

恢复气象信号

根据托普利兹矩阵的结构，，根据以下步骤恢复气象信号：取各托普利兹矩阵Z_T的第一 列及第一行形成M个L×1维的向量[z₁(m)，z₂(m)，...，z_L(m)]^Tm＝1，2，...，M，然后将上述M个L×1 维的向量组成一个L×M维的矩阵

取Z_L×M中的第i行向量作 为恢复出的第i个距离单元的气象信号

本对比算法中，IALM算法中截断点为气象信号S_L×M降序排列第5个奇异值处。

为了清晰的了解本发明算法对WTC以及噪声抑制效果，图2给出了信噪比(SNR)为30db 时的气象信号多普勒功率谱，图中仿真结果显示，噪声和缺陷信号重构误差引起的基底噪 声在矩阵补全抑制WTC之前波动较大；矩阵补全前后气象信号中心频率约为360Hz，接近中 心频率时，噪声干扰的影响很微弱，气象信号恢复精度更高。本发明的WIALM算法可以进一 步降低随机噪声的影响，恢复的功率谱可以在IALM算法的基础上再降低5～15dB的噪声功率， 比IALM算法更接近气象信号功率谱真值。

图3给出了气象信号幅值图，由图3可知，本发明的WIALM算法恢复的幅值比IALM算法更接近气象信号幅度真值。

为了定量分析IALM算法和本发明的WIALM算法的恢复性能，定义如下均方根误差(RMSE) 作为性能指标：

为了分析输入信号信噪比对本算法性能的影 响，进行100次蒙特卡洛实验，将不同信噪比下的信号恢复误差绘制出如图4所示的曲线 图，其中信噪比变化范围为0dB～30dB。从图中可以看出，对缺失数据进行矩阵补全的均方 根误差随着输入信号的信噪比的增加而降低，输入信号的信噪比越高，进行奇异值分解时， 噪声因素对奇异值分解的影响越小，矩阵补全恢复精度越高。由图可知本发明的WIALM算 法的RMSE相较于IALM算法降低的更快，当SNR大于7dB后，本发明的WIALM算法的RMSE 整体低于IALM算法的RMSE，故WIALM算法的恢复精度更高。

该发明针对气象信号观测矩阵的先验知识，通过对核范数中的每个奇异值赋予不同的 权值，避免了用同一常数对所有的奇异值进行阈值化，能精确的获取原矩阵在低秩上的最 优逼近矩阵。仿真实验结果表明，本发明有效提高了气象信号的恢复精度，更好的抑制风 力涡轮机杂波和噪声，适合工程应用。

以上所述仅为本申请的优选实施例而已，并不用于限制本申请，对于本领域的技术人 员来说，本申请可以有各种更改和变化。凡在本申请的精神和原则之内，所作的任何修改、 等同替换、改进等，均应包含在本申请的保护范围之内。

Claims (2)

1.一种基于先验知识加权的气象雷达风力涡轮机杂波抑制方法，其特征在于，包括如下步骤：

步骤一、气象雷达回波信号的建模，具体为：输入气象雷达回波信号，假定第i个距离单元同时包含WTC信号，第i个距离单元第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为相干积累脉冲数；

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

在第i个距离单元两侧各取

个距离单元，其中L为距离单元数，并将第i个距离单元中的回波信号[x_i(1),x_i(2),...,x_i(M)]置零，可得观测矩阵X_L×M：

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

则回波信号x_i(m)在第i个距离单元置零后构建的低秩托普利兹矩阵T为:

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

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

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

WTC信号w_i(m)在第个i距离单元置零后构建的低秩托普利兹矩阵W为零矩阵；

步骤三、通过加权奇异值矩阵补全模型抑制WTC信号：

s.t.P_Ω(T)＝P_Ω(S+N+W+C)；

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

||S||_*,ω是气象信号Toeplitz矩阵S加权核范数：

σ_i(S)表示矩阵S的第i个奇异值，r表示矩阵S的奇异值个数，ω＝[ω₁,ω₂,…,ω_r]为非负加权矢量满足：

其中，ε是为了避免出现奇异值为0的情况而添加的正极小值，v是调整参数满足

步骤四：利用加权非精确增广拉格朗日乘子法WIALM求解矩阵补全模型，逐脉冲输出抑制WTC后矩阵补全后的气象信号

步骤四中拉格朗日函数为：

其中，Y＝Y₀+μ(T-S-N-C-W)为拉格朗日乘子矩阵，Y₀为拉格朗日乘子矩阵初值，取值为0；μ＞0表示惩罚因子，λ为正则化参数，设置为

||·||_*,ω表示加权核范数，||·||₁表示矩阵的l₁范数，

表示取复数的实部，||·||_F表示F范数，

tr(·)表示取矩阵的迹，＜X,Y＞＝tr(X^HY)表示矩阵的内积；

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

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

2)更新

首先利用公式(U,Σ,V^H)＝svd(T-C_k-N_k-W_k+μ_k ^-1Y_k)更新加权矢量ω，利用收缩运算符S_W(∑)＝max(∑_i,i-ω_i,0)i＝1,2,…,r求解S_k+1＝US_W(∑)V^H，其中r为矩阵(T-C_k-N_k-W_k+μ_k ^-1Y_k)奇异值数量；其中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)更新WTC矩阵W：

4)更新噪声矩阵N：

5)更新地杂波矩阵C：

6)更新Y_k：Y_k+1＝Y_k+μ_k(T-S_k+1-N_k+1-W_k+1-C_k+1)；

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

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

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

9)结束循环，输出：

2.根据权利要求1所述的基于先验知识加权的气象雷达风力涡轮机杂波抑制方法，其特征在于，利用WIALM迭代算法逐次输出补全后恢复出的各脉冲下回波信号的Toeplitz矩阵为：

提取矩阵

的第一列与第一行，将其构成M个L×1维的列向量

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

取矩阵

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

Images