一种基于矩阵补全算法的五叶片风力发电机的气象雷达杂波抑制方法

1.本发明属于通信领域，具体涉及一种基于矩阵补全算法的五叶片风力发电机的气象雷达杂波抑制方法。


背景技术：

2.近年来，通过详细分析了气象雷达不同工作模式下风力发电机杂波与气象回波的时、频域分布特性，提出了多重二次插值、高维拟合插值、样条插值和统计不变插值、距离—多普勒谱回归、递归稀疏重构等风力涡轮机杂波抑制算法，但由于大气目标的多样性与复杂性，插值与回归模型的参数设置如今尚缺乏定量研究。而且受到风电场规模、风机转速、气象雷达工作模式等实际条件的限制，上述算法均无法同时兼顾风力涡轮机杂波抑制与气象信息无损恢复。
3.现今用于气象雷达杂波抑制的方法主要是从频域出发，对杂波进行抑制，如mtd、 mti等。上述气象雷达杂波抑制算法都是基于杂波信号与气象回波在多普勒域或二维距离—多普勒域可分离条件，采用滤波或分离算法对杂波进行有效抑制。但风力发电机杂波由于大型叶轮高速旋转，其多普勒频谱展宽极大，发电机杂波与大气目标在时、频域均严重混叠，从而气象回波淹没在风力发电机杂波中，导致时、频滤波或分离算法均无法在抑制杂波信号的同时实现对气象信号的无损保留。因此上述滤波算法等方法均失效。


技术实现要素：

4.发明目的：为解决现有滤波算法均无法在抑制杂波信号的同时实现对气象信号的无损保留的问题，本发明提出了一种基于矩阵补全算法的五叶片风力发电机的气象雷达杂波抑制方法，将稀疏优化理论引入气象雷达风力发电机杂波抑制中，重点研究小型五叶片风轮机基于矩阵补全理论的气象回波二维高精度稀疏恢复方法，该方法针对mc理论，构建满足随机采样的气象回波数据低秩恢复矩阵，并设计相应的高效迭代求解算法。通过该算法实现对气象雷达五叶片风力发电机杂波的抑制，精确恢复出要检测的气象信号。
5.技术方案：一种基于矩阵补全算法的风力发电机的气象雷达杂波抑制方法，包括以下步骤：
6.步骤1：构建气象雷达回波信号的数学模型：x(t)＝s(t)+c(t)+w(t)+n(t)，其中，t 代表时间，s(t)表示气象目标回波信号，c(t)表示地杂波信号，w(t)表示风力发电机杂波信号，n(t)表示噪声信号，x(t)表示雷达回波信号；分别对s(t)、c(t)、w(t)、n(t)矩阵化，得s、c、w、n，并得到快拍矩阵形式的雷达回波信号矩阵：x＝s+c+w+n，雷达回波信号矩阵中的每个元素均为复数值；
7.步骤2：利用地杂波抑制方法将回波信号矩阵中的地杂波信号滤除，得滤除地杂波信号后的雷达回波信号矩阵x＝s+w+n；
8.步骤3：基于矩阵补全模型依据矩阵元素，利用非精确增广拉格朗日乘子法解矩阵
补全模型的最优化问题，恢复气象数据。
9.进一步的，步骤3中，所述的非精确增广拉格朗日乘子法对应的拉格朗日函数l(s,w,n,y,μ)，表示为：
[0010][0011]
其中，y＝y0+μ(x-s-n-w)表示拉格朗日乘子，初值令其为0，μ＞0表示惩罚因子，||
·
||表示核范数，r《
·
》表示取复数实部，《x,y》＝tr(xhy)表示矩阵的内积。
[0012]
进一步的，所述的非精确增广拉格朗日乘子法，包括以下步骤：
[0013]
s1：初始化参数：拉格朗日乘子初值y0＝0，风力发电机杂波信号初值w0＝0，噪声信号n＝0，ρ＞1为一常数，迭代次数k＝0，噪声系数η；
[0014]
s2：根据下式，求解得到s
k+1
＝argminl(sk,wk,nk,yk,μk)；
[0015][0016][0017]
其中，svd表示矩阵的奇异值分解；
[0018]
s3：更新wtc矩阵：其中，ω为采样元素的集合，为正交投影算子，表示投影到仅在指标集ω非零的稀疏矩阵子空间上的映射；
[0019]
s4：根据下式，求解得到n＝argminl(s
k+1
,wk,n,yk,μk)；
[0020][0021]
s5：y
k+1
＝yk+μk(x-s
k+1-n
k+1-w
k+1
)；
[0022]
s6：μ
k+1
＝ρμk；
[0023]
s7：判断||s
k-s
k-1
||f/||sk||f≤η是否成立，若不成立，则令k
←
k+1，转到s2，否则结束循环，输出
[0024]
有益效果：本发明通过矩阵补全，能够利用空时信息高精度恢复气象信号，而多重二次插值仅利用了时域信息，造成了先验知识的浪费，降低了恢复能力。因此，本发明与多重二次插值相比，本发明的矩阵补全模型的选取更加简单方便，且恢复出气象回波不容易存在分辨率的损失。
附图说明
[0025]
图1为风轮机与雷达的几何关系图；
[0026]
图2为气象信号自身的功率谱、受噪声干扰的功率谱、同时受噪声和瞬态干扰污染的功率谱以及通过低秩快拍矩阵恢复后的功率谱的对比图。
具体实施方式
[0027]
现结合附图和实施例对本发明的技术方案做进一步说明。
[0028]
本发明的一种基于矩阵补全算法的风力发电机的气象雷达杂波抑制方法，包括以
下步骤：
[0029]
步骤1：建立气象雷达回波信号的数学模型，对气象信号进行分析。具体为：
[0030]
在气象雷达系统中，假设气象雷达回波由四部分构成，包括：目标回波信号、地杂波信号、风力发电机杂波信号以及噪声信号。因此,可理想化的将气象雷达回波描述为如下形式：x(t)＝s(t)+c(t)+w(t)+n(t)；其中，t代表时间，s(t)表示气象目标回波信号， c(t)表示地杂波信号，w(t)表示风力发电机杂波信号，n(t)表示噪声信号，x(t)表示气象雷达回波信号。
[0031]
当气象雷达发射功率为p
t
，雷达到气象目标的距离为r时，给定的天线增益为g，那么入射能流密度si为：
[0032][0033]
设目标的后向散射截面积rcs为σ，那么从目标返回天线的散射能流密度为：
[0034][0035]
设天线的有效截面积为其中λ表示光的波长，则天线所接收的功率，即气象目标的雷达方程为：
[0036][0037]
在理想情况下，默认散射粒子只进行一次散射，那么对第l个距离单元内所有散射粒子回波进行相干叠加，即为气象信号回波。第l个距离环(即第l个距离单元内所有散射粒子)，第m个脉冲采样可以表示为：
[0038][0039]
其中，u为第l个距离单元内散射粒子的总数，ao为第o个散射粒子的幅度，由雷达方程决定ω
t
为时域角频率，其中，fd为多普勒频率， fr为脉冲重复频率，在已知径向速度前提下，多普勒频率近似默认为恒定值。
[0040]
在不影响分析的前提下，提取回波信号中的高频因子，得到目标气象信号第m个脉冲采样下的回波表达式如下：
[0041][0042]
上述公式成立的条件为：有l个距离单元以及m个脉冲采样。
[0043]
步骤1中的初值设定为：雷达载频fe＝5.5
×
109hz，发射功p
t
＝1250mhz，脉冲数 m＝128，距离单元数l＝100。
[0044]
步骤2：建立地杂波分布模型。具体为：
[0045]
雷达系统在运行时发射电磁波束，遇到障碍目标后会形成回波，即雷达杂波。除了
少数固定的目标外，大多数地物杂波都有着极为复杂的建模过程，且都是随机分布的。因此，往往需要采用一些合理的数学模型来对地杂波的分布特性进行分析。因此，本发明采用由幅度概率密度以及杂波序列的相关函数来描述杂波的统计模型。
[0046]
为了更好的构造杂波序列的相关函数，可选用莱斯分布模型作为杂波分析的选择，其概率密度函数为：
[0047][0048]
其中，i0(
·
)为修正的第一类贝索函数，ε为地杂波的标准差，若将地杂波视为主信号与服从瑞利分布的多径信号分量的和，则f为地杂波主信号幅度的峰值，r为服从瑞利分布的多径信号，其表达式为：
[0049]
r＝xicosωct+jyisinωct xi,yi∈n(0,σ2)
ꢀꢀ
(7)
[0050]
其中，ωc为杂波载频，t为时间，xi和yi分别代表实部和虚部的系数，ri为地杂波瑞利分布分量的包络振幅：
[0051][0052]
又有杂波的标准差其中表示雷达的工作波长，εv表示风速的速度标准差，所以：
[0053][0054]
步骤2中的初值设定为：雷达高度h＝1000m，风机转速f
rot
＝12r/min，雷达-轴心间距d＝36km，叶片长度j＝40m，叶片个数n＝5，α＝90
°
，β＝90
°
，φ＝60
°
。
[0055]
步骤3：基于快拍矩阵补全，对风力发电机杂波进行抑制。具体为：
[0056]
随着风力发电的快速发展，风力发电机的装机量急速攀升，由风电叶片高速旋转所引起的风力涡轮机杂波抑制已经成为了杂波抑制的关键性问题。
[0057]
利用散射点叠加理论，对风力涡轮机杂波建模，图1为风轮机与雷达的几何关系图。
[0058]
多数条件下风力涡轮机的桅杆和轮机舱部位不会运动，它们生成的回波较为稳定，其数据很容易获得，因此只需研究风轮机叶片的回波模型。
[0059]
将叶片轴心设置成中心o，把同旋转面垂直的方位设成如图1所示的空间坐标系。图中，р为风机叶片上的任意点，点b为气象雷达，α、β分别表示雷达的方位角与俯仰角，φ为雷达的探测视线los与风机叶片的夹角，li为叶片上任一点p到o的距离，r 为气象雷达与轴心o的径向距离。
[0060]
首先对风机叶片与桅杆划分点目标，以五个叶片的风力涡轮机为例，假设每个风机叶片长度为j，桅杆的高度为h，设风机叶片、桅杆点目标的点间距为h米，得到风机叶片与桅杆的点目标数p_b与p_t：
[0061][0062]
风机叶片上任意一点р到雷达的距离为r
p
，理想情况下(li/r)2→
0，风机叶片上任
一点到雷达的距离为r
l
，由图1可知：
[0063][0064]
其中，为风机叶片和雷达波束的夹角。
[0065]
设θ
t
为叶片与y轴的夹角随时间的变化量，满足：
[0066][0067][0068]
其中，f
rot
表示风机转速，θ0表示叶片起始夹角。
[0069]
桅杆上任一点到雷达的距离为r2，由图1可知：
[0070][0071]
其中，ki、kj为桅杆上任意一点到中心o的距离。
[0072]
风力涡轮机杂波的信号为桅杆与叶片的信号叠加：
[0073][0074]
式中，rcs1与rcs2分别是叶片与桅杆的横截面积。
[0075]
已知雷达回波信号的观测矩阵为距离-脉冲二维矩阵，列为脉冲数，行为距离单元个数，分别用m，l表示：
[0076][0077]
式中，x为雷达回波信号的距离-脉冲二位矩阵中的元素，列为脉冲数，行为距离单元个数，分别用m、l表示
[0078]
矩阵补全不仅对矩阵性质有要求，还对观测矩阵中数据元素的采样方式存在一定限制，即矩阵补全的两个强非相关的约束条件为：
[0079]
(1)矩阵补全中的原始矩阵满足低秩(或者近似低秩)条件；
[0080]
(2)矩阵补全中观测矩阵的元素采样集必须满足均匀随机采样条件。
[0081]
为了满足以上条件，同时根据气象信号导向矢量的定义，可以发现信号在同一脉冲下矢量信号满足等比特性，因此针对上式的观测矩阵，依次提取各脉冲下的信号矢量，根据信号的等比特点，分别将其矩阵化，重构出m个距离-脉冲二维域的低秩观测矩阵，矩阵第p行，第q列的元素用k
p,q
表示，即x
m1xm2
＝(k
p,q
)
m1xm2
：
[0082][0083]
气象信号的单脉冲下的矢量本是一个l
×
1维的向量，在上式中，将其重新写为一个等效的m1×
m2维快拍矩阵，满足：
[0084]
l＝m1×
m2ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(16)
[0085][0086]
因为气象回波信号的表达式为：x(t)＝s(t)+c(t)+w(t)+n(t)，矩阵的秩一般远远小于气象雷达的脉冲数m，因此，总能重构出一个快拍矩阵x，使其满足 rank(x)＜＜min(m1,m2)。
[0087]
首先，分别对s(t)，c(t)，w(t)，n(t)矩阵化，得s,c,w,n，下式的回波信号即表示为快拍矩阵形式：x＝s+c+w+n。
[0088]
其次，对回波信号进行预处理，利用现有的杂波抑制方法将地杂波滤除，消除地杂波对矩阵补全的影响，得x＝s+w+n。
[0089]
因为上述气象信号矩阵元素均为复数值，为此需要将非精确增广的拉格朗日乘子法 (ialm)分别应用到复矩阵中各个单元的实部和虚部，然后再将两者重新合并为复信号，对应的拉格朗日函数l(s,w,n,y,μ)，表示为：
[0090][0091]
其中，y＝y0+μ(x-s-n-w)表示拉格朗日乘子，初值令其为0，μ＞0表示惩罚因子，||
·
||表示核范数，r《
·
》表示取复数实部，《x,y》＝tr(xhy)表示矩阵的内积。
[0092]
在许多实际问题中，如何将稀疏向量空间过渡到低秩矩阵空间是压缩感知衍生出矩阵补全一个关键问题。对于某一矩阵，假设只采样获取矩阵的部分样本元素，而另一部分元素丢失后无法获取，根据压缩感知理论，可以通过最小秩约束优化模型，合理准确地填充这些缺失的元素：minrank(x)s.t.p
ω
(m)＝p
ω
(x)；
[0093]
其中，ω为采样元素的集合，p
ω
为正交投影算子，表示投影到仅在指标集ω非零的稀疏矩阵子空间上的映射，当(i，j)∈ω时，m
ij
表示采样元素。
[0094]
一般情况下，气象雷达的接收信号x(t)同时包含了气象信号s(t)、杂波信号c(t)、噪声信号n(t)，在脉冲-多普勒二维域，将其构造为二维复合矩阵。由于wtc的外界干扰，需要将wtc置零去除，因而只观测到该矩阵的d个元素，且d远小于数据矩阵总元素个数m
×
n。当回波矩阵x满足低秩条件时，则有可能从d个已知元素中恢复矩阵x的全部元素。
[0095]
已知指标集中的元素s
ij
∈r，其中ω为x 的d个已采样元素的集合，可通过求解如下秩最小化问题实现对未知元素的补全： minrank(s)s.t.p
ω
(s)＝p
ω
(x)；其中x和s都是m
×
n维的矩阵，p
ω
表示为：
[0096][0097]
然而，由于矩阵具有非连续、非凸的秩函数，因而直接求解秩最小化问题比较困难。为了得到一个更加容易计算的凸优化问题，用凸包核范数‖x‖.来表示上中的函数 rank(x)，得到了最小化核范数优化模型：
[0098]
min||s||s.t.p
ω
(s)＝p
ω
(x)
[0099]
上式所示核范数凸优化问题，通过引入m
×
n维实矩阵变量n，增广拉格朗日乘子法高效的解决该核范数凸优化问题，将上式转化为对如下优化问题的求解：
[0100]
min||s||s.t.x＝s+n p
ω
(x)＝0
[0101]
凸优化问题利用ialm算法在更新稀疏噪声分量的约束p
ω
(n)＝0时最小化 l(s,w,n,y,μ)，其算法步骤如下所示：
[0102]
s1：输入观测样本x
ij
，(i,j)∈ω，惩罚因子μ0，噪声系数η；
[0103]
s2：初始化拉格朗日乘子初值y0＝0，风力发电机杂波信号初值w0＝0，噪声信号 n＝0，ρ＞1为一常数，迭代次数k＝0；
[0104]
s3：当未收敛时，从s4和s5中解出s
k+1
＝argminl(sk,wk,nk,yk,μk)；
[0105]
s4：svd表示矩阵的奇异值分解，奇异值分解理论
[0106][0107]
其中，∑r＝diag(γ1,γ2,γ3,
…
,γr)，γr为矩阵x的全部奇异值，u和v包含相应的左右奇异向量；
[0108]
s5：
[0109]
s6：更新wtc矩阵：其中，ω为采样元素的集合，为正交投影算子，表示投影到仅在指标集ω非零的稀疏矩阵子空间上的映射；
[0110]
s7：s8解出n＝argminl(s
k+1
,wk,n,yk,μk)；
[0111]
s8：
[0112]
s9：更新yk：y
k+1
＝yk+μk(x-s
k+1-n
k+1-w
k+1
)；
[0113]
s10：更新μk至μ
k+1
：μ
k+1
＝ρμk；
[0114]
s11：若下式不成立，则算法未收敛，令k
←
k+1，转到s3，否则转到s12，；
[0115]
s12：结束循环，输出
[0116]
步骤3中的初值设定为：η＝10-3
，μ0＞0。
[0117]
现对64个脉冲下的阵列回波进行仿真分析。假设由100个距离单元构成了等距线阵，检测出wtc存在于第30个距离单元中，对其置零处理后，以快拍形式重构补全矩阵。为了分析不同信噪比下的信号恢复能力，将信噪比设置为0～35db，考虑到测试的准确性，每个信噪比下均做100次蒙特卡罗实验。
[0118]
图2展示了气象信号自身的功率谱、受噪声干扰的功率谱、同时受噪声和瞬态干扰污染的功率谱以及通过低秩快拍矩阵恢复后的功率谱的对比，从图2中可以看出，受噪声及瞬态干扰污染叠加干扰的气象信号功率谱已失真，气象信号受噪声干扰的功率谱仅在350hz附近一小段频率段与原信号功率谱发生重叠，而采用低秩快拍矩阵恢复出来的信号的功率谱在很大程度上与原信号功率谱发生重叠，可见低秩快拍矩阵补全在一定程度上能够较好的抑制噪声的干扰，尽可能的恢复气象信号。
[0119]
在本发明中，使用低秩快拍矩阵补全算法代替工程中广泛应用的多重二次插值算法来抑制涡轮风机产生的杂波，有利于提高信号的信噪比，在恢复气象信号的同时，可以更有效地抑制噪声，提高恢复精度。与多重二次插值相比，矩阵补全模型的选取更加简单方便，且恢复出气象回波不容易存在分辨率的损失。此外，矩阵补全能够利用空时信息高精度恢复气象信号，而多重二次插值仅利用了时域信息，造成了先验知识的浪费，降低了恢复能力，因此本发明所提算法优于传统的多重二次插值算法。