合成孔径声纳(synthetic aperture sonar，SAS)是一种高分辨水下成像声纳，可以获得湖/海底全距离恒定的二维高分辨图像[1-2]，在水下地形测绘、考古搜救、沉物打捞、猎雷及重要海域的高分辨侦察等领域具有广阔的应用前景．运动补偿(motion compensation，MOCO)是SAS成像过程中极为关键的步骤[3-4]，其本质是从运动传感器或回波数据中估计声纳运动误差，通过MOCO提高二维声图质量．声纳运动误差[5-7]主要表现为沿成像坐标系3个轴的线运动，即纵荡、侧摆和升沉，以及沿3个轴的角运动，即横滚、俯仰和偏航，合称为六自由度(six-degrees-of-freedom，SDOF)运动误差．目前常用的MOCO方法主要分为两类．一类被称为基于回波数据的MOCO方法，又可分为基于等效相位中心阵元的MOCO方法[8-9]和基于偏移子图像自聚焦的MOCO方法[10]，该类方法从声纳回波数据中估计声纳基阵的侧摆、偏航等误差，通过运动补偿提高SAS图像质量．一般地，该类方法无法获得全部SDOF运动误差，因而无法彻底补偿．另一类被称为基于运动传感器的MOCO方法[11-13]，该类方法根据SAS平台安装的惯性导航系统(inertial navigation system，INS)提供的姿态、位置等数据，以及INS与声纳基阵之间的相对位置关系，解算声纳基阵的SDOF运动误差，通过运动补偿提高SAS图像质量，其优势是可以获得声纳基阵的实际位置和全部SDOF误差．基于此，容易利用逐点MOCO(point-by-point MOCO，PPMOCO)算法实现SDOF补偿[12]．PPMOCO算法本质上属于二维时域算法，因而运算效率极低．文献[13]对PPMOCO算法进行多核并行化设计，有效提高了运算效率，但是并没有改变PPMOCO算法的本质．为解决补偿方位空不变SDOF运动误差补偿及PPMOCO算法效率低的问题，提出一种多子阵SAS运动补偿与成像快速算法，首先建立方位空不变SDOF运动误差存在情形下的多子阵SAS双根号形式距离历程，利用泰勒级数展开对双根号形式距离历程进行近似并保留至方位时间的四阶项，利用级数反演方法推导出严格解析的点目标响应二维谱，通过对传统距离多普勒(range-Doppler，RD)算法[14]进行改进，实现方位空不变SDOF运动误差的快速补偿与高分辨成像．仿真试验和实测数据成像结果验证了所提算法的有效性．1 多子阵SAS距离历程理想情况下，多子阵SAS沿理想航迹以速度v作匀速直线运动．由于SDOF运动误差的存在，声纳基阵偏离理想航迹，沿实际航迹运动，如图1所示．在方位时间0时刻，以发射阵中心为原点o，以理想航迹为x轴(向前为正)建立用于成像的xyz三维空间直角坐标系，即成像坐标系．其中，z轴垂直于湖/海底平面，向上为正；y轴垂直于xoz平面，向左为正．点o'为点o在湖/海底平面的投影，Rmin，Rmax和W分别表示最近测绘距离、最远测绘距离和有效测绘带宽．点P为成像场景中任一点目标，坐标用(0，yn，zn)表示．在信号发射t时刻，发射阵坐标用(vt，0，0)表示，在t+τi时刻，第i个接收子阵接收到目标反射回波，此时子阵i的坐标用(vt+Δdi+vτi，0，0)表示，其中：Δdi为发射阵与子阵i中心之间的距离；τi为子阵i接收信号的双程传播时间，i=1，2，…，M，M表示用于成像的子阵个数．10.13245/j.hust.240309.F001图1实际情形多子阵SAS距离历程几何模型由于SDOF运动误差的影响，收发阵元在信号收发时刻的位置会偏离理想坐标．利用空间直角坐标转换原理，可以根据t时刻声纳基阵的运动误差及发射阵的理想坐标，解算出发射阵的实际坐标，即为(vt+lx，ly，lz)，式中lx，ly和lz分别表示声纳基阵纵荡、侧摆和升沉等线运动误差的空不变部分．注意，lx，ly和lz可以根据INS输出的位置、姿态等数据，以及INS与声纳基阵之间的相对位置，利用空间直角坐标转换原理推算得到．此时，发射阵与点目标P之间的距离变化为RTr(t;r)=[(vt+lx)2+(ly-yn)2+(lz-zn)2]1/2，(1)式中r=yn2+zn2为目标与声纳的最近距离．对于同一发射信号和同一目标，不同子阵接收目标反射回波的时刻不同．当方位空不变SDOF运动误差存在时，可以根据运动误差与信号接收时刻子阵i的理想坐标，解算出信号接收时刻子阵i的实际坐标，即(vt+vτi*+Δdik1+lx,Δdik2+ly，Δdik3+lz)，其中：k1=cosθycosθz；k2= -cosθy∙sinθz；k3=-sinθy；θy和θz分别为声纳基阵俯仰和偏航的空不变部分，同样可以根据INS输出的位置、姿态等数据及INS与声纳基阵之间的相对位置关系推算得到．τi*为待求变量，表示方位空不变SDOF运动误差存在时子阵i接收信号的双程传播时间．因此，接收信号时刻子阵i与点目标P之间的距离为RRr,i(t;r)=[(vt+vτi*+Δdik1+lx)2+(Δdik2+ly-yn)2+(Δdik3+lz-zn)2]1/2. (2)结合式(1)和(2)可知，当方位空不变SDOF运动误差存在时，子阵i接收信号的双程距离历程为Rs,i(t;r)=RTr(t;r)+RRr,i(t;r)．(3)根据τi*的定义可知τi*=c-1Rs,i(t;r)，其中c为水下声速．将式(1)~(3)代入，参考理想情况下τi精确表达式的求解方法[15]，可以解得τi*=[q2+q22+q1q3]/q1，(4)式中：q1=c2-v2；q2=cRTr(t;r)+v(vt+Δdik1+lx)；q3=(Δdi)2+2Δdi[k1(vt+lx)+k2(ly-yn)+k3(lz-zn)]．将式(4)代入式(2)和(3)即可得Rs,i(t;r)的解析表达式．2 距离历程近似及误差分析如式(1)~(3)所示，由于SDOF运动误差的存在，Rs,i(t;r)的表达式为典型的双根号形式，这导致很难利用传统的驻定相位原理[16]推导严格解析的点目标响应二维谱，进而影响运动补偿及成像快速算法的设计．为解决该问题，通常对式(3)进行近似．由式(4)可知：τi*是决定式(3)复杂程度的重要因素．如果将其直接代入式(2)和(3)，会大幅增加式(3)的复杂度，因此一般先对τi*进行简化．考虑到τi*的距离空变和方位空变特性，这里近似为τi*≈2r/c+v2t2/(rc)．将其代入式(3)，并将式(3)在t=0处进行泰勒级数展开且保留至t的N阶项，即RA,i(t;r)≈δ0,i+δ1,it+δ2,it2+⋯+δN,itN，(5)式中δn,i表示泰勒级数展开过程中关于t的各阶项系数，具体为δn,i=1n!dnRs,i(t;r)dtnt=0，(6)其中n=0,1,⋯,N．注意，式(5)将作为后续推导点目标响应二维谱解析表达式的基础．显然，对于式(3)的近似会引入一定的相位误差，误差较大时可能影响成像质量．为避免该问题，通常须将相位误差控制在π/8范围之内．为了验证本文距离历程近似的有效性，同时确定N的取值，这里对不同距离历程近似引入的相位误差进行仿真分析．仿真过程中，距离历程真值由式(3)计算，其中τi*的取值由式(4)计算，以确保距离历程真值的准确性．通常收发阵元间隔Δdi越大，距离历程近似引入的相位误差越大．为了确保误差分析的阵元适应性，误差分析中Δdi用其最大值代替．仿真参数如下：水下声速为1 500 m/s，收发阵元最大间隔为0.84 m，距离向起止范围为40~540 m，方位向起止范围为-50~50 m，纵荡、侧摆和升沉均为0.3 m，俯仰和偏航均为2°．图2绘制了方位空不变SDOF运动误差存在时不同距离历程近似方法引入的相位误差，图中：r为距离向；x为方位向；Δϕ为相位误差(单位为π)．其中：图2(a)为理想距离历程模型引入的相位误差；图2(b)和(c)分别为本文距离历程近似方法在N分别取3和4时引入的相位误差；图2(d)对比了N分别取4，5和6时本文距离历程近似不同距离上的最大相位误差．10.13245/j.hust.240309.F002图2不同距离历程近似方法引入的相位误差如图2(a)所示，当方位空不变SDOF运动误差存在时，理想距离历程模型引入的最大相位误差达到58.57π，远大于π/8，难以满足运动补偿与成像算法设计要求．本文距离历程近似方法考虑方位空不变SDOF运动误差，相位误差大幅降低．然而，当N=3时，最大相位误差仍达2.1π，无法满足要求．但是，当N取4，5和6时，相位误差最大值均不超过0.03π，如图2(c)和(d)所示，均可以满足运动补偿与成像算法设计要求．由于N取值越大，式(5)的表达式越复杂，以其为基础的运动补偿与成像算法越复杂，运算量也会相应增大．为了降低公式复杂度，减小运算量，取N=4．3 算法实现严格解析的点目标响应二维谱是建立多子阵SAS运动补偿与成像快速算法的关键．3.1　点目标响应二维谱以子阵i为例，将其回波信号解调至基带，其点目标响应的二维时域形式可以表示为si(τ,t;r)=Aω(t)p[τ-RA,i(t;r)/c]×exp[-j2πf0RA,i(t;r)/c],式中：τ为距离向快时间；A为信号幅度；ω(⋅)为基阵指向性参数；p(τ)=rect(τ/Tr)exp(jπμτ2)为线性调频信号的复包络，其中，rect(⋅)为门函数，Tr为脉冲宽度，μ为调频斜率；f0为信号载频．为简化分析，后续推导忽略与成像质量无关的幅度项．依据级数反演理论以及傅里叶变换(Fourier transform，FT)的性质，可以获得子阵i的点目标响应二维谱[14]，其相位形式[15]表示为ϕi(fr,fa;r)=ϕrc(fr)+ϕac,i(fa;r)+ϕrcmc,i(fr,fa;r)+ϕsrc,i(fr,fa;r)+ϕres,i(r),式中：ϕrc(fr)=-πfr2/μ为距离调制项，是距离向脉压(range compression，RC)的来源项；ϕac,i(fa;r)为方位调制项，仅与方位频率fa有关，是方位向脉压(azimuth compression，AC)的来源项，具体为ϕac,i(fa;r)=b1(r)f0[h12(fa;r)-δ1i,c2]+b2(r)f0[h13(fa;r)-δ1i,c3]+b3(r)f0[h14(fa;r)-δ1i,c4], (7)其中，b1(r)=πc/[2δ2,i(r)]，b2(r)=πc2δ3,i(r)/ [4δ2,i3(r)]，b3(r)=πc3[9δ3,i2(r)-4δ2,i(r)δ4,i(r)]/ [32δ2,i5(r)]，分别表示由δ2,i(r)~δ4,i(r)决定的系数项，h1(fa;r)=fa0+δ1i,c，fa0=fa/f0，δ1i,c=δ1,i(r)/ c；ϕrcmc,i(fr,fa;r)为fr的1阶项与fa的耦合项，是距离徙动校正(range cell migration correction，RCMC)的来源项，具体为ϕrcmc,i(fr,fa;r)=fr[-2πδ0i,c+b1(r)∙h1(fa;r)(δ1i,c-fa0)+b2(r)h12(fa;r)∙(δ1i,c-2fa0)+b3(r)h13(fa;r)(δ1i,c-3fa0)], (8)其中δ0i,c=δ0,i(r)/c，ϕsrc,i(fr,fa;r)为fr的2~4阶项与fa的耦合项，是二次距离压缩(secondary range compression，SRC)的来源项，具体为ϕsrc,i(fr,fa;r)=fr2fa02/f0[b1(r)+3b2(r)h1(fa;r)+6b3(r)h12(fa;r)]-fr3fa02/f0[b1(r)+3δ1i,c)+b2(r)(4fa0+2b3(r)h1(fa;r)(5fa0+3δ1i,c)]+fr4fa02/f0[b1(r)+b2(r)(5fa0+3δ1i,c)+b3(r)(15fa02+20δ1i,cfa0+6δ1i,c2)], (9)ϕres,i(r)=f0[-2πδ0i,c+b1(r)δ1i,c2+b2(r)δ1i,c3+b3(r)δ1i,c4], (10)与fr和fa均无关，是剩余相位校正(redundant phase correction，RPC)的来源项．观察式(7)~(10)可知：AC，RCMC和SRC的来源项均与式(6)中δ0,i~δ4,i密切相关，这是本文距离历程近似中取N=4的结果．若取N为5或6，显然会增加点目标响应二维谱的复杂度，进而增加后续运动补偿与成像算法的复杂度，导致运算量增加．3.2　运动补偿与成像快速算法由式(7)~(10)可知，AC，RCMC，SRC和RPC的来源项均与子阵序号i有关，因此这些操作须要对不同阵元的回波数据单独处理，这里将其统称为阵元级回波数据处理．注意，虽然上节已经推导了子阵i点目标响应二维谱，但是由于多子阵SAS成像原理决定了该二维谱在方位向是欠采样的，因此通常须要将子阵i的二维谱沿方位向扩展，以满足方位向采样要求．图3给出了多子阵SAS方位空不变SDOF运动误差补偿与成像快速算法的结构图．首先对多子阵SAS回波数据进行RC处理，然后进行阵元级回波数据处理，最后沿方位向叠加，即可实现方位空不变SDOF运动误差补偿与高分辨成像．10.13245/j.hust.240309.F003图3本文运动补偿与成像快速算法结构图RC处理可以在距离频域通过与参考函数HRC,i相乘实现，其中HRC,i(fr)=exp[-jϕrc(fr)]．阵元级回波数据处理是算法的核心．在阵元级回波数据处理之前，须从多子阵SAS回波数据中分离出各阵元的回波数据．对于子阵i的回波数据，首先利用二维FT变换到二维频域，然后通过方位谱扩展使其满足方位向采样要求．由于二维频域不能同时表征距离和距离频率，SRC处理可以通过与一个参考距离rref处的参考函数HSRC,i相乘实现，其中HSRC,i(  fr, fa;rref )=exp[ -jϕsrc,i  (  fr, fa;rref ) ]．当测绘带较宽时，可以利用距离向分块方法保证不同距离处的SRC处理精度．RCMC处理通常在RD域利用sinc插值实现，因此须利用距离向逆FT(inverse FT，IFT)将SRC之后的二维频域数据变换回RD域，须校正的距离徙动量为Δri=-cϕrcmc,i(fr,fa;r)/(4πfr)．由于方位空不变SDOF运动误差的影响，Δri的最大值可能并不出现在最远距离处，在实际处理中须重点考虑．AC和RPC处理一般在RD域通过与参考函数HAC,i相乘实现，即HAC,i(fa;r)=exp[-jϕac,i(fa;r)]∙exp[-jϕres,i(r)].在阵元级回波数据处理之后，将所有阵元的回波数据沿方位向叠加，即可完成方位空不变SDOF运动误差补偿和高分辨成像．4 仿真试验及实测数据成像为验证本文算法的有效性，进行仿真及实测数据成像试验．仿真试验分别利用文献[14]所提基于四阶多项式模型(fourth order polynomial model，FOPM)的RD算法、文献[13]所提PPMOCO算法及本文算法对3种情形下的点目标回波数据进行成像．图像质量指标[16]一般用分辨率(I)、峰值旁瓣比(P)和积分旁瓣比(L)衡量．仿真参数如下：信号载频为80 kHz，带宽为20 kHz，脉冲宽度为20 ms，平台航速为2.5 m/s，声纳基阵纵荡、侧摆和升沉误差仍为0.3 m，俯仰和偏航误差为2°，点目标坐标为(0 m，806 m，-40 m)，即方位向为0 m、距离向为200 m．实测数据成像利用3种算法对ChInSAS系统某航次海试中人造梯台形目标回波数据进行成像，系统参数与仿真试验一致，方位空不变SDOF运动误差根据INS实测数据推算得到．4.1　仿真试验情形1 仅存在纵荡、侧摆、升沉等平动误差．图4(a)~(c)分别绘制了FOPM RD算法、PPMOCO算法和本文算法的成像结果，图4(d)和(e)分别对比了3种算法的距离剖面图和方位剖面图，图像质量指标如表1情形1所示．从图像质量上看，虽然没有考虑平动误差，但是FOPM RD算法在距离向和方位向上均聚焦良好；相比之下，本文算法和PPMOCO算法同样聚焦良好．从定位精度上看，FOPM RD算法在距离向和方位向均存在一定偏移，分别约为-0.23和-0.3 m；本文算法和PPMOCO算法则定位准确，如图4(d)和(e)所示．10.13245/j.hust.240309.F004图4仅存在平动误差时3种算法成像结果对比10.13245/j.hust.240309.T001表1图像质量指标对比类别距离剖面图方位剖面图情形算法I/mP/dBL/dBI/mP/dBL/dB1FOPM RD0.038-13.3-10.10.079-14.2-11.6PPMOCO0.037-13.2-10.40.08-14.3-11.4本文0.037-13.3-10.40.08-14.1-11.52FOPM RD0.085-10.6-8.90.08-14.216.7PPMOCO0.038-13.8-10.40.08-14.2-11.5本文0.037-13.3-10.40.08-14.3-11.53FOPM RD0.083-11.4-9.20.08-14.115.6PPMOCO0.037-13.0-10.30.08-14.2-11.5本文0.037-13.3-10.40.08-14.1-11.5情形2 仅存在俯仰、偏航等角度误差．图5(a)~(c)分别绘制了3种算法的成像结果，图5(d)和(e)分别对比了3种算法关于点目标本身的距离剖面图和方位剖面图，图像质量指标如表1情形2所示．从图像质量上看：由于未考虑角度误差，因此FOPM RD算法在方位向出现大量的虚假目标，部分虚假目标的强度甚至已经超过设定目标，如图5(a)所示．经过测量，与成像场景中最强点相比，设定点目标的强度约为-5 dB；同时，由于大量虚假目标的出现，其方位剖面图的积分旁瓣比高达16.7 dB，远超理论值-11 dB．实际工作中，高强度虚假目标的存在，可能导致真实目标被掩盖．如图5(d)所示，FOPM RD算法在距离向存在明显的双峰效应，由于仿真数据中仅包含单个点目标，因此这里认为双峰为同一目标．经过量测，双峰对应的距离分辨率约为0.085 m，远差于理论值0.037 m．注意，经过量测，图5(d)中两个单峰对应的分辨率分别约为0.029和0.032 m，均优于理论值，因此认为其中较高峰值对应于点目标主瓣，较低峰值对应于第一旁瓣是不合理的．10.13245/j.hust.240309.F005图5仅存在角度误差时3种算法成像结果对比本文算法和PPMOCO算法均聚焦良好，如图5(b)~(e)所示．从定位精度上看，在不考虑虚假目标的情况下，FOPM RD算法在方位向定位准确，如图5(e)所示；在距离向，若将双峰的中间位置定义为目标的距离向坐标，则存在不大于一个距离向分辨单元的定位误差，如图5(d)所示．在实际工作中，真实目标可能已经被虚假目标所掩盖而形成较大的定位误差．相比之下，本文算法和PPMOCO算法均有效补偿了角度误差，定位准确，如图5(d)~(e)所示．情形3 同时存在平动误差和角度误差．图6(a)~(c)分别绘制了3种算法的成像结果，图6(d)~(e)分别对比了3种算法关于点目标本身的距离剖面图和方位剖面图，图像质量指标如表1情形3所示．从图像质量上看，FOPM RD算法在方位向出现大量的虚假目标，部分虚假目标强度超过设定点目标强度，与成像场景中最强点相比，设定目标强度约为-5.2 dB，积分旁瓣比达到15.6 dB，如图6(a)和(e)所示；在距离向上，同样出现明显的双峰效应，距离分辨率约为0.083 m，远差于理论值，如图6(d)所示．这些特征与仅存在角度误差时的特征类似，故方位向虚假目标和距离向双峰效应的出现可能主要是由于角度误差未被补偿．从定位精度上看，FOPM RD算法在距离向和方位向上均存在一定偏移，分别约为-0.22和-0.3 m，如图6(d)和(e)所示，该偏移量与仅存在平动误差时的偏移量基本一致，因此二维定位误差可能主要是由于平动误差未被补偿．相比之下，本文算法和PPMOCO算法不仅聚焦良好，而且定位准确，如图6(b)~(e)所示．上述分析充分说明了本文算法的有效性．10.13245/j.hust.240309.F006图6同时存在平动和角度误差时3种算法成像结果对比4.2　实测数据成像ChInSAS系统某航次海试数据中目标类型为高2 m、顶宽0.4 m、底宽4.4 m、斜面倾角45°的人造梯台形体目标．根据INS输出的位置和姿态数据，推算出方位空不变纵荡、侧摆和升沉误差分别为-0.03，0.02和0.03 m，俯仰和偏航误差分别为-0.9°和2.24°．在实际工作中，目标类型大多为面目标或体目标，一般用目标图像的边缘锐化（清晰）评价实测数据图像质量的优劣．图7(a)~(c)分别给出了FOPM RD算法、PPMOCO算法和本文算法对梯台形目标回波数据的成像结果，目标边缘如图中A，B和C箭头所指区域．与PPMOCO算法相比，FOPM RD算法在3个边缘区域明显较为模糊，如图7(a)和(b)所示．这是由于该组数据存在较小的方位空不变SDOF运动误差，而FOPM RD算法并未补偿，从而导致其所成图像在方位向出现虚假目标，在距离向分辨率下降，共同导致边缘模糊现象发生．与FOPM RD算法相比，本文算法补偿了方位空不变SDOF运动误差，所成图像目标边缘处更为清晰，如图7(c)所示，其程度与PPMOCO算法相当，这表明了本文算法在实测数据处理中的有效性．10.13245/j.hust.240309.F007图7实测数据运动补偿成像结果对比5 结语仿真试验和实测数据成像结果不仅表明了所提算法的有效性，而且揭示了方位空不变SDOF运动误差对成像质量的影响．由于本文算法仅考虑方位空不变SDOF运动误差，对运动传感器的测量精度要求较低，因此可以广泛应用于运动传感器精度不高时的快速运动补偿和高分辨成像．