TDM-MIMO FCMW 雷达的相位补偿
TDM-MIMO 下的 FCMW 雷达的相位补偿
像是 2T4R TDM MIMO 雷达，两个 Tx 在时间上轮流发射 chirp，虚拟出 1T8R 的天线阵列，以获得更高的角度分辨率。然而，Tx2 形成的四根虚拟天线相对于 Tx1 形成的四根虚拟天线将会产生一个普勒相偏 \(\Delta \varphi\)。角度估计前须进行相位矫正，才能保证角度估计的正确性。
多普勒相偏
速度为 \(V_r\) 的目标产生的多普勒频偏 \(\Delta f\) 和相偏 \(\Delta \varphi\)：
\[\Delta f=\frac{2V_rf_c}{c}=\frac{2V_r}{\lambda}
\]
\[\Delta \varphi=2\pi\cdot \Delta f \cdot T_c=\frac{4\pi V_r f_c T_c}{c}=\frac{4\pi V_r T_c}{\lambda}
\]
- \(f_c\)，雷达中心频率
- \(\lambda\)，雷达波长
- \(T_c\)，chirp 时长
根据目标速度进行相位补偿
经过 Range-Doppler 处理后，速度信息变为了离散的数字。假设速度维度已由 fftshift() 处理，即多普勒单元索引 \(n\) 的范围在 \((-\frac{N}{2},\frac{N}{2})\)，其中 \(N\) 是 Doppler-FFT 的长度，则 \(V_r\) 的转换如下：
\[\begin{split}
\Delta \varphi
&=\frac{4\pi V_r T_c}{\lambda}\\
&=\frac{4\pi \cdot n \Delta V_r\cdot T_c}{\lambda}\\
&=\frac{4\pi T_c}{\lambda}\cdot n \frac{\lambda}{2NT_c}\\
&=\frac{2\pi n}{N}\\
\end{split}
\]
假设在 \(M\times N\) 根虚天线上对应目标的 2D-FFT 结果序列为 \(X(m,n)\)，其中 \(m\) 表示这个 2D-FFT 结果对应的发射天线编号（从 0 编号），\(n\) 表示这个 2D-FFT 结果对应的接收天线编号，那么相偏补偿：
\[X'(m,n)=X(m,n)\times e^{-jm\Delta \varphi \cdot \frac{M}{2}}
\]
换句话说，从第二根发射天线开始就要进行补偿
Python 代码与实验
假设该雷达系统是 2T8R， Range-Doppler 处理后的数据 range_doppler 维度为 [N_rx*N_tx, doppler_bin, range_bin]，且 range_doppler[:8] 是第一个 Tx 的数据、range_doppler[8:] 是第二个 Tx 的数据，则补偿代码如下：
# 获取 dop_bin 长度
dop_bin = np.arange(range_doppler.shape[-2])
# 计算相位差
delta_phi = np.exp(-2j * np.pi * (dop_bin - len(dop_bin) / 2) / len(dop_bin)) / N_tx
# 进行补偿
range_doppler[8:, :, :] *= delta_phi[None, :, None]
仿真的结果说明方法很有效。
完整实验代码示例：使用 4T4R 系统
使用 radarsimpy 进行雷达系统仿真。以下代码写得有些乱，供参考。
创建雷达系统
import numpy as np
from radarsimpy import Radar, Transmitter, Receiver
from scipy.constants import speed_of_light
wavelength = speed_of_light / 60.5e9
N_tx = 4
N_rx = 4
tx = []
tx.append(dict(location=(0, -N_rx / 2 * wavelength * 3, 0),delay=0))
tx.append(dict(location=(0, -N_rx / 2 * wavelength * 2, 0),delay=20e-6))
tx.append(dict(location=(0, -N_rx / 2 * wavelength, 0),delay=40e-6))
tx.append(dict(location=(0, 0, 0),delay=60e-6))
# 间距为 λ/2
rx = [dict(location=(0, wavelength / 2 * idx, 0)) for idx in range(N_rx)]
tx = Transmitter(f=[61e9, 60e9], t=16e-6, tx_power=15, prp=80e-6, pulses=512, channels=tx)
rx = Receiver(fs=20e6, noise_figure=8, rf_gain=20, load_resistor=500, baseband_gain=30, channels=rx)
radar = Radar(transmitter=tx, receiver=rx)
显示阵列示意图。
import plotly.graph_objs as go
fig = go.Figure()
fig.add_trace(
go.Scatter(
x=radar.radar_prop["transmitter"].txchannel_prop["locations"][:, 1] / wavelength,
y=radar.radar_prop["transmitter"].txchannel_prop["locations"][:, 2] / wavelength,
mode='markers',
name='Transmitter',
opacity=0.7,
marker=dict(size=10),
)
)
fig.add_trace(
go.Scatter(
x=radar.radar_prop["receiver"].rxchannel_prop["locations"][:, 1] / wavelength,
y=radar.radar_prop["receiver"].rxchannel_prop["locations"][:, 2] / wavelength,
mode='markers',
opacity=1,
name='Receiver',
)
)
fig.update_layout(
title='Array configuration',
height=400,
width=700,
xaxis=dict(title='y (λ)'),
yaxis=dict(title='z (λ)', scaleanchor="x", scaleratio=1),
)
fig.show()
设置目标，产生回波数据
true_theta = [30]
speed = [12]
ranges = [30]
targets = []
for theta, v, r in zip(true_theta, speed, ranges):
target = dict(
location=(r * np.cos(np.radians(theta)), r * np.sin(np.radians(theta)), 0),
speed=(v * np.cos(np.radians(theta)), v * np.sin(np.radians(theta)), 0),
rcs=10,
phase=0,
)
targets.append(target)
from radarsimpy.simulator import simc
data = simc(radar, targets, noise=False)
timestamp = data['timestamp']
baseband = data['baseband']
print(baseband.shape) # (16, 512, 320)
Range-Doppler
首先编写所需函数。
import numpy as np
def range_fft(X, n=None, window=None):
if window is not None:
shape = [1] * X.ndim
shape[-1] = X.shape[-1]
X = X * window.reshape(shape)
output = np.fft.fft(X, n, axis=-1)
return output
def doppler_fft(X, n=None, window=None):
if window is not None:
shape = [1] * X.ndim
shape[-2] = X.shape[-2]
X = X * window.reshape(shape)
output = np.fft.fftshift(np.fft.fft(X, n, axis=-2), axes=-2)
return output
def angle_fft(X, n=None, window=None):
if window is not None:
shape = [1] * X.ndim
shape[0] = X.shape[0]
X = X * window.reshape(shape)
output = np.fft.fftshift(np.fft.fft(X, n, axis=0), axes=0)
return output
然后进行 Range-Doppler。
from scipy import signal
range_window = signal.windows.chebwin(radar.sample_prop["samples_per_pulse"], at=80)
dop_window = signal.windows.chebwin(radar.radar_prop["transmitter"].waveform_prop["pulses"], at=60)
ranged = range_fft(baseband.conj(), window=range_window)
range_doppler = doppler_fft(ranged, window=dop_window)
绘制图像。
from scipy.constants import speed_of_light
max_range = ( # 计算最大探测距离
speed_of_light # 光速
* radar.radar_prop["receiver"].bb_prop["fs"] # 采样率
* radar.radar_prop["transmitter"].waveform_prop["pulse_length"] # 脉冲时长 Tc
/ radar.radar_prop["transmitter"].waveform_prop["bandwidth"]
/ 2
)
range_axis = np.linspace(0, max_range, radar.sample_prop["samples_per_pulse"], endpoint=False)
doppler_axis = (
np.linspace(
-1 / radar.radar_prop["transmitter"].waveform_prop["prp"][0] / 2,
1 / radar.radar_prop["transmitter"].waveform_prop["prp"][0] / 2,
radar.radar_prop["transmitter"].waveform_prop["pulses"],
endpoint=False,
)
* speed_of_light
/ 60.5e9
/ 2
)
fig = go.Figure()
fig.add_trace(
go.Heatmap(
z=20 * np.log10(np.mean(np.abs(range_doppler), axis=0)).T,
x=doppler_axis,
y=range_axis,
colorscale='Viridis',
zmin=-40,
zmax=0,
)
)
fig.update_layout(
title='Range-Doppler map',
height=400,
width=700,
xaxis=dict(title='Doppler (m/s)'),
yaxis=dict(title='Range (m)'),
)
fig.show()
相位补偿与 Angle-FFT
进行相位补偿。
speed_bin = np.arange(range_doppler.shape[-2])
delta_phi = -2j * np.pi * (speed_bin - len(speed_bin) / 2) / len(speed_bin) / N_tx
range_doppler_2 = range_doppler.copy().reshape(
N_tx, N_rx, range_doppler.shape[-2], range_doppler.shape[-1]
)
range_doppler_2[1, :, :, :] *= np.exp(delta_phi[None, :, None])
range_doppler_2[2, :, :, :] *= np.exp(delta_phi[None, :, None] * 2)
range_doppler_2[3, :, :, :] *= np.exp(delta_phi[None, :, None] * 3)
range_doppler_2 = range_doppler_2.reshape(N_tx * N_rx, range_doppler.shape[-2], range_doppler.shape[-1])
进行 Angle-FFT，绘制对比图像。
angled_2 = angle_fft(range_doppler_2, n=128)
det_idx = [np.argmax(np.mean(np.abs(range_doppler[:, 0, :]), axis=0))]
dop_idx = [np.argmax(np.mean(np.abs(range_doppler[:, :, det_idx[0]]), axis=0))]
##############################################
target_zero = dict(
location=(
ranges[0] * np.cos(np.radians(true_theta[0])),
ranges[0] * np.sin(np.radians(true_theta[0])),
0,
),
speed=(0, 0, 0),
rcs=10,
phase=0,
)
baseband_zero = simc(radar, [target_zero], noise=False)['baseband']
ranged_zero = range_fft(baseband_zero.conj(), window=range_window)
range_doppler_zero = doppler_fft(ranged_zero, window=dop_window)
angled_zero = angle_fft(range_doppler_zero, n=128)
##############################################
angled = angle_fft(range_doppler, n=128)
##############################################
azimuth = np.linspace(-1, 1, 128)
azimuth = np.arcsin(azimuth) / np.pi * 180
fig = go.Figure()
angle_data_zero = angled_zero[:, 256, det_idx[0]]
angle_data_without = angled[:, dop_idx[0], det_idx[0]]
angle_data_2 = angled_2[:, dop_idx[0], det_idx[0]]
angle_data_zero = 20 * np.log10(np.abs(angle_data_zero) + 1e-6)
angle_data_without = 20 * np.log10(np.abs(angle_data_without) + 1e-6)
angle_data_2 = 20 * np.log10(np.abs(angle_data_2) + 1e-6)
fig.add_trace(go.Scatter(x=azimuth, y=angle_data_zero, name='速度为 0'))
fig.add_trace(go.Scatter(x=azimuth, y=angle_data_without, name='有速度但不补偿'))
fig.add_trace(go.Scatter(x=azimuth, y=angle_data_2, name='有速度且补偿'))
fig.update_layout(
height=400,
width=700,
title='某一个 range 上的角度估计',
)
print(f"速度为 0：{azimuth[np.argmax(angle_data_zero)]}")
print(f"有速度但不补偿：{azimuth[np.argmax(angle_data_without)]}")
print(f"有速度且补偿：{azimuth[np.argmax(angle_data_2)]}")
fig.show()
参考来源
- “MIMO Radar”，https://www.ti.com/lit/an/swra554a/swra554a.pdf
- Ray，“干货|利用MATLAB实现FMCW MIMO雷达的超分辨测角”，https://zhuanlan.zhihu.com/p/44359309
- 调皮连续波，“雷达仿真 | FMCW TDMA-MIMO毫米波雷达信号处理仿真”，https://zhuanlan.zhihu.com/p/576353487
- “基于 AWR1642 汽车雷达的速度扩展算法研究”，https://www.ti.com.cn/cn/lit/an/zhca901/zhca901.pdf
- Zach，“DoA Estimation”，https://radarsimx.com/2022/12/12/doa-estimation/
浙公网安备 33010602011771号