Passive ranging method for shallow water based on dual depth Warping transform information fusion
-
摘要:
针对浅海环境中利用单深度接收信号自相关函数的Warping变换方法测距时易出现的多值性问题, 提出一种双深度Warping变换融合测距方法。该方法将双深度测量场的Warping变换提取结果与具有深度差异的拷贝场本征函数信息融合, 判别简正波干涉项模态, 提高测距的准确度。此外, 针对实际海洋环境中测量场Warping变换特征峰位置难以准确提取的问题, 将Warping变换的频率不变性与图像形态学方法结合以增强特征频率结构, 提高后续测距的稳健性。海试数据处理结果显示, 测距结果与实际距离符合较好, 表明所提方法具有实际可行性。
Abstract:In shallow-water environments, to address the issue of multi-value in estimating the distance of a sound source using the autocorrelation function Warping transform of single-depth receiver signal, a method of dual-depth Warping transform information fusion is proposed for distance measurement. This method integrates the Warping transform extraction results from the dual-depth measurement field with depth-differentiated eigenfunction information derived from the copy field, enabling the identification of the order of the normal modal interference terms, consequently enhancing the accuracy of distance measurement. Additionally, there is another need in accurately extracting the characteristic peak positions of the Warping transform in actual marine environments. In this study, the frequency invariance of Warping transform is jointly applied with image morphology to strengthen the characteristic frequency structure, thereby enhancing the robustness of ranging algorithm. The real-data experiments show that the ranging results align well with the ground truth, indicating that the proposed passive ranging method is feasible for practical application.
-
引言
水下无源测距隐蔽性强、安全性高, 相关研究备受关注, 但在无源测距的实际应用中, 接收信号常被环境噪声或平台噪声污染, 同时海洋环境参数也难以准确获知, 在海洋环境噪声干扰和参数不确知条件下如何实现有效测距成为探索无源测距新方法的重要研究方向。
根据简正波理论[1], 水声信号在波导中的传播呈现频散特性, 一些学者提出通过消频散变换[2-3]或采用时频分析方法[4]来克服模态的非线性群延迟。近年来, 在简正波模态提取[5-6]、参数反演[7-9]以及无源定位[10-12]时, 常将Warping变换与简正波的频散特性相结合[13], 对接收信号进行时域重采样, 使其转化为一系列脉冲信号。针对接收信号Warping变换, Touze等提出了Pekeris信道环境下的Warping变换算子[14], Bonnel等利用Warping变换估计各阶简正波时延差实现无源测距[15], Lopatka等利用Warping变换进行特征提取实现了宽带声源定位[16]。为拓展Warping变换的适用范围, 牛海强等修正了非理想波导条件下的Warping算子[17], Brown等提出适用于分层声学环境的时域Warping算子并在仿真条件下的中纬度深海环境中证明了其有效性[18]。由于接收信号Warping变换存在信号传播时间未知的问题, 戚聿波等推导了信号自相关函数Warping变换表达式并应用于单水听器无源测距[19], 并且引入波导不变量提出β-Warping变换实现了浅海脉冲声源的测距[20]。王冬等进一步讨论了负跃层环境下自相关函数Warping变换与β-Warping的适用条件[21], 提出对能量密度函数进行Warping变换并应用于无源测距[22]。Warping变换也用于获取声源或海洋环境信息, 刘志韬等利用自相关函数Warping变换实现浅海的声源深度判决[23], 董阁等对自相关函数进行时域Warping变换得到距离特征量, 并将其作为新的观测信息应用于目标运动分析中, 估计目标运动信息[24], 高伟等将β-Warping变换与奇异值分解相结合, 使用一维宽带声强度估计波导不变量[25]。
对接收信号进行Warping变换的效果依赖于对信号到达时间的准确估计, 而在无源测距中这一参数是未知的, 戚聿波等指出了自相关函数Warping变换简正波相干项的频率不变特征[26]。王冬等分析了海洋环境参数失配对自相关函数Warping变换测距结果的影响[27], 其中海深失配影响较大, 声速剖面的小幅度偏移和海底底质参数失配影响较小, 说明自相关函数Warping变换方法对海洋环境参数失配具有一定的宽容性。但实际应用中, 简正波干涉项阶数判断往往需要结合引导声源或其他方法, 孟瑞洁等提出通过水平阵端射方向获取水平波数差结合拷贝场对干涉项进行模态判决[28]。此外, 基于Warping变换的测距方法依赖于特征频率的峰值位置提取, 易受环境噪声的干扰。因此, 在利用自相关函数Warping变换进行无源测距时, 需要找到一种能够有效判别简正波干涉项阶数, 且在有背景噪声干扰时仍可成功测距的方法。
为解决无源测距的多值性问题, 本文提出一种通过垂直双阵元接收信号联合实现浅海无源测距的方法, 将双接收深度的本征函数差异与双接收深度的接收信号自相关函数Warping变换结果差异相匹配, 判别简正波相干项的阶数并实现测距。该方法可以减弱或消除使用单接收深度Warping变换进行测距时因简正波号数未知而产生的测距模糊。低信噪比下, Warping变换性能下降、特征峰值难以提取, 可利用Warping变换的频率不变性, 通过图像处理方法获取特征频率, 以进一步提高无源测距稳健性。
1. 浅海波导自相关函数Warping变换
根据简正波理论, 假设波导水平不变, 远场条件下浅海水中声场可以表示为一系列简正波叠加的形式:
p(r,f,z)=S(f)iρ(zs)√8π re−iπ/4M∑m=1ψm(zs)ψm(z)eikrm(f)r√krm(f), (1) 其中,
S(f) 为声源频谱,r 为声源与接收点之间的水平距离,zs 为声源深度,z 为接收深度,ψm(z) 是第m阶简正波在深度z处的本征函数值,krm=βm+iαm 为第m阶简正波的水平波数。不同阶简正波同一频率下的水平波数krm 不同, 同一阶简正波不同频率下的水平波数也不相同, 导致简正波在传播过程中出现频散现象。令
Am(r,f,z)=ie−iπ/4ρ(zs)√8πkrm(f)rψm(zs)ψm(z) 表示第m阶简正波的幅度, 声压表达式可简化为
p(r,f,z)=S(f)M∑m=1Am(r,f,z)eikrmr, (2) 则接收信号自相关函数的频域表达式(即功率谱密度)为
ξ(r,f,z)=|S(f)|2⋅(M∑m=1Am(z)A∗m(z)+M∑m=1M∑n≠mAm(z)A∗n(z)ei(krm−krn)r). (3) 根据信道传播特性, 简正波的幅度随频率是缓变的, 在有限带宽内, 式(3)中第一项的变化可以忽略, 自相关函数频域表达式的第二项为干涉起伏项, 可简化为
˜ξ(r,f,z)=|S(f)|2M∑m=1M∑n≠mAm(z)A∗n(z)ei(krm−krn)r. (4) 假设声源谱级在一定带宽内不变或缓变, 对式(4)在有限带宽内进行逆傅里叶变换, 并将结果在时间上延迟
tr , 可得自相关函数的时域表达式:˜ξ(rs,t,z)=∫+∞−∞|S(f)|2(M∑m=1M∑n≠mAm(z)A∗n(z)⋅ei[(krm−krn)r−2π ftr])ei2π ftdf, (5) 其中,
tr = rs/c0 为延迟时间,c0 为海水的平均声速,rs 为声源与接收器的水平距离。无源测距, 即估计接收距离rs (未知量)。在理想波导条件下, 第m阶简正波的水平波数可表示为
krm=2π c√f2−f2cm, (6) 其中,
fcm 为第m阶简正波的截止频率。将式(6)代入式(5), 通过泰勒一阶近似与稳相法近似[29], 可得˜ξ(rs,t,z)=M∑m=1M∑n≠mAm(tμmn√t2−t2r)A∗n(tμmn√t2−t2r)⋅tr√μmn(t2−t2r)34ei(2π μmn√t2−t2r) , (7) 其中,
μmn=√f2cn−f2cm 表示第m阶与第n阶简正波干涉项的特征频率。令
h(t)=√t2+t2r 表示时域的Warping算子, 对自相关函数进行Warping变换, 可得(Wh˜ξ)(rs,t,z)=M∑m=1M∑n≠mAm(√t2+t2rμmnt)A∗n(√t2+t2rμmnt)⋅tr√μmnt32ei(2π μmnt) . (8) 经傅里叶变换后,
(Wh˜ξ)(rs,t,z) 转换为频域上的一系列脉冲。由式(8)可知, 当
tr 已知时, Warping变换后第m, n阶简正波的特征频率只与特征频率μmn 有关。令rs 表示真实的声源距离, 设定一个参考声源距离r0 ,˜μmn 表示参考声源距离下得到的第m, n阶简正波的特征频率, 则特征频率μmn 与˜μmn 存在如下关系:˜μmn=√rsr0μmn. (9) 因此, 在无源测距时, Warping变换的特征频率受简正波干涉项特征频率
μmn 和参考声源距离r0 这两个因素的影响。对接收信号的自相关函数进行Warping变换, 若提取出简正波测量场干涉项的特征频率˜μmn , 且已知其对应的简正波干涉项阶数, 便可通过声场模型计算理论特征频率μmn , 并结合式(9)实现目标的无源测距。2. Warping变换测距方法
2.1 单接收深度Warping变换测距分析
当浅海波导只存在一阶简正波时, 无法根据自相关函数Warping变换结果得到特征峰, 本文只考虑传播简正波阶数大于1的情况。
自相关函数的Warping变换存在干涉项特征频率相似的问题, 特征频率相近会影响简正波干涉项的阶数判断。因此在得到Warping变换理论特征频率后, 需要利用层次聚类将相近的特征频率归为一类, 直至每一类的类内特征频率差值均小于某一阈值, 并用
{{\boldsymbol{F}}_c}\left( k \right) 表示第k类的聚类结果,{f_c}\left( k \right) 表示类内特征频率的均值, 图1为层次聚类的示意图。对单接收信号进行自相关Warping变换后, 可根据幅值提取特征频率峰。由于不能有效地判断提取的特征频率对应的简正波号数, 可将每个峰值与聚类后的每一类特征频率的中心频率通过式(9)计算得到一系列的声源距离估计值后再进行筛选。如表1所示,
{{{\rm X}}_p} 为提取的第p个特征峰对应的频率,{r_{p - q}} 为假设第p个特征峰对应于第q簇理论特征频率时计算得到的声源距离估计值。表 1 不同特征频率对应的声源距离估计值假设距离{a_r} {f_c}\left( {\text{1}} \right) … {f_c}\left( K \right) {X_1} {r_{1 - 1}} … {r_{1 - K}} {X_2} {r_{2 - 1}} … {r_{2 - K}} {X_3} {r_{3 - 1}} … {r_{3 - K}} 若提取的特征峰的数量为H, 将提取出的一组峰值按照频率从大到小排列。由于提取的每个特征峰频率归属于不同的特征频率簇, 并且每个特征峰频率对应的实际声源距离估计值同属于一个距离范围, 因此可以根据这两个条件对声源估计值进行筛选, 即每个特征峰对应的声源距离估计值满足式(10),则认为其同属于一个距离范围。除此之外, 这些簇的中心频率的大小顺序和间隔关系理论上都与提取的特征峰的频率序列相对应。结合以上信息筛选出特征峰序列对应的聚类标签和距离, 匹配过程如图2所示。
\frac{{\max \left( {\left| {{r_{m,i}} - {r_{n,j}}} \right|} \right)}}{{\max \left( {r{ _{m,i}}} \right)}} < {d_{{\text{thresh}}}} , (10) 其中, m和n为提取的特征峰编号,
m \in [1, \cdots ,H] ,n \in [1, \cdots ,H] , i, j为理论特征频率对应的簇的编号i \in [1, \cdots ,N] ,j \in [1, \cdots ,N] , 且m \ne n ,i \ne j ,{d_{{\text{thresh}}}} 的值设置为0.1。在提取自相关函数Warping变换的特征峰的幅值和位置时, 设定幅值最大值的0.7倍为峰值筛选的阈值, 并且设置特征峰提取数量上限为3, 若高于阈值的局部极大值数量超过3, 则将局部极大值的幅值从高至低排序, 选取前3个幅值最高的局部极大值, 记录其对应的频率以及相应的幅值。采用水体等声速的Pekeris波导仿真条件, 海深为80 m, 水体声速1500 m/s, 海底声速为1600 m/s, 海底密度为1.9 g/cm3, 声源深度为5 m, 频率范围为100~200 Hz。计算150 Hz处的本征函数值, 虚部幅值明显增大, 因此水体中主要存在前5阶简正波。设单接收阵元的接收深度为30 m, 分别仿真水平距离为11.0 km和12.5 km的接收信号, 参考声源距离为8.0 km, 进行自相关函数Warping变换后结果如图3所示, 设定层次聚类的频率间隔阈值为1 Hz, 红色圆圈表示筛选出的特征频率峰。
水平距离的变化范围为11~16 km时, Warping变换频谱幅值如图4(a)所示, 单深度接收信号的测距结果如图4(b)所示, 其中圆圈为测距结果, 虚线为实际距离, 由于特征频率峰对应的模态阶数信息的缺失, 测距结果存在着多值性, 无法直接判断出正确结果, 因此标记出所有可能的距离值。由图4可以看出, 若可提取的特征频率峰值数量较多时, 仅通过本节所提的峰值提取与筛选也能得到正确的测距结果。图3(a)中根据Warping变换后峰值的强度关系仅提取出一个特征频率, 在仅有接收深度信息以及海洋环境信息的情况下, 只能根据该特征频率以及理论特征频率, 通过式(9)计算声源距离估计值。由于特征峰对应简正波干涉项模态数相关信息的缺失, 测距的结果存在多值性, 若无法知道特征频率对应的简正波干涉项阶数, 测距结果与真实距离之间会存在较大偏差。
2.2 双深度信息融合无源测距方法
自相关函数的频谱幅值可表示为
B\left( {f,{{\textit z}_s},{\textit z}} \right) = \frac{{{\psi _m}\left( {{{\textit z}_s}} \right){\psi _m}\left( {\textit z} \right){\psi _n}\left( {{{\textit z}_s}} \right){\psi _n}\left( {\textit z} \right)}}{{8\pi r{\rho ^2}({{\textit z}_s})\sqrt {{k_{rm}}\left( f \right){k_{rn}}\left( f \right)} }} . (11) 在同一接收距离、不同接收深度下, 特征频率的强度受本征函数值
{\psi _m}\left( {\textit z} \right) 影响, 即特征频率幅值具有深度依赖性。在水平不变波导下,{\psi _m}\left( {\textit z} \right) 在一定的频率范围内随频率是缓变的, 在较窄的频率范围内可以近似将其视为一个常量, 在不同接收深度条件下, 特征频率的强度主要取决于\left| {{\psi _m}\left( {{{\textit z}_s}} \right){\psi _m}\left( {\textit z} \right){\psi _n}\left( {{{\textit z}_s}} \right){\psi _n}\left( {\textit z} \right)} \right| 的值。而对于同一对简正波干涉项, 同一声源在两个接收深度点处的{\psi _m}\left( {{{\textit z}_s}} \right){\psi _n}\left( {{{\textit z}_s}} \right) 相同, 只需要考虑两个接收深度之间{\psi _m}\left( {\textit z} \right){\psi _n}\left( {\textit z} \right) 的本征函数信息差异即可。若能将两个接收深度的本征函数信息差异与Warping变换结果差异实现有效对应, 则可以对测距过程中的多值现象进行消除。针对测量场的接收信号, 将两个深度的特征频率峰值根据频率由大到小进行组合和排序, 如图5所示, 只保留组合后符合频率间隔关系的单深度筛选结果,从而得到双深度接收条件下的筛选组合。
为了简化对比, 每个深度的特征频率幅值如表2所示, 其中
{{\varLambda }_{q,p}} 表示第q个深度接收信号特征频率的第p个峰值的幅值。根据{{\varLambda }_{1,p}} 和{{\varLambda }_{2,p}} 的大小关系组成一个仅包含+1和−1的序列C:表 2 双深度特征频率幅值序列接收深度 ① 号特征频率
峰幅值② 号特征频率
峰幅值③ 号特征频率
峰幅值深度1 {{\varLambda }_{1,1}} {{\varLambda }_{1,2}} {{\varLambda }_{1,3}} 深度2 {{\varLambda }_{2,1}} {{\varLambda }_{2,2}} {{\varLambda }_{2,3}} {C}\left( p \right) = \left\{ {\begin{array}{*{20}{c}} {{\text{ }}1,{\text{ }}{{\varLambda }_{1,p}} \geqslant {{\varLambda }_{{\text{2}},p}},} \\ { - 1,{\text{ }}{{\varLambda }_{1,p}} < {{\varLambda }_{{\text{2}},p}}.} \end{array}} \right. (12) 为了将Warping变换的深度差异信息与拷贝场的本征函数信息进行对应, 首先计算拷贝场的本征函数值, 得到如下信息:
\left[ {\begin{array}{*{20}{c}} 0&{{\psi _1}\left( {{{\textit z}_q}} \right){\psi _2}\left( {{{\textit z}_q}} \right)}& \cdots &{{\psi _1}\left( {{{\textit z}_q}} \right){\psi _M}\left( {{{\textit z}_q}} \right)} \\ 0&0& \ddots & \vdots \\ 0&0&0&{{\psi _{M - 1}}\left( {{{\textit z}_q}} \right){\psi _M}\left( {{{\textit z}_q}} \right)} \\ 0&0&0&0 \end{array}} \right]\left( {q \in [1,2]} \right) , (13) 其中, q表示接收阵元的编号, 靠近海面的阵元
q{\text{ = 1}} , 靠近海底的阵元q{\text{ = 2}} , M为水体中存在的简正波的数量。对式(13)中的上三角部分重新排列, 得到
{{\boldsymbol{\varPsi }}_R}\left( {\textit z} \right) = \left[ {{{\varGamma }_1},{{\varGamma }_2}, \cdots,{{\varGamma }_l} ,\cdots,{{\varGamma }_L}} \right], (14) 其中,
{ \varGamma}_{l} = {\psi }_{m}\left({\textit z}\right){\psi }_{n}\left({\textit z}\right),~{l} = {m} + \sum_{i=0}^{n-2}{i},~ 1 \leq m < n \leq M 。由于特征频率有重复和近似的可能, 按照特征频率层次聚类后的分组对本征函数矩阵
{{\boldsymbol{\varPsi }}_R} 进行分类和叠加:{{{\widehat \varPsi }}_R}\left( {k,{\textit z}} \right) = \sum\limits_{{f_l} \in {{\boldsymbol{F}}_c}\left( k \right)} {{{\varGamma }_l}} . (15) 由于在浅海波导中不同阶简正波随距离的衰减程度不同, 因此在考虑每一簇本征函数的作用时, 不能忽略简正波的衰减差异直接进行叠加。
根据简正波微扰理论, 当频率较低时, 海水的声衰减通常忽略不计, 导致衰减的主要因素是海底的吸收, 利用微扰法得到非理想波导下第m号简正波的衰减可表示为
{\alpha _m}\left( f \right) = \frac{1}{{{k_{rm}}}}\int_0^{D_{\rm{Sea}} }{\frac{{\alpha \left( {\textit z} \right)2{\text{π}}f\psi _m^2\left( {\textit z} \right)}}{{c\left( {\textit z} \right)\rho \left( {\textit z} \right)}}} {\text{d}}{\textit z} + {\text{δ}}{a_m} , (16) 其中,
{\text{δ}}{a_m} 表示由海底引起的吸收,f 为信号的频率,\alpha 表示水体的声吸收系数, DSea为海深。海底的本征函数可表示为[1]
{\psi _m}\left( {\textit z} \right) = {\psi _m}\left( {D_{\rm{Sea}} }\right){{\text{e}}^{ - {\gamma _m}\left( {{\textit z} - {D}_{\rm{Sea}}} \right)}} , (17) 其中,
{\gamma _m} = \sqrt {k_{rm}^2 - {{( {\omega /{c_b}} )}^2}} , cb为海底声速。由海底引起的吸收可表示为{\text{δ}}{\alpha _m} = \frac{{\psi _m^2\left( {D}_{\rm{Sea}} \right){\alpha _b}\omega }}{{2{k_{rm}}{c_b}{\rho _b}\sqrt {k_{rm}^2 - {{\left( {\omega /{c_b}} \right)}^2}} }} . (18) 为表述方便, 后面用
{\alpha _{bm}} 表示第m阶简正波的海底衰减。将简正波衰减因子引入简正波幅值公式, 可得{\widetilde A_m}\left( {r,f,{\textit z}} \right) = \frac{{{\text{i}}{{\text{e}}^{ - {{{\mathrm{i}}\pi }}/4}}}}{{\rho ({{\textit z}_s})\sqrt {8{\text{π }}{k_{rm}}\left( f \right)r} }}{\psi _m}\left( {{{\textit z}_s}} \right){\psi _m}\left( {\textit z} \right){{\text{e}}^{ - {\alpha _{bm}}r}} . (19) 自相关函数的频谱幅值重新表示为
\widetilde B\left( {r,f,{{\textit z}_s},{\textit z}} \right) = \frac{{{\psi _m}\left( {{{\textit z}_s}} \right){\psi _m}\left( {\textit z} \right){\psi _n}\left( {{{\textit z}_s}} \right){\psi _n}\left( {\textit z} \right)}}{{8{\text{π }}r{\rho ^2}({{\textit z}_s})\sqrt {{k_{rm}}\left( f \right){k_{rn}}\left( f \right)} }}{{\text{e}}^{ - \left( {{\alpha _{bm}} + {\alpha _{bn}}} \right)r}} . (20) 聚类后的本征函数矩阵元素也重新表示为
{{{\widetilde \varPsi }}_R}\left( {k,{\textit z}} \right) = \sum\limits_{{f_l} \in {F_c}\left( k \right)} {{{{\widetilde \varGamma }}_l}} , (21) 其中,
{\widetilde{\varGamma}_l} = {\varGamma }_{l}\left({\textit z}\right){\text{e}}^{-\left({\alpha }_{bm} + {\alpha }_{bn}\right)r} 。由此可以得到拷贝场下的本征函数信息
{{\boldsymbol{\widetilde \varPsi }}_R} , 根据式(22)将双深度的筛选结果的每一个组合生成一个序列D:D_{k}\left( p \right) = \left\{ \begin{array}{ll} 1, & {{{\widetilde {\varPsi }}}_R}\left({m_{kp}} {,{{\textit z}_1}} \right) \geqslant {{{\widetilde {\varPsi } }}_R}\left(m_{kp} ,{{\textit z}_2} \right), \\ - 1, & {{{\widetilde {\varPsi } }}_R}\left( m_{kp}{,{{\textit z}_1}} \right) < {{{\widetilde{\varPsi } }}_R}\left( m_{kp}{,{{\textit z}_2}} \right). \end{array} \right. (22) 将序列C与序列D结果进行序列匹配, 若满足如下条件:
{D_k} \cdot C = \left\| {{D_k}} \right\| , (23) 则认为此时
{D_k} 序列双深度筛选结果中对应的第k个测距结果{r_k} 为期望的测距结果, 若经过双深度序列匹配测距仍存在测距结果的多值可能, 根据简正波的衰减规律, 可判定对应理论特征频率最小的一组为测距结果。由于非理想波导的Warping变换结果呈现一定的窄带特性, 并且实际环境下易受到背景噪声的干扰, 若聚类后类间特征频率间隔较小, 例如随着海深的增加或频带升高时, 海水中存在的简正波阶数也会随之增加, 导致简正波干涉项的类间特征频率间隔过小, 影响干涉项模态阶数判别和测距方法的性能, 因此本文背景下选择的频率范围满足类间频率间隔大于2 Hz。
将双接收点分别放置于深度30 m, 40 m处, 信号频带为100~200 Hz, 参考声源距离为8 km, 当声源与接收器水平距离为11 km时, 对两个接收深度的接收信号分别进行自相关函数Warping变换, 结果如图6所示, 其中实线为深度1的Warping变换频谱, 虚线为深度2的Warping变换频谱, 深度1提取的特征频率峰用“+”表示, 深度2提取的特征频率峰用“○”表示。
由图6可知, 当存在两个垂直深度不同的接收点时, 由于本征函数的差异导致两个接收深度的Warping变换特征频率峰值出现的位置以及同一位置处出现的特征频率幅值都存在差异。根据两个接收深度处的本征函数以及理论衰减因子得到的序列D如图7所示。利用双深度信息融合方法进行无源测距, 可判断出图6中两个深度的特征频率分别对应于图7中的框线部分, 编号①与编号②的特征峰分别对应于第1, 3阶简正波干涉项与第2, 4阶简正波干涉项。当声源与接收阵的水平距离变化范围为11~16 km时, 得到的测距结果如图8所示。由图可知, 本节提出的双深度信息联合测距方法可以有效消除测距的多值现象。
3. 仿真分析
在声源频谱为随机连续谱并且存在背景噪声的条件下进行仿真, 验证第2节所提双深度联合测距方法对舰船辐射噪声连续谱信号测距的有效性。
3.1 Warping变换特征频率图像增强预处理
在水下接收信号中, 舰船频谱主要由线谱和连续谱组成。利用Warping变换进行舰船辐射噪声连续谱无源测距, 其性能会受到连续谱幅值的随机起伏变化和背景噪声影响。在2.1节仿真条件下, 当接收深度为30 m时, 在均匀声源谱(纯信号、不包含背景噪声)、随机声源谱(纯信号、不包含背景噪声)和随机声源谱包含0 dB随机噪声的条件下进行接收信号自相关函数Warping变换, 结果如图9所示。
虽然连续谱的随机性不完全满足Warping变换稳相近似条件, 但是在一定的信噪比条件下, 目标声信号干涉项的特征频率依然能够呈现出一定的条纹特性。由于环境噪声的时空随机特性, 当信噪比较低时, 其Warping变换结果会出现不同分布的亮点或亮斑, 可利用图像处理方法来增强条纹的结构清晰度: 首先通过合并J个相邻帧的Warping变换数据形成图像矩阵, 并通过维纳滤波自适应二值化[30]; 对二值化后的图像进行腐蚀[31]操作, 在去除杂点的同时尽量保留有结构的图像, 图像腐蚀所用的模板为
\left[ {\begin{array}{*{20}{c}} 0&1&0 \\ 1&1&1 \\ 0&1&0 \end{array}} \right]. 随后利用相同的模板对图像进行膨胀, 还原其中的条纹结构部分; 将此矩阵与原图像矩阵相点乘, 以频率轴为横轴, 距离轴为纵轴, 再将整个矩阵进行90°范围内的小角度积分, 得到最终增强结果, 图10为条纹增强示意图。
3.2 双深度联合测距仿真分析
设置声源深度为5 m, 声源与接收阵的水平距离为4~20 km, 接收阵的两个阵元深度分别为30 m和40 m, 并且加入0 dB高斯随机背景噪声。假设参考声源距离为8 km, 对接收信号自相关函数进行Warping变换, 如图11所示。使用3.1节提出的预处理方法对接收信号自相关函数Warping变换结果进行特征频率条纹结构增强, 图像增强帧数设置为J=30, 增强后的结果如图12所示。对自相关函数Warping变换结果矩阵和图像增强结果矩阵进行归一化, 提取特征频率条纹结构, 并计算条纹部分平均能量与背景部分平均能量的比值:
{P_e} = \frac{{{{\overline p}_s}}}{{{{\overline p}_b}}}, (24) 其中,
{\overline p_s} 为条纹部分的能量均值,{\overline p_b} 为背景部分的能量均值。当接收深度为30 m时, 图像增强前后的{P_e} 值分别为4.1和120.7; 当接收深度为40 m时, 图像增强前后的{P_e} 值分别为3.9和168.8。经过图像增强处理后, 两个接收深度的条纹与背景的平均能量比值分别提升了29倍和43倍。通过单深度Warping变换测距方法分别对两个阵元的接收信号进行测距, 结果如图13(a)(b)所示, 联合双接收深度信息的无源测距结果如图13(c)所示, 红色虚线为声源与接收器之间的实际距离。在单深度测距结果中, 同一时刻存在一个或多个测距结果, 而双深度信息融合无源测距的测距结果具有单值性, 在包含背景噪声时依然可以实现无源测距。双深度联合测距的相对误差如图14所示, 当声源与接收器之间的水平距离较小(小于7 km)时, Warping变换的特征峰易出现重叠和畸变, 并且在背景噪声的干扰下, 特征频率条纹更难提取, 导致无源测距方法的性能下降, 测距误差达10%。当声源与接收器之间的水平距离大于7 km时, 测距误差为4%。
令背景噪声从−20 dB变化至20 dB, 利用平均相对误差(相对测距误差的统计平均值)对本文所提测距方法在低信噪比时的测距有效性进行分析, 如图15所示。单深度阵元工作带宽内信噪比高于−2 dB时, 平均相对误差在10%以下, 但随着信噪比降低, 测距误差也随之增加至50%。
4. 海试数据分析
中国科学院声学研究所1999年3月在黄海海域进行了一次声传播测量实验, 实验海域海深约37 m, 海水声速剖面近似等声速, 水体声速为1477 m/s。仿真计算的水体密度为1.0
{\text{g/c}}{{\text{m}}^{\text{3}}} , 基底声速为1605 m/s, 基底密度为1.95{\text{g/c}}{{\text{m}}^{\text{3}}} , 基底衰减系数为0.164{{{\mathrm{dB}}/\lambda }} 。接收阵为等间距16阵元垂直阵, 信号处理的工作频带为200~300 Hz。接收阵元深度为5~31 m, 声源的距离变化范围为9~14 km。首先通过Kraken计算拷贝场, 在250 Hz处, 第6阶及更高阶数的简正波在所选频带内的水平波数虚部较大, 因此随距离衰减较大, 同时考虑到海底对高阶简正波的衰减作用, 此时实验海域中传播的主要为前5阶简正波。计算拷贝场的本征函数值, 选取接收深度为12 m和19 m的两个阵元进行分析, 对接收信号进行自相关Warping变换, 结果如图16所示。由于背景噪声和声源谱的联合影响, 导致Warping变换后的结果中存在很多杂点, 一些杂点的强度甚至与干涉的特征频率条纹的强度相近, 特征频率条纹图像处理增强后的结果如图17所示。当接收深度为12 m时, 图像增强前后的
{P_e} 值分别为2.5和161.5; 当接收深度为19 m时, 图像增强前后的{P_e} 值分别为1.7和23.7, 经过图像增强后条纹部分的平均能量与背景部分的平均能量分别提升了65倍和14倍。由图16可以看出, 在12 min内, 特征频率从50 Hz变化至40 Hz, 而图17中特征频率从80 Hz变化至65 Hz, 特征频率条纹所在的频率位置和强度均存在差异。下面进行无源测距。首先利用单深度测距方法, 分别对两个深度的条纹增强后的结果进行测距, 测距结果如图18(a)(b)所示, 尽管利用条纹增强手段提取出了特征频率, 但是两个阵元的单深度测距结果均存在测距模糊现象。结合两个深度的拷贝场的本征函数与实测场的特征频率位置和幅值差异进行双深度测距, 结果如图18(c)所示。测距结果的相对误差整体在5%以内, 如图19所示, 可见信息融合方法对舰船辐射噪声的连续谱无源测距具有明显的效果。
5. 结论
针对单深度阵元接收信号自相关函数的Warping变换测距方法的多值性问题, 提出了双深度信息融合的浅海Warping变换无源测距方法。在海洋环境参数以及双阵元接收深度已知的条件下, 该方法首先通过层次聚类和特征峰值筛选方法得到初步的单深度测距结果, 然后联合两个深度实测场的特征频率峰值强度差异和拷贝场的本征函数信息差异, 通过序列匹配方式确定特征频率峰的简正波干涉项模态阶数, 消除测距结果的多值性, 实现无源测距。为提高算法的抗噪性能, 还通过图像形态学方法增强Warping变换的特征频率条纹, 提高了测距方法的稳健性。
所提方法在海水声速随海深不变或缓变的情况下, 单深度接收信噪比高于−2 dB时的测距效果较好。海试数据测距结果的相对误差整体在5%以内, 验证了双深度信息融合方法的有效性。对于负跃层等声速随深度变化的浅海波导中的浅深度声源, 由于低阶简正波的相速度大于水体声速, 导致低阶简正波不激发或部分激发, 使自相关函数Warping变换失效, 这种情况下的测距方法还需进一步研究。
-
表 1 不同特征频率对应的声源距离估计值
假设距离{a_r} {f_c}\left( {\text{1}} \right) … {f_c}\left( K \right) {X_1} {r_{1 - 1}} … {r_{1 - K}} {X_2} {r_{2 - 1}} … {r_{2 - K}} {X_3} {r_{3 - 1}} … {r_{3 - K}} 表 2 双深度特征频率幅值序列
接收深度 ① 号特征频率
峰幅值② 号特征频率
峰幅值③ 号特征频率
峰幅值深度1 {{\varLambda }_{1,1}} {{\varLambda }_{1,2}} {{\varLambda }_{1,3}} 深度2 {{\varLambda }_{2,1}} {{\varLambda }_{2,2}} {{\varLambda }_{2,3}} -
[1] Jensen F B, Kuperman W A, Porter M B, et al. Computational ocean acoustics. Comput. Phys., 1995; 9(1): 55−56 DOI: 10.1063/1.4823373
[2] Gao D Z, Wang N, Wang H Z. A dedispersion tranform for sound propagation in shallow water waveguide. J. Comput. Acoust., 2010; 18(3): 245−257 DOI: 10.1142/S0218396X10004188
[3] 王宁, 高大治, 王好忠. 频散、声场干涉结构、波导不变量与消频散变换. 哈尔滨工程大学学报, 2010; 31(7): 825−831 DOI: 10.3969/j.issn.1006-7043.2010.07.002 [4] Bonnel J, Le Touze G, Nicolas B, et al. Physics-based time-frequency representations for underwater acoustics: Power class utilization with waveguide-invariant approximation. IEEE. Signal Proc. Mag, 2013; 30(6): 120−129 DOI: 10.1109/MSP.2013.2267651
[5] Tan T, Godin O A, Lefebvre A, et al. Characterizing the seabed by using noise interferometry and time warping. Proc. Mtgs. Acoust., 2018; 35(1): 070001 DOI: 10.1121/20000955
[6] Bonnel J, Caporale S, Thode A. Waveguide mode amplitude estimation using warping and phase compensation. J. Acoust. Soc. Am., 2017; 141(3): 2243−2255 DOI: 10.1121/1.4979057
[7] Bonnel J, Dosso S E, Chapman N R. Bayesian geoacoustic inversion of single hydrophone light bulb data using warping dispersion analysis. J. Acoust. Soc. Am., 2013; 134(1): 120−130 DOI: 10.1121/1.4809678
[8] Bonnel J, Gervaise C, Nicolas B, et al. Single-receiver geoacoustic inversion using modal reversal. J. Acoust. Soc. Am., 2012; 131(1): 119−128 DOI: 10.1121/1.3664083
[9] Bonnel J. Geoacoustic inversion in a dispersive waveguide using warping operators. J. Acoust. Soc. Am, 2011; 130(2): 101−107 DOI: 10.1121/1.3611395
[10] 李晓曼, 朴胜春, 张明辉, 等. 一种基于单水听器的浅海水下声源被动测距方法. 物理学报, 2017; 66(18): 114−127 DOI: 10.7498/aps.66.184301 [11] Bonnel J, Thode A. Using warping processing to range bowhead whale sounds from a single receiver. J. Acoust. Soc. Am., 2017; 19(1): 070066 DOI: 10.1121/1.4800509
[12] 姚美娟, 鹿力成, 马力, 等. 一种基于双引导声源和warping变换的拷贝声场计算方法. 声学学报, 2016; 41(1): 73−80 DOI: 10.15949/j.cnki.0371-0025.2016.01.009 [13] Baraniuk R G, Jones D L. Unitary equivalence: New twist on signal processing. IEEE Trans. Signal Process., 1995; 43(10): 2269−2282 DOI: 10.1109/78.469861
[14] Touzé G L, Nicolas B, Mars J I, et al. Matched representations and filters for guided waves. IEEE Trans. Signal Process., 2009; 57(5): 1783−1795 DOI: 10.1109/TSP.2009.2013907
[15] Bonnel J, Touzé G L, Nicolas B, et al. Automatic and passive whale localization in shallow water using gunshots. OCEANS 2008, Canada, IEEE, 2008
[16] Lopatka M, Touzé G L, Nicolas B, et al. Underwater broadband source localization based on modal filtering and features extraction. EURASIP J. Adv. Signal Process., 2010; 2010: 304103 DOI: 10.1155/2010/304103
[17] Niu H Q, Zhang R H, Li Z L. Theoretical analysis of warping operators for non-ideal shallow water waveguides. J. Acoust. Soc. Am., 2014; 136(1): 53−65 DOI: 10.1121/1.4883370
[18] Brown M G. Time-warping in underwater acoustic waveguides. J. Acoust. Soc. Am., 2020; 147(2): 898−910 DOI: 10.1121/10.0000693
[19] Qi Y B, Zhou S H, Ren Y, et al. Passive source range estimation with a single receiver in shallow water. Chinese Journal of Acoustics, 2015; 34(1): 1−14 DOI: 10.15949/j.cnki.0217-9776.2015.01.001
[20] 戚聿波, 周士弘, 张仁和, 等. 一种基于 β-warping 变换算子的被动声源距离估计方法. 物理学报, 2015; 64(7): 241−246 DOI: 10.7498/aps.64.074301 [21] 王冬, 郭良浩, 刘建军. 浅海负跃层条件下的warping变换宽带无源测距方法. 声学学报, 2019; 44(2): 145−154 DOI: 10.15949/j.cnki.0371-0025.2019.02.001 [22] 王冬, 郭良浩, 刘建军, 等. 一种基于warping变换的浅海脉冲声源被动测距方法. 物理学报, 2016; 65(10): 142−150 DOI: 10.7498/aps.65.104302 [23] 刘志韬, 郭良浩, 闫超. 利用自相关函数warping变换的浅海声源深度判别. 声学学报, 2019; 44(1): 28−38 DOI: 10.15949/j.cnki.0371-0025.2019.01.004 [24] 董阁, 郭良浩, 徐鹏, 等. 利用信号自相关函数warping变换的浅海水下移动观测平台机动优化方法. 声学学报, 2020; 45(6): 811−823 DOI: 10.15949/j.cnki.0371-0025.2020.06.003 [25] Gao W, Wang S D, Li X L. FPM-β: A method for waveguide invariant estimation using one-dimensional broadband acoustic intensity. JASA Express Lett., 2021; 1(8): 084802 DOI: 10.1121/10.0005842
[26] 戚聿波, 周士弘, 张仁和, 等. 水平变化浅海声波导中模态特征频率与声源距离被动估计. 物理学报, 2014; 63(4): 187−195 DOI: 10.7498/aps.63.044303 [27] 王冬. 水平阵被动测距及其应用研究. 博士学位论文, 北京: 中国科学院大学, 2018: 31−37 [28] 孟瑞洁, 周士弘, 李风华, 等. 浅海波导中低频声场干涉简正模态的判别. 物理学报, 2019; 68(13): 134304 DOI: 10.7498/aps.68.20190221 [29] 李家春, 周显初. 数学物理中的渐近方法. 北京: 科学出版社, 1998 [30] 周光, 罗元, 张毅, 等. 基于维纳滤波的自适应二值化技术的研究. 微计算机信息, 2008; 24(21): 231−233 DOI: 10.3969/j.issn.1008-0570.2008.21.094 [31] Gonzalez R, Woods R E. 数字图像处理. 第3版. 阮秋琦, 阮宇智译. 北京: 电子工业出版社, 2017: 404−406