相关测速声呐模型参数矩阵化求解方法.pdf

第 42 卷第 4 期2023 年 8 月Vol.42,No.4Aug.,2023声学技术Technical Acoustics相关测速声呐模型参数矩阵化求解方法崔晓康1,2，薛凤杰1,3，赵二亮1,3，王长红1,3(1.中国科学院声学研究所海洋声学技术中心，北京 100190；2.中国科学院大学，北京 100190；3.北京市海洋声学装备工程技术研究中心，北京 100190)摘要：相关测速声呐的基本原理是“波形不变性”。该原理指出：发射具有一定时间间隔的信号，接收回波相关性最大的两个基元之间的间距与载体移动速度成正比。相关测速算法首先通过回波数据构建数据相关矩阵，然后利用速度估计算法将理论相关模型和回波数据构建的相关矩阵进行拟合得到最优参数信息，进而解算出载体速度信息。文章以简化最优参数求解流程为目标，提出了基于矩阵推导的参数拟合运算方法。该方法通过将底混响时空相关函数模型中的目标参量进行分离，进而将函数转变为矩阵形式，随后通过矩阵QR分解等手段，根据局部最小二乘法推导目标参数的求解公式，并直接求解出目标参数的解析解。仿真和湖上试验数据验证了该方法可以作为一种速度解算方法的可行性，并且具有运算复杂度低的优点。关键词：相关测速声呐；最小二乘；矩阵化中图分类号：O427.9 文献标志码：A 文章编号：1000-3630(2023)-04-0541-06Parametric matrix solution method of correlation velocimetry sonar modelCUI Xiaokang1,2,XUE Fengjie1,3,ZHAO Erliang1,3,WANG Changhong1,3(1.Ocean Acoustics Technology Center,Institute of Acoustics,Chinese Academy of Sciences,Beijing 100190,China;2.University of Chinese Academy of Sciences,Beijing 100190,China;3.Beijing Engineering Research Center for Ocean Acoustic Equipment,Beijing 100190,China)Abstract:The basic principle of correlation velocimetry sonar is waveform invariance,which is described as follows:when a carrier transmits signals with a certain time interval,the distance between the two elements with the greatest correlation between the received echoes is proportional to the carriers moving speed.The correlation velocimetry algorithm firstly constructs the data correlation matrix through the echo data,then uses the velocity estimation algorithm to fit the theoretical correlation model with the correlation matrix constructed by the echo data to get the optimal parameter information,and then calculates the carrier velocity information.In order to simplify the process of solving optimal parameters,a parameter fitting algorithm based on matrix derivation is proposed in this paper.In this method,the target parameters in the space-time correlation function model of bottom reverberation are separated,and then the function is expressed in matrix form.Finally,by matrix QR decomposition and local least squares,the analytical solution formula of the target parameters is obtained.The feasibility of this scheme as a velocity solution method is verified by simulation and on-lake test data,and it has the advantages of low computational complexity.Key words:correlation velocimetry sonar;least squares;matrix0引 言随着对海洋资源探索的不断深入，潜航设备对测速导航功能的需求和设备所能提供的测速精度都在不断提升。由于传播问题和存在累积误差，全球定位系统和惯性导航都不能满足水下长时间的使用。目前主要的水下测速导航设备是测速声呐。测速声呐包含多普勒测速声呐和相关测速声呐两种，按照不同的原理和功能进行划分，可对应4种不同类型的设备：分别为声相关计程仪(Acoustic Correlation Log,ACL)、声相关流速剖面仪(Acoustic Correlation Current Profiler,ACCP)、声多普勒计程仪(Acoustic Doppler Log,ADL)和声多普勒流速剖面仪(Acoustic Doppler Current Profiler,ADCP)。多普引用格式：崔晓康,薛凤杰,赵二亮,等.相关测速声呐模型参数矩阵化求解方法J.声学技术,2023,42(4):541-546.CUI Xiaokang,XUE Fengjie,ZHAO Erliang,et al.Parametric matrix solution method of correlation velocity measurement sonar modelJ.Technical Acoustics,2023,42(4):541-546.DOI:10.16300/ki.1000-3630.2023.04.019收稿日期:2022-03-01;修回日期:2022-05-29基金项目:中科院声学所青年英才基金(QNYC201807)。作者简介:崔晓康(1993),男,河南平顶山人,硕士研究生,研究方向为相关测速声呐。通信作者:王长红,E-mail:2023 年声学技术勒测速声呐主要依靠发射声波和散射体之间的多普勒效应来解算载体速度信息1；相关测速声呐依靠多个接收阵元接收信号的延时相关性来测量载体速度信息2。相对于多普勒测速声呐相关测速声呐具有以下两点优势：(1)相同尺寸的换能器，发射声波频率更低，作用距离更远；(2)声波垂直向下发射，受载体姿态影响较小。但是目前相关测速声呐相对于多普勒测速声呐测速精度较低。一直以来研究人员都在为提升测速精度不断努力。冯雷博士提出了简化的理论相关函数模型，使得参数解算更加简便，便于实时运算3；蒯多杰根据修正后的理论相关函数值对最小二乘的目标函数计算公式进行了自适应加权4。易卉琴等5基于真实数据的相关模型将理论相关函数添加了方向幅值调节因子和旋转因子，使得理论相关函数俯视图形态由圆变为椭圆，更符合实际数据的相关形态。王映春等考虑了阵元坐标误差给测速结果带来的影响，利用相关延时为0的数据和理论相关函数来估计阵元声中心，进而提高测速精度6。文献7-8中分析了载体纵横摇对设备测速的影响，并提出了一种目标函数添加权重调节项来提升速度估计精度的方法。本文基于局部最小二乘算法，利用矩阵推导直接求解出目标的解析解，具有无须迭代和运算复杂度低的优点。1基本工作原理波形不变性是声相关技术的基本工作原理。波形不变性原理在常见的流体散射和海底散射等各种声波散射过程中都适用9-10。本文以最简单的一发两收的声呐阵型来阐述相关测速声呐的测速原理。相关测速声呐采用的是宽波束发射换能器，声波以一定的波束角垂直向下发射，以一定间隔连续发射两次信号，对应的发射声波海底散射区域图如图1所示(图1中位置1 3个阵元从右往左分别是发射阵元T、接收阵元R1和R2)。当以时间间隔发射两次声波，经过相似的传播区域后以时间间隔被两个接收阵元接收，此时两个阵元接收的信号具有极大的相似性。根据几何知识可以得出式(1)所示的关系式：V=-d122(1)式(1)表明，载体速度V可以通过接收阵元R1和R2构成的矢量d12和两次发射信号的时间间隔来确定。波形不变性原理为载体测速提供了一种求解思路。该原理主要涉及载体相对于散射体的速度、接收阵元构成的空间相关矢量和发射信号相关延时3个参数。通过式(1)将求速度问题转换成了求矢量d12和相关延时的问题，矢量和相关延时两个参数如何获得取决于相关测速声呐的设计。有3种比较通用的方法3：(1)时间相关方法。固定矢量d12，通过相关函数搜索相关延时，然后得到速度信息。(2)空间相关方法。固定相关延时时，在较多的空间矢量中搜索出适合的矢量d12，进而得到速度信息。(3)时空相关方法。通过多种相关延时与阵型设计的多个矢量d12相结合，进而得到速度信息。通常采用的速度解算方法是空间相关方法，即利用有限的阵元进行特定的布局，利用回波信号构建空间相关矩阵，借助理论相关函数来拟合数据相关矩阵的数据来获取目标函数，最终求取极值来获取水平速度信息。ACL数据处理流程如图2所示。理论相关函数的准确性对速度估计结果的精度起着决定性作用。在确定了理论相关函数后，模型图1 声波在海底散射区域的示意图Fig.1 Schematic diagram of the sound scattering area on the sea floor图2 ACL数据处理流程图Fig.2 Flow chart of ACL data processing542第 4 期崔晓康等：相关测速声呐模型参数矩阵化求解方法拟合的精确程度决定了速度估计结果的精度。理论相关函数也被称作混响时空相关模型，最早由国外的 Dickey 推导给出2，国内最早由 Zhu等11根据 Faure，Olshevskii 和 Middleton 提出的FOM模型推导建立，随后经过冯雷博士、蒯多杰博士以及易卉琴等研究人员的不断改进后3-5，目前采用的理论相关函数为RT(dV)=()exp-j0(+2Vz+zc)Rmin+(1-Rmin)exp-14()ke2()2Vx+x cosc+()2Vy+y sinc2+ac()2Vy+y cosc-()2Vx+x sinc2(2)其中：是相关时延，V是三维速度，()是发射信号相关函数，k=0/c是波数。在确定了底回波的理论相关函数后，速度估计方法就是决定性能和精度的关键。目前采用的速度解算方法是以局域最小二乘法为核心的参数分布估计法。该方法采用了水平速度和垂向速度分步估计的思路，利用局域最小二乘法估计水平速度，利用解析方法求解垂向速度，是目前最有效的速度解算方法。针对以局域最小二乘法为核心的参数估计方法，当局域范围确定后，以最小二乘确定目标函数，搜索使目标函数最小的参数变量。最优参数求解流程如图3所示。搜索方向和步长采用拟牛顿法，用近似黑塞(Hessian)矩阵来避免Hessian矩阵求取时的复杂运算。近似 Hessian 矩阵采用的是 1970 年前后由Broyden、Fletcher、Goldfarb 以 及 Shanno 提 出BFGS算法，具有不用求逆和不用计算二阶导数等优点，在此用Bk代替，Bk的迭代公式为Bk+1=Bk+ykyTkyTksk-()Bksk()BkskTsTkBksk(3)式中：sk=k-k-1，yk=gk-gk-1。sk代表目标参数k的差值，yk代表一阶梯度gk的差值。根据简化后的近似Hessian矩阵即根据式(3)，不断迭代来获得一定相对偏差Dgk下的目标参数，进而解算出速度信息。k+1-k=B-1()k+1Dgk(4)此算法在参数的求解过程中需要多次迭代，参数多，过程复杂，而且需要考虑众多的边界条件限制。为了简化求解流程，降低算法复杂度，本文经过系统的理论推导，提出了基于矩阵推导的模型参数求解方法。2基于矩阵推导的模型参数求解方法基本思路：分离变量使得模型函数可以通过矩阵表示，然后基于最小二乘进行理论推导，最终根据参与运算的数据点直接推导出可以使误差最小的参数，整个过程不需要迭代运算。基于矩阵的推导方法的目的同迭代方法一样，也是获得最优的模型参数。因此分析还是要基于底回波混响相关模型。目前所采用的底回波模型是经过添加椭圆系数以及旋转角后的模型，模型如式(5)所示。首要目标是把底回波模型中的参数分离，转化为矩阵形式：R12(d12V)=()exp -j0+()2Vz+z12c exp-()ke2()2Vx+x12cos+()2Vy+y12sin2a()2Vy+y12cos-()2Vx+x12sin2(5)图3 最优参数求解流程图Fig.3 Flow chart of solving optimal parameters5432023 年声学技术对式(5)进行归一化处理后，可得：R12(d12V)=exp-()2Vx+x12cos+()2Vy+y12sin2+a()2Vy+y12cos-()2Vx+x12sin2(6)为了运算清晰，令x0=2Vx，y0=2Vy，替换后如式(7)所示：R=exp-()x0+x cos+()y0+y sin2+a()y0+y cos-()x0+x sin2(7)下面基于式(7)从理论上分析求解最优参数的可行性。首先对式(7)取对数，可得lnR=ln-()x0+x cos+()y0+y sin2+a()y0+y cos-()x0+x sin2(8)为了将运算转变为矩阵运算，对式(8)进行变量分离，结果为lnR=-()cos2+asin2 x2+()sin2+acos2 y2+2sincos()1-a xy+2()cos2+asin2 x0+2sincos()1-a y0 x+2()acos2+sin2 y0+2sincos()1-a x0y+()cos2+asin2 x20+()sin2+acos2 y20+2sincos()1-a x0y0-ln (9)式(9)已经完成了变量的分离。下一步将基于式(9)推导最优参数。假如参与拟合的矢量点有N个，则将这N个矢量数据点写成矩阵的形式为A=BC，其中：A为N1的列向量，其元素为ai=rilnri(10)B为N6型的矩阵，其行向量为bi=ri rixi riyi rixiyi rixi2 riyi2，i=1,N(11)C为61型的列向量：C=c1 c2 c3 c4 c5 c6T(12)式中：c1=-()cos2+asin2 x20+()sin2+acos2 y20+2sincos()1-a x0y0-lnc2=-2()cos2+asin2 x0+2sincos()1-a y0c3=-2()acos2+sin2 y0+2sincos()1-a x0c4=-2sincos()1-ac5=-()cos2+asin2c6=-()sin2+acos2(13)N点数据的误差列向量为：E=A-BC，采用最小二乘拟合，使 N 个数据点的均方误差(Mean Square Erro,MSE)最小，即：Ems=1NE22=1NETE=1N(A-BC)T(A-BC)(14)为了简化计算，将矩阵B进行QR分解，即B=QR，分解后Q是一个NN型的正交矩阵，R是一个N6型的上三角矩阵，对E=A-BC两边同时乘以QT，推导可得：QTE=QTA-QTBC=QTA-QTQRC=QTA-RC(15)误差项为E22=QTE22=QTA-RC22(16)令：QTA=STR=R10(17)其中，S为6维列向量，T为N6维列向量，R1为66型的上三角矩阵，则式(14)可以进一步转化为EMS=1NE22=1NQTA-RC22=1N(S-R1C22+T22)(18)从式(18)很容易可以看出，当S=R1C时，式(14)可以取得最小值，因此只需要求得最小值时的C即可，由S=R1C可得：C=R-11S(19)根据式(13)以及式(19)得到的结果解线性方程组，即可求得参数x0和y0：x0=-2c2c6+c4c34c5c6-c24y0=-2c5c3+c2c44c5c6-c24(20)已知x0=2Vx、y0=2Vy，再根据式(20)即可解算出速度信息。3数据验证3.1仿真数据验证通过一个二维旋转高斯模型产生样本数据，将样本数据输入新的模型参数求解方法，根据结果进行验证。采用的二维高斯旋转模型如式(21)所示。分析544第 4 期崔晓康等：相关测速声呐模型参数矩阵化求解方法结果如表1所示。f=exp-cos()x-x0+sin()y-y02-acos()y-y0-sin()x-x02(21)验证方案为固定模型参数，通过改变样本点数，针对不同样本点数得到的结果验证求解方案的可行性。仿真参数设定x0=0.5 y0=1=30，a=0.8，从中心点向外随机取点，分别在选取样本点为1010、88、55、33的情况下进行了测试，运算结果如表1所示。从表1可以看出，在不同采样点下，仿真测试结果都能正确地解算出目标参数，在没有噪声的情况下偏差基本为0。测试结果验证了基于矩阵推导求解最优参数的可行性。3.2实验数据验证采用湖上试验数据进行结果验证，试验设备为75 kHz声相关计程仪(ACL)，比测设备为GPS(型号为ProPak-G2)，湖况两级左右，测试水域水深65 m左右，各航速情况下的试验数据都在4 200帧左右。矩阵化处理程序和原程序处理结果如表2、3所示。由于矩阵化运算中舍弃了幅值调节因子，因此所使用的测速范围会有一定程度的下降。速度在7 kn(1 kn=1.852 kmh-1)以上测速偏差明显增大，表2中只列出7 kn及以下航速下的运算结果。对比表2、表3中不同航速下湖试数据的矩阵运算和原程序运算的结果可以看出，虽然矩阵运算结果在航速5 kn和7 kn时相对偏差比原程序大，但是偏差方向一致，而且整体测速结果稳定，测速偏差范围在0.2%以内。原程序运算结果在不同航速下偏差有正有负。在测速结果整体一致性上矩阵运算优于原程序求解。另外，以上运算所使用的样本点数(接收阵元两两之间构成相关矢量，对应的相关值被称为样本点)是基于大量数据验证，是原程序求解方法最优的选择。考虑到矩阵求解法舍弃了幅值调节因子，造成一部分样本点的模型失配，对此采取了减少样本点和增加样本点的尝试。减少样本点是尽量利用满足模型的样本点进行运算。增加样本点是在样本点不足的情况下，样本点新增和带来的模型失配之间会有一个最佳样本点数，使得测速取得最优结果。通过大量的数据验证，本文获得了以下经验结果：(1)基于式(1)可以看出，测速范围与相关延时和相关矢量相关，阵元布局决定了相关矢量，当相关延时加倍时，测速范围就会减半，也就是说相同的相关矢量数量，所能表示的测速范围降低为原来的一半，这样会使得同样速度范围内可以利用的样本点数增多，因此通过将相关延时更改为2来使相关运算的衰减速度降低，增加可以利用的相关点数；(2)原程序是在原模型下进行运算，在12个样本点时结果最优。而本文基于矩阵运算的方法，对使用的模型进行了一定调整，因此需要通过调整参与运算的样本点数来搜索更适合矩阵运算的样本点数。通过对比不同点数下的测速结果，发现2延时下采用15点做速度解算得到的测速结果最优(速度小于7 kn)，而且测速均方差也降低到原来的70%左右。另一方面，对航速710 kn情况也针对不同点数做了进一步分析。发现当参与相关运算的点数为13时测速结果稳定。不同速度下的矩阵运算测速结果如表4所示。整体偏差方向一致，不同速度下的相对偏差在0.3%以内。以上结果进一步验证了矩阵化求解算法的可行性。表2不同航速下湖试数据的矩阵运算结果Table 2Matrix calculation results of lake-trial data at dif-ferent speeds航速/kn5678910偏差/%0.380.450.49均方差34.235.239.5注：1 kn=1 852 kmh-1。表3不同航速下湖试数据的原程序处理结果Table 3Original program processing results of lake-trial da-ta at different speeds航速/kn5678910偏差/%0.040.590.210.230.190.18均方差22.821.920.120.919.517.4表1仿真数据验证结果Table 1Verification results of simulation data样本点数101088554433x0.50.50.50.50.5y11111cos20.750.750.750.750.750 001a0.80.80.80.80.85452023 年声学技术4结 论本文通过理论分析得到了基于矩阵推导的最优参数的解析解，并通过仿真验证和湖上试验数据，验证了该方案作为一种速度解算方法的可行性。另外通过对时间复杂度分析，原程序由于需要求逆运算和多次迭代，时间复杂度为MN3量级，其中M代表迭代运算次数。基于矩阵推导的运算方法虽然也涉及矩阵求逆，但是为上三角矩阵求逆，整体时间复杂度为N2量级，其中乘法运算和加法运算次数相当。因此，基于矩阵推导的最优参数求解方法相较于基于拟牛顿法的迭代参数求解，在时间复杂度上具有一定的优势。本文由于无法应用幅值调节因子，因此带来了一定的模型失配问题，导致文中所提测速方案没有达到最优的测速结果。后续研究将考虑从提升可以利用的样本点数和提高矩阵推导算法的容错性两方面进一步深入研究，以获得更加精确的估计参数。参考文献1王长红.宽带多普勒测速原理与技术研究D.北京:中国科学院声学研究所,1996.WANG Changhong.Research on principle and technology of broadband Doppler velocimetryD.Beijing:Institute of Acoustics,Chinese Academy of Sciences,1996.2DICKEY F R Jr.Velocity measuring correlation sonarJ.The Journal of the Acoustical Society of America,1981,70(3):914.3冯雷.相关测速声呐理论、信号处理研究及海试验证D.北京:中国科学院大学,2005.4蒯多杰,王长红,冯雷,等.海底回波空间相关特性研究J.声学学报(中文版),2009,34(5):385-395.KUAI Duojie,WANG Changhong,FENG Lei,et al.Spatial correlation analysis of sea-bottom backscatteringJ.Acta Acustica,2009,34(5):385-395.5易卉芹,陈龙,王长红.相关测速声呐时空相关函数模型改进研究J.声学学报,2011,36(5):484-488.YI Huiqin,CHEN Long,WANG Changhong.The research of temporal-spatial correlation function model improvement of correlation sonarJ.Acta Acustica,2011,36(5):484-488.6王映春,王长红,陈龙,等.相关测速声呐阵元坐标估计算法研究J.声学学报,2013,38(1):29-34.WANG Yingchun,WANG Changhong,CHEN Long,et al.The research of element coordinates estimating algorithm of correlation sonarJ.Acta Acustica,2013,38(1):29-34.7闫善勇.相关测速声呐时空相关特性及信号处理算法D.北京:中国科学院大学,2019.YAN Shanyong.Spatial-temporal correlation characteristics and signal processing algorithm of correlated velocimeter so-narD.Beijing:University of Chinese Academy of Sciences,2019.8闫善勇,赵二亮,邱薇,等.基于加权最小二乘的声相关计程仪速度精确估计J.声学技术,2019,38(5):502-507.YAN Shanyong,ZHAO Erliang,QIU Wei,et al.Accurate ve-locity estimation of acoustic correlation log based on weighted least squaresJ.Technical Acoustics,2019,38(5):502-507.9BRADLEY S E,DEINES K L,ROWE F D.Acoustic correla-tion current profilerC/IEEE Journal of Oceanic Engineer-ing.IEEE,2002:408-414.10 DENBIGH P N.A design study for a correlation log to mea-sure speed at seaJ.Journal of Navigation,1982,35(1):160-184.11 ZHU W Q,WANG C H,HUANG Y,et al.Theoretical model and error analysis of acoustic correlation current profiler(AC-CP)C/Oceans 99.MTS/IEEE.Riding the Crest into the 21st Century.Conference and Exhibition.Conference Proceed-ings(IEEE Cat.No.99CH37008).Seattle,WA,USA.IEEE,2002:382-385vol.1.12 FENG L,WANG C H,GONG L H,et al.Deep scattering layer observation by acoustic correlation current profilerC/OCEANS 2008-MTS/IEEE Kobe Techno-Ocean.Kobe,Ja-pan.IEEE,2008:1-4.13 朱正宜,邱薇,陈龙,等.相关测速声纳垂向速度测量与底深的关系J.海洋技术,2010,29(3):44-46,58.ZHU Zhengyi,QIU Wei,CHEN Long,et al.The relationship between vertical velocity measurement of acoustic correlation sonar and marine depthJ.Ocean Technology,2010,29(3):44-46,58.表4模型适当调整后不同航速下湖试数据的矩阵运算结果Table 4Matrix calculation results of lake-trial data at dif-ferent speeds after appropriate model adjustments航速/kn5678910偏差/%0.350.410.470.580.600.39均方差24.725.125.827.325.023.6546