CN110174650A - 基于两维联合插值的气象雷达风电场杂波抑制方法 - Google Patents

Abstract

本发明公开了一种基于两维联合插值的气象雷达风电场杂波抑制方法，该方法针对低信噪比情况下，首先对同时包含风力涡轮机杂波(WTC)和气象信号的距离单元置零，以该距离单元为中心，两端各选取2个相邻距离单元，通过相干积累提高信噪比。再分别计算出距离域一维插值气象信号估计值和多普勒域一维插值气象信号估计值，最后根据最小二乘法计算出的两维权系数进行两维联合插值。仿真实验结果表明，本方法在低信噪比情况下，可有效提高插值精度，实现气象信号的精确恢复，且运算量小，适合工程应用。

Description

基于两维联合插值的气象雷达风电场杂波抑制方法

技术领域

本发明涉及一种基于两维联合插值的气象雷达风电场杂波抑制方法，特别涉及一种基于距离—多普勒域两维联合插值的气象雷达风电场杂波抑制方法，属于气象雷达杂波抑制领域。

背景技术

石油价格的不断上涨，以及全球气候变化的加剧，为许多替代能源——可再生清洁能源，铺平了道路，其中增长最快的是风力发电。风力发电作为一种清洁可再生能源，已得到世界各国广泛关注。由于风电场数量和规模呈指数式增长，风力涡轮机叶片高速旋转引起的时变运动杂波已成为制约气象雷达探测性能的主要限制因素。

传统杂波抑制技术都是基于杂波信号与气象回波在多普勒域或距离－多普勒域可分离条件，采用滤波或分离算法对杂波进行有效抑制。但风力涡轮机杂波由于大型叶片高速旋转，其多普勒频谱严重展宽甚至产生混叠，因此气象信号淹没在风力涡轮机杂波中，导致时、频滤波或分离算法均无法在抑制杂波信号的同时实现对气象信号的无损保留。因此时域滤波、频域滤波等方法均失效。对此，国内外科学家就插值算法进行了深入研究，常用的插值算法有拉格朗日插值、样条插值、统计不变插值(ISI)等。传统插值算法在低信噪比情况下，无法恢复出气象信号。且传统插值算法只是利用了气象数据的空间信息，而气象数据是典型的时空信号。面向空间的插值算法直接应用到时空过程可能导致有价值信息在时间维的丢失，影响插值精度。

发明内容

本发明所要解决的技术问题是提供一种基于两维联合插值的气象雷达风电场杂波抑制方法，改善了低信噪比情况下，传统插值无法恢复出气象信号的问题，并且将一维插值扩展到二维插值，有效地提高了插值性能，实现了气象信号的精确恢复。

本发明一种基于距离—多普勒域两维联合插值的气象雷达风电场杂波抑制方法，在低信噪比情况下，首先进行相干积累提高信噪比，将气象雷达回波信号转换至距离—多普勒域后，同时联合距离域和多普勒域的信息，提高插值精度，实现气象信号的精确恢复。

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

本发明提供一种基于两维联合插值的气象雷达风电场杂波抑制方法，，该方法的具体步骤如下：

步骤1，对同时包含风力涡轮机杂波WTC和气象信号的距离单元置零，以该距离单元为中心，两端各选取2个相邻距离单元，通过相干积累提高信噪比；

步骤2，对雷达扫描回波信号进行预处理，构建雷达插值初始矩阵；

步骤3，分别计算距离域一维插值气象信号估计值和多普勒域一维插值气象信号估计值；

步骤4，根据最小二乘法计算出的两维权系数，进行两维联合插值，得到两维联合插值气象信号估计值；

步骤5，输出抑制WTC后的信号。

作为本发明的进一步技术方案，步骤1具体为：

对同时包含WTC和气象信号的第a个距离单元置零，对第a±i个距离单元进行相干积累，通过快速傅立叶变换实现：

式中，i＝1,2，x_a±i(n)、分别为第a±i个距离单元雷达当前扫描与前一次扫描回波信号，X_a±i(k)、分别为x_a±i(n)、经过傅立叶变换后的结果，K为脉冲数。

作为本发明的进一步技术方案，步骤2具体为：

对第a±i个距离单元雷达当前扫描回波信号进行频谱分析，分别估计出其气象信号频谱峰值对应的多普勒单元；

设各单元气象信号频谱峰值所在的最小多普勒单元序号为m_min、最大多普勒单元序号为m_max，选取雷达两次扫描回波的第m_min-10～m_max+10个多普勒单元为插值样本，将其他多普勒单元置零，以该结果针对对应的X_a±i(k)、进行更新，并构建雷达插值初始矩阵Z、Z₀：

作为本发明的进一步技术方案，多普勒峰值有如下两种特殊情况：

(1)多普勒峰值出现在零频附近，此时m_min-10＜0，则取第1～m_min、m_max～m_max+10、K-10+m_min～K个多普勒单元为插值样本；

(2)多普勒峰值出现在脉冲重复频率附近，此时m_max+10＞K，则取第1～10-K+m_max、m_min-10～m_min、m_max～K个多普勒单元为插值样本。

作为本发明的进一步技术方案，步骤3中距离域一维插值气象信号估计值为：

式中，为第a个距离单元的的第k个多普勒单元距离域一维插值气象信号估计值，α_a±i(k)为第a±i个距离单元在第k个多普勒单元的权重因素，径向基函数R_a为第a个距离单元的距离，R_a±i为第a±i个距离单元的距离，b为可调参数，||·||表示范数，

作为本发明的进一步技术方案，步骤3中多普勒域一维插值气象信号估计值为：

式中，为第a个距离单元的的第k个多普勒单元多普勒域一维插值气象信号估计值，N为被选中距离单元的总数，雷达当前扫描回波信号与前一次扫描回波信号的互相关矩阵

作为本发明的进一步技术方案，步骤4中两维联合插值气象信号估计值为：

式中，为第a个距离单元的第k个多普勒单元两维联合插值气象信号估计值，

作为本发明的进一步技术方案，步骤5中输出的抑制WTC后信号为：

式中，为第a个距离单元的第n个脉冲抑制WTC后的气象信号高精度估计值。

本发明采用以上技术方案与现有技术相比，具有以下技术效果：

1.基于距离—多普勒域两维联合插值的气象雷达风电场杂波抑制这一方法，运算量低，有利于工程实现；

2.仿真实验结果表明，在低信噪比情况下，基于距离—多普勒域两维联合插值的气象雷达风电场杂波抑制方法可有效提高插值精度，实现气象信号的精确恢复。

附图说明

图1是风力涡轮机与雷达关系图；

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

图3是气象雷达回波信号距离—多普勒图；

图4是三种插值方法的误差分析对比图。

图5是三种插值方法恢复出的气象信号速度估计值的均值。

具体实施方式

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

本发明基于距离—多普勒域两维联合插值的气象雷达风电场杂波抑制方法，如图2所示，其主要步骤如下：

步骤一、雷达回波信号仿真，具体为：

在本发明中，雷达回波信号包括：地杂波信号、气象信号、风力涡轮机信号，气象雷达回波信号距离—多普勒图如图3所示，三种回波信号仿真如下：

1.1地杂波信号：采用距离单元划分法，每个距离单元的宽度为f_s为采样频率，c为光速。将一个距离单元划分成180个杂波散射单元，每个距离单元的杂波数据即为这180个杂波散射单元的叠加。第l个距离单元，第n个脉冲地杂波回波信号为：

其中，ω_t＝2πf_d/f_r为时域角频率，R_l为第l个距离单元的距离，F为天线增益，f_d为多普勒频率，f_r为脉冲重复频率，K表示脉冲数，L表示气象雷达与风力涡轮机之间回波信号移动路径上距离单元的总数。

1.2气象信号：气象目标是分布式目标，其回波信号可由多个散射单元的回波叠加形成。对第l个距离单元内所有散射单元的回波进行相干叠加，即为气象信号回波。第l个距离单元，气象信号在第n个脉冲采样可以表示为:

其中，U为第l个距离单元内散射单元的总数，A_u为气象目标粒子u的幅度，由雷达方程决定。ω_t为时域角频率，由下式可得：

其中，v_u为气象目标粒子u的径向速度，λ为雷达波长。

1.3风力涡轮机信号：风力涡轮机由叶片、桅杆和轮机舱组成，风力涡轮机杂波主要由叶片和桅杆杂波组成。叶片是转动的，桅杆是静止的，采用散射点叠加模型。

叶片信号：风力涡轮机有三个叶片，它们的间隔是120度，所以两两之间在角度上有一个120度的滞后，首先考虑一个叶片，那么其它两个叶片只是相差120度。将每个叶片看成由多个散射点组成，一个叶片的杂波就是这些散射点的叠加，三个叶片散射点的叠加就是风力涡轮机叶片杂波。第n个脉冲，风力涡轮机叶片一个散射点p的回波信号：

其中，f_c为雷达频率，为时延，T_r为脉冲重复周期，R_p为散射点p到雷达的距离：

其中，R_g为风力涡轮机到雷达的距离，d_p为叶片散射点p到叶片轴心的距离，为参考叶片与雷达视线的夹角：

其中，α,β分别为方位角和俯仰角，θ为参考叶片与Y轴正方向的夹角。

去掉载波，可得基带信号为：第n个脉冲，所有叶片总的回波信号为(风力涡轮机叶片数一般为3个)：

其中，P为叶片散射粒子的总数。

桅杆信号：桅杆是静止的，桅杆杂波也看成是多个散射点的叠加。第n个脉冲，桅杆散射点q的回波信号为：

其中，R_q为桅杆散射点q到雷达的距离：

其中，d_q为桅杆散射点q到叶片轴心的距离。

去掉载波，可得基带信号：第n个脉冲，整个桅杆的回波信号为：

其中，M为桅杆散射点的总数。

综上所述，第n个脉冲，风力涡轮机杂波信号为：

步骤二、相干积累，具体为：

对同时包含WTC和气象信号的距离单元置零，设该距离单元序号为a，对第a±i(i＝1,2)个距离单元进行相干积累，可通过快速傅立叶变换实现：

式中，x_a±i(n)、分别为雷达当前扫描与前一次扫描回波信号，X_a±i(k)、分别为x_a±i(n)、经过傅立叶变换后的结果，K为脉冲数。经过相干积累，信噪比提高了K倍。

步骤三、雷达扫描回波信号预处理，具体为：

对训练单元(第a±i个距离单元)雷达当前扫描回波信号进行频谱分析，分别估计出其气象信号频谱峰值对应的多普勒单元。设各训练单元气象信号频谱峰值所在的最小多普勒单元序号为m_min，最大多普勒单元序号为m_max。考虑雷达天线副瓣的影响，选取雷达两次扫描回波的第m_min-10～m_max+10个多普勒单元为插值样本，将其他多普勒单元置零，以该结果针对对应的X_a±i(k)进行更新(因为k是从1～K的，此处指的是逐个k以此进行更新)，并构建雷达插值初始矩阵Z₀,Z：

由于气象信号径向速度的未知性，其多普勒峰值可能出现两种特殊情况：(1)多普勒峰值出现在零频附近，此时m_min-10＜0；(2)多普勒峰值出现在脉冲重复频率附近，此时m_max+10＞K；由于频域信号的周期性，针对以上两种特殊情况，采用如下方法选取插值样本：(1)当m_min-10＜0时，取第1～m_min、m_max～m_max+10、K-10+m_min～K个多普勒单元为插值样本；(2)当m_max+10＞K时，取第1～10-K+m_max、m_min-10～m_min、m_max～K个多普勒单元为插值样本。

步骤四、距离域一维插值，具体为：

距离域一维插值逐多普勒单元分别利用径向基函数来近似缺失数据，径向基函数定义如下：

其中，R_a为第a个距离单元的距离，R_a±i为第a±i个距离单元的距离，b为可调参数，||·||表示范数。

由下式可解出距离域一维插值的权重因素α_a±i(k)：Z＝Q_rα

其中，Z为步骤三中构建出的雷达插值初始矩阵。可得：即解出α_a±i(k)，则可由插值方程得出距离域一维插值气象信号估计值：

其中，为第k个多普勒单元距离域一维插值气象信号估计值，系数α_a±i(k)为第a±i个距离单元在第k个多普勒单元的权重因素。

步骤五、多普勒域一维插值，具体为：

计算雷达当前扫描回波信号与前一次扫描回波信号的互相关矩阵

其中，N为被选中距离单元的总数，本发明中N＝4，[·]^T表示转置。

可得互相关矩阵的第一阶矩：其中E[·]表示期望。根据互相关矩阵的第一阶矩，平移气象雷达前一次扫描回波信号，并求平均，得到第k个多普勒单元多普勒域一维插值气象信号估计值：

步骤六、两维联合插值，具体为：

设为第k个多普勒单元两维联合插值气象信号估计值，令：

其中，w为两维权系数。引入残差平方和函数D：

式中，

为了使插值精度尽可能高，实现气象信号的精确恢复，由最小二乘法可知，应使D(w)尽可能小。当时，D(w)取最小值，记作：

通过对D(w)进行微分求最值，可以得到：因此可解出

根据上式求得的两维权系数可得出两维联合插值气象信号估计值：

其中，为第k个多普勒单元两维联合插值气象信号估计值。

步骤七、输出抑制WTC后的信号，具体为：

对步骤六中求得的两维联合插值气象信号估计值进行逆傅立叶变换，可得第n个脉冲抑制WTC后的气象信号高精度估计值：

其中，为第a个距离单元第n个脉冲抑制WTC后的气象信号高精度估计值。

下面通过计算机仿真验证本发明方法的有效性。雷达系统仿真参数如表1所示，风力涡轮机系统仿真参数如表2所示。在本仿真实验中，以第25个距离单元为WTC和气象信号共存的距离单元。表3为第25个距离单元信号速度估计值的误差，由表3可以看出，第25个距离单元信号速度估计值的误差为4m/s左右，表明该距离单元已被WTC污染，无法直接从该距离单元准确估计出气象信号的速度。本仿真实验选取第23、24、26、27个距离单元的回波信号对第25个距离单元进行气象信号的恢复。采用交叉验证法来评价插值算法的性能，本发明选取平均绝对误差(MAE)来描述插值精度。为了说明本发明方法的有效性，选取多重二次插值进行对比，分别分析多重二次插值与基于距离—多普勒域两维联合插值的性能。

表1雷达系统仿真参数

载频(GHz)	5.5
		脉冲重复频率(Hz)	1000
脉冲个数	128
		雷达高度(m)	1000

表2风力涡轮机仿真参数

叶片长度(m)	26
		桅杆高度(m)	70
轴心到雷达的距离(km)	36
		转叶转速(r/min)	15
旋转面夹角(°)	73.5
		参考叶片初始夹角(°)	90

表3第25个距离单元信号速度估计值的误差

SNR/dB	-3	-2	-1	0	1	2	3	4	5
										MAE/(m/s)	4.153	4.188	4.171	4.187	4.182	4.142	4.179	4.184	4.177

本仿真实验分析了不同信噪比情况下的插值性能，首先由插值方法恢复出被WTC污染距离单元的气象信号，然后估计其恢复出的气象信号速度，计算其平均绝对误差、均值与方差。以下仿真实验结果均依据100次独立蒙特卡洛实验。图4为三种插值方法的误差分析对比图，从图中可以看出，在SNR＝-3dB时，时域多重二次插值的平均绝对误差为2.2m/s，频域多重二次插值的平均绝对误差为0.2m/s，距离—多普勒域两维联合插值的平均绝对误差为0.1m/s。由图4曲线趋可以看出，在低信噪比情况下，距离—多普勒域两维联合插值的平均绝对误差最小，插值精度最高，插值性能得到极大提升，能够实现气象信号的精确恢复。

表4三种插值方法恢复出的气象信号速度估计值的方差

图5为三种插值方法恢复出的气象信号速度估计值的均值，由图5可以看出，时域多重二次插值恢复出的气象信号速度估计值受噪声的影响较大，无法精确估计出气象信号的速度。频域多重二次插值恢复出的气象信号速度估计值与真值误差较大，最终无法收敛于真值。本发明方法恢复出的气象信号速度估计值与真值误差最小，最终收敛于真值。

表4为三种插值方法恢复出的气象信号速度估计值的方差。从表4可以看出，在低信噪比情况下，时域多重二次插值恢复出的气象信号速度估计值的方差最大，其速度估计波动较大，本发明方法恢复出的气象信号速度估计值的方差最小，表明其速度估计波动较小，速度估计值最为精确。

本发明公开了一种基于距离—多普勒域两维联合插值的气象雷达风电场杂波抑制方法，该发明针对低信噪比情况下，首先对同时包含风力涡轮机杂波(WTC)和气象信号的距离单元置零，以该距离单元为中心，两端各选取2个相邻距离单元，通过相干积累提高信噪比。再分别计算出距离域一维插值气象信号估计值和多普勒域一维插值气象信号估计值，最后根据最小二乘法计算出的两维权系数进行两维联合插值。仿真实验结果表明，本方法的插值性能明显优于多重二次插值，可有效提高插值精度，实现气象信号的精确恢复，且运算量小，适合工程应用。

以上所述，仅为本发明中的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉该技术的人在本发明所揭露的技术范围内，可理解想到的变换或替换，都应涵盖在本发明的包含范围之内，因此，本发明的保护范围应该以权利要求书的保护范围为准。

Claims (8)

1.基于两维联合插值的气象雷达风电场杂波抑制方法，其特征在于，该方法的具体步骤如下：

步骤1，对同时包含风力涡轮机杂波WTC和气象信号的距离单元置零，以该距离单元为中心，两端各选取2个相邻距离单元，通过相干积累提高信噪比；

步骤2，对雷达扫描回波信号进行预处理，构建雷达插值初始矩阵；

步骤3，分别计算距离域一维插值气象信号估计值和多普勒域一维插值气象信号估计值；

步骤4，根据最小二乘法计算出的两维权系数，进行两维联合插值，得到两维联合插值气象信号估计值；

步骤5，输出抑制WTC后的信号。

2.根据如权利要求1所述的基于两维联合插值的气象雷达风电场杂波抑制方法，其特征在于，步骤1具体为：

对同时包含WTC和气象信号的第a个距离单元置零，对第a±i个距离单元进行相干积累，通过快速傅立叶变换实现：

式中，i＝1,2，x_a±i(n)、分别为第a±i个距离单元雷达当前扫描与前一次扫描回波信号，X_a±i(k)、分别为x_a±i(n)、经过傅立叶变换后的结果，K为脉冲数。

3.根据如权利要求2所述的基于两维联合插值的气象雷达风电场杂波抑制方法，其特征在于，步骤2具体为：

对第a±i个距离单元雷达当前扫描回波信号进行频谱分析，分别估计出其气象信号频谱峰值对应的多普勒单元；

设各单元气象信号频谱峰值所在的最小多普勒单元序号为m_min、最大多普勒单元序号为m_max，选取雷达两次扫描回波的第m_min-10～m_max+10个多普勒单元为插值样本，将其他多普勒单元置零，以该结果针对对应的X_a±i(k)、进行更新，并构建雷达插值初始矩阵Z、Z₀：

4.根据如权利要求3所述的基于两维联合插值的气象雷达风电场杂波抑制方法，其特征在于，多普勒峰值有如下两种特殊情况：

(1)多普勒峰值出现在零频附近，此时m_min-10＜0，则取第1～m_min、m_max～m_max+10、K-10+m_min～K个多普勒单元为插值样本；

(2)多普勒峰值出现在脉冲重复频率附近，此时m_max+10＞K，则取第1～10-K+m_max、m_min-10～m_min、m_max～K个多普勒单元为插值样本。

5.根据如权利要求3所述的基于两维联合插值的气象雷达风电场杂波抑制方法，其特征在于，步骤3中距离域一维插值气象信号估计值为：

式中，为第a个距离单元的的第k个多普勒单元距离域一维插值气象信号估计值，α_a±i(k)为第a±i个距离单元在第k个多普勒单元的权重因素，径向基函数R_a为第a个距离单元的距离，R_a±i为第a±i个距离单元的距离，b为可调参数，||·||表示范数，

6.根据如权利要求5所述的基于两维联合插值的气象雷达风电场杂波抑制方法，其特征在于，步骤3中多普勒域一维插值气象信号估计值为：

式中，为第a个距离单元的的第k个多普勒单元多普勒域一维插值气象信号估计值，N为被选中距离单元的总数，雷达当前扫描回波信号与前一次扫描回波信号的互相关矩阵

7.根据如权利要求6所述的基于两维联合插值的气象雷达风电场杂波抑制方法，其特征在于，步骤4中两维联合插值气象信号估计值为：

式中，为第a个距离单元的第k个多普勒单元两维联合插值气象信号估计值，

8.根据如权利要求7所述的基于两维联合插值的气象雷达风电场杂波抑制方法，其特征在于，步骤5中输出的抑制WTC后信号为：

式中，为第a个距离单元的第n个脉冲抑制WTC后的气象信号高精度估计值。