[DSP案例]HRRP应用与从C++到MindSpore的移植实践

一、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 核心算法

步骤分解:

  1. 相位补偿:首先对输入数据(复数形式)进行相位补偿,乘以预计算好的相位补偿因子(用于校正二次相位误差)。

  2. 多模糊数处理:计算需要搜索的模糊数总数(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中。

  3. 选择最优结果:将所有模糊数对应的结果堆叠起来(使用ops.cat),然后找到熵最小的那个索引(xuhao = ops.argmin(shang)),并取出该索引对应的结果(min_result = results[xuhao])。

  4. 更新输出:将最优结果的第一行(即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赋能奠定了坚实基础。