2022年华为杯数学建模竞赛A题论文和代码

1. 任务背景和问题重述

1 . 1 任务背景

自动驾驶可有效提高交通效率和交通安全，已成为智能交通发展的重要趋势[1]。自动驾驶车辆上路行驶的前提是安全可靠，其首要基础任务是在移动环境中对车辆周围交通目标精确测量定位。目前，自动驾驶移动场景定位主要采用可见光、激光、毫米波雷达等传感器。与其他传感器相比，毫米波雷达能够不受太阳光、雨雪雾等复杂天气条件的影响，可对交通环境进行全天时全天候的测量感知，是实现全自动驾驶最为重要的传感器之一[2]。

毫米波雷达通常采用发射调频连续波(Frequency modulated continuous wave, FMCW)信号实现移动场景目标的定位。具体而言，雷达按照一定周期重复发射FMCW Chrip信号，对目标回波在快时间域内进行傅立叶变换得到目标的径向距离位置，同时对相邻回波在慢时间域内进行相参处理得到目标的相对移动速度。在目标的方位角度定位方面，雷达采用时分复用的多输入多输出天线发射与接收信号，形成一个大孔径的虚拟天线阵列。

以图1-1为例，间距为 2λ 的发射两个阵元依次发射信号，再由间距为 0.5λ的四个阵元接收，进而形成间距为0.5λ的八个阵元虚拟天线阵列，其中λ为雷达发射波的中心波长。在以雷达为中心原点的极坐标系中，假设某一目标的坐标为(R, θ),虚拟天线阵元的间距为d,阵元总数量为Na,阵列总长度为L。通过几何关系可知，对于同一个目标，其第 n 个阵元接收信号的方位相位变化项为ej2π λ ndsinθ。实际上，方位空间存在多个目标，假设目标总数量为M,则FMCW雷

达的第n个阵元接收信号可表达为

i=1

式中t为一个Chirp周期内的快时间，c为电磁波传播的速度，ai 与θi 分别为是第i个目标的强度与方位角度，Ri 为第i个目标的径向距离。

          由上式可知，雷达目标定位是对目标距离和角度估计的问题。目标的距离可在快时间域内进行傅立叶变换得到，距离分辨率与信号的带宽B成反比，为 2B。          c由于毫米波 FMCW 雷达信号带宽大，因此距离分辨率较高 [3]。在方位角度上，雷达虽然通过了虚拟天线阵列对目标的方位角度进行了多次空间采样，但其角度分辨率为 Lcosθ。当天线阵元间隔为 λ 2 及目标方位角度为零时，雷达检测目标具有最佳方位角度分辨率为 Na 弧度[4]。

在自动驾驶等移动场景中，提升目标的定位分辨率对自动驾驶安全具有重要意义。本题所指的目标超分辨是指，在多物体信号叠加时的最大化物体定位精度。基于毫米波FMCW雷达对移动场景超分辨率定位，需要设计鲁棒的低运算复杂度在线算法，以实时超分辨定位到物体。

目前FMCW雷达为了减小强目标信号的旁瓣对微弱信号的影响，现有方法通过加对称窗函数[5]，再对距离和方位向进行两维快速傅立叶变换（Fast Fourier Transformation，FFT)来进行目标的距离和角度估计。该方法运算简单、复杂度低，但加窗后会带来目标能量的主瓣宽度扩展，进而导致目标的分辨率下降，不能实现目标的超分辨。传统的谱估计算法有最大似然估计方法、最小方差无失真响应（Minimum variance distortionless response，MVDR)）、多重信号分类(Multiple signal classification，MUSIC）算法、ESPRIT（Estimation of Signal Parameters via Rotational Invariance Technique）、基于最大熵方法等[5]。这些谱估计方法运算量大，比如多重信号分类MUSIC算法通过信号空间和噪声空间的估计，通过两者的正交特性实现目标的超分辨定位。此外，作为稀疏信号处理中的一类，现有的压缩感知算法利用了空间物体分布的稀疏性，可以有效提升分辨率，但压缩感知贪婪的求解算法不利于低复杂度算法的开发[6]。

          综上所述，面对日益增长的目标超分辨定位需求，如何通过建模并设计相应的算法提高分辨率是一个亟待解决的问题。

1 . 2 问题重述

本文试图解决移动场景下目标超分辨定位问题，对题目提供对数据进行处理，建立精确的定位模型，同时通过推导开发数学模型方法自适应噪声环境，最终在噪声污染条件下对目标超分辨定位和动态轨迹估计。具体地，本文需要完成以下四个任务：

任务一、无噪条件下的目标超分辨定位。

          针对提供的无噪声仿真数据，建立定位模型，计算出物体相对位置，并以二维极坐标图（横坐标表示距离，纵坐标表示角度）展示。

2. 模型假设和符号说明

2 . 1 模型假设

假设一、雷达天线与目标之间的相对距离满足天线远场测量条件。

          假设二、雷达系统参数满足实验条件需求，对所测的目标不会产生距离模糊和速度模糊。

假设三、雷达天线阵元之间不存在串扰情况，隔离度良好。问题二与问题四的天线噪声是由自身通道或者天线老化等产生的,从而导致天线阵元数据受到噪声污染。

2 . 2 符号说明

7

进一步采用分离变量法，上式可等效为：

u-x= 0

通过增扩拉格朗日方法，可建立以下目标优化函数模型

上式中µ,η1 与η2 为模型的参数。

上式可通过引入一个中间向量d进行交替求解,求解表达式为

xk+1 =dk +vk

上式中soft为软域值操作运算，其定义为

          由求解过程公式(11)可知，该算法只有一步存在DFT变换，因此运算量小。另外软域值操作运算可保持信号的相位信息，有利用阵列信号相位重建恢复。

最终可得到超分辨增强后的雷达方位角度信号

3 . 4 模型与算法结果

首先考虑模型在迭代过程收敛方面，因为L1正则化的最小二乘法模型处理的信号为复正弦信号，因此该算法可快速收敛。图3-2为处理本任务提供的数据收敛过程。

为了验证提出方法的有效性，首先模拟仿真了 100 个阵元的一维天线阵列信号对比实验。利用点数为256的DFT对雷达的方位信号进行构建DFT变换矩阵，结果如图3-3所示。由对比结果可知，提出的方法有效增强了目标方位分辨率。

接下来展示对提供数据的处理结果。提供的数据的天线阵元数量为 86 个，模型算法处理过程中DFT点数设置为512个。图3-4为模型与算法结果对比，左边一列是整体视野范围内雷达数据处理后的图，右边一列为左边目标的细节放大。由图可知，提出的两种方法都可把原本无法分辨的目标分离，基于L1正则化的最小二乘法模型结果优于基于DFT插值的超分辨处理方法。

az_sig_dft=rafData(:,idx);
figure;plot(imag(az_sig_dft));xlim([0 86]);xlabel('方位阵元采样点');ylabel('信号强度');

az_sig_dft=fftshift(fft(az_sig_dft,AzNfft));az_sig_dft=az_sig_dft/max(abs(az_sig_dft));

% Define functions (Matlab function handles) Trans_num=AzNfft;
H = @(x) InvTransformation(x,Na,Trans_num);HT = @(x) Transformation(x,Na,Trans_num);

p = Trans_num; % 变换点数

Nit = 20;                       % 迭代次数

mu =1e-1;                % 模型参数

%方法2：基于L1正则最小二乘法的方位向超分辨算法处理

y=rafData(:,idx);
% [coef, cost] = mimo_anttena_salsa(y, H, HT, p, mu, Nit);

lambda = 7;
[coef, cost] = mimo_anttena_denoising_salsa(y, H, HT, p, lambda,   mu, Nit);

err = y - InvTransformation(coef,Na,Trans_num);err_max = max(abs(err));

lsz=2;

figure;plot(cost,'linewidth',lsz);

xlabel('迭代次数');ylabel('目标函数值');

ImageFormat(16);

az_sig_ls=fftshift(abs(coef));

az_sig_ls=az_sig_ls/max(abs(az_sig_ls));

figure; plot(azz,20*log10(abs(az_sig_dft)),'linewidth',lsz);hold on;plot(azz,20*log10(abs(az_sig_ls)),':','linewidth',lsz);xlabel('方位角度 (°)');ylabel('相对强度 (dB)');
legend('DFT 插值', '基于L1正则化的最小二乘法');

xlim([-20 20]);ylim([-40 1]);

ImageFormat(16);

data2=data*0+eps;

data2(:,119)=fftshift(coef);

a_data=(abs(fftshift(coef)));

figure;plot(a_data);

% 方位目标峰值检测

num_train=27.4;

num_guard=0;

rate_fa=1e-1;

num_train_half = round(num_train / 2);

num_guard_half = round(num_guard / 2);
num_side = num_train_half + num_guard_half;
alpha = num_train*(rate_fa.^(-1/num_train) - 1);

peak_idx = 0;

pk=0;

for k=num_side : AzNfft - num_side

tem=a_data(k-num_side+1:k+num_side);

[M,Ix] = max(tem,[],1);

if (k ~= k-num_side+Ix)

sum1 = sum(a_data);
sum2 =sum(a_data(k-num_guard_half+1:k+num_guard_half));p_noise = (sum1 - sum2) / num_train ;

threshold = alpha * p_noise;

if a_data(k)> threshold

pk=pk+1;

peak_idx(pk)=k;

end

end

end

for k=1:pk

theta(k)=azz(peak_idx(k))/180*pi;

rho(k)=Ridx(idx);

end

figure;%雷达目标信息极坐标显示

apo = polaraxes;
polarplot(theta,rho,'o','MarkerSize',3,'MarkerFaceColor','b');thetalim([-20 20]);ImageFormat(16);

figure;%雷达超分辨处理后数据极坐标显示

PloarPlot(Ridx,azz,data2);ImageFormat(16);

%程序名称：Q2_main.m

%任务二：含噪雷达数据MIMO定位模型
%该程序根据含噪雷达数据MIMO定位模型，计算出物体相对位置，结果以二维极坐标显示