水下回波处理中分数阶傅里叶变换的带通采样实现方法*

詹 飞 马晓川 吴永清

(1 中国科学院声学研究所 中国科学院水下航行器信息技术重点实验室 北京 100190)

(2 中国科学院大学 北京 100049)

0 引言

主动声呐系统常使用线性调频(Linear frequency modulated，LFM)信号[1-2]对水下目标进行探测与参数估计。实际应用中，根据LFM信号的宽多普勒容限特性，通常选取零多普勒速度的单个副本信号和接收回波进行匹配滤波。为提升主动声呐系统的距离和速度分辨力，通常会增大发射信号的时间带宽积。LFM 信号时间带宽积的增加导致信号多普勒容限变化[3]，多普勒效应带来的影响不可忽略，单副本匹配滤波处理方法会降低系统检测性能并带来距离估计偏差，而且也无法获得目标速度信息。分数阶傅里叶变换(Fractional Fourier transform，FrFT)作为一种线性时频分析工具，在处理LFM 信号时具有独特优势[4]。目前已有大量研究将FrFT 方法及改进方法应用于处理LFM 回波[5-10]，从而获得目标参数的估计。研究结果表明，在低信噪比[11]和强混响背景[12-13]下利用FrFT处理LFM回波能提高系统检测性能，实现对目标回波时延和目标速度的有效估计。

主动声呐发射信号的脉宽增加时，若按窗处理接收回波，则为了完整捕获目标回波，处理窗宽度需相应增大。此外，脉宽增加带来的速度分辨力提高使得系统可以在设定速度范围内实现更精细的搜索。系统性能的提升却带来了回波处理算法计算复杂度的显著增加，尤其是奈奎斯特采样下FrFT的回波处理算法。理论计算表明FrFT 方法的计算复杂度与FrFT 变换阶数数量(等效为多普勒处理通道数量)呈线性变化趋势，与处理数据的长度呈对数变化趋势。在9 个波束方向上以0.1 m/s 间隔搜索-4～4 m/s速度范围内的水下目标时，奈奎斯特采样下处理一窗32768 点数据时FrFT 方法的计算复杂度高达100 GFloats 以上。高的计算复杂度使FrFT 方法难以应用于具有实时性要求的主动声呐系统，尤其是功耗和体积受限的无人水下航行器(Unmanned underwater vehicle，UUV)平台。

本文基于主动声呐回波后处理框架，结合LFM信号提出FrFT 的带通采样实现方法，通过对LFM信号时频特性直线在分数阶域的投影进行修正，使利用FrFT 方法处理带通采样回波数据时，可获得正确的目标参数估计。在FrFT 方法计算复杂度理论分析基础上，计算机仿真数据和UUV 平台湖试数据的处理结果验证了FrFT 的带通采样实现方法的正确性，数据处理时间能够满足UUV 平台处理的实时性要求，从而实现了FrFT 方法的工程化实时应用。

1 系统描述

1.1 分数阶傅里叶变换

FrFT 是传统傅里叶变换的广义形式，可以理解为信号在时频二维平面上进行逆时针旋转。信号x(t)的p阶连续FrFT线性积分形式定义[4]为

其中，积分核为

FrFT 的实现过程主要分为7 个步骤：(1)HT：对一帧时域采样数据执行Hilbert变换，得到时域数据的解析形式；(2)INTERP2：对解析数据进行两倍插值，可转换到频域实现；(3)C_MUL1：插值后的解析数据与LFM 信号1 执行时域乘法运算；(4)C_CONV：乘法运算后的数据与LFM 信号2 执行时域卷积运算，可转换到频域实现；(5)C_MUL2：卷积运算后的数据与LFM 信号1 以及复幅度因子执行时域乘法运算；(6)EXTRACT2：对乘法运算后的数据进行两倍抽取，得到与输入原始数据长度一致的FrFT输出结果；(7)NORMA：对FrFT输出结果执行归一化运算，获得归一化幅度信息。

1.2 系统处理流程

基于FrFT 的主动声呐LFM 回波处理流程如图1所示，主动声呐系统发射LFM 脉冲信号后，声呐基阵开始接收回波信号。接收信号经过调理和带通采样后，送入声呐系统信号处理单元。信号处理单元首先对采集数据执行多个波束方向的波束形成运算[14]，得波束加权时域输出数据。公式(1)显示FrFT运算主要是基于解析信号的，因此执行FrFT 运算前，先对波束加权输出实数据进行Hilbert变换转换为复数，再以变换阶数p为变量，对复数数据进行不同变换阶数的FrFT 运算，输出结果构成变换阶数—分数阶域二维平面((p，u)平面)的形式。对FrFT 输出结果进行归一化运算后，从(p，u)平面中搜索最大值，并与设定门限进行比较。若最大值超过判决门限，则根据最大值坐标计算出LFM 信号调频斜率和中心频率参数，进而得到目标距离与目标速度估计信息。

1.3 计算复杂度分析

算法各步骤计算复杂度公式如表1所示。表中M表示声呐基阵阵元数，N表示一窗数据的长度(设置为2 的整数次幂有利于运算)；NB表示执行波束加权的波束方向数量，NBW表示执行基于频域离散傅里叶变换(Discrete Fourier transform，DFT)波束形成的频域子带数量，NP表示FrFT变换阶数p的数量；Om= 6 表示复数乘法所需浮点运算次数，Oa=2表示复数加法所需浮点运算次数。对比算法各步骤计算复杂度公式可知，C_CONV 步骤的计算复杂度最大。

图1 基于FrFT 的主动声呐LFM 回波处理流程Fig.1 Processing diagram of LFM echo for active sonar based on FrFT

表1 算法各步骤计算复杂度公式Table1 Computation complexity formula of each step of the method

工程应用时，奈奎斯特采样条件下系统采样率通常设置为4～10倍的信号上限频率。对于具有较高中心频率的系统，若接收信号为带通信号，则常采用带通采样降低系统设计难度。对于下限频率为fL、上限频率为fH、带宽为B的带通信号，带通采样[15]下采样率需满足2fH/m≤fs≤2fL/(m-1)，其中，符号为向下取整。

利用FrFT 分别处理两种采样条件下计算机仿真的LFM 回波数据。仿真条件如下：主动声呐发射LFM 信号的中心频率26 kHz，带宽3 kHz，脉宽120 ms。主动声呐系统与目标相对速度为-2 m/s，相对距离750 m。考虑噪声限制情况，目标回波信噪比为0 dB。系统采集2 s 回波数据，利用FrFT按窗对回波数据进行处理。为使单个处理窗覆盖完整目标回波，处理窗宽设置为两倍信号脉宽，并且处理窗按照50%重叠。奈奎斯特采样时设置系统采样率fs= 136.533 kHz，则处理窗点数Nn= 32768。带通采样时设置系统采样率fs= 34.133 kHz，则处理窗点数Nb= 8192。奈奎斯特采样下处理窗点数为带通采样下处理窗点数的4 倍。通过计算，两种采样条件下第9 个处理窗覆盖完整目标回波，第8 个和第10 个处理窗覆盖部分目标回波。仿真数据的处理结果如图2所示，两种采样条件下，第8～10 窗FrFT输出归一化幅度明显大于其他窗，第9窗输出归一化幅度最大。估计出的目标信息如表2所示，带通采样时利用修正的回波时延估计公式计算目标回波到达时间。处理结果表明，两种采样条件下都能获得较准确的目标参数。

图2 仿真数据处理结果Fig.2 The processing results for simulated data

表2 处理仿真数据得到的目标信息Table2 Target parameters obtained by processing simulated data

在此仿真条件下，LFM 回波处理算法计算复杂度如图3所示。理论计算结果显示算法中C_CONV 步骤的计算复杂度占算法总计算复杂度的95%以上，因此算法总计算复杂度的变化趋势与C_CONV 步骤的计算复杂度变化趋势一致。图3(a)和图3(b)显示计算复杂度与波束数量NB和变换阶数NP呈线性关系。当NB= 3、NP= 81 时，带通采样下算法计算复杂度约为8.5 GFloats，奈奎斯特采样下算法计算复杂度约为37.8 GFloats，是带通采样的4.45 倍。图3(c)显示计算复杂度与数据长度N(处理窗点数)呈对数关系，随着N的增加，计算复杂度明显降低。图3表明带通采样可显著降低回波处理算法的计算复杂度。

图3 算法的计算复杂度变化趋势Fig.3 Change trend of computation complexity of the method

2 基于FrFT的目标参数估计

2.1 目标回波时延估计

在工程应用中，通常对数据按窗进行处理，为保证处理窗能够覆盖完整目标回波，要求处理窗宽(即观测时长t0)大于信号脉宽。假设脉宽为Δt、带宽为Δf的LFM 信号在数据处理窗内的时延偏移量为τ。文献[12—13]中给出了基于FrFT 的目标回波时延估计公式，该公式建立了分数阶域u与时延ut的联系(由于采样率发生变化，将τ改写为ut)。当u轴旋转到最佳位置，最佳旋转角与变换阶数-分数阶域平面峰值坐标对应，则时延ut为

式(4)成立的基础是量纲归一化时用处理窗宽t0代替信号脉宽Δt，用系统采样率fs代替信号带宽Δf得到尺度因子S，并且系统采样率fs满足奈奎斯特采样定理。此时信号的能量限定在以原点为中心、以Δx/2 为半径的圆内[17]，其中当系统采样率fs不满足奈奎斯特采样时，即带通采样下无法利用该公式计算出正确的目标回波时延。为解决此问题，本文对该公式进行了修正。带通采样条件下，经过量纲归一化后信号的能量被限制在时频平面上以为中心、以Δx/2为半径的圆内。图4显示了m= 2 时FrFT 示意图，信号的瞬时频率函数与新坐标系频域轴的交点在旋转了-α角度的u轴上的投影为

图4 带通采样条件下FrFT 示意图Fig.4 Diagram of FrFT with bandpass sampling

同样将τ改写为ut，当u轴旋转到最佳角度时，时延ut与分数阶域u的关系为

则处理窗内LFM回波时延偏移量为

以上分析利用了时频平面上LFM 信号时频特性直线在分数阶域坐标轴的投影，得到目标回波时延与LFM信号参数的关系。根据对FrFT计算过程的分析得到计算目标回波时延的另一种思路：当处理窗宽t0大于信号脉宽Δt时，利用FrFT 对一窗数据进行处理得到的中心频率是将整窗数据看作LFM 信号时的中心频率(即t0/2 时刻的瞬时频率)，并非处理窗内LFM 信号真实中心频率fc。则根据图5所示的3种情况，得到处理窗内回波时延为

式(8)与公式(7)一致，根据记录的处理窗序号，以及处理窗内回波时延偏移量ut0，即可计算出目标回波到达时刻，进而计算出目标距离。

图5 回波时延估计示意图Fig.5 Diagram of the time delay estimation

2.2 目标速度估计

利用FrFT 估计LFM 信号参数的基本思路是以变换阶数p为变量，对接收回波进行0～2 阶的FrFT 运算，从得到的变换阶数-分数阶域平面上的能量二维分布中搜索峰值点，并计算旋转角度。为避免产生较大估计误差，仿真测试变换阶数p的间隔需要至少为10-5量级，这导致搜索算法的运算量过于庞大。事实上可直接利用待测目标范围计算出理论最佳旋转角α的范围，进而仅在此范围内完成对接收数据的FrFT 运算即可。向量形式的最佳旋转角α公式为

式(9)中，η ≈1+(2vd)/c为多普勒压缩因子，vd为主动声呐系统与目标的相对运动速度。经过量纲归一化处理，处理窗内LFM信号的调频斜率和中心频率分别为

再根据主动声呐系统的速度计算出目标速度。若在公式(9)计算出的理论最佳旋转角α的范围内对接收数据进行处理，则可直接根据索引值得到相对速度估计值，进而计算出目标速度。

3 数据分析

3.1 数据处理结果

以UUV 平台为对象，验证FrFT 的带通采样实现方法的正确性及实时性。测试数据选用2015年冬季UUV 平台的千岛湖实验数据。UUV 平台搭载主动声呐基阵，以固定速度航行，航行过程中间隔发射脉宽为120 ms 的LFM 脉冲信号。实验选用应答器作为目标，与UUV 平台处于同一深度，应答器可模拟具有不同散射强度和速度的点目标散射回波。相对于UUV 平台，应答器初始方位为θb= (-7°，0°)，距离为r0=1790 m，速度为vt=-2 m/s。

处理UUV 平台湖试数据时，设置50%重叠处理窗，并且处理窗宽为两倍信号脉宽，对应处理点数N= 8192。对单次航行任务连续接收的29 ping LFM回波数据(序号为0～28)进行处理，得到的目标方位、距离和速度信息如图6所示。第0～2 ping数据对应目标方位为(-7°，0°)，UUV通过前3 ping数据处理结果估计出目标参数后，会及时调整其航向，因此其他ping数据对应目标方位为(0°，0°)。目标距离估计值为900～1790 m 范围，由于UUV 发射LFM 脉冲的时间间隔并不相同，距离并非线性变化。由于第2、第9、第16和第17 ping接收数据中目标回波出现严重畸变，估计出的目标速度与设定值存在较大偏差，其他ping 数据对应的目标速度在-2 m/s附近变化。

图6 所有ping 数据目标信息估计Fig.6 Estimation results of target parameters for all ping data

分别选取第0 ping 和第23 ping 数据进行处理，计算得到第0 ping 目标回波信噪比约4.5 dB，第23 ping目标回波信噪比约2.6 dB。每窗数据处理结果如图7所示，将目标出现数据窗的FrFT结果转换为模糊图的形式，如图8所示。图7(a)显示第0 ping数据中第20、第21 和第22 窗数据处理结果超过判别门限，第21 窗数据完全覆盖目标回波，归一化结果最大。图7(b)显示第23 ping 数据中第12、第13和第14 窗数据处理结果超过判别门限，第13 窗数据完全覆盖目标回波，归一化结果最大，估计出的目标信息如表3所示。采用50%重叠窗处理方式，最大值前一窗和后一窗数据恰好都能够覆盖部分目标回波，通过FrFT 处理可获得近似的目标回波时延估计。实际应用时，可利用FrFT的该特性对目标距离进行多次确认，从而得到可靠的目标距离估计。

图7 每窗数据处理结果Fig.7 Processig results for data of each window

表3 处理湖试数据得到的目标信息Table3 Target parameters obtained by processing experimental data

图8 目标出现数据窗的模糊图Fig.8 Ambiguity graphs of the data when target exists

3.2 性能分析

在具有实时性处理要求的主动声呐系统中，通常选择按窗处理回波数据。为了保证处理窗能够捕捉到完整目标回波，要求处理窗宽大于信号发射脉宽，并且处理窗间具有一定的重叠性。处理一窗数据的时间小于数据的更新时间时能满足实时性要求。根据UUV 平台发射LFM 信号参数与系统采样率fs，处理一窗8192 点数据时，处理时间小于120 ms时可满足实时性。根据1.3 节分析，制约FrFT 方法实时应用的瓶颈在于C_CONV 步骤高的计算复杂度，而C_CONV步骤主要是FFT/IFFT运算。计算得NB= 3、NP= 81 时，带通采样下算法的计算复杂度约8.5 GFloats，其中C_CONV 步骤中FFT/IFFT 运算的计算复杂度约8.3 GFloats。工程实现时，C_CONV 步骤中用于时域卷积的LFM 信号可预先生成，则该步骤计算复杂度降低为5.5 GFloats。UUV 平台信号处理单元为嵌入式图形处理器(Graphics processing unit，GPU)平台，执行16N点批处理FFT 运算的计算吞吐量约80 GFlops，则预估该步骤执行时间约69 ms。在此条件下，分别对第0 ping 共35 窗数据和第23 ping共16 窗进行处理，处理时间如图9所示。在嵌入式GPU 平台上，所有窗数据的处理时间约90 ms，始终小于120 ms，能够满足UUV 数据处理的实时性要求。利用带通采样下FrFT 的LFM 回波处理方法，UUV 平台能在3 个波束方向上实时搜索速度在-4～4 m/s 范围内的水下目标。

图9 每窗数据处理时间Fig.9 Processig times for data of each window

4 结论

本文基于主动声呐回波后处理框架，结合LFM信号提出FrFT 的带通采样实现方法，通过对LFM信号时频特性直线在分数阶域的投影进行修正，使利用FrFT 方法处理带通采样回波数据时，可获得正确的目标参数估计。对FrFT 方法计算复杂度的理论分析结果表明，处理相同脉宽回波数据时，若奈奎斯特采样下处理窗宽为带通采样下处理窗宽的4倍，则带通采样下算法计算复杂度可降低至奈奎斯特采样下算法计算复杂度的22%。计算机仿真数据和UUV 平台湖试数据的处理结果验证了FrFT 的带通采样实现方法的正确性，数据处理时间能够满足UUV 平台处理的实时性要求，从而实现了FrFT方法的工程化实时应用。

下一步将深入研究FrFT 的带通采样实现方法的性能，以及如何进一步降低FrFT 方法的计算复杂度，以满足主动声呐系统中更广泛的应用需求。