一、HRRP技术概述与应用背景
高分辨率距离像(High Resolution Range Profile, HRRP)是雷达信号处理中的核心技术之一,它通过对雷达回波信号进行脉冲压缩和相干处理,获取目标在距离维上的散射点分布特征。HRRP具有以下重要特点:
- 高分辨率 :能够区分相距较近的散射中心
- 特征稳定性 :为目标识别提供稳定的物理特征
- 计算效率 :相比二维成像,计算复杂度较低
在军事和民用领域,HRRP广泛应用于:
- 目标识别与分类
- 导弹制导
- 自动目标识别(ATR)系统
- 雷达导引头信号处理
二、C++实现的HRRP处理流程分析
原始C++代码(Keystone.cpp)实现了一个完整的HRRP处理链,主要包含以下核心模块:
2.1 Keystone变换模块
void Keystone_interp(float **data_r,float **data_i,int fwn,int jln, float f0,float fs,float B,int moshuzhouqishu)
- 距离向FFT变换
- 方位向插值补偿(解决距离徙动问题)
- 距离向IFFT恢复时域信号
2.2 插值算法实现
void stoltsun(float *a_in, float *a_value, int a_size, float *b_in, int b_size, float *b_value)
- 三次样条插值算法
- 边界处理和无效点处理
- 内存动态分配管理
2.3 多模糊数处理
void Gethrrp_v2_multimohu_v2(float **dr,float **di,int hangmc,int liejl,float f0,float fs,float B,int mohuhalfshu, float PRF,float vdengxiao,float mubiaojuli)
- 遍历多个模糊数候选
- 基于熵最小化原则选择最优结果
- 相位补偿和运动补偿
三、MindSpore移植架构设计
将原有代码重构为基于计算图的MindSpore Cell,实现雷达信号处理流程模块化。
3.1 数据读取
在Python中使用MindSpore时,我们可以使用NumPy的fromfile接口读二进制格式文件,因此在Python代码开头需要导入NumPy库。
示例代码:
def read_hrrp_data(filename, maichongshu, jln, rfftn):
"""
读取HRRP二进制数据文件并处理
"""
with open(filename, 'rb') as fp:
expected_size = maichongshu * jln * 2
data = np.fromfile(fp, dtype=np.float32)
if len(data) != expected_size:
raise ValueError(f"数据长度错误:期望{expected_size}个点,实际读取{len(data)}个点")
data = data.reshape(maichongshu, jln, 2)
datar = data[:, :, 0]
datai = data[:, :, 1]
# 创建填充0的完整数组
datar_full = np.zeros((maichongshu, rfftn), dtype=np.float32)
datai_full = np.zeros((maichongshu, rfftn), dtype=np.float32)
datar_full[:, :jln] = datar
datai_full[:, :jln] = datai
return datar_full, datai_full
读到的数据需要转换为复数形式,以便后续处理。
datar, datai = read_hrrp_data("hrrp_400_512.dat", maichongshu, jln, rfftn)
input_data = np.zeros(datar.shape, dtype=np.complex64)
input_data = datar + 1j * datai
3.2 数据预处理
从算法整体分析,数据读取后到核心计算之前的步骤,主要是对相位补偿因子计算,Keystone变换参数预计算,泰勒窗函数预计算,循环移位索引预计算。这部分都是核心计算的前期准备,建议将这部分代码封装在__init__函数中完成,不放在construct函数中可以避免额外的开销,当实例化一个类时自动触发一次__init__函数。 示例代码:
class Hrrp(nn.Cell):
def __init__(self, hangmc, liejl, f0, fs, B, mohuhalfshu, PRF, vdengxiao, mubiaojuli):
super(Hrrp, self).__init__()
self.hangmc = hangmc
self.liejl = liejl
self.f0 = f0
self.fs = fs
self.B = B
self.mohuhalfshu = mohuhalfshu
self.PRF = PRF
self.vdengxiao = vdengxiao
self.mubiaojuli = mubiaojuli
self.fft = mr.FFT()
self.ifft = mr.IFFT()
self.fftshift = mr.FFTShift()
self.abs = mr.ComplexAbs()
self.stoltsun = mr.Stoltsun(dim=0)
self.mul = ops.Mul()
self.exp = ops.Exp()
self.reduce_sum = ops.ReduceSum(keep_dims=True)
self.gather = ops.Gather()
wl = LC / self.f0
Ka = 2.0 * self.vdengxiao**2 / (wl * self.mubiaojuli)
mctime = (np.arange(self.hangmc) - self.hangmc * 0.5) / self.PRF
xiangwei = np.pi * Ka * mctime**2
self.phase_comp = np.exp(1j * xiangwei).astype(np.complex64)
self.phase_comp = self.phase_comp.reshape(self.hangmc, 1)
self.phase_comp = ms.Tensor(self.phase_comp)
fwn = self.hangmc
self.fwn = fwn
rnfft = Get2intm(self.liejl)
self.rnfft = rnfft
fwn_8 = (fwn + 7) // 8 * 8 # 8字节对齐
self.fwn_8 = fwn_8
self.xold = (np.arange(fwn_8, dtype=np.float32) - fwn / 2.0)[:, np.newaxis]
self.xold = np.broadcast_to(self.xold, (fwn_8, rnfft))
i_arr = np.arange(rnfft, dtype=np.float32)
sigma_arr = (f0 + (i_arr - rnfft * 0.5) * fs / rnfft) / f0
self.xnew = self.xold * sigma_arr
self.xold = ms.Tensor(self.xold)
self.xnew = ms.Tensor(self.xnew)
self.mchtr = self.TaiLeWindow(self.hangmc, 50)
window_2d = self.mchtr[:, np.newaxis] # 转换为列向量便于广播
self.window_2d = ms.Tensor(window_2d, dtype=ms.float32)
self.j_constant = ms.Parameter(ms.Tensor(1.0j, dtype=ms.complex64), name="j_constant")
self._precompute_cirshift_indices()
def _precompute_cirshift_indices(self):
rnfft = self.rnfft
shift = self.rnfft // 2
indices = np.arange(rnfft, dtype=np.int32)
indices_shifted = (indices - shift) % rnfft
self.cirshift_indices = ms.Tensor(indices_shifted, dtype=ms.int32)
def construct(self, input_data):
...
3.3 核心算法
步骤分解:
-
相位补偿:首先对输入数据(复数形式)进行相位补偿,乘以预计算好的相位补偿因子(用于校正二次相位误差)。
-
多模糊数处理:计算需要搜索的模糊数总数(mohushu = 2 * mohuhalfshu + 1),然后对每一个模糊数进行以下操作:
a. 计算当前模糊数对应的偏移量(offset = mh - mohuhalfshu),这个偏移量将用于Keystone变换中的相位补偿。
b. 进行Keystone变换(Keystone_interp),这是HRRP算法的核心,用于校正距离徙动。
c. 对Keystone变换后的数据,进行HRRP生成(Gethrrp),得到每一距离单元的幅度和相位信息,同时返回处理后的数据(result)和幅度数据(hrrp)。
d. 计算该HRRP的熵(Gethrrp_v2shang),熵值越小表示能量越集中,成像质量越好。
e. 将当前模糊数对应的结果(result)和熵值(shang_value)分别保存到列表results和shang中。
-
选择最优结果:将所有模糊数对应的结果堆叠起来(使用ops.cat),然后找到熵最小的那个索引(xuhao = ops.argmin(shang)),并取出该索引对应的结果(min_result = results[xuhao])。
-
更新输出:将最优结果的第一行(即HRRP处理后的第一行,包含了最优的相位和幅度信息)与原始输入数据的第二行到最后一行拼接起来,形成最终的输出。这里之所以只替换第一行,是因为在HRRP处理中,我们只关心每个距离单元的最强散射点,而原始数据的第一行被用来存储这个信息。
示例代码如下:
class Hrrp(nn.Cell):
def __init__(self, hangmc, liejl, f0, fs, B, mohuhalfshu, PRF, vdengxiao, mubiaojuli):
super(Hrrp, self).__init__()
self.hangmc = hangmc
self.f0 = f0
self.fs = fs
self.B = B
self.mohuhalfshu = mohuhalfshu
self.PRF = PRF
self.vdengxiao = vdengxiao
self.mubiaojuli = mubiaojuli
...
wl = LC / self.f0
Ka = 2.0 * self.vdengxiao**2 / (wl * self.mubiaojuli)
mctime = (np.arange(self.hangmc) - self.hangmc * 0.5) / self.PRF
xiangwei = np.pi * Ka * mctime**2
self.phase_comp = np.exp(1j * xiangwei).astype(np.complex64)
self.phase_comp = self.phase_comp.reshape(self.hangmc, 1)
self.phase_comp = ms.Tensor(self.phase_comp)
...
def construct(self, input_data):
out = input_data * self.phase_comp
mohushu = 2 * self.mohuhalfshu + 1
results = []
shang = []
for mh in range(mohushu):
offset = mh - self.mohuhalfshu
tmp = self.Keystone_interp(out, self.f0, self.fs, self.B, offset)
hrrp, result = self.Gethrrp(tmp, self.hangmc, self.liejl)
shang_value = self.Gethrrp_v2shang(hrrp)
shang.append(shang_value)
result = ops.expand_dims(result, 0)
results.append(result)
results = ops.cat(results, axis=0)
shang = ops.cat(shang)
xuhao = ops.argmin(shang)
min_result = results[xuhao]
out = ops.cat([min_result[0:1, :], out[1:, :]], axis=0)
return out
其中Keystone_interp是核心代码,变换原理 Keystone变换是一种距离徙动校正算法,用于解决以下问题:
变换公式: 原始坐标:t (方位时间), f (距离频率) 变换后:t’ = t × (f₀/(f₀+f))
代码实现如下:
# 1. 距离向FFT(转到距离频率域)
# 2. 相位补偿(补偿距离徙动),方位向插值(Keystone变换核心)
# 3. 距离向IFFT(转回时域)
def Keystone_interp(self, input_data, f0, fs, B, moshuzhouqishu):
# 获取维度
fwn = self.fwn
fwn_8 = self.fwn_8
# 计算距离FFT长度
rnfft = self.rnfft
# 模块1: 距离向FFT
out = self.FFT4f(input_data, 1)
nyquist_idx = rnfft//2
new_nyquist_values = (out[:, nyquist_idx-1] + out[:, nyquist_idx+1]) / 2.0
# 使用concat操作替换特定列
left_part = out[:, :nyquist_idx]
nyquist_part = new_nyquist_values.expand_dims(1) # 增加维度
right_part = out[:, nyquist_idx+1:]
# 重新拼接
out = ops.concat([left_part, nyquist_part, right_part], axis=1)
# 计算无小点数量
wuxiaodiansujl = int(np.floor(0.5 * (fs - B) * rnfft / fs))
wuxiaodiansujl = (wuxiaodiansujl >> 3) << 3 # 8字节对齐
wuxiaodiansujl = max(0, wuxiaodiansujl)
# 模块2: 方位向插值
pindian = (ops.arange(0, rnfft, dtype=ms.float32) - 0.5 * rnfft) * fs / rnfft
j_indices = ops.arange(0, fwn, dtype=ms.float32).reshape(-1, 1)
angle = -2 * PI * j_indices * pindian / f0 * moshuzhouqishu
complex_angle = angle.astype(ms.complex64)
j_times_angle = complex_angle * self.j_constant
phase_comp_all = self.exp(j_times_angle)
temp_complex_top = out * phase_comp_all
temp_complex = ops.Pad(((0, fwn_8 - fwn), (0, 0)))(temp_complex_top)
temp_complex_valid = temp_complex[:, wuxiaodiansujl:rnfft - wuxiaodiansujl]
xnew_valid = self.xnew[:, wuxiaodiansujl:rnfft - wuxiaodiansujl]
xold_valid = self.xold[:, wuxiaodiansujl:rnfft - wuxiaodiansujl]
tmp = self.stoltsun(temp_complex_valid, xnew_valid, xold_valid)
update_part = tmp
left_part = out[:fwn, :wuxiaodiansujl]
right_part = out[:fwn, rnfft - wuxiaodiansujl:]
# 在列方向拼接三部分
out = ops.concat([left_part, update_part, right_part], axis=1)
# 模块3: 距离向IFFT
out = self.cirshift_optimized(out)
out = self.IFFT4f(out, 0)
return out
其中Keystone变换核心是方位向插值,这里的插值函数是使用SciPy的CubicSpline进行插值。
3.4 数据后处理
这里将结果导出成二进制数据的dat文件,用于和正确结果比较。确认结果无误后,可通过 mindspore.export 导出 MINDIR 模型,便于在 MindSpore Lite 端部署(板卡侧运行)。
完整代码参考链接: HRRP(高分辨率距离像算法)
四、实验结果验证
通过对比C++和MindSpore输出结果,确保:
- HRRP幅度特征一致
- 相位信息准确
- 熵最小化选择相同
五、总结与展望
本次从C++到MindSpore的移植不仅实现了功能迁移,更通过深度学习框架的特性:
- 为后续的AI增强提供了基础框架
- 建立了传统信号处理与深度学习融合的范例
注 :HRRP作为雷达信号处理的核心技术,其高效实现对于实时目标识别至关重要。通过将传统C++算法迁移到MindSpore框架,我们构建了一个可扩展、易优化的现代化信号处理流水线,为后续的AI赋能奠定了坚实基础。
