近些年来，世界各国对风电产业的投资日益扩大，全球风力发电累计装机容量呈指数上升趋势^[1]。然而，风电场会对空中交通管制系统、防空预警监视系统、导航系统、气象观测系统等雷达设备产生严重的影响^[2-6]。风轮机杂波(Wind Tubrine Clutter, WTC)不仅会造成虚警问题，还可能掩盖目标信号，造成雷达系统的探测盲区。另外，风轮机杂波还会导致风电场附近目标检测概率的下降，因此风轮机杂波被当作一种杂波。现阶段空管雷达多为低分辨率两坐标雷达，难以通过高度信息或距离信息对目标和风轮机杂波进行区分，因此风电场上方或附近的目标杂波均会受到WTC的影响，而空管雷达需要通过对飞行目标持续检测以实现对目标的连续跟踪和引航，因此对消除杂波和噪声的要求非常高，任何虚假目标都可造成航迹的不连续甚至偏航、误航。

近年来，相关科研组织和学者已经提出了一些风电场杂波抑制方案。这些方案措施可分为两大类：一类是减少风轮机杂波进入雷达系统，另一类是对进入雷达系统中的WTC进行抑制处理。文献[7]提出加涂吸波材料减小风轮机的雷达散射截面(RCS)以达到抑制WTC产生的效果，但该措施不仅增大了风轮机叶片重量，降低了电能转化率，而且提高了成本，因此无法大面积推广。文献[8]提出通过部署补盲雷达，利用多部雷达的信息融合来抑制风轮机杂波的干扰。文献[9]提出通过优化风轮机的部署来减小WTC对雷达系统的影响，但是该方法依赖于地形因素且不适用于已建成的风电场。上述方法都是通过相关措施减少WTC进入雷达系统以达到杂波抑制效果，难以有效解决已部署风电场的干扰问题。而在实际应用中常利用信号处理手段对进入雷达系统中的WTC进行处理。由于风轮机叶片在频域上产生连续变化的多普勒频率，存在频谱展宽效应^[10]，导致WTC分散到多个非零频的多普勒滤波器组，因此传统的动目标显示(MTI)、动目标检测(MTD)效果显然很差。文献[11]利用谱估计和补偿的办法，通过剔除WTC所在距离单元以达到抑制的目的，但是该方法只适用于目标和WTC不在同一个距离单元，适用范围窄。Uysal团队^[12-13]提出了一种基于形态学成分分析(MCA)的抑制方法，在不同原子基上对风轮机和目标信号进行重构和变换，根据各成分间的稀疏性差异实现了对WTC的分离。文献[14]提出对风轮机的杂波进行重构，通过把重构数据作为滤波器将回波信号中的WTC滤除。但是，其重构过程为了保证精确度需要生成数目庞大的原子字典，运算量大、字典冗余度高且重构精度相对较低。文献[15]在目标参数估计领域提出通过对字典库进行分解，构建规模更小的动态字典库来解决原子数过多带来的计算和存储负担问题。

因此，本文为解决WTC引起的目标遮蔽和虚假目标问题，对进入雷达系统的WTC提出了一种动态化稀疏重构抑制的方法。该方法首先根据WTC时频特征对风轮机叶片转速和初始相位进行粗估计，以缩小字典构建参数范围，再利用正交匹配追踪(OMP)算法在字典中匹配最佳原子来逐级重构，最后通过滤除WTC以达到抑制的目的。理论分析和仿真结果均表明，该方法可实现对杂波信号的精准重构，达到对WTC的有效抑制。同时，针对实际应用场景，仿真实验表明，该稀疏重构方法对WTC也有明显的抑制效果，具有一定的工程价值。

1 WTC建模分析

由于风轮机桅杆、引擎舱等产生的静态杂波可通过MTI滤波器滤除，因此本文只针对由风轮机叶片产生的动态杂波进行研究。以风轮机旋转中心为原点O，建立如图 1坐标系。其中：P为雷达站位置；r_{n, i}为第n个叶片上第i个散射点的雷达视线(LOS)距离；R为雷达中心距离风轮机旋转中心距离；l_i为风轮机上散射点距离叶片旋转中心O的长度，l_i≪R；θ为叶片旋转角；α为雷达方位角；φ为叶片与雷达视线间的夹角。

由于l_i≪R，则可近似为^[16]

$ \begin{array}{l} {r_{n,i}}\left( t \right) = {\left( {{R^2} + l_i^2 - 2R{l_i}\cos \varphi \left( t \right)} \right)^{1/2}} \approx \\ \;\;\;\;\;R - {l_i}\cos \varphi \left( t \right) \end{array} $	（1）

根据空间几何关系可知

$ \begin{array}{l} \cos \varphi (t) = \sin \alpha \cos \theta (t) = \\ \;\;\;\;\;\;\sin \alpha \cos \left( {{\theta _0} + 2{\rm{ \mathit{ π} }}{f_{{\rm{rot}}}}t + \frac{{2{\rm{ \mathit{ π} }}}}{N}\left( {n - 1} \right)} \right) \end{array} $	（2）

式中：θ₀为该叶片初始旋转角；f_rot为叶片转速, 单位为r/s。

假设雷达为线性调频号体制，其发射信号表达式为^[17]

$ \begin{array}{l} {s_{\rm{t}}}\left( {\hat t,{t_m}} \right) = A{\rm rect}\left( {\frac{{\hat t}}{{{T_{\rm{P}}}}}} \right)\exp \left( {{\rm{j}}2{\rm{ \mathit{ π} }}\left( {{f_{\rm{c}}}t + \frac{1}{2}k{{\hat t}^2}} \right)} \right)\\ \;\;\;\;\; - \frac{{{T_{\rm{P}}}}}{2} < \hat t < \frac{{{T_{\rm{P}}}}}{2} \end{array} $	（3）

式中：f_c为载频；T_P为脉冲宽度；A为天线增益；k为调频率；$\hat{t} $为快时间；t_m为慢时间，t_m=mT_r，T_r为脉冲重复时间；t为总时间，且$t=\hat{t}+t_{m} $。

因此，风轮机叶片的回波信号可表示为

$ \begin{array}{l} {s_{\rm{w}}}\left( {\hat t,{t_m}} \right) = {A_{\rm{w}}}\sum\limits_{n = 1}^N {\sum\limits_{i = 1}^M {{\sigma _{n,i}}} } rect\left( {\frac{{\hat t - 2{r_{n,i}}\left( {{t_m}} \right)/c}}{{{T_{\rm{p}}}}}} \right) \cdot \\ \;\;\;\;\;\exp \left\{ {{\rm{j}}2{\rm{ \mathit{ π} }}\left[ {{f_{\rm{c}}}\left( {t - \frac{{2{r_{n,i}}\left( {{t_m}} \right)}}{c}} \right) + \frac{1}{2}k \cdot } \right.} \right.\\ \;\;\;\;\;\left. {\left. {{{\left( {\hat t - \frac{{2{r_{n,i}}\left( {{t_m}} \right)}}{c}} \right)}^2}} \right]} \right\} \end{array} $

（4）

式中：N为风轮机叶片数；每个叶片上有M个散射点；σ_{n, i}为第n个叶片上第i个散射点的散射系数；A_w为风轮机方向的雷达天线增益。

假设风电场周围有一匀速运动目标，其径向速度为v_tar，则目标的回波信号为^[18]

$ \begin{array}{l} {s_{{\rm{tar}}}}\left( {\hat t,{t_m}} \right) = {A_{{\rm{tar}}}}{\rm rect}\left( {\frac{{\hat t - 2{r_{{\rm{tar}}}}\left( {{t_m}} \right)/c}}{{{T_{\rm{p}}}}}} \right) \cdot \\ \;\;\;\;\;\;\exp \left\{ {{\rm{j}}2{\rm{ \mathit{ π} }}\left[ {{f_{\rm{c}}}\left( {t - \frac{{2{r_{{\rm{tar}}}}\left( {{t_m}} \right)}}{c}} \right) + } \right.} \right.\\ \;\;\;\;\;\;\left. {\left. {\frac{1}{2}k{{\left( {\hat t - \frac{{2{r_{{\rm{tar}}}}\left( {{t_m}} \right)}}{c}} \right)}^2}} \right]} \right\} \end{array} $

（5）

式中：r_tar(t_m)=R_tar－v_tart；σ_tar为目标散射系数；R_tar为目标初始距离；A_tar为运动目标方向的雷达天线增益。

因此，混合回波信号去载频并作脉压处理后可表示为

$ \begin{array}{l} s\left( {\hat t,{t_m}} \right) = \\ \;\;\;\;\;\;\;{A_{\rm{w}}}{\sigma _{\rm{w}}}\sum\limits_{n = 1}^N {\sum\limits_{i = 1}^M {{T_{\rm{P}}}} } {\rm sinc}\left[ {B\left( {\hat t - \frac{{2{r_{n,i}}\left( {{t_m}} \right)}}{c}} \right)} \right] \cdot \\ \;\;\;\;\;\;\;\;\exp \left( {\frac{{{\rm{j}}4{\rm{ \mathit{ π} }}R}}{\lambda }} \right)\exp \left( {\frac{{ - {\rm{j}}4{\rm{ \mathit{ π} }}}}{\lambda }{l_i}\sin \alpha \cos \theta (t)} \right) + \\ \;\;\;\;\;\;\;\;{A_{{\rm{tar}}}}{\sigma _{{\rm{tar}}}}{T_{\rm{P}}} {\rm sinc} \left[ {{B_{\rm{f}}}\left( {\hat t - \frac{{2{r_{{\rm{tar}}}}\left( {{t_m}} \right)}}{c}} \right)} \right] \cdot \\ \;\;\;\;\;\;\;\;\exp \left( {\frac{{{\rm{j}}4{\rm{ \mathit{ π} }}{R_{{\rm{tar}}}}}}{\lambda }} \right)\exp \left( { - {\rm{j}}2{\rm{ \mathit{ π} }}{f_{{\rm{d}} - {\rm{tar}}}}t} \right) \end{array} $

（6）

式中：f_d-tar=2v_tar/λ；B为雷达带宽；exp(j4πR/λ)和exp(j4πR_tar/λ)为常数相位项；exp(－j4πl_icos θ(t)/λ)和exp(－j2πf_d-tart)为多普勒相位项，分别包含了风轮机和目标的多普勒信息；σ_w为整个风轮机叶片的等效散射系数。根据信号相位与频率间的转化关系可知，第i个散射点的多普勒频率为

$ {f_{{\rm{d}}i}} = \frac{{4{\rm{ \mathit{ π} }}{f_{{\rm{rot}}}}{l_i}\sin \alpha }}{\lambda }\sin \theta (t) $	（7）

因此WTC的最大多普勒频率来自叶尖，其理论最大值为

$ {f_{{\rm{dmax}}}} = \frac{{4{\rm{ \mathit{ π} }}{f_{{\rm{rot}}}}L}}{\lambda } $	（8）

而方位角α的存在相当于对散射点距离乘以一个固定系数sin α，此时，WTC的最大多普勒频率为

$ f_{{\rm{dmax}}}^ * = \frac{{4{\rm{ \mathit{ π} }}{f_{{\rm{rot}}}}L}}{\lambda }\sin \alpha $	（9）

其中，由回波时频图可以获得回波的最大多普勒频率f_dmax^*，而由式(6)可以看出，回波的周期性主要由相位中的cos θ(t)函数项决定，而方位角α并不会影响到cos θ(t)的周期，所以也并未影响到回波的周期大小, 风轮机叶片转速f_rot可由2.2节微动参数粗估计方法获得。因此可以根据风轮机叶片转速f_rot、叶片长度L以及回波的最大的多普勒频率f_dmax^*求得风轮机与雷达间的方位角：

$ \alpha = \arcsin \frac{{f_{{\rm{dmax}}}^*}}{{{f_{{\rm{dmax}}}}}} \le \arcsin \frac{{f_{{\rm{dmax}}}^*\lambda }}{{4{\rm{ \mathit{ π} }}{f_{{\rm{rot}}}}L}} $	（10）

假设某风轮机转速为0.5 r/s, 叶片转动初相为0，叶片长度为20 m，方位角为90°, 其时域WTC如图 2(a)所示，对其做短时傅里叶变换如图 2(b)所示。

由图 2(a)可以看出，由于叶片转动不同脉冲间的脉冲幅度呈现辛格函数变化规律，当叶片与LOS垂直时达到最大值，且其周期与叶片旋转周期相同，这被称作时域“闪烁”现象，该时刻为闪烁时刻。由图 2(b)可以看出WTC的时频域特征主要由叶尖引起的正弦曲线包络和旋转中心引起的零频，以及闪烁时刻产生的频带3部分组成。由于风轮机杂波存在以上特征，而一般目标不具备该特征，因此存在通过信号形态将风轮机杂波与目标分离以达到杂波抑制的可能性。

2 WTC抑制原理

本文所提的动态字典稀疏重构风轮机杂波的抑制原理，首先是基于WTC时频特征对重构字典的参数进行粗估计得到字典的生成范围；其次，在粗估计风轮机微动参数的基础上动态构建过完备字典，其动态构建过程是利用OMP算法匹配最佳原子并以此为基础更新过完备字典，再次匹配；最后，利用动态字典实现对WTC精准稀疏重构，通过滤去重构的风轮机杂波信号，进而得到目标信号，完成对风轮机杂波的抑制。其抑制流程如图 3所示。

2.1 原子构造

由于风轮机回波信号具有明显的微多普勒特征，其瞬时频率呈现正弦函数变化趋势，而目标信号不具备这种时频特征，而Chen给出了旋翼目标的近似数学表达式^[19]，因此字典原子可表示为

$ \begin{array}{l} d\left( t \right) = {\rm sinc}\left( {\frac{{2L}}{\lambda }\cos \left( {vt + \mu } \right)} \right) \cdot \\ \;\;\;\;\;\;\exp \left( {{\rm{j}}\frac{{2L}}{\lambda }{\rm{ \mathit{ π} }}\cos \left( {ut + \mu } \right)} \right) \end{array} $	（11）

式中：υ为旋转因子系数，与风轮机转速f_rot有关; μ=θ₀+2π(n－1)/N为相位因子，与风轮机叶片初始旋转角相关，其中风轮机叶片长度可知，当风轮机与雷达不在同一水平面时，用转换等效长度$ \hat{L}=L \sin \alpha$代替L，即可将风轮机与雷达不在同一水平面内。因此，字典参数主要为叶片转速f_rot和叶片初始相位θ₀。假设某原子相关参数与图 2中风轮机参数相同，该字典原子的时频特性如图 4所示。

对比图 2和图 4可以看出，原子具有与WTC模型十分相似的正弦曲线包络、零频线以及频带的时频特性，只是在幅度上存在差异，因此该原子所构建的字典可作为WTC的重构字典。

2.2 基于时频特征的微动参数粗估计

对WTC生成重构字典，在未有任何先验信息条件下，需要在很大范围内生成过完备字典，字典原子数目越多，其对运行内存要求也越高，而字典也具有很大的冗余度。为了减少字典的冗余度，提出了一种基于WTC时频特征的字典参数粗估计方法，通过对叶片转速和叶片初相进行估计，缩小字典原子的生成范围，之后在估计的微动参数的基础上生成字典，极大地减小了字典原子个数，提高了算法的效率。

如图 2(b)所示，WTC的时频特征主要具备以下特点：

1) WTC时频特征由零频带、正弦包络和时频“闪烁”频带3部分组成，其中正弦包络周期与风轮机叶片旋转周期保持一致。

2) 时频“闪烁”的频带范围在[－f_dmax, f_dmax]之间，具有回波中最宽的频率带，且回波能量远大于零频带和正弦包络。

3) 当且仅当风轮机叶片与LOS垂直时，时频“闪烁”出现，其频率达到最大值。如图 5(a)所示，当风轮机叶片转向雷达方向时，“闪烁”频带为正频带，当风轮机叶片背离雷达方向时，“闪烁”频带为负频带。因此在一个叶片旋转周期内，一个叶片将出现一正一负两个“闪烁”频带。

因此WTC时频的“闪烁”频带相关信息可作为微动参数的提取依据。对于一个三叶片的风轮机而言，由于风轮机的构造的对称性，其相邻2个叶片间的角度间隔约为120°，则在2个相邻正向闪烁时刻内叶片共转过了1/3转，因此可根据同向闪烁时刻间的间隔平均值Δt，可求得风轮机的转速估计值为

$ {{\hat f}_{{\rm{rot}}}} = 1/3\Delta t $	（12）

同时，假设初始时刻t₀的叶片初始相位为θ₀，当叶片达到第一个闪烁时刻t₁时，叶片转过的角度为θ₁，如图 5(b)所示。叶片相位的变化实际上是时频图的平移，利用几何知识可得叶片的初始相位估计值为

$ {{\hat \theta }_0} = \frac{{\rm{ \mathit{ π} }}}{2} - {\theta _1} = \frac{{\rm{ \mathit{ π} }}}{2} - 2{\rm{ \mathit{ π} }}{{\hat f}_{{\rm{rot}}}}{t_1} $	（13）

因此，对初始相位的求解实际上转化为对$\hat{f}_{\mathrm{rot}} $和t₁的求解，式(12)已给出$\hat{f}_{\mathrm{rot}} $的求解方式，只需根据时频图提取第1个闪烁时刻t₁的值，即可求得该叶片的初始相位。单独以第1个闪烁时刻t₁作为风轮机叶片初始相位的求解依据会因为测量误差等造成估计精度下降，存在不确定性。但可以通过所有闪烁时刻求得第1个闪烁时刻的平均值，降低测量意外的发生。如图 5(c)所示，每2个相邻同向闪烁时刻之间相差1/3f_rot，因此平均值t₁的求解过程表示为

$ \begin{array}{l} {{\bar t}_1} = \\ \;\;\;\;\frac{{{t_1} + \left( {{t_2} - 1/3{f_{{\rm{rot}}}}} \right) + \cdots + \left[ {{t_n} - (n - 1)/3{f_{{\rm{rot}}}}} \right]}}{n} \end{array} $	（14）

将t₁和$ \hat{f}_{\mathrm{rot}}$代入式(13)即可求得该叶片初始相位，而其他叶片的初始相位可在$\hat{\theta}_{0} $基础上加上相应的叶片角度间隔即可。

基于WTC时频特征的微动参数粗估计实质上是利用时频图的图像特征对$ \hat{f}_{\mathrm{rot}}$和t₁进行求解，因此可以利用图像处理的平滑处理和二值化处理^[20-21]，只保留WTC闪烁频带。首先，对回波时频图进行平滑处理，剔除一些突变的数据点，如过大的噪声等。其次，由于风轮机回波时频“闪烁”频带占据了回波的大部分能量，闪烁时刻的频带信号能量也要明显强于噪声，因此可以设置相应的能量门限η，对回波时频图进行二值化处理，剔除与时频“闪烁”无关的时频单元。之后对二值化处理后的时频数据通过寻峰函数提取相关时间信息。由于短时傅里叶变换时间分辨率问题，得到的频带实际占据了多个时间采样点，从而可能会造成估计偏差。但是该步骤只是对微动参数进行粗估计，以达到缩小动态字典的生成范围的目的，因此可以把二值化结果的频带中心的时间点作为闪烁时刻的时间点，近似估计出风轮机的转速和叶片初相。

2.3 动态字典构建

构建过完备字典D=[d_k]，主要是对原子参数{f_rot, θ}进行离散化表示。参数间隔选取越小，其原子与WTC的匹配度越高，但参数间隔越小则原子数目越多，字典也越大，匹配算法的收敛速度也越慢。为了进一步提高重构信号的精度，本文在2.2节粗估计微动参数的基础上，提出通过构建动态字典以利用较少的字典原子达到更好的抑制效果。为了满足稀疏分解的条件，原子d_k的长度与信号长度一致，且需要对字典作归一化处理。

该字典动态构建方法对参数均匀离散化{jΔf, kΔθ}(1≤j≤M₁, 1≤k≤M₂)离散的原子间步长为Δd_λ{Δf, Δθ}，初级字典的原子步长可由粗估计参数的误差决定。基于OMP算法可匹配出本级字典中与风轮机杂波匹配度最高的原子d_r，以该原子的基础可在参数范围[d_r－Δd_λ, d_r+Δd_λ]上更新生成新一级的字典D_new，则新一级各原子间步长更新为$\Delta \widetilde{d}_{\lambda}\left\{\Delta f / M_{1}, \Delta \theta / M_{2}\right\} $，动态字典生成过程如图 6所示。因此字典一直处于动态筛选、动态更新的过程中，且字典范围在逐级缩小，原子间步长$\Delta \widetilde{d}_{\lambda} $在逐级减小，即原子的准确度在逐级提升。

采用动态字典构造可以将原子数一直控制在K=M₁M₂个，其构建实质上是对目标参数的逐级迭代匹配，在不增加运行内存的前提下实现了核原子更精准的构建，从而实现了对风轮机杂波信号s_w更加精准的稀疏重构，可以达到更好的风轮机杂波抑制效果。

2.4 基于OMP算法的匹配原理

在动态字典构建过程中，需要利用OMP算法筛选出最匹配原子进而缩小原子构建的范围。OMP算法利用过完备字典原子D对信号s(t)其进行线性组合匹配。由于OMP算法对选出对的原子进行了正交化处理，因此被匹配过的原子不会被再次匹配，这使得OMP算法的收敛速度比匹配追踪(MP)算法更快^[22]。因此，信号s(t)可以由每次迭代匹配的原子以及残差来线性表示：

$ s\left( t \right) = \sum\limits_r {{\alpha _r}{d_r}} + {\rm{res}} = {s_{{\rm{wr}}}} + {\rm{res}} $	（15）

式中：α_r为稀疏系数；d_r为OMP选择的匹配原子；s_wr为重构的WTC；res为残差。

OMP算法步骤为

步骤1  初始化残差res₀=s。

步骤2  在字典D中选择与残差res₀最匹配的原子，即两者内积最大，该原子表示为d_r1。

$ {d_{r1}} = \arg \left\{ {\sup \left| {\left\langle {s,{d_k}} \right\rangle } \right|} \right\} $	（16）

步骤3  将所选择的原子组成匹配原子矩阵Φ，定义Φ所在空间的正交投影算子为

$ P = \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}{\left( {{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}^{\rm{T}}}\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}} \right)^{ - 1}}{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}^{\rm{T}}} $	（17）

在信号中减去R₀在此空间上的正交投影，得到新的残差res₁

$ {\rm{re}}{{\rm{s}}_1} = {\rm{re}}{{\rm{s}}_0} - {P_0}{\rm{re}}{{\rm{s}}_0} $	（18）

步骤4  更新残差res，并循环执行步骤2和步骤3

$ {\rm{re}}{{\rm{s}}_{m + 1}} = {\rm{re}}{{\rm{s}}_m} - {P_m}{\rm{re}}{{\rm{s}}_m} $	（19）

步骤5  当循环次数大于N，跳出循环，计算此时的稀疏系数矩阵。

因为只有混合回波信号内的WTC部分能够在该过完备字典上稀疏分解，所以目标信号和噪声信号作为残差被保留下来。

2.5 应用分析

前文是在观测时间足够长的基础上进行的理论研究，但是在实际应用中，空管雷达天线大多数情况下以4~15 r/min的扫描频率f_scan对空域进行扫描，由于天线和风轮机均处于不同的运动中，因此在一个天线扫描周期内得到WTC的相干脉冲十分有限，且风轮机处于慢速旋转状态，很难在雷达一帧回波数据内充分表现WTC的特征。但是雷达天线扫描模式下的回波信号可以看作是由非扫描模式下回波数据按f_scan采样频率采样获得。假设雷达天线波束等效宽度为β，单位为rad，在一个雷达天线波束宽度内，天线可看作固定不动，雷达在每个方位角的波束驻留时间t_w可表示为

$ {t_{\rm{w}}} = \frac{\beta }{{2{\rm{ \mathit{ π} }}{f_{{\rm{rot}}}}}} $	（20）

则在天线转动情况下每帧回波数据中WTC持续时间为t_w。

对于本文所提稀疏重构算法而言，天线转动只是影响到了需要重构信号s_wr的长度，由于其回波可看作是按f_scan采样频率采样获得，相应地，也需要对原子字典进行同频采样，这样即可得到天线转动情况下的原子字典。所以本质上该基于重构的抑制方法也可适用于天线转动下情况。由于t_w持续时间较短，而风轮机旋转周期大于波束驻留时间，因此在一帧数据中风轮机回波的周期特性不能完整体现，难以基于时频特征提取到风轮机微动参数。但是风轮机的时频周期特征可采取对风轮机所在区域进行短暂凝视的策略，或借鉴雷达杂波图思想利用长期数据积累获得，通过时频周期特征估测风轮机有关微动参数，实现对重构字典范围的进一步缩小，减少字典的冗余度，提高算法的运行效率。

需要指出的是，基于微动特征的参数估计作为一种先验信息，只是为了给出重构原子参数的初始值，微动参数估计精度越高，该重构算法占用内存越少、运行时间越快。在工程实践中，该初始值也可以根据经验值给出，比如利用现代空管雷达的气象监测功能估测风速大小并结合风轮机参数可换算出叶片旋转频率的估计值，或通过统计手段利用该地风轮机转速的平均值作为重构字典的初始值，但经验起始值的低精确度会在一定程度上增加算法的运算量。

3 仿真结果及分析

为了验证本文所提方法对风轮机杂波信号的抑制作用，本文采用2组实验进行分析。实验1采用理想仿真场景，仿真分析本文方法的抑制性能和稳定性，以验证该方法的最优性能，见3.1节；实验2考虑实际应用场景，即在天线转动、观测时间有限情况下，验证本文方法在工程应用中的有效性。仿真参数设置如表 1所示，雷达天线波束形状为余割平方，叶片散射点间隔为λ/10。

3.1 仿真实验1 3.1.1 WTC信号仿真及影响分析

假设风电场内有3台同类型风轮机，风轮机尺寸相同且不存在互相干扰问题，风轮机具体参数设置如表 2所示。假设风轮机2上空5 km高度有一匀速飞行目标，其进入风电场的初速度v₀=50 m/s，风轮机2和目标处在同一雷达波束内，回波信号建模仿真如图 7所示。

从图 7可以看出，回波在18、20、28 km回波幅度均随脉冲数出现了明显的周期性起伏变化，在不同的脉冲时刻WTC幅度也不同。

图 8为对雷达回波数据进行脉冲积累处理结果。由于受到WTC微动特征的影响，回波在18、20、28 km处产生了3条频带。风轮机1产生的频带为[-442.63, 442.63]Hz, 风轮机2产生的频带为[-416.89, 416.89]Hz, 风轮机3产生的频带为[-420.35, 420.35]Hz。其中，目标的多普勒频率被风轮机所产生的频带完全淹没，脉冲积累难以提取到目标的相关信息。

为了定性分析WTC对目标检测的影响，每隔250个脉冲选取一组数据进行CA-CFAR检测，这8组数据可以代表不同脉冲时刻WTC对目标检测的影响。结果如图 9所示，图中虚线代表CFAR门限。风轮机1~3分别在第240、267、374距离单元，由图 9可以看出，风轮机1在第250、第1 000、第1 500和第2 000个脉冲信号幅度较强，均被当作目标检测出来；风轮机3在第500、第1 250、第1 500、第1 750和第2 000个脉冲被当作目标检测出来，因此由于WTC幅度变化会造成雷达目标检测的虚警问题。而风轮机2与目标在同一距离单元，WTC与目标信号交织在一起，无法对其进行有效区分，造成目标信息提取的困难。因此风轮机的存在不仅会造成对目标检测的遮蔽问题，由于WTC幅度的起伏变化，还会造成检测的虚警问题。

3.1.2 WTC微动参数粗估计仿真

对回波进行微动参数粗估计，以目标所在距离单元信号为例，其回波时频图和频率图如图 10所示。根据时频图信号强度分布，本文选取时频单元最大信号能量的80%这一经验值作为分割门限η，对其时频图经过平滑、二值化图像处理后如图 11所示。其中黑色代表提取出的闪烁时刻频带，白色为背景。

时频单元的最大信号能量为30 dB, 因此设门限为η=24 dB，由于WTC频带杂波大于20 dB，因此只有闪烁频带被提取出来，噪声、目标信号均被滤去。此时根据图可清晰地计算出风轮机2的转速和叶片初相。同理可分别求得风轮机1和风轮机2的相关参数，其微动参数粗估计结果如表 3所示。

由表 3可以看出基于WTC时频特征的微动参数提取方法具有一定精确度，其中叶片初相估计误差明显大于转速误差，这是因为叶片初相需要利用转速估计值求解而造成了误差积累。同时，粗估计误差主要由于图像的分辨率以及STFT的时间分辨率造成的，因此无法对微动参数做到精准提取，但是本文只是利用粗估计结果可以缩小字典的生成范围，减少字典的冗余度。

3.1.3 杂波抑制性能仿真与分析

对3.1.1节中仿真回波进行杂波抑制。利用3.1.2中粗估计结果，可将字典参数生成范围缩小在$ \left[\hat{f}_{\text {rot }}-5 \% \hat{f}_{\text {ret }}, \hat{f}_{\text {rot }}+5 \% \hat{f}_{\text {rot }}\right]、\left[\hat{\theta}_{0}-5 \% \hat{\theta}_{0}\right.\left.\hat{\theta}_{0}+5 \% \hat{\theta}_{0}\right]$之间，设置5%为浮动范围是为了尽可能大地减小粗估计误差带来的影响。假设动态字典大小为K=50×50个原子，动态字典生成级数为3级。

对每个距离单元数据进行动态字典稀疏重构的杂波抑制后，同样每隔250个脉冲选择数据对其进行CA-CFAR检测，结果如图 12所示。对比图 9可以看出，由于风轮机1和风轮机3的杂波均被有效抑制，因此未被CFAR检测出来，所以本文所提方法能够有效地抑制WTC，可以解决目标检测的虚警问题。同时，风轮机2杂波也被有效抑制，目标信号功率趋于20 dB稳定值，所以该方法对目标与风轮机杂波在同一距离单元时也可同样适用。因此该滤波方法既抑制了WTC, 又对目标回波功率影响很小。

以目标与风轮机2为例，单独取出目标所在距离单元信号，其杂波抑制后如图 13所示。图 13(a)、图 13(b)是以粗估计结果构造匹配原子对WTC抑制后的时频图和频率图；图 13(c)和图 13(d)是在粗估计参数误差范围内利用动态字典杂波抑制后的时频图和频率图；图 13(e)和图 13(f)是杂波抑制后经过三脉冲MTI滤波器的时频图和频谱图。

由图 13(a)和图 13(b)可以看出，该杂波稀疏重构方法对原子参数精度较为敏感，尽管粗估计微动参数结果具有一定精确度，但利用其稀疏重构的杂波与原杂波信号存在失配现象，在零频附近抑制效果较好，在抑制性能较差，因此无法实现目标与杂波的有效分离。因此利用粗估计微动参数重构杂波对WTC并未得到有效地抑制。

由图 13(c)和图 13(d)可以看出，在粗估计微动参数的基础上，进行进一步动态字典的OMP算法，可以对WTC进行有效地抑制。但需要指出的是，原子模型与仿真信号在零频附近存在失配现象，这是由于模型间差异性造成的，因此在图 13(c)和图 13(d)中可以看出WTC在零频仍有小部分残余，由于这部分残余仅在零频附近，因此可当作固定杂波，通过MTI滤波器可对其进行滤除，如图 13(e)和图 13(f)所示。所以本文方法可以实现对风轮机杂波信号的有效抑制。

为了量化分析WTC抑制效果，本文引入杂波改善因子I的概念^[23]，其定义为杂波抑制滤波器输出端信杂比与输入端信杂比的比值：

$ I = \frac{{{{\left( {{s_{{\rm{tar}}}}/{s_{\rm{w}}}} \right)}_{\rm{o}}}}}{{{{\left( {{s_{{\rm{tar}}}}/{s_{\rm{w}}}} \right)}_{\rm{i}}}}} $	（21）

图 14仿真了基于微动参数粗估计结果的动态字典和一般动态字典2种方法下的杂波改善因子I随字典级数的变化结果，其中2种方法每级字典的大小均固定为K=50×50。

由图 14可知，利用时频特征的微动参数粗估计结果可以在动态字典的第2级对WTC进行了比较好的抑制，杂波改善因子最终稳定在30.21 dB，这是由于微动参数粗估计结果缩小动态字典的生成范围，进而提高动态字典的方法效率。而对于未采用粗估计结果的动态字典，WTC在字典的前3级的抑制效果相对较差，在第4级才基本达到与“粗估计+动态字典”方法的抑制效果，并且最终杂波改善因子稳定在30.18 dB。对于动态字典，其每级要在K=50×50个原子间进行OMP匹配，多一级字典就增加了一定量的计算量。因此相比较而言，基于微动参数粗估计结果的动态字典具有更高的效率和更好的抑制效果, 也充分说明了利用微动参数粗估计结果缩小字典参数范围的必要性。

3.1.4 杂波抑制稳定性仿真与分析

该仿真是动态字典原子数保持为K=50×50个下进行的。用平均相对误差来衡量稀疏重构杂波与WTC之间准确度的指标，其求解公式为

$ {\rm mean}\_{\rm mse} = E\left( {\frac{{{{\left| {{s_{\rm{w}}} - {s_{{\rm{wr}}}}} \right|}^2}}}{{{{\left| {{s_{\rm{w}}}} \right|}^2}}}} \right) $	（22）

在杂噪比CNR=50 dB前提下，雷达相关参数保持不变，以风轮机1相关参数为例，通过蒙特卡罗实验仿真了动态字典级数与重构相对误差之间的变化关系，如图 15所示，蒙特卡罗次数为100次。其中0级代表粗估计参数重构原子误差。

由图 15可知，粗估计重构原子的平均相对误差相对较大，经过两级动态字典稀疏重构后，重构相对误差大幅度下降，当动态字典生成3级以上，其重构相对误差随动态字典级数变化较小，误差趋于收敛状态。因此，该方法一般对动态字典生成3级即可得到较为准确的重构WTC。

同时，该稀疏重构方法对原子微动参数精度较为敏感，尽管当粗估计微动参数结果较高，但是重构误差却相对较大，杂波改善因子也很小。在CNR=50 dB前提下，通过调整重构相对误差大小，使误差在10^-3~10⁰量级均有分布，图 16仿真了重构相对误差与改善因子间的关系。

由图 16可知，重构相对误差与杂波改善因子整体呈反比关系，但是对于不同重构相对误差量级，其反比增速不一样。当重构相对误差在大于10^-1量级时，改善因子小于10 dB，而粗估计重构相对误差正在这个范围之内，因此其粗估计重构的抑制性能相对较差。当重构相对误差达到10^-2量级时，其杂波改善因子达到明显提升。当重构相对误差达到10^-3量级时，杂波改善因子可达到30 dB以上。由图 15可知，本文所采用动态稀疏重构方法在动态字典达到3级以上时，重构相对误差稳定在10^-3量级，因此本文所提方法杂波改善因子可达30 dB以上。

图 17通过蒙特卡罗实验仿真了不同杂噪比(CNR)条件下，该抑制方法的WTC杂波改善因子I变化情况。由仿真结果可知在CNR为[10, 30]dB之间时，杂噪比越高，改善因子也越大。当杂噪比大于30 dB之后，该算法对WTC的杂波改善因子趋于稳定。杂噪比的大小在对算法对WTC抑制会有一定的影响，但是在实际应用中，风轮机叶片往往具有很大的RCS，因此其回波功率较大，CNR要大于SNR，因此这种影响实际上是微小的，另一方面也可说明该抑制方法较为稳定且对噪声有一定的容忍度。

3.2 仿真实验2

假设雷达天线波束等效宽度为β=0.052 rad，雷达天线的扫描频率f_scan=5 r/min，其他雷达仿真参数如表 1所示，风轮机参数与风轮机2参数相同。由于在一个叶片旋转周期内，当叶片与雷达视线垂直时，WTC幅度和多普勒频率范围达到最大，此时对目标检测的影响程度也达到最大，因此该部分仿真以风轮机与雷达视线垂直情况下的回波数据为研究对象。

由式(20)可知，波束驻留时间t_w=0.1 s。根据3.1节分析可知，WTC幅度随脉冲数呈现周期起伏变化特性，而在该天线扫描模式下，回波数据只是相当于完整周期中的一段，因此其回波并未表现出旋转周期特性，图 18(a)为天线扫描情况下回波时频特征，图 18(b)为回波频率特征。

由仿真结果可知，由于该扫描时刻风轮机叶片垂直于LOS, 此时WTC多普勒频率范围达到最大，WTC产生“时频闪烁”现象，由图 18(a)可知在[-416.89, 0]Hz具有较大的幅度，而目标的多普勒频率为334 Hz，因此风轮机产生的多普勒频率与目标多普勒频率重叠在一起，造成目标检测的困难，如图 18(b)所示。同时，风轮机回波的零频带特征依然存在，“时频闪烁”在回波中所占比例较大，更增加了目标检测的难度。

利用本文所提稀疏重构抑制方法对该段扫描数据处理，并经过MTI滤波后，其结果如图 19所示，图 19(a)为利用本文所提方法抑制后的时频图，图 19(b)为抑制后的频谱图，其中动态字典大小仍为个原子，字典参数初始值为$ \hat{f}_{\text {rot }}=0.49 \mathrm{r} / \mathrm{s}$，$\hat{\theta}_{0}=1.21 \mathrm{rad} $，动态字典生成级数为3级。对比图 18(a)和图 19(a)可知，“时频闪烁”被有效地抑制，经过MTI处理后，零频附近的杂波被完全抑制，而目标回波信息被保留下来。通过图 19(b)可知，目标回波频率信息被完整显露出来，剩余噪声和杂波相较目标频率很小，并未影响到目标频率信息的提取。因此由仿真结果可知，本文所提抑制方法在天线扫描模式下同样适用。

4 结论

1) 实验结果显示，基于时频特征的WTC微动参数粗估计具有一定精度，但WTC稀疏重构对原子参数较为敏感，粗估计重构结果不足以有效抑制WTC。

2) 本文在粗估计基础上利用动态字典不仅实现对风轮机杂波的精准重构，还降低了字典冗余度。仿真实验证明了该方法可以有效地抑制WTC，解决风电场引起的虚警问题。同时该方法的稳定性高，其杂波改善因子在高CNR条件下可以保持良好的性能。

3) 通过仿真实验表明该稀疏重构抑制方法在雷达天线扫描模式下同样适用。