CN113009419B - 一种基于频域互相关匹配的目标深度估计方法 - Google Patents

一种基于频域互相关匹配的目标深度估计方法 Download PDF

Info

Publication number
CN113009419B
CN113009419B CN202110210715.XA CN202110210715A CN113009419B CN 113009419 B CN113009419 B CN 113009419B CN 202110210715 A CN202110210715 A CN 202110210715A CN 113009419 B CN113009419 B CN 113009419B
Authority
CN
China
Prior art keywords
target
acoustic array
depth
acoustic
frequency domain
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202110210715.XA
Other languages
English (en)
Other versions
CN113009419A (zh
Inventor
张波
李风华
张燕君
余炎欣
朱飞龙
李文
翟铎
刘若芸
杨习山
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Institute of Acoustics CAS
Original Assignee
Institute of Acoustics CAS
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Institute of Acoustics CAS filed Critical Institute of Acoustics CAS
Priority to CN202110210715.XA priority Critical patent/CN113009419B/zh
Publication of CN113009419A publication Critical patent/CN113009419A/zh
Application granted granted Critical
Publication of CN113009419B publication Critical patent/CN113009419B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S5/00—Position-fixing by co-ordinating two or more direction or position line determinations; Position-fixing by co-ordinating two or more distance determinations
    • G01S5/18—Position-fixing by co-ordinating two or more direction or position line determinations; Position-fixing by co-ordinating two or more distance determinations using ultrasonic, sonic, or infrasonic waves
    • G01S5/28—Position-fixing by co-ordinating two or more direction or position line determinations; Position-fixing by co-ordinating two or more distance determinations using ultrasonic, sonic, or infrasonic waves by co-ordinating position lines of different shape, e.g. hyperbolic, circular, elliptical or radial
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/14—Fourier, Walsh or analogous domain transformations, e.g. Laplace, Hilbert, Karhunen-Loeve, transforms
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/15—Correlation function computation including computation of convolution operations
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Data Mining & Analysis (AREA)
  • Theoretical Computer Science (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Algebra (AREA)
  • General Engineering & Computer Science (AREA)
  • Software Systems (AREA)
  • Databases & Information Systems (AREA)
  • Computing Systems (AREA)
  • Remote Sensing (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Evolutionary Biology (AREA)
  • Operations Research (AREA)
  • Probability & Statistics with Applications (AREA)
  • Measurement Of Velocity Or Position Using Acoustic Or Ultrasonic Waves (AREA)

Abstract

本发明公布了一种基于频域互相关匹配的目标深度估计方法，包括：步骤1)利用两条近似平行的声学阵列从被测目标辐射噪声数据中提取归一化平均频域互相关函数，即频域伪格林函数；步骤2)在深度和距离上分别划分网格，根据声学阵列布放海域的水文环境，计算各网格点上对应的频域伪格林函数；步骤3)将步骤1)提取的频域伪格林函数与各网格点上计算得到的频域伪格林函数进行匹配处理，获得水中目标的深度估计结果。本发明的方法将目标深度估计和目标距离估计分离，可以在未知目标距离且不对目标距离进行估计的情况下获得目标深度估计结果，克服传统匹配场处理方法需要同时估计目标距离和目标深度或者先估计出目标距离然后再估计目标深度的限制。

Description

一种基于频域互相关匹配的目标深度估计方法
技术领域
本发明属于海洋声学领域，具体涉及一种基于频域互相关匹配的目标深度估计方法。
背景技术
目标的深度信息是非常重要的目标参数。根据声学传感器接收信号给出可靠的目标深度估计结果，是解决水中目标深度分辨的有效途径。匹配场处理方法是对目标进行深度估计的方法之一。这种方法充分利用了声场的干涉特性，采用声传播模型来计算拷贝场向量，然后将拷贝场与测量场进行匹配处理，从而获得未知变量的匹配处理结果。自从1976年匹配场(MFP)定位概念提出以来，国内外已发表了数百篇相关论文，进行了大量的海上实验，取得了不少进展。但是，海洋环境条件十分复杂，环境参数往往是时变、空变的，匹配场处理方法面临海洋环境参数无法准确测量的难题，从而拷贝场计算不准确，导致匹配处理结果与真实结果偏差太大，甚至无法给出置信度超过设定阈值的处理结果。也有学者提出利用声能流有功分量和无功分量正负号的变化来判别水上水下目标的方法，该方法要求将声学传感器布设在特定深度上，此深度需要根据声场分析频率和水深、底质、声速剖面等海洋环境计算得到，而实际应用中很难精确计算此深度，限制了该方法的使用。另外，可以利用脉冲信号的频散特征对目标声源进行深度估计，这种方法要求声源发出的是时域上很短的脉冲信号，而实际应用中所接收到的往往是连续的辐射噪声信号，各类船只很少会辐射这类持续时间很短但强度很大的瞬态脉冲信号。
发明内容
针对声学目标深度分辨这一需求，本发明根据浅海低频宽带声场的水平纵向互相关函数的物理特征，提出了一种基于频域互相关匹配的目标深度估计方法，适用于浅海海域的低频宽带被动目标深度估计，为提高浅海海域的目标深度分辨能力提供技术手段。该方法根据浅海低频宽带声场的水平纵向互相关函数的物理特征，利用两条近似平行的声学阵列提取归一化平均频域互相关函数，即频域伪格林函数，通过将实际提取的伪格林函数与在不同深度上仿真计算的伪格林函数进行匹配处理，给出水中目标的深度估计结果。
为了实现上述目的，本发明提供了一种基于频域互相关匹配的目标深度估计方法，所述方法包括：
步骤1)利用两条近似平行的声学阵列从被测目标辐射噪声数据中提取归一化平均频域互相关函数，即频域伪格林函数；
步骤2)在深度和距离上分别划分网格，根据声学阵列布放海域的水文环境，计算各网格点上对应的频域伪格林函数；
步骤3)将步骤1)提取的频域伪格林函数与各网格点上计算得到的频域伪格林函数进行匹配处理，获得水中目标的深度估计结果。
作为上述方法的一种改进，所述步骤1)具体包括：
步骤1-1)对声学阵列A一段时间内的多拍数据进行处理，获得声学阵列A的宽带波束输出结果；获得某一目标相对于声学阵列A的随时间t变化的方位历程θA(t)；
步骤1-2)对声学阵列B一段时间内的多拍数据进行处理，获得声学阵列B的宽带波束输出结果；获得某一目标相对于声学阵列B的随时间t变化的方位历程θB(t)；
步骤1-3)若|θA(t)-θB(t)|≤α，α为阈值；转入步骤4)；否则改变时间段，转入步骤1)，直到在宽带波束输出结果中有满足|θA(t)-θB(t)|≤α的目标出现；
步骤1-4)对两条声学阵列在目标方位的波束形成输出结果进行预处理；对预处理后的波束数据进行互相关运算，并在一定时间内积累，求取其平均值，即为伪格林函数。
作为上述方法的一种改进，所述步骤1-1)具体包括：
步骤1-1-1)将声学阵列A各阵元接收到的声信号分成多拍，每拍时间长度为ΔT，前后相邻两拍重叠P％，并对每一拍信号做傅里叶变换得到频域信号；声学阵列A的频域信号为sA(f；t)＝[sA1(f；t),sA2(f；t),...,sAM(f；t)]T；f是信号频率，t是当前拍数据的起始时刻；M是声学阵列A的阵元个数；(·)T表示转置；
步骤1-1-2)以声学阵列A的几何中心为参考点，计算声学阵列A各频点的在波束形成扫描方位θ的加权向量wA(f,θ)＝[wA0(f,θ),wA1(f,θ),...,wAM-1(f,θ)]T，其中
Figure BDA0002952167050000021
其中，τAm(θ)表示声学阵列A的第m+1个阵元相对于参考点的时延：
Figure BDA0002952167050000031
v(θ)＝-[cosθ,sinθ]T，
pAm＝[pxm,pym]T，
v(θ)是入射信号的单位向量，pAm为第m+1个阵元相对于参考点的二维坐标；
步骤1-1-3)计算声学阵列A指向方位θ的波束形成结果bA(f,θ；t)：
bA(f,θ；t)＝wA H(f,θ)sA(f；t)
其中，(·)H表示共轭转置；θ的取值范围为[θ0-45°，θ0+45°]，θ0表示声学阵列A和声学阵列B的参考阵元连线的角度；
步骤1-1-4)对声学阵列A的一段时间[t1，t2]内的多拍数据进行步骤1-1-1)至步骤1-1-3)的处理，获得声学阵列A的宽带波束输出结果HA(θ,t)：
Figure BDA0002952167050000032
其中，符号“Σ”表示求和；
步骤1-1-5)根据声学阵列A的宽带波束输出结果，获得某一目标相对于声学阵列A的随时间t变化的方位历程θA(t)。
作为上述方法的一种改进，所述步骤1-2)具体包括：
步骤1-2-1)将声学阵列B各阵元接收到的声信号分成多拍，每拍时间长度为ΔT，前后相邻两拍重叠P％，并对每一拍信号做傅里叶变换得到频域信号；声学阵列B的频域信号为sB(f；t)＝[sB1(f；t),sB2(f；t),...,sBN(f；t)]T；f是信号频率，选取的范围为[fL,fH]；t是当前拍数据的起始时刻；N是声学阵列B的阵元个数；
步骤1-2-2)以声学阵列B的几何中心为参考点，计算声学阵列B各频点的在波束形成扫描方位θ的加权向量wB(f,θ)＝[wB0(f,θ),wB1(f,θ),...,wB,N-1(f,θ)]T，其中
Figure BDA0002952167050000033
其中τBn(θ)表示声学阵列B的第n+1个阵元相对于参考点的时延：
Figure BDA0002952167050000034
v(θ)＝-[cosθ,sinθ]T，
pBn＝[pxn,pyn]T，
v(θ)是入射信号的单位向量，pBn为第n+1个阵元相对于参考点的二维坐标；
步骤1-2-3)计算声学阵列B指向方位θ的波束形成结果
bB(f,θ；t)＝wB H(f,θ)sB(f；t)
其中，(·)H表示共轭转置；
步骤1-2-4)对声学阵列B的一段时间[t1，t2]内的多拍数据进行步骤1-2-1)至步骤1-2-3)的处理，计算声学阵列B的宽带波束输出结果HB(θ,t)：
Figure BDA0002952167050000041
步骤1-2-5)根据声学阵列B的宽带波束输出结果获得同一目标相对于声学阵列B的随时间t变化的方位历程θB(t)。
作为上述方法的一种改进，所述步骤1-4)具体包括：
步骤1-4-1)对两条声学阵列在目标方位的波束形成输出结果进行预处理，预处理后声学阵列A、B在目标方位的波束形成输出
Figure BDA0002952167050000042
分别为：
Figure BDA0002952167050000043
Figure BDA0002952167050000044
其中，k为波数，k＝2πf/c0；L为声学阵列A、B几何中心水平间距；
Figure BDA0002952167050000045
为分析时间段内起始时刻所选目标相对于声学阵列B的水平方位；
IA(f)为所选取目标相对于声学阵列A的平均声强谱：
Figure BDA0002952167050000046
IB(f)为所选取目标相对于声学阵列B的平均声强谱：
Figure BDA0002952167050000047
其中，Nt为总拍数；Nt的取值范围为[600/(ΔT(1-P％))，1800/(ΔT(1-P％))]；ΔT为单拍数据的时间长度，单位s；P％为前后相邻两拍重叠率；
步骤1-4-2)对预处理后的波束数据进行互相关运算，并在一定时间内积累，求取其平均值，获得频域伪格林函数
Figure BDA0002952167050000051
Figure BDA0002952167050000052
(·)*表示共轭。
作为上述方法的一种改进，所述步骤2)具体包括：
步骤2-1)在深度上按照[0,4H/5]间隔划分为Qz个网格，H为水深；在距离上按照[Rmin,Rmax]间隔划分为QR个网格，Rmin为满足声学阵列远场要求所需的最小距离，Rmax为声学阵列对目标的最大理论探测距离；
步骤2-2))根据水深、声速剖面和海底参数这些海洋环境参数，采用声场计算方法计算简正波信息：第l号简正波的水平波数kl＝μl+iηl，μl、ηl分别为kl的实部、虚部；第l号简正波的特征函数Ψl(z)；其中，l＝1,2，…，I；I为简正波总号数；
步骤2-3)根据简正波信息计算各网格点的频域伪格林函数：
第q个网格点的频域伪格林函数为：
Figure BDA0002952167050000053
其中，B(z1,z2,zq)为归一化因子：
Figure BDA0002952167050000054
z1,z2表示两个水听器的深度，其水平间距为L，两个水听与声源的水平间距分别为r1,r2；zq、rq分别为第q个网格点对应的声源深度和与两条声学阵列几何中心连线中点的水平距离，q＝1,2,…Q，Q为网格总数，Q＝Qz×QR。
作为上述方法的一种改进，深度网格数Qz的取值范围为[8H/5λ,40H/5λ]，H为水深，λ为分析带宽内中心频率所对应的声波波长；距离网格数QR的取值范围为[(Rmax-Rmin)/10,(Rmax-Rmin)]。
作为上述方法的一种改进，所述步骤3)具体包括：
步骤3-1)第q个网格点的频域伪格林函数与步骤1)提取的频域伪格林函数
Figure BDA0002952167050000061
的相关系数W(zq,rq)：
Figure BDA0002952167050000062
(·)*表示共轭；
步骤3-2)对于所有的W(zq,rq)，其最大值其对应网格点的深度即为被测目标的深度估计结果。
本发明的优势在于：
1、本发明的方法将目标深度估计和目标距离估计分离，可以在未知目标距离且不对目标距离进行估计的情况下获得目标深度估计结果,克服了传统匹配场处理方法需要同时估计目标距离和目标深度或者先估计出目标距离然后再估计目标深度的限制；
2、本发明的方法对海洋环境参数精确度要求较低，因此提高了方法的宽容性。
附图说明
图1为本发明中水听器及声源的相对位置示意图；
图2是本发明的基于频域互相关匹配的目标深度估计方法的流程图；
图3为本发明一个实施例中海水平均声速剖面；水深95m，海水平均声速剖面如图2所示，海面声速1544m/s,海面至65m深度之间声速随深度逐渐减小至1522m/s,65m至95m深度近似为等温层；海底声速1608m/s；
图4为本发明一个实施例中两条声学阵列相对位置图；右上的点线表示声学阵列A，左下的点线表示声学阵列B；声学阵列A、B均包含41个阵元；声学阵列A、B几何中心水平距离为720m；两条声学阵列(声学阵列A、B)布放于海底表面；
图5是本发明一个实施例中声学阵列A的目标方位历程图；灰度图为声学阵列A波束形成输出结果，分析带宽为100Hz至110Hz；虚线为根据灰度图获得的目标方位随时间的变化曲线；横轴为水平方位角，单位度，纵轴为时间；
图6(a)是本发明一个实施例中获得的频域伪格林函数实部随频率和时间变化的二维图像，横轴为频率，单位Hz，纵轴为绝对时间；
图6(b)为04:10时刻频域伪格林函数实部随频率的变化曲线；
图7为02:30时刻相关系数W(zq,rq)二维图，横轴为水平距离rq，纵轴为深度zq；该图中，W(zq,rq)的局部极大值位于5m深度处，如图中“+”所示，据此判断被测目标的深度估计值为5m；
图8为目标深度估计结果随时间的变化曲线；纵轴为深度，单位m；横轴为时间；该时间段内目标深度估计最大为12.5m，最小为2.5m，均值为4.24m，该测深结果与根据船舶自动识别系统的该目标为水面商船的信息相符；
图9为根据船舶自动识别系统信息获得的被测目标相对于声学阵列的水平距离随时间的变化曲线；根据船舶自动识别系统信息，可判断被测目标为水面目标，且该水面目标对于声学阵列的水平距离从70km均匀变化至10km；
图10为水下人工声源相关系数W(zq,rq)二维图，横轴为水平距离rq，纵轴为深度zq；该图中，W(zq,rq)的局部极大值位于40m深度处，如图中“+”所示，据此判断被测目标的深度估计值为40m；该人工声源实际吊放深度为50m。
具体实施方式
针对声学目标深度分辨这一需求，本发明根据浅海低频宽带声场的水平纵向互相关函数的物理特征，提出了一种基于频域互相关匹配的目标深度估计方法，适用于浅海海域的低频宽带被动目标深度估计，为提高浅海海域的目标深度分辨能力提供技术手段。
本发明所采用的技术方案是：
根据简正波理论，深度为zs的点声源在水平距离r、深度z处的产生的声场(格林函数)可表示为
Figure BDA0002952167050000071
其中，km为第m号简正波的水平波数，km＝μm(f)+jηm(f)，μm、ηm分别为km的实部、虚部，均为频率f的函数。j为虚数符号；e为自然对数的基底，e＝2.718……。
如图1所示，两个圆点表示水听器，其水平间距为L，深度分别为z1、z2，与声源(五角星所示)水平间距分别为r1、r2；一个圆点位于两个水听器的中心位置，与声源水平间距为R；声源相对于两个水听器连线的方位为θ。
对于方位为θ的声源，利用图1所示两个水听器，可获得频域的归一化平均互相关函数
Figure BDA0002952167050000081
其中B(z1,z2,zq)为归一化因子
Figure BDA0002952167050000082
符号“Σ”为求和算子，π为圆周率。
鉴于公式(2)所表示的归一化平均互相关函数与公式(1)表示的格林函数的相似性，我们称之为频域“伪格林函数”。
在深度和距离上划分网格，然后根据海洋环境参数和公式(2)计算得到的各网格点的伪格林函数。将实际提取的伪格林函数与各网格点上计算得到的伪格林函数进行匹配处理，获得水中目标的深度估计结果。
需要指出的是，以上论述是针对两个水听器的情形的。但是一般情况下，为了提高处理增益，实际应用中通常用两个近似平行的声学阵列代替上述的两个水听器，以获得空间增益，提高信噪比。
如图2所示，本发明提出了一种基于频域互相关匹配的目标深度估计方法，包括如下步骤：
S1，处理实际数据，提取频域伪格林函数。
S1-1，对布设在水中的两条近似平行的声学阵列(声学阵列A、B)所接收的声学数据分别进行波束形成，其波束形成结果记为bA(f,θ；t)和bB(f,θ；t)。其中，f是信号频率，选取的范围为[fL,fH]；t是当前拍数据的起始时刻；θ为预成波束方位。
优选的，fL、fH的典型取值范围为[2c0/H,20c0/H]。c0为参考声速，通常取为水体声速的均值。H为水深。
优选的，θ的典型取值范围为[θ0-45°，θ0+45°]。θ0表示声学阵列A和声学阵列B的参考阵元连线的角度。
S1-2，根据上一步结果，选取声学阵列正横方位附近的水面船只作为辐射噪声的声源，其相对于声学阵列A、B的随时间t变化的方位历程分别记为θA(t)、θB(t)。
如果所分析的时间段内有多个目标出现，应尽量选取声学阵列正横方位附近的目标。
θA(t)、θB(t)应满足|θA(t)-θB(t)|≤2°。
S1-3，对两条声学阵列在目标方位的波束形成输出结果进行预处理。预处理后声学阵列A、B在目标方位的波束形成输出分别为
Figure BDA0002952167050000091
Figure BDA0002952167050000092
其中，k为波数，k＝2πf/c0；L为声学阵列A、B几何中心水平间距；
Figure BDA0002952167050000093
为分析时间段内起始时刻所选目标相对于声学阵列B的水平方位；
Figure BDA0002952167050000094
Figure BDA0002952167050000095
Nt为总拍数。
优选的，Nt的典型取值范围为[600/(ΔT(1-P％))，1800/(ΔT(1-P％))]。ΔT为单拍数据的时间长度，单位s；P％为前后相邻两拍重叠率。
S1-4，对预处理后的波束数据进行互相关运算，并在一定时间内积累，求取其平均值，获得频域伪格林函数
Figure BDA0002952167050000096
Figure BDA0002952167050000097
(·)*表示共轭。
S2，根据声学阵列布放海域的水文环境，仿真计算频域伪格林函数。
S2-1，在深度和距离上划分网格。将深度在[0,4H/5]范围内划分为QZ个网格，H为水深。将距离在[Rmin,Rmax]范围内划分为QR个网格，Rmin为满足声学阵列远场要求所需的最小距离，Rmax为声学阵列对目标的最大理论探测距离。
设网格总数为Q＝QR×QZ，QR为距离网格数，QZ为深度网格数。
优选的，QZ的典型取值范围为[8H/5λ,40H/5λ]。H为水深，λ为分析带宽内中心频率所对应的声波波长。
优选的，QR的典型取值范围为[(Rmax-Rmin)/10,(Rmax-Rmin)]。Rmin、Rmax的单位为km。
S2-2，根据水深、声速剖面、海底参数等海洋环境参数，采用声场计算程序(比如Kraken或KrakenC等)计算简正波信息：第m号简正波的水平波数km＝μm+jηm，μm、ηm分别为km的实部、虚部；第m号简正波的特征函数Ψm(z)。其中m＝1,2，…，M。M为简正波总号数。
S2-3，根据简正波信息，计算各网格点的频域伪格林函数。对于第q个网格点，其频域伪格林函数为
Figure BDA0002952167050000101
其中B(z1,z2,zq)为归一化因子
Figure BDA0002952167050000102
zq、rq分别为第q个网格点对应的声源深度和至两条声学阵列几何中心连线中点的水平距离；
Figure BDA0002952167050000103
为步骤S1中获得的目标方位。
S3，匹配处理，获得目标深度。
S3-1，计算每个网格点的频域伪格林函数与实际数据提取的频域伪格林函数
Figure BDA0002952167050000104
的相关系数。对于第q个网格点，其对应的相关系数为
Figure BDA0002952167050000105
(·)*表示共轭。
S3-2，对于所有的W(zq,rq)，q＝1,2,3，…Q，取局部最大值，其对应网格点的深度即为被测目标的深度估计结果。
由于上述算法对目标距离非常不敏感，因此局部最大值对应网格点的距离与目标真实距离可能有较大差异，因此不能作为目标距离的估计结果。
如果目标距离是已知的，记为rs，那么仅对目标距离附近范围内的网格点进行局部最大值搜索即可。距离搜索的典型范围为[0.8rs,1.2rs]。
本发明根据浅海低频宽带声场的水平纵向互相关函数的物理特征，提出了一种基于伪格林函数的目标深度频域匹配估计算法，适用于浅海海域的低频宽带被动目标深度估计，为提高浅海海域的目标深度分辨能力提供技术手段。
现通过附图和具体实施例来阐述本发明的技术方案。
2020年5月在南海北部海域进行了一次海上试验。水深95m，海水平均声速剖面如图3所示，海面声速1544m/s,海面至65m深度之间声速随深度逐渐减小至1522m/s，65m至95m深度近似为等温层。海底声速1608m/s。两条声学阵列(声学阵列A、B)布放于海底表面，阵形布置如图4所示。
分析2020年5月19日0点至5点时间段内的数据。该时间段内，将数据分成每拍10s的多拍数据，相邻拍重叠40％。声学阵列A的波束形成输出HA(θ,t)如图5所示，在47°方位有目标，其方位随时间的变化曲线如图5中虚线所示。分别获得声学阵列A、B在目标方位的波束输出并进行归一化处理，得
Figure BDA0002952167050000111
和
Figure BDA0002952167050000112
在频域做互相关运算，并在一定时间段内取平均，得频域伪格林函数
Figure BDA0002952167050000113
本实施例中，Nt取100，对应时间段长度为10分钟；频率分析范围为100Hz至150Hz。图6(a)显示了
Figure BDA0002952167050000114
的实部随频率和时间变化的二维图像，横轴为频率，单位Hz，纵轴为时间；图6(b)为04:10时刻
Figure BDA0002952167050000115
的实部随频率变化曲线。
将深度在[2.5m,80m]范围内划分为36个网格，步长2.5m。将距离在[10km,70km]范围内划分为71个网格，步长1km。根据图3所示的海域环境参数，采用Kraken程序计算简正波信息，并据此计算各网格点的频域伪格林函数。对于第q个网格点，其频域伪格林函数为
Figure BDA0002952167050000116
其中B(z1,z2,zq)为归一化因子：
Figure BDA0002952167050000121
zq、rq分别为第q个网格点对应的声源深度和至两条声学阵列几何中心连线中点的水平距离；
Figure BDA0002952167050000122
为目标方位。
本实施例中，声学阵列A、B深度z1＝z2＝95m，声学阵列A、B几何中心的水平间距L＝720m，目标方位
Figure BDA0002952167050000123
将每个网格点的频域伪格林函数
Figure BDA0002952167050000124
与实际数据提取的频域伪格林函数
Figure BDA0002952167050000125
进行相关运算，其相关系数为
Figure BDA0002952167050000126
(·)*表示共轭。对于所有的W(zq,rq)，q＝1,2,3，…Q，取局部最大值，其对应网格点的深度即为被测目标的深度估计结果。
02:30时刻，相关系数W(zq,rq)的二维图如图7所示，横轴为水平距离rq，纵轴为深度zq。该图中，W(zq,rq)的局部极大值位于5m深度处，如图中“+”所示，据此判断被测目标的深度估计值为5m。
图8为目标深度估计结果随时间的变化曲线。纵轴为深度，单位m；横轴为时间。该时间段内目标深度估计最大为12.5m，最小为2.5m，均值为4.24m，该测深结果与根据船舶自动识别系统的该目标为水面商船的信息相符。根据船舶自动识别系统信息，该水面目标相对于声学阵列的水平距离随时间的变化曲线如图9所示，在本实施例所对应时间段内，该水面目标对于声学阵列的水平距离从70km均匀变化至10km。
另外，采用上述方法，对声学阵列接收的人工吊放声源发射的宽度声学信号进行了处理，处理频段为100Hz-150Hz。图10为其相关系数W(zq,rq)二维图，横轴为水平距离rq，纵轴为深度zq。该图中，W(zq,rq)的局部极大值位于40m深度处，如图中“+”所示，据此判断被测目标的深度估计值为40m，与声源实际吊放深度(50m)基本相符。
上述结果验证了方法的有效性。
最后所应说明的是，以上实施例仅用以说明本发明的技术方案而非限制。尽管参照实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，对本发明的技术方案进行修改或者等同替换，都不脱离本发明技术方案的精神和范围，其均应涵盖在本发明的权利要求范围当中。

Claims (7)

1.一种基于频域互相关匹配的目标深度估计方法，所述方法包括：
步骤1)利用两条近似平行的声学阵列从被测目标辐射噪声数据中提取归一化平均频域互相关函数，即频域伪格林函数；
步骤2)在深度和距离上分别划分网格，根据声学阵列布放海域的水文环境，计算各网格点上对应的频域伪格林函数；
步骤3)将步骤1)提取的频域伪格林函数与各网格点上计算得到的频域伪格林函数进行匹配处理，获得水中目标的深度估计结果；
所述步骤1)具体包括：
步骤1-1)对声学阵列A一段时间内的多拍数据进行处理，获得声学阵列A的宽带波束输出结果；获得某一目标相对于声学阵列A的随时间t变化的方位历程θA(t)；
步骤1-2)对声学阵列B一段时间内的多拍数据进行处理，获得声学阵列B的宽带波束输出结果；获得某一目标相对于声学阵列B的随时间t变化的方位历程θB(t)；
步骤1-3)若|θA(t)-θB(t)|≤α，α为阈值；转入步骤4)；否则改变时间段，转入步骤1)，直到在宽带波束输出结果中有满足|θA(t)-θB(t)|≤α的目标出现；
步骤1-4)对两条声学阵列在目标方位的波束形成输出结果进行预处理；对预处理后的波束数据进行互相关运算，并在一定时间内积累，求取其平均值，即为伪格林函数。
2.根据权利要求1所述的基于频域互相关匹配的目标深度估计方法，其特征在于，所述步骤1-1)具体包括：
步骤1-1-1)将声学阵列A各阵元接收到的声信号分成多拍，每拍时间长度为ΔT，前后相邻两拍重叠P％，并对每一拍信号做傅里叶变换得到频域信号；声学阵列A的频域信号为sA(f；t)＝[sA1(f；t),sA2(f；t),...,sAM(f；t)]T；f是信号频率，t是当前拍数据的起始时刻；M是声学阵列A的阵元个数；(·)T表示转置；
步骤1-1-2)以声学阵列A的几何中心为参考点，计算声学阵列A各频点的在波束形成扫描方位θ的加权向量wA(f,θ)＝[wA0(f,θ),wA1(f,θ),...,wAM-1(f,θ)]T，其中
Figure FDA0003208689190000011
其中，τAm(θ)表示声学阵列A的第m+1个阵元相对于参考点的时延：
Figure FDA0003208689190000021
v(θ)＝-[cosθ,sinθ]T，
pAm＝[pxm,pym]T，
v(θ)是入射信号的单位向量，pAm为第m+1个阵元相对于参考点的二维坐标，c0为参考声速；
步骤1-1-3)计算声学阵列A指向方位θ的波束形成结果bA(f,θ；t)：
bA(f,θ；t)＝wA H(f,θ)sA(f；t)
其中，(·)H表示共轭转置；θ的取值范围为[θ0-45°，θ0+45°]，θ0表示声学阵列A和声学阵列B的参考阵元连线的角度；
步骤1-1-4)对声学阵列A的一段时间[t1，t2]内的多拍数据进行步骤1-1-1)至步骤1-1-3)的处理，获得声学阵列A的宽带波束输出结果HA(θ,t)：
Figure FDA0003208689190000022
其中，符号“Σ”表示求和；
步骤1-1-5)根据声学阵列A的宽带波束输出结果，获得某一目标相对于声学阵列A的随时间t变化的方位历程θA(t)。
3.根据权利要求2所述的基于频域互相关匹配的目标深度估计方法，其特征在于，所述步骤1-2)具体包括：
步骤1-2-1)将声学阵列B各阵元接收到的声信号分成多拍，每拍时间长度为ΔT，前后相邻两拍重叠P％，并对每一拍信号做傅里叶变换得到频域信号；声学阵列B的频域信号为sB(f；t)＝[sB1(f；t),sB2(f；t),...,sBN(f；t)]T；f是信号频率，选取的范围为[fL,fH]；t是当前拍数据的起始时刻；N是声学阵列B的阵元个数；
步骤1-2-2)以声学阵列B的几何中心为参考点，计算声学阵列B各频点的在波束形成扫描方位θ的加权向量wB(f,θ)＝[wB0(f,θ),wB1(f,θ),...,wB,N-1(f,θ)]T，其中
Figure FDA0003208689190000031
其中τBn(θ)表示声学阵列B的第n+1个阵元相对于参考点的时延：
Figure FDA0003208689190000032
v(θ)＝-[cosθ,sinθ]T，
pBn＝[pxn,pyn]T，
v(θ)是入射信号的单位向量，pBn为第n+1个阵元相对于参考点的二维坐标；
步骤1-2-3)计算声学阵列B指向方位θ的波束形成结果
bB(f,θ；t)＝wB H(f,θ)sB(f；t)
其中，(·)H表示共轭转置；
步骤1-2-4)对声学阵列B的一段时间[t1，t2]内的多拍数据进行步骤1-2-1)至步骤1-2-3)的处理，计算声学阵列B的宽带波束输出结果HB(θ,t)：
Figure FDA0003208689190000033
步骤1-2-5)根据声学阵列B的宽带波束输出结果获得同一目标相对于声学阵列B的随时间t变化的方位历程θB(t)。
4.根据权利要求3所述的基于频域互相关匹配的目标深度估计方法，其特征在于，所述步骤1-4)具体包括：
步骤1-4-1)对两条声学阵列在目标方位的波束形成输出结果进行预处理，预处理后声学阵列A、B在目标方位的波束形成输出
Figure FDA0003208689190000034
分别为：
Figure FDA0003208689190000035
Figure FDA0003208689190000036
其中，k为波数，k＝2πf/c0；L为声学阵列A、B几何中心水平间距；
Figure FDA0003208689190000037
为分析时间段内起始时刻所选目标相对于声学阵列B的水平方位；
IA(f)为所选取目标相对于声学阵列A的平均声强谱：
Figure FDA0003208689190000041
IB(f)为所选取目标相对于声学阵列B的平均声强谱：
Figure FDA0003208689190000042
其中，Nt为总拍数；Nt的取值范围为[600/(ΔT(1-P％))，1800/(ΔT(1-P％))]；ΔT为单拍数据的时间长度，单位s；P％为前后相邻两拍重叠率；
步骤1-4-2)对预处理后的波束数据进行互相关运算，并在一定时间内积累，求取其平均值，获得频域伪格林函数
Figure FDA0003208689190000043
Figure FDA0003208689190000044
(·)*表示共轭。
5.根据权利要求4所述的基于频域互相关匹配的目标深度估计方法，其特征在于，所述步骤2)具体包括：
步骤2-1)在深度上按照[0,4H/5]间隔划分为Qz个网格，H为水深；在距离上按照[Rmin,Rmax]间隔划分为QR个网格，Rmin为满足声学阵列远场要求所需的最小距离，Rmax为声学阵列对目标的最大理论探测距离；
步骤2-2))根据水深、声速剖面和海底参数这些海洋环境参数，采用声场计算方法计算简正波信息：第l号简正波的水平波数kl＝μl+iηl，μl、ηl分别为kl的实部、虚部；第l号简正波的特征函数Ψl(z)；其中，l＝1,2，…，I；I为简正波总号数；
步骤2-3)根据简正波信息计算各网格点的频域伪格林函数：
第q个网格点的频域伪格林函数为：
Figure FDA0003208689190000045
其中，B(z1,z2,zq)为归一化因子：
Figure FDA0003208689190000051
z1,z2表示两个水听器的深度，其水平间距为L，两个水听与声源的水平间距分别为r1,r2；zq、rq分别为第q个网格点对应的声源深度和与两条声学阵列几何中心连线中点的水平距离，q＝1,2,…Q，Q为网格总数，Q＝Qz×QR，Qz为深度网格数，QR为距离网格数。
6.根据权利要求5所述的基于频域互相关匹配的目标深度估计方法，其特征在于，深度网格数Qz的取值范围为[8H/5λ,40H/5λ]，H为水深，λ为分析带宽内中心频率所对应的声波波长；距离网格数QR的取值范围为[(Rmax-Rmin)/10,(Rmax-Rmin)]。
7.根据权利要求6所述的基于频域互相关匹配的目标深度估计方法，其特征在于，所述步骤3)具体包括：
步骤3-1)第q个网格点的频域伪格林函数与步骤1)提取的频域伪格林函数
Figure FDA0003208689190000052
的相关系数W(zq,rq)：
Figure FDA0003208689190000053
(·)*表示共轭；
步骤3-2)对于所有的W(zq,rq)，其最大值其对应网格点的深度即为被测目标的深度估计结果。
CN202110210715.XA 2021-02-25 2021-02-25 一种基于频域互相关匹配的目标深度估计方法 Active CN113009419B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110210715.XA CN113009419B (zh) 2021-02-25 2021-02-25 一种基于频域互相关匹配的目标深度估计方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110210715.XA CN113009419B (zh) 2021-02-25 2021-02-25 一种基于频域互相关匹配的目标深度估计方法

Publications (2)

Publication Number Publication Date
CN113009419A CN113009419A (zh) 2021-06-22
CN113009419B true CN113009419B (zh) 2021-11-09

Family

ID=76386192

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110210715.XA Active CN113009419B (zh) 2021-02-25 2021-02-25 一种基于频域互相关匹配的目标深度估计方法

Country Status (1)

Country Link
CN (1) CN113009419B (zh)

Citations (14)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2011056830A1 (en) * 2009-11-03 2011-05-12 Michael Bruno Passive acoustic underwater intruder detection system
CN103513249A (zh) * 2012-06-20 2014-01-15 中国科学院声学研究所 一种宽带相干模基信号处理方法及系统
CN104678384A (zh) * 2013-11-28 2015-06-03 中国科学院声学研究所 一种波束域的声压差互相关谱分析水下目标速度估计方法
CN104820218A (zh) * 2015-05-26 2015-08-05 中国科学院声学研究所东海研究站 一种基于频域自相关的浅海海底单参数反演方法
CN107085216A (zh) * 2017-04-20 2017-08-22 国家海洋局第三海洋研究所 一种基于单水听器的深海水声被动测距测深方法
CN107329114A (zh) * 2017-06-21 2017-11-07 歌尔股份有限公司 声源定位方法和装置
CN109443516A (zh) * 2018-12-25 2019-03-08 西北工业大学 一种基于噪声场垂直振速信号的海底声速被动获取方法
CN109489799A (zh) * 2018-12-25 2019-03-19 西北工业大学 一种基于双矢量水听器的海底声速分步反演方法
CN109815942A (zh) * 2019-03-18 2019-05-28 西北工业大学 基于海洋环境噪声信号的简正波特征提取方法
CN110146846A (zh) * 2019-06-06 2019-08-20 青岛理工大学 一种声源位置估计方法、可读存储介质及计算机设备
CN110554357A (zh) * 2019-09-12 2019-12-10 苏州思必驰信息科技有限公司 声源定位方法和装置
CN110824429A (zh) * 2019-10-28 2020-02-21 西北工业大学 深海环境下利用非同步垂直阵的宽带声源被动定位方法
CN111323824A (zh) * 2020-03-04 2020-06-23 中国科学院声学研究所 一种模块式阵列声波测井仪及井外界面探测方法
CN111580048A (zh) * 2020-05-09 2020-08-25 中国科学院声学研究所 一种利用单矢量水听器的宽带声源深度估计方法

Family Cites Families (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7315488B2 (en) * 2006-06-06 2008-01-01 Raytheon Company Methods and systems for passive range and depth localization
EP2263097B1 (en) * 2008-03-12 2013-06-19 Raytheon Company Autonomous sonar system and method
CN105467364B (zh) * 2015-11-20 2019-03-29 百度在线网络技术（北京）有限公司 一种定位目标声源的方法和装置
CN112034441B (zh) * 2020-09-14 2023-11-24 中国科学院声学研究所 基于格林函数二维解卷积的阵不变量被动定位方法及系统

Patent Citations (14)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2011056830A1 (en) * 2009-11-03 2011-05-12 Michael Bruno Passive acoustic underwater intruder detection system
CN103513249A (zh) * 2012-06-20 2014-01-15 中国科学院声学研究所 一种宽带相干模基信号处理方法及系统
CN104678384A (zh) * 2013-11-28 2015-06-03 中国科学院声学研究所 一种波束域的声压差互相关谱分析水下目标速度估计方法
CN104820218A (zh) * 2015-05-26 2015-08-05 中国科学院声学研究所东海研究站 一种基于频域自相关的浅海海底单参数反演方法
CN107085216A (zh) * 2017-04-20 2017-08-22 国家海洋局第三海洋研究所 一种基于单水听器的深海水声被动测距测深方法
CN107329114A (zh) * 2017-06-21 2017-11-07 歌尔股份有限公司 声源定位方法和装置
CN109443516A (zh) * 2018-12-25 2019-03-08 西北工业大学 一种基于噪声场垂直振速信号的海底声速被动获取方法
CN109489799A (zh) * 2018-12-25 2019-03-19 西北工业大学 一种基于双矢量水听器的海底声速分步反演方法
CN109815942A (zh) * 2019-03-18 2019-05-28 西北工业大学 基于海洋环境噪声信号的简正波特征提取方法
CN110146846A (zh) * 2019-06-06 2019-08-20 青岛理工大学 一种声源位置估计方法、可读存储介质及计算机设备
CN110554357A (zh) * 2019-09-12 2019-12-10 苏州思必驰信息科技有限公司 声源定位方法和装置
CN110824429A (zh) * 2019-10-28 2020-02-21 西北工业大学 深海环境下利用非同步垂直阵的宽带声源被动定位方法
CN111323824A (zh) * 2020-03-04 2020-06-23 中国科学院声学研究所 一种模块式阵列声波测井仪及井外界面探测方法
CN111580048A (zh) * 2020-05-09 2020-08-25 中国科学院声学研究所 一种利用单矢量水听器的宽带声源深度估计方法

Non-Patent Citations (14)

* Cited by examiner, † Cited by third party
Title
Arrival-time structure of the time-averiaged ambient noise cross-correlation function in an oceanic waveguide;Sabra, KG，等;《JOURNAL OF THE ACOUSTICAL SOCIETY OF AMERICA》;20000128;第164-174页 *
Geoacoustic inversion of tow-ship noise via near-field-matched-field processing;D.J. Battle，等;《IEEE Journal of Oceanic Engineering》;20031027;第28卷(第3期);第454-467页 *
Improving resolution and SNR of correlation function with the increase in bandwidth of recorded noise fields during estimation of bottom profile of ocean;Md. Jahangir Alam，等;《OCEANS"10 IEEE SYDNEY》;20100527;第1-7页 *
Performance Evaluation of Acoustic Model-Based Blind Channel Estimation in Ocean Waveguides;Wei Feng，等;《IEEE Access》;20180423(第6期);第27239-27250页 *
Source Localization With Multiple Hydrophone Arrays via Matched-Field Processing;Dag Tollefsen，等;《EEE Journal of Oceanic Engineering》;20130731;第42卷(第3期);第654-662页 *
Using ocean ambient noise for array self-localization and self-synchronization;K.G. Sabra，等;《IEEE Journal of Oceanic Engineering 》;20051024;第30卷(第2期);第338-347页 *
一种基于warping变换的浅海脉冲声源被动测距方法;王冬, 等;《物理学报》;20160413(第10期);第1-9页 *
一种基于模态域波束形成的水平阵被动目标深度估计;李鹏，等;《物理学报》;20170401;第66卷(第8期);第1-13页 *
利用矢量海洋环境噪声提取声场格林函数;周建波; 等;《声学学报》;20190515(第3期);第337-344页 *
基于EMD的噪声互相关提取时域格林函数研究;王奉宝;《电子技术与软件工程》;20140604(第10期);第82-84页 *
基于路径选择的深海水下运动目标被动深度估计;刘炎堃;等;《应用声学》;20200522;第1-8页 *
基于频域互相关的LED晶粒阵列快速视觉检测;李鹤喜，等;《计量学报》;20180722(第4期);第30-34页 *
水下宽带脉冲声源的单水听器被动测距方法;李晓曼，等;《中国声学学会水声学分会2019年学术会议论文集》;20190525;全文 *
海面波浪对海洋环境噪声垂直空间相关性的影响;张乾初, 等;《声学学报》;20190315(第2期);第189-200页 *

Also Published As

Publication number Publication date
CN113009419A (zh) 2021-06-22

Similar Documents

Publication Publication Date Title
CN113011006B (zh) 一种基于互相关函数脉冲波形匹配的目标深度估计方法
CN108226933B (zh) 一种基于条纹干涉结构的深海宽带目标深度估计方法
CN107179535A (zh) 一种基于畸变拖曳阵的保真增强波束形成的方法
CN112987004B (zh) 一种浅海环境下基于水平阵列的水面水下目标分类方法
CN111580048B (zh) 一种利用单矢量水听器的宽带声源深度估计方法
CN110824429B (zh) 深海环境下利用非同步垂直阵的宽带声源被动定位方法
CN108828522A (zh) 一种利用垂直阵lcmv波束形成的水下目标辐射噪声测量方法
CN101915922A (zh) 拖曳线列阵被动测距方法
CN108845325A (zh) 拖曳线列阵声纳子阵误差失配估计方法
CN111025273B (zh) 一种畸变拖曳阵线谱特征增强方法及系统
CN111537982B (zh) 一种畸变拖曳阵线谱特征增强方法及系统
CN114280541B (zh) 一种基于深海分布式垂直线列阵的目标被动定位方法
CN109100711B (zh) 一种深海环境下单基地主动声纳低运算量三维定位方法
CN104793212A (zh) 利用声波海底反射实现主动声纳远程探测的方法
CN109061654B (zh) 一种深海环境下单圆环阵主动三维定位方法
CN104714235A (zh) 一种低频双矢量水听器阵列的测距方法及系统
CN112098938B (zh) 一种基于六元锥矢量阵的水声目标降维匹配声场定位方法
CN115656994B (zh) 双基地有源探测拖曳阵阵形实时校准方法
CN113009419B (zh) 一种基于频域互相关匹配的目标深度估计方法
CN111679248A (zh) 一种基于海底水平l型阵列的目标方位和距离联合稀疏重构定位方法
CN113009418B (zh) 一种基于伪格林函数脉冲时延的目标深度估计方法
AU2016228476B2 (en) Method for calculating a confidence echo signal that is exempt from multipath propagation effects and for determining a distance and/or a direction to an echo source and device and vehicle
CN113009417B (zh) 利用声场干涉特性的海底声学阵列阵形估计方法
CN108761465A (zh) 一种基于距离维波束形成的三元阵被动测距方法
CN113138365B (zh) 一种基于深度学习的单矢量水听器方位估计方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant