自定义函数s_blade
上篇博客 海上风电场对雷达性能的影响——绕射损耗我们讨论了海上风电场对雷达波绕射损耗的影响。在本篇博客中，我们主要对风轮机雷达回波进行建模与仿真。
基于风轮机散射点叠加的模型，研究了风轮机雷达回波信号的仿真问题，考虑了雷达入射波的方位角和俯仰角对回波信号的影响，建立了任意观测点处雷达回波信号模型及其回波仿真曲线，对雷达回波信号的时域、频域和时频域特性进行了分析。
1 风轮机回波模型
在远场情况下，叶片和桅杆可等效为细长的圆柱体，并可将圆柱体沿轴向分割成一系列薄圆片，风轮机回波可等效为一系列薄圆片回波的合成。由于桅杆和轮机舱属于静止目标，其回波易于分析和计算，故主要介绍叶片的回波模型。以叶片轴心为中心，垂直于旋转面的方向为X轴建立如图1所示的坐标系。
图 1 风轮机与雷达几何关系图
由于在博客中输入公式比较繁琐，建模过程可以参考相关文献，或者参考我自己写的word文档
最后得到了风轮机回波表达式：
其中
是风轮机的叶片回波，而
是风轮机的桅杆回波。
M为桅杆的高度；
β为雷达LOS与沿桅杆轴线向上方向(即Z轴方向)的夹角；
L为叶片长度；
N为风轮机的叶片个数；与雷达 LOS 相应的夹角分别为
φ1、
φ2、...φN；λ为雷达的波长。
2 风轮机回波信号的仿真
2.1 对风轮机回波信号的时域分析
对风轮机回波进行时域仿真并分析。根据本课题中给出风轮机特性表和雷达特性表，设置叶片长度L=66m，桅杆高度M=90m，风轮机叶片轴心距雷达的距离为r=10km，叶片数目N=3；设置雷达中心频率
=9GHZ，脉冲宽度τ=1us，脉冲重复频率PRF=1000HZ。另外针对风轮机的额定功率P=5MW以及额定风速w=11.5m/s，推断出风机在额定风速工作下的风机转速范围11—12r/min，在本文取风机转速
=15r/min。另外，为了便于分析，令雷达的方位角α=pi/2，雷达的俯仰角β=pi/2，叶片1和Y轴的夹角θ1=0。MATLAB代码如下，得到的风轮机回波时域波形见图2。
close all
clear all
% 《风轮机雷达回波的仿真与分析》
% 作者：张尺&Sangtro
% 日期：2022.9.20
% -----------雷达和风机的属性----------- %
L = 66; % 风轮机叶片长度 66 m
M = 90; % 桅杆高度 90 m
r = 10000; % 设风轮机叶片轴心距雷达的距离为 10 km
f_rot = 15; % 风机的转速 f_rot = 15 r/ min
num_offan_blades = 3; % 风机叶片数目
c = 3e8; % 光速
f_c = 9000000000; % 雷达频率 f_c = 9 GHz, X波段
tao = 1e-6; % 脉冲宽度 tao = 1 us
PRF = 1000; % 脉冲重复频率 PRF = 1000 Hz
time_windows = 3; % 时窗 time_windows = 3s
angle_a = pi/2; % 雷达的方位角 α
angle_b = pi/2; % 雷达的俯仰角 β
lambda = c/f_c; % 发射信号波长 λ = c / f_c
fs = 30000; % 采样频率
t = 1/fs:1/fs:60/f_rot; % 时间变量 t
% -----------叶片和Y轴的夹角θ，这里采用弧度制----------- %
angle1 = f_rot/60*2*pi*t; % Y轴与叶片1夹角θ1
angle2 = f_rot/60*2*pi*t+2*pi/num_offan_blades; % Y轴与叶片2夹角θ2
angle3 = f_rot/60*2*pi*t+2*pi/num_offan_blades*2; % Y轴与叶片3夹角θ3
% -----------雷达和叶片的夹角φ，计算cosφ----------- %
% 由空间两直线的夹角公式可得，叶片与雷达 LOS 的夹角的余弦值为:
% cosφ = sinβsinαcosθ+ cosβsinα
cos_angle_c1 = sin(angle_b)*sin(angle_a)*cos(angle1)+cos(angle_b)*sin(angle1);
cos_angle_c2 = sin(angle_b)*sin(angle_a)*cos(angle2)+cos(angle_b)*sin(angle2);
cos_angle_c3 = sin(angle_b)*sin(angle_a)*cos(angle3)+cos(angle_b)*sin(angle3);
% -----------雷达的回波信号：叶片的回波信号和桅杆的回波信号----------- %
% 叶片的回波信号sblade
sblade = s_blade(L,t,f_c,r,lambda,cos_angle_c1,cos_angle_c2,cos_angle_c3);
smast = s_mast(M,t,f_c,r,lambda,angle_b); % 桅杆的回波信号smast
swt = sblade + smast; % 风轮机回波信号:swt(t)=sblade(t)+smast(t)
% -----------1.雷达的回波信号的时域分析----------- %
amplitude_swt = abs(swt); % 风轮机回波信号的振幅
% 画图1——风轮机回波信号的时域波形图
figure(1)
plot(t,amplitude_swt,'b')
grid
xlabel ('时间 /s');
ylabel ('振幅 /V');
title('风轮机回波信号时域波形图')自定义函数 s_blade() 的MATLAB代码：
function sblade = s_blade(L,t,f_c,r,lambda,cos_angle_c1,cos_angle_c2,cos_angle_c3)
% 《风轮机雷达回波的仿真与分析》
% 作者：张尺&Sangtro
% 日期：2022.9.20
% 风机叶片的总回波信号
% 变量1：imaginary1 = 2Lcos(φk)/λ
imaginary1 = 2*L*cos_angle_c1/lambda;
imaginary2 = 2*L*cos_angle_c2/lambda;
imaginary3 = 2*L*cos_angle_c3/lambda;
% 变量2：variable_imaginary1 = 2πjLcos(φk)/λ
variable_imaginary1 = pi*imaginary1*sqrt(-1);
variable_imaginary2 = pi*imaginary2*sqrt(-1);
variable_imaginary3 = pi*imaginary3*sqrt(-1);
sblade1 = exp(variable_imaginary1).*sinc(imaginary1)*L;
sblade2 = exp(variable_imaginary2).*sinc(imaginary2)*L;
sblade3 = exp(variable_imaginary3).*sinc(imaginary3)*L;
% % 变量3：载波 carrier_wave 2πj f_c t
% carrier_wave = exp(2*pi*sqrt(-1)*f_c*t);
%
% % 变量4：恒定相位 constant_phase -4πjr/λ
% constant_phase = exp(-4*pi*sqrt(-1)*r/lambda);
%
% sblade1 = sblade1.*carrier_wave*constant_phase;
% sblade2 = sblade2.*carrier_wave*constant_phase;
% sblade3 = sblade3.*carrier_wave*constant_phase;
sblade = sblade1 + sblade2 + sblade3;
return自定义函数 s_mast() 的MATLAB代码：
function smast = s_mast(M,t,f_c,r,lambda,angle_b)
%------桅杆和轮机舱可等效为静止目标，轮机舱的回波较小，忽略不计。
% 《风轮机雷达回波的仿真与分析》
% 作者：张尺&Sangtro
% 日期：2022.9.20
% 风机叶片的总回波信号
imaginary = -2*M*cos(angle_b)/lambda;
variable_imaginary = pi*imaginary*sqrt(-1);
smast = exp(variable_imaginary)*sinc(imaginary)*M;
% % 变量1：载波 carrier_wave 2πj f_c t
% carrier_wave = exp(2*pi*sqrt(-1)*f_c*t);
%
% % 变量2：恒定相位 constant_phase -4πjr/λ
% constant_phase = exp(-4*pi*sqrt(-1)*r/lambda);
%
% smast = smast.*carrier_wave*constant_phase;
return
图2 风轮机回波时域波形图
从图2 可以得到，当雷达波束垂直照射叶片的一个侧边时，回波达到最强，出现一个闪烁，称之为波形峰。随着叶片的转动，波形峰周期性出现。在本次仿真中，把风轮机回波的时间长度定义为叶片旋转一周所用的时间，即的范围是0—60/
(s)。在此时间段，波形峰一共出现了6次，也就意味着雷达波束6次垂直照射风机叶片。为了更清晰的说明问题，做出了风轮机回波和风机叶片位置对照图，如图3 所示。
图3 风轮机回波和风机叶片位置对照图
从图3 中可以看出，在一个旋转周期中，每个叶片都经历了2次被雷达波束垂直照射，而峰值在150~160V之间。例如在波形峰1中，根据风机叶片的位置推算θ1=pi/6、θ2=5pi/6和θ3=3pi/2，求此时回波信号的振幅：
在风轮机回波时域波形图中可清楚地看到风轮机叶片产生的波形峰，任意2个波形峰之间的时间间隔相同，该间隔与叶片的旋转速度直接相关，波形峰的宽度与叶片长度、发射信号波长和叶片旋转速度相关。因而由波形峰之间的时间间隔及波形峰宽度可以推导出叶片旋转速度、叶片长度等相关信息。但是，叶片的旋转特性仅从时域观察是不够的。如果叶片个数未知，仅从风轮机回波时域波形图中无法判断叶片数目的奇偶，也无法判断波形峰信号多普勒谱的正负。为此，需要分析风轮机回波信号的频域特性。
2.2 对风轮机回波信号的频域分析
MATLAB代码如下:
% 《风轮机雷达回波的仿真与分析》
% 作者：张尺&Sangtro
% 日期：2022.9.20
% -----------2.雷达的回波信号的频域分析----------- %
% 对叶片的回波信号进行傅里叶变换
% 这里用的是MATLAB内部函数-快速傅里叶变换
f_sblade = fftshift(fft(sblade,length(sblade)));
% 这里用的是自定义函数，如果采样频率太高，导致内存不读
% f_sblade = dft(sblade,length(sblade));
f = (-length(f_sblade)/2+1:length(f_sblade)/2)*fs/length(f_sblade);
amplitude_sblade = abs(f_sblade);
norm_sblade = (amplitude_sblade - min(amplitude_sblade))/(max(amplitude_sblade)-min(amplitude_sblade));
% 画图2——叶片的回波信号的频域波形图
figure(2)
plot(f,norm_sblade,'b')
grid
xlabel ('频率 /HZ');
ylabel ('幅度 /V');
title('叶片的回波信号的傅里叶变换')自定义函数 dft() 的MATLAB代码：
function [Xk]=dft(xn,N)
% 这里是双边谱
% 作者：张尺&Sangtro
% 日期：2022.9.20
n=-N/2+1:N/2;
k=-N/2+1:N/2;
Xk=xn*(exp(-sqrt(-1)*2*pi/N).^(n'*k)); % n'为n的转置, w = 2π/ N单个叶片产生的频谱将从0频到最大线速度对应的最大多普勒频率，整个连成一片。当叶片迎着雷达射线方向旋转时，产生正的多普勒谱区；当叶片背着雷达射线方向旋转时，将产生负的多普勒谱区。图4 是图2 的回波信号频谱图，其中去除了桅杆的频谱分量，这是因为桅杆静止不动，其频谱分量集中在0频处，且幅度远大于叶片频谱分量的幅度，为了便于分析由叶片旋转产生的多普勒频移，故去除了桅杆的频谱分量。从图4 中可以看到，频谱以0频为中心向左右扩展，呈一定的宽度，对应旋转叶片产生的正多普勒谱区和负多普勒谱区，注意此时纵轴显示的是归一化后的幅度值。
图4 风轮机回波信号频域波形图
由于转速
=20r/min时，叶片的最大线速度vmax=
=103.62m/s，那么最大的多普勒频率为
，从图4 可以看出，仿真结果与理论计算结果一致。但是，从风轮机回波信号的频域特性还是无法得知每个波形峰的频谱，以及正负多普勒是否同时出现等特性。因此需要进行风轮机回波信号的时频域特性分析。
2.3 对风轮机回波信号的时频分析
分析时域可以得到信号随时间变化的信息，分析频域可以得到信号随频率变化的信息，这两者都只能分析时域或频域，而不能同时观察时域和频域。时频分析法提供了时间域和频率域的联合分布信息，清楚地描述了信号频率随时间变化的关系。短时傅里叶变换定义为：
其中 swt(t) 为输入的回波信号；w(m)是窗函数，在本此仿真中利用Hanning窗作为窗函数，X(n,w)是时间和频率的二维函数。MATLAB代码如下, 仿真结果如图5 所示。
% 《风轮机雷达回波的仿真与分析》
% 作者：张尺&Sangtro
% 日期：2022.9.20
% -----------3.雷达的回波信号的短时傅里叶变换----------- %
wlen = 128; % 窗口大小，推荐取 2 的幂次
hop = wlen/4; % hop size 即移动步长，一般要取一个小于 wlen 的数，推荐取 2 的幂次
nfft = 128; % FFT 点数，理论上应该不小于wlen，推荐取 2 的幂次
win = hanning(wlen, 'periodic');
freqrange = 'centered';
[S,F,T] = spectrogram(swt, win, wlen - hop, nfft, fs, freqrange);
% 画图3——雷达的回波信号的短时傅里叶变换
figure(3)
PlotSTFT(T,F,S,win);
title('风轮机回波信号时频波形图');
图5 风轮机回波信号时频分析图
从图5 可以看出，时频分析后会出现多普勒频率闪烁的现象，闪烁出现的周期是叶片转速和叶片个数的函数。在叶片旋转面内，当叶片与雷达LOS垂直时，叶片上所有散射点回波具有相同的相位，叠加后幅度增加，能量急剧上升；当叶片与雷达LOS不垂直时，叶片上各散射点回波的相位不一样，叠加后部分相消，能量下降。所以图5 中多普勒频率闪烁之间的能量微弱很多。3个叶片的叶尖相当于做圆周运动，图5 中3组能量微弱的正弦形状的多普勒频率对应的是叶片叶尖散射点回波。零多普勒频率对应的是风轮机静止部分(桅杆)的回波。多普勒扩展带宽是叶片旋转速度、叶片长度、叶片与雷达LOS夹角的函数。
3 雷达方位角和俯仰角对回波信号影响
风轮机的旋转面会随着风向的改变而自动旋转，旋转到与风向垂直的位置，从而雷达相对于风轮机的方位角和俯仰角也会随着改变。另外，在不同的位置观察风轮机回波信号时，雷达的方位角和俯仰角也将不同。由于叶片的旋转，不同方位角和俯仰角处风轮机叶片相对于雷达的旋转速度不同，故回波信号会发生相应的变化。下面分别分析雷达的方位角和俯仰角对风轮机叶片回波信号的影响。
3.1 方位角对回波信号的影响
令叶片1和Y轴的夹角 θ1=0，当雷达的俯仰角 β=pi/2 时，雷达的方位角α=0、pi/6、pi/3、pi/2时的回波信号时域波形图和频域波形图见图6 , 此时频谱的纵轴依旧显示的是归一化后的幅度值。
(1) α=0
(2) α=pi/6
(3) α=pi/3
(4) α=pi/2
图6 风轮机回波随变化的时域波形图和频域波形图
从图6(a)中可以看出，当俯仰角β=pi/2、方位角α=0时，即波束与叶片旋转面垂直，在叶片旋转过程中，叶片上的所有散射点与雷达之间的距离保持不变，散射点相对于雷达的速度大小为0。此时，所有散射点回波具有相同的相位，叠加后得到如图6(a)所示的结果。当俯仰角β=pi/2，方位角α从0~pi/2变化时，叶片上散射点相对于雷达的速度越来越大，故产生的多普勒谱越来越宽。由此可得，雷达相对于风轮机的方位角影响回波的频谱分布。
3.2 俯仰角对回波信号的影响
当雷达的方位角α=pi/2时，雷达的俯仰角β=0、pi/6、pi/3、pi/2时的回波信号时域波形图和频域波形图见图7, 此时频谱的纵轴依旧显示的是归一化后的幅度值。
(1) β=0
(2) β=pi/6
(3) β=pi/3
(4) β=pi/2
从图7 中可以看出，在雷达的方位角 α=pi/2 的情况下，当俯仰角β从0~pi/2变化时，风轮机回波的频谱图保持不变。在时域图中，波形峰的幅度和宽度不变，但位置发生了改变。这是因为当俯仰角β改变时，叶片上的散射点在整个旋转周期内产生的多普勒谱不变，然而叶片边缘与波束垂直时叶片的位置将发生改变。因此，波形峰出现的位置发生了改变。
参考文献：
[1] 近海风电场障碍下海事雷达图像模拟研究_王树武
[2] 海上风电场对岸基雷达的效能影响_杨守峰
[3] 风机对岸基雷达和船载雷达的影响测试
[4] 风轮机雷达回波的仿真与分析_何炜琨
