CN112379380B - 基于均值法再处理截断矩阵补全的风电场杂波抑制方法 - Google Patents

Abstract

本发明公开了一种基于均值法再处理截断矩阵补全的风电场杂波抑制方法，利用观测矩阵的先验知识，首先将被风力涡轮机杂波(WTC)污染的距离单元置零，并逐脉冲将距离向量重构为随机采样的低秩托普利兹矩阵，然后针对气象信号S_L×M奇异值序列呈快速衰减的特点，在奇异值收缩过程中选取合适的截断点进行截断，最后利用矩阵补全恢复缺失元素并对补全出的目标信号的托普利兹矩阵进行均值处理，从而输出气象信号。仿真结果表明，本发明能同时抑制WTC和噪声干扰，高精度恢复气象信号，具有很好的工程应用前景。

Description

基于均值法再处理截断矩阵补全的风电场杂波抑制方法

技术领域

本发明涉及气象雷达杂波抑制领域，尤其涉及一种基于均值法再处理截断矩阵补全的风电场杂波抑制方法。

背景技术

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

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

发明内容

本发明所要解决的技术问题是提供一种基于均值法再处理截断矩阵补全的风电场杂波抑制方法，将均值法再处理截断矩阵理论引入气象雷达WTC抑制，研究基于矩阵补全的气象雷达小型风电场杂波抑制方法，利用在观测矩阵中气象信号S_L×M奇异值收缩过程中选取合适的截断点进行截断，利用矩阵补全恢复缺失元素并对补全出的目标信号的托普利兹矩阵进行均值处理，从而输出气象信号，能够规避上述几种方法中造成气象信息缺失，高精度矩阵补全被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)及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)结束循环，输出：

及

进一步地，根据托普利兹矩阵的结构，逐脉冲提取出目标信号的托普利兹矩阵

的第l条对角线的元素diag(Z_Tk，l)，对其均值处理： a_l＝mean(diag(Z_Tk，l))，l∈{1，..i..，m₁+m₂-1}，得到M个均值处理后的托普利兹矩阵

取各托普利兹矩阵

的第一列及第一行形成M个L×1维的向量

然后将上述M个L×1维的向量组成一个L×M维矩阵

取矩阵

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

本发明的有益效果在于：利用在观测矩阵中气象信号S_L×M奇异值收缩过程中选取合适的截断点进行截断，利用矩阵补全恢复缺失元素并对补全出的目标信号的托普利兹矩阵进行均值处理，从而输出气象信号，抑制了风力涡轮机杂波及噪声信号，提高了气象信号矩阵补全的精度，运算量低，实用性强，具有很好的工程应用前景。

附图说明

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

图2为选取气象信号S_L×M不同数量奇异值后恢复的气象信号功率谱对比图；

图3为噪声干扰下截断矩阵补全算法及本发明的均值法再处理截断矩阵补全算法恢复的气象信号功率谱对比图；

图4为传统矩阵补全算法、截断矩阵补全算法及本发明的均值法再处理截断矩阵补全算法恢复的气象信号速度值对比图；

图5为传统矩阵补全算法、截断矩阵补全算法及本发明的均值法再处理截断矩阵补全算法恢复的气象信号误差对比图；

图6为传统矩阵补全算法(NNM算法)、截断矩阵补全算法(TNNR算法)及本发明的均值法再处理截断矩阵补全算法(TNNR算法)恢复的气象数据信噪比对比图。

具体实施方式

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

本发明主要研究一种基于均值法再处理截断矩阵补全的风电场杂波抑制方法，图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)及n_i(m)分别为气象信号、地杂波信号、风力涡轮机杂波WTC信号和噪声信号，M为脉冲数，s_i(m)与n_i(m)为目标信号，记为z_i(m)＝s_i(m)+n_i(m)，可取M＝64；

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

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

个距离单元，其中L为距离单元数，可取L＝65，即在第33个距离单元两边各取32个距离单元，并将第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，m₁＝m₂＝33。令托普利兹矩阵的第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个距离单元，即观测矩阵中WTC信号W_L*M：

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

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

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

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

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

本实施例中选取截断位置为气象信号S_L×M第6个奇异值处，即0＝6，

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

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

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

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

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

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

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

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

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

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

3)(U，∑，V^H)＝svd(X_T-C_Tk-N_Tk-W_Tk+μ_k ^-1Y_Tk)，对奇异值在O处截断得到截断后的奇异值矩阵

4)

5)更新W_T：

其中

表示Ω以外的指标集；

6)更新N_T：

7)更新C_T：

8)更新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))；

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

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

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

11)结束循环，输出：

及

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

及

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

恢复气象信号

根据托普利兹矩阵的结构，逐脉冲提取出目标信号的托普利兹矩阵

的第l 条对角线的元素diag(Z_Tk，l)，对其均值处理：a_l＝mean(diag(Z_Tk，l))，l∈{1，..i..，m₁+m₂-1}，得到M个均值处理后的托普利兹矩阵

其中

当l取1时，此时R_l第m₁*1位置元素为1，其他均为0，对应到对目标信号中s₁(m)+n₁(m)元素求取均值，当l取2时，此时R_l第(m₁-1)*1、m₁*2位置元素为1，其他均为0，对应到对目标信号中s₂(m)+n₂(m)元素求取均值，以此类推，当l取m₁时，此时R_l第i*i(i＝1，2，...，m₁)位置元素为1，其他均为0，对应到目标信号中

元素求取均值，当l取m₁+m₂-1时，此时R_l第1*m₂位置元素为1，其他均为0，对应到目标信号中s_L(m)+n_L(m)元素求取均值，从而得到单个脉冲下均值处理后的拓普利兹矩阵

取各托普利兹矩阵

的第一列及第一行形成M个L×1维的向量

然后将上述M个L×1维的向量组成一个L×M维矩阵

取矩阵

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

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

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

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

需要说明：本实施例中对比的截断矩阵补全算法是没有均值预处理目标信号的托普利兹矩阵

在本发明的步骤三后，根据以下步骤恢复气象信号：取各托普利兹矩阵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个距离单元的气象信号

如图2显示了根据观测矩阵中气象信号S_L×M不同个数奇异值恢复出气象信号的图像，假设观测矩阵中气象信号矩阵S_L×M的秩为R，当选取奇异值个数远远小于R时，例如前1％奇异值数量进行重构信号时，此时恢复出的气象信号图像远离气象信号的原始图像，当选取前10％奇异值进行重构信号时，此时恢复出的气象信号图像离近气象信号的原始图像越来越近，由于气象信号矩阵S_L×M奇异值在O处呈现快速衰落的现象即前O个奇异值包含了气象的大部分信息，随着选取的奇异值个数增加到O时，重构后的图像与原始图越来越接近，因此可知基于截断核范数的矩阵补全算法可以得到待补全矩阵在秩为O上的最优的近似。

为了清晰的了解本发明算法对WTC以及噪声抑制效果，图3给出了噪声干扰下截断矩阵补全算法及本发明的均值法再处理截断矩阵补全算法恢复的气象信号功率谱对比图，此时截断点为6，图中仿真结果显示，截断矩阵算补全算法能够抑制WTC，但在在抑制噪声信号方面效果不太理想；本发明均值法再处理截断矩阵补全算法在信噪比SNR＝30dB的情况下，多普勒频率在300-420Hz之间时，噪声干扰较弱，气象信号恢复精度较高，峰值旁瓣对高干扰噪声有较明显的抑制作用，在300HZ左右噪声功率降低约20dB，旁瓣处噪声的波动起伏降低，提高了信噪比，对气象信号恢复精度更高，对噪声抑制效果更好。

图4为传统矩阵补全算法、截断矩阵补全算法及本发明的均值法再处理截断矩阵补全算法恢复的气象信号速度值，恢复后的气象信号速度为：

其中∠为取相位，

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

为恢复后的气象信号的共轭转置。径向速度估计在低信噪比下波动较为明显，与真值存在较大偏差，但随着信噪比的增加，径向速度估计的误差逐渐减小，最终收敛到真值。对比三种算法不难发现，在相同性噪比的情况下，本发明提出的均值法再处理截断矩阵补全算法恢复更平稳，与原值相差更小，收敛性更高，性能更佳。

图5为传统矩阵补全算法、截断矩阵补全算法及本发明的均值法再处理截断矩阵补全算法恢复的气象信号误差对比图，气象信号的估计值与真实值的均方根误差定义为

从100个独立蒙特卡罗实验的平均结果可以看出，由图5可知随着输入信号信噪比的提高，传统矩阵补全算法、截断矩阵补全算法及本发明的均值法冉处理截断矩阵补全算法恢复的气象信号误差均明显减小在相同信噪比时，本发明的均值法冉处理截断矩阵补全算法恢复的气象信号误差更小，性能更优。

图6为传统矩阵补全算法(NNM算法)、截断矩阵补全算法(TNNR算法)及本发明的均值法冉处理截断矩阵补全算法(TNNR算法)恢复的气象数据信噪比对比图，从图中可以三种算法在恢复缺失数据时输入信噪比均小于输出信噪比，即恢复气象信号的信号子空间中的噪声扰动量都随之降低。对比图中三种算法的性能曲线可知，本发明的均值法冉处理截断矩阵补全算法(TNNR算法)降噪效果明显优于截断矩阵补全算法(TNNR算法)及NNM算法，且当输入信噪比在10～15dB时，信号子空间与噪声子空间中噪声扰动量的比例随之变化，噪声子空间中比例增大，降噪性能更明显。

本发明首先该发明利用气象信号空间相关性，首先检测出回波信号中同时含有气象信号和风电场杂波的距离单元，将其置零，然后依次将单个脉冲下的向量重构为满足零元素随机分布的低秩托普利兹矩阵，利用在观测矩阵中气象信号S_L×M奇异值收缩过程中选取合适的截断点进行截断，利用IALM迭代算法恢复缺失元素并对补全出的目标信号的托普利兹矩阵进行均值处理，从而输出气象信号，抑制了风力涡轮机杂波及噪声信号，提高了气象信号矩阵补全的精度，运算量低，实用性强，具有很好的工程应用前景。

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

Claims (1)

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)及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个到第o个奇异值相加之和：

o对应的位置为气象信号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)结束循环，输出：

及

根据托普利兹矩阵的结构，逐脉冲提取出目标信号的托普利兹矩阵

的第l条对角线的元素diag(Z_Tk,l)，对其均值处理：a_l＝mean(diag(Z_Tk,l)),l∈{1,..i..,m₁+m₂-1}，得到M个均值处理后的托普利兹矩阵

取各托普利兹矩阵

的第一列及第一行形成M个L×1维的向量

然后将上述M个L×1维的向量组成一个L×M维矩阵

取矩阵

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

Images