基于Matlab模拟用于海况海洋学研究的 X 波段雷达系统（附源码）

目录

一、定义雷达系统参数

二、对海面进行建模

 三、配置雷达收发器

四、生成数据多维数据集

五、处理海面回波

六、总结

七、程序

---

海事雷达系统在充满挑战的动态环境中运行。为了改进对感兴趣目标的检测并评估系统性能，必须了解海面返回的性质。

在本例中，将模拟用于海况海洋学研究的 X 波段雷达系统。雷达系统是一个固定的海上平台。将首先使用 Elfouhaily 光谱模型在海况 4 处生成移动的海面。然后，将从海面生成 IQ 回波，并研究模拟海面信号的统计数据和时频行为。

一、定义雷达系统参数

定义一个距离分辨率约为 30 米的 X 波段雷达系统。使用函数验证范围分辨率是否符合预期。

rng(2021)       % Initialize random number generator

% Radar parameters
freq = 9.39e9;  % Center frequency (Hz)
prf = 500;      % PRF (Hz)
tpd = 100e-6;   % Pulse width (s)
azBw = 2;       % Azimuth beamwidth (deg)
depang = 30;    % Depression angle (deg)
antGain = 45.7; % Antenna gain (dB)
fs = 10e6;      % Sampling frequency (Hz)
bw = 5e6;       % Waveform bandwidth (Hz)
bw2rangeres(bw)

接下来，将海面定义为风速为 4 节、分辨率为 4 米的海况 2。将海面长度设置为 512 米。

% Sea surface
seaState = 4;   % Sea state number 
vw = 19;        % Wind speed (knots)
L = 512;        % Sea surface length (m) 
resSurface = 2; % Resolution of sea surface (m)

% Calculate wavelength and get speed of light
[lambda,c] = freq2wavelen(freq);

% Setup sensor trajectory and simulation times
rdrht = 100;                                   % Radar platform height (m) 
rdrpos = [-L/2 0 rdrht];                       % Radar position [x y z] (m)
numPulses = 1024;                              % Number of pulses 
scenarioUpdateTime = 1/prf;                    % Scenario update time (sec) 
scenarioUpdateRate = prf;                      % Scenario update rate (Hz)
simTime = scenarioUpdateTime.*(numPulses - 1); % Overall simulation time (sec) 

二、对海面进行建模

定义方案。接下来，添加海动。保持属性默认值以使用Elfouhaily频谱。将分辨率设置为 2 米。Elfouhaily模型由海光谱和角扩散函数组成。海谱模型获取物理属性，如风速和获取。

还可以通过将属性设置为来导入自己的自定义模型。文献中的替代海光谱模型包括JONSWAP，Donelan-Pierson和Pierson-Moskovitz光谱。

% Create scenario
scene = radarScenario('UpdateRate',scenarioUpdateRate, ...
    'IsEarthCentered',false,'StopTime',simTime);

% Define Elfouhaily sea spectrum 
seaSpec = seaSpectrum('Resolution',resSurface) % See Elfouhaily reference

现在，选择一个反射率模型。雷达工具箱™为海面提供 9 种不同的反射率模型，涵盖广泛的频率、掠角和海况。星号表示默认模型。在命令窗口中键入有关每个型号的使用和适用掠角的更多信息。海水反射率模型如下。

• APL：半经验模型支持在 1 至 100 GHz 范围内的频率上具有低掠角。 适用于海况 1 至 6。

• GIT：半经验模型支持在 1 至 100 GHz 范围内的频率上具有低掠角。 适用于海况 1 至 6。

• Hybrid：半经验模型支持在 500 MHz 至 35 GHz 频率范围内进行中低掠角。 适用于海况 0 至 5。

• Masuko：支持在 8 至 40 GHz 频率范围内具有低掠角的经验模型。 适用于海况 1 至 6。

• Nathanson：支持在 300 MHz 至 36 GHz 范围内的频率上具有低掠角的经验模型。 适用于海况 0 到 6。

• NRL*：支持在 500 MHz 至 35 GHz 范围内频率范围内具有低掠角的经验模型。 适用于海况 0 至 6。

• RRE：支持在 9 至 10 GHz 频率范围内实现低掠角的数学模型。 适用于海况 1 至 6。

• Sittrop：支持在 8 至 12 GHz 频率范围内具有低掠角的经验模型。 适用于海况 0 至 7。

• TSC：支持在 500 MHz 至 35 GHz 范围内的频率上具有低掠角的经验模型。 适用于海况 0 到 5。

对于此示例，将反射率设置为 GIT（佐治亚理工学院）模型。

% Define reflectivity model
pol = 'V'; % Polarization 
reflectModel = surfaceReflectivity('Sea','Model','GIT','SeaState',seaState,'Polarization',pol)

使用海向雷达方案添加海面。将先前定义的海光谱和反射率模型分配给海面。，和属性用于生成Elfouhaily海光谱。通过对这些特性的深思熟虑的选择和海面的分辨率，可以产生不同波龄的重力和毛细波。

毛细管波是高频、小波长的波，其嵴为圆形，V形槽小于2厘米。重力波是更长的低频波。该物体无法模拟破波，破波是达到临界振幅的波，导致波能转化为湍流动能。风速越大，传递到海面的能量就越浪越大。

取水是风畅通无阻地吹过的水的长度。Elfouhaily模型将获取与反波时代联系起来 

 海面的高程值可以通过使用海面上的方法获得。使用helperSeaSurfacePlot绘制雷达和海面。

% Configure sea surface
knots2mps = 0.514444; % Knots to meters/sec
vw = vw*knots2mps; % Wind speed (m/s)
seaSurf = seaSurface(scene,'SpectralModel',seaSpec,'RadarReflectivity',reflectModel, ...
    'WindSpeed',vw,'WindDirection',0,'Boundary',[-L/2 L/2; -L/2 L/2])

% Plot sea surface and radar 
x = -L/2:resSurface:(L/2 - 1); 
y = -L/2:resSurface:(L/2 - 1); 
[xGrid,yGrid] = meshgrid(x,y);
z = height(seaSurf,[xGrid(:).'; yGrid(:).'],scene.SimulationTime); 
helperSeaSurfacePlot(x,y,z,rdrpos)

 调查海面高程值的统计数据。首先，使用辅助海面CDF计算并绘制累积分布函数。可以使用此图来确定观测值小于或等于特定高程值的概率。例如，90% 的高程将小于或等于约 1 米。95%的海拔将小于或等于约1.5米。如果要更改生成的海面的统计数据以增加海拔，则可以增加风速以向海面输送更多能量。

% CDF
helperSeaSurfaceCDF(z)

接下来，使用助手估计显著波高估计有效波高。估计值是通过取波高的最高 1/3 的平均值获得的，其中波高从波谷到波峰定义，如图所示。

% Significant wave height
actSigHgt = helperEstimateSignificantWaveHeight(x,y,z)

 海况 4 的预期波高在 1.25 到 2.5 米之间。请注意，模拟的海面在预期值的范围内。

expectedSigHgt = [1.25 2.5]; % Sea state 4
actSigHgt >= expectedSigHgt(1) && actSigHgt <= expectedSigHgt(2)

最后，绘制海面高度随时间变化的子集，以可视化由于Elfouhaily海光谱引起的运动。请注意，雷达数据收集的模拟时间仅运行 2 秒，但绘制 10 秒的时间段以更好地了解运动。取消注释帮助程序海面运动绘图以绘制移动的海面。

% Plot sea surface motion 
plotTime = 0:0.5:10; % Plot time (sec)
% helperSeaSurfaceMotionPlot(x,y,seaSurf,plotTime); % Uncomment to plot motion

 三、配置雷达收发器

在本节中，配置雷达系统属性。定义天线和发射的线性调频 （LFM） 波形。将雷达传感器分配给雷达平台。

% Create a radar platform using the trajectory information
rdrplat = platform(scene,'Position',rdrpos);

% Create a radar sensor looking to the East
rdr = radarTransceiver('MountingAngles',[0 depang 0],'NumRepetitions',1);

% Configure the LFM signal of the radar
rdr.Waveform = phased.LinearFMWaveform('SampleRate',fs,'PulseWidth',tpd, ...
    'PRF',prf,'SweepBandwidth',bw);

% Set receiver sample rate and noise figure
rdr.Receiver.SampleRate = fs;
rdr.Receiver.NoiseFigure = 10; 

% Define transmit and receive antenna and corresponding parameters
ant = phased.SincAntennaElement('Beamwidth',azBw);
rdr.TransmitAntenna.OperatingFrequency = freq;
rdr.ReceiveAntenna.OperatingFrequency = freq;
rdr.Transmitter.Gain = antGain;
rdr.Receiver.Gain = antGain;
rdr.TransmitAntenna.Sensor = ant;
rdr.ReceiveAntenna.Sensor = ant;

% Add radar to radar platform
rdrplat.Sensors = rdr;

四、生成数据多维数据集

现在场景和雷达系统已经定义，使用该方法从海面生成返回，并使用该方法收集 IQ 数据，两者都驻留在对象中。默认情况下，模拟主瓣中的杂波返回。

使用帮助程序帮助程序 PlotRawIQ和帮助程序更新绘图原始 IQ绘制未处理的原始 IQ 数据。

% Collect clutter returns with the clutterGenerator
clutterGenerator(scene,rdr); 

% Run the scenario
numSamples = 1/prf*fs;
maxRange = 20e3; 
Trng = (0:1/fs:(numSamples-1)/fs);
rngGrid = [0 time2range(Trng(2:end),c)];
[~,idxTruncate] = min(abs(rngGrid - maxRange));
iqsig = zeros(idxTruncate,numPulses);
ii = 1; 
hRaw = helperPlotRawIQ(iqsig);
while advance(scene)
    tmp = receive(scene); % nsamp x 1
    iqsig(:,ii) = tmp{1}(1:idxTruncate); 
    if mod(ii,128) == 0
        helperUpdatePlotRawIQ(hRaw,iqsig);
    end
    ii = ii + 1;
end
helperUpdatePlotRawIQ(hRaw,iqsig);

五、处理海面回波

用于脉冲压缩返回的 IQ 数据。使用帮助程序范围时间图可视化范围时间行为。如果海面是静态的，会在绘图中看到笔直的水平线，但海面返回表现出活力。在给定的范围内，幅度随时间变化。请注意，由于几何形状（雷达高度、波束宽度和俯角），回波发生在一小部分范围内。

% Pulse compress
matchingCoeff = getMatchedFilter(rdr.Waveform);
rngresp = phased.RangeResponse('RangeMethod', 'Matched filter', ...
    'PropagationSpeed',c,'SampleRate',fs);
[pcResp,rngGrid] = rngresp(iqsig,matchingCoeff); 

% Plot
pcsigMagdB = mag2db(abs(pcResp));
[maxVal,maxIdx] = max(pcsigMagdB(:));
pcsigMagdB = pcsigMagdB - maxVal;
helperRangeTimePlot(rngGrid,prf,pcsigMagdB);

 使用最大值的范围指数，并使用帮助程序MagTimePlot可视化此范围箱的幅度与时间的关系。

% Plot magnitude versus time
[idxRange,~] = ind2sub(size(pcsigMagdB),maxIdx); 
helperMagTimePlot(pcsigMagdB(idxRange,:),numPulses,prf,rngGrid(idxRange));

 接下来，使用helperSTFTPlot生成相同范围箱的时频图。使用短时傅里叶变换函数生成绘图。

% STFT
[S,F,T] = stft(pcResp(idxRange,:),scenarioUpdateRate);
helperSTFTPlot(S,F,T,lambda,rngGrid(idxRange));

 请注意，由于海面的运动，范围箱的检测速度会随时间而变化。在较长的仿真时间内，可以检测到周期性。 最后，使用helperHistPlot，从 180 到 210 米范围内的数据量级值形成直方图。请注意，直方图的形状类似于威布尔分布。如果您有统计和机器学习工具箱™，则可以使用该函数从数据量级中获取 Weibull 分布的参数。

% Look at a subset of data in range and convert to decibel scale
[~,idxMin] = min(abs(rngGrid - 180)); 
[~,idxMax] = min(abs(rngGrid - 210)); 
pcsigMagdB = mag2db(abs(pcResp(idxMin:idxMax,:)));

% Remove any inf values
pcsigMagdB(isinf(pcsigMagdB(:))) = []; 

% Shift values to be positive
pcsigMagdB = pcsigMagdB(:) - min(pcsigMagdB(:)) + eps; 

% Weibull parameters
% Note: These values can be obtained if you use fitdist(pcsigMagdB,'weibull')
scale = 37.8589;
shape = 7.80291;

% Plot histogram with theoretical PDF overlay
helperHistPlot(pcsigMagdB,scale,shape);

 六、总结

在此示例中，学习了如何：

• 构建雷达场景，

• 模拟移动的海面，

• 调查海面统计数据和行为，

• 对模拟的 IQ 返回执行时频处理。

这个例子可以很容易地扩展到包括目标，以支持海上监视和雷达成像等应用。

七、程序

使用Matlab R2022b版本，点击打开。（版本过低，运行该程序可能会报错）

打开下面的“RadarReturnsFromMovingSeaSurfacesExample.mlx”文件，点击运行，就可以看到上述效果。

关注下面公众号，点击文章《基于Matlab模拟用于海况海洋学研究的X波段雷达系统》​​​​​​​，获取下载链接。