波导不变量最先由学者Chuprov提出并用于描述浅海波导简正波频散特性^[1]，它定义为简正波相慢度差与群慢度差之比，可以描述声场干涉条纹斜率与声源距离、频率之间的关系，在反演、定位等方面有广泛应用，因此，估计波导不变量具有重要意义。目前常见的波导不变量估计方法主要分为三类。第一类是利用声场仿真软件计算简正波的相群速度，根据波导不变量定义，获取简正波相干项的波导不变量。然而该方法需要事先获取波导环境参数信息，以保证估计结果的准确性。第二类是利用人工投掷爆炸声源或其他类似的脉冲声源，对水听器接收信号声强进行频移补偿实现波导不变量估计。任云等^[2]利用一定水平间距的双水听器接收宽带脉冲声源信号，计算双水听器接收信号的空间相关系数，提取相关系数峰值对应的波导不变量作为估计结果。该类方法一般要求信号来自双水听器的端射方向附近。第三类是利用舰船辐射噪声源运动形成低频分析记录谱(low frequency analysis and recording, LOFAR)中条纹状干涉图样，通过提取条纹信息估计波导不变量。Rouseff等^[3]将距离-频率域干涉谱变换到波数-时延域，通过极坐标变换估计能应脊斜率实现波导不变量估计；Gall等^[4]对干涉谱进行距离域傅里叶变换得到相对频散曲线，其斜率对应波导不变量的倒数；Turgut等^[5]采用Hough变换参数映射方法，得到最符合声强干涉条纹特征的波导不变量值。然而，在一些场景下，机会声源(商船、货轮等)辐射噪声LOFAR谱中宽带连续谱干涉条纹非常微弱，上述估计方法应用效果欠佳。鉴于机会声源辐射噪声中线谱成分比连续谱成分平均高10~20 dB，相较连续谱而言能量高的线谱更易观测^[6]。因此，研究利用线谱信息实现波导不变量估计具有可行性。

与连续谱干涉条纹物理机理类似，线谱声强随距离起伏变化也是由简正波两两干涉引起的。与连续谱干涉条纹不同的是，线谱声强干涉起伏体现在距离维度。Harms和Young等^[7-8]发现在距离域上经过尺度伸缩变换的不同频率线谱声强干涉起伏之间存在正比关系，并且基于正比关系包含波导不变量这一事实，通过构建似然函数实现波导不变量估计；随后研究人员在此基础上开展了一系列线谱声源参数估计研究^[9-11]。然而，上述波导不变量估计方法需要将线谱附近(不包含线谱)的接收辐射噪声视作瑞利分布，用来估计似然函数未知参数——噪声方差；而线谱附近的辐射噪声是由服从对数正态分布的连续谱和服从瑞利分布的噪声共同构成，这会导致噪声方差估计结果存在误差进而影响波导不变量估计结果。此外，通常水听器记录声场以时间作为自变量，而在上述估计方法中线谱声强以距离作为自变量，需要通过声源径向速度将随时间变化的声强映射到随距离变化的声强，其前提是声源径向运动速度恒定且已知。而实际中，由于存在声源相对水听器的最近通过距离，声源径向运动速度会随时间发生变化。

针对已有研究的不足，本文提出一种利用线谱声强干涉起伏相关系数估计波导不变量方法。考虑声源沿水听器径向运动和非径向运动2种模型，对线谱声强干涉起伏进行重采样以及利用包含波导不变量的尺度伸缩变换关系，得到满足正比关系的线谱声强干涉起伏。在假设区间范围内对波导不变量进行搜索，计算参考线谱声强干涉起伏与尺度伸缩变换后的线谱声强干涉起伏之间的相关系数，相关系数峰值对应的波导不变量搜索值作为估计结果。由于线谱声强干涉起伏的准周期结构，相关系数变化曲线会出现旁瓣，通过多线谱声强处理可以有效降低旁瓣影响。所提方法直接对水听器记录的以时间为自变量的声场进行处理，无需声源径向运动速度恒定且已知。

1 浅海波导线谱声强干涉起伏 1.1 接收信号模型

以自身旋转部件的机械噪声和螺旋桨推进噪声为主的噪声源可以建模为K个低频线谱分量的叠加^[8]

(1)

式中: k表示线谱分量的号数; A_k和ω_k分别表示第k个线谱的幅度和角频率; φ_k表示线谱的相位。本文主要关注声源与水听器传播距离变化引起的幅度起伏, 暂不考虑声源线谱幅度起伏, 假设声源在线谱频率处辐射能量是稳定的。

根据简正波理论, 波导中位于(0, z_s)的声源到单水听器(r, z_r)的信道传递函数G(r, ω)可以表示为^[12]

(2)

式中: k_rm表示第m阶简正波的水平波数; ϕ_m表示第m阶简正波的模态深度函数; B_m表示声压的幅度项。水听器接收信号频谱表示为

(3)

式中, S(ω)为声源s(t)的频谱。

1.2 线谱声强干涉起伏

定义第k个线谱声强为^[12]

(4)

式中, Δk_rmn为第m阶与第n阶简正波的水平波数差, 即Δk_rmn=k_rm-k_rn。(4)式省略了随距离和频率慢变的非相干项, 保留了简正波两两干涉叠加的相干项, 该项决定了线谱声强随距离的起伏变化^[13]。

通常在多号简正波频散特性非常一致的浅海波导中, 波导不变量与相互干涉的简正波号数m, n无关。宽带连续谱辐射噪声的接收信号声强在距离-频率平面上的恒定强度条纹轨迹满足^[12]

(5)

式中, r₀和ω₀分别表示参考距离和参考频率, β表示波导不变量。根据(5)式和文献[7], 轨迹上任意两频率点ω₁和ω₂处声强I₁(r₁)和I₂(r₂)满足如下关系

(6)

式中: A₁和A₂分别表示频率ω₁和ω₂的声源频谱幅度; r₁和r₂为轨迹上两频率点声强对应的距离。

由于简正波干涉效应, 在一段距离区间内线谱声强呈现起伏变化, 把变化的声强称为线谱声强干涉起伏。(6)式表明, 参考频率ω₁的线谱声强干涉起伏与经过距离域尺度伸缩变换后频率ω₂的线谱声强干涉起伏满足正比关系, 且该正比关系包含波导不变量。这是本文估计波导不变量的基础。然而, 以上的正比关系是相对距离而言, 在大多数情况下, 水听器记录声场是以时间作为自变量, 此时正比关系并不直接满足。因此, 有必要研究时间域接收信号线谱声强干涉起伏估计波导不变量方法。

2 波导不变量估计

声源相对水听器的运动方式可以分为径向运动和非径向运动。不同运动方式下的时间域接收信号线谱声强干涉起伏不同。下面分别对这2种情况下的波导不变量估计方法进行研究。

2.1 径向运动情形

声源以绝对速度v₀沿水听器径向匀速运动, 任意时刻t的声源与水听器之间的距离表示为

(7)

式中, t_CPA表示最近通过时刻, 即声源运动到相对水听器的最近通过位置(the closest position of approach, CPA)处的时刻。对于浅海波导中单水听器接收的机会声源线谱而言, 可基于文献[14]事先获知最近通过时刻t_CPA。τ表示以t_CPA为参考的相对时间(下文称为相对时间)。

将(7)式代入(6)式, 由于频率ω₁和ω₂声强之间的正比关系与声源绝对速度v₀无关, 于是有

(8)

式中, τ₁和τ₂为轨迹上两频率点声强对应的相对时间。

为利用(8)式估计波导不变量, 引入相关系数概念, 并利用相关系数刻画参考频率线谱声强干涉起伏与经过时间域尺度伸缩变换后的频率线谱声强干涉起伏的相似程度, 实现波导不变量估计, 具体过程如下。

首先, 利用(5)式和(7)式, 可得尺度伸缩关系

(9)

给定一段时间区间[τ_L, τ_H]内的参考频率ω₁线谱声强干涉起伏I₁(下文称为参考线谱声强干涉起伏), 基于尺度伸缩关系对频率ω₂线谱声强干涉起伏I₂进行插值变换, 得到(下文称为经过尺度伸缩变换后的线谱声强干涉起伏), 即

(10)

式中，τ_L和τ_H分别表示起始时刻和结束时刻。定义I₁与之间的相关系数为

(11)

式中，I₁(τ)表示对时间区间[τ_L, τ_H]内的I₁(τ)进行中心化处理, 即I₁(τ)=I₁(τ)-E[I₁(τ)], E[I₁(τ)]表示I₁(τ)的均值。表示对进行中心化处理, β_s表示波导不变量搜索量。

由(11)式可知, 相关系数随着波导不变量的变化而变化, 当波导不变量是真值时, 相关系数达到最大值, 表明I₁与之间的相似度最高，满足正比关系。因此, 本文在[β₁, β₂] 区间范围内对不同波导不变量计算相关系数, 相关系数峰值对应的β便是参数估计结果。

2.2 非径向运动情形

图 1中CPA处与水听器之间存在一段距离, 称为最近通过距离, 记为r_CPA; 运动声源以v₀匀速驶向CPA处, 任意时刻t运动声源与水听器之间的距离为

(12)

显然, 距离和时间之间是非线性关系。若声强在时间域是均匀采样, 则在距离域是非均匀采样。

(12) 式两端同时除以v₀, 得到新的时间^[11]

(13)

将原始采样时间t代换为新的采样时间t′, 对线谱声强干涉起伏I_k(t)进行重采样, 得到重采样后的线谱声强干涉起伏

(14)

式中, 表示声源的最近通过距离与速度之比, 简称距速比。根据(6)式, 声强在距离域均匀采样是正比关系存在的基础, 由于t′与r之间满足正比关系并且I′_k(t′)是在新的时间域t′均匀采样, 所以条纹轨迹上任意两频率点ω₁和ω₂处声强I′₁(t′₁)和I′₂(t′₂)满足

(15)

式中, t′₁和t′₂为轨迹上两频率点声强对应的时间。

利用相关系数刻画在一段新的时间区间内参考频率重采样线谱声强干涉起伏与经过尺度伸缩变换后的重采样线谱声强干涉起伏之间的相似程度, 实现距速比和波导不变量联合估计, 具体过程如下。

首先, 利用(5)式和(12)式, 可得尺度伸缩关系

(16)

给定一段时间区间[t′_L, t′_H]内的参考频率ω₁重采样线谱声强干涉起伏I′₁(下文称为重采样参考线谱声强干涉起伏), 基于尺度伸缩关系对频率ω₂重采样线谱声强干涉起伏I′₂进行插值变换, 得到(下文称为经过尺度伸缩变换后的重采样线谱声强干涉起伏), 即

(17)

式中, t′_L和t′_H分别表示起始时刻和结束时刻。定义I′₁与之间的相关系数为

(18)

式中：b_s表示距速比搜索量；β_s表示波导不变量搜索量。

由(18)式可知, 相关系数随着距速比和波导不变量的变化而变化, 当距速比及波导不变量都是真值时, 相关系数达到最大值, 表明I′₁与之间的相似度最高, 满足正比关系。因此, 在区间[b₁, b₂] 和[β₁, β₂] 对不同距速比和波导不变量计算相关系数, 相关系数峰值对应的b, β便是参数估计结果。

3 仿真实验与分析

为利用数值仿真验证前述波导不变量估计方法的有效性, 假设波导环境为：水深为100 m, 水体声速为1 500 m/s, 半空间声速为1 750 m/s, 密度为1 750 kg/m³, 吸收系数为0.2 dB/λ, 其中λ为波长。声源深度为40 m, 水听器深度为70 m。仿真利用Kraken声场仿真软件。

3.1 径向运动情形

假设声源以2 m/s匀速运动, 时间区间为0~4 200 s, 运动声源在CPA处的时刻为3 800 s。信噪比(signal-to-noise ratio, SNR)定义为带宽内信号功率与高斯白噪声功率的比值。水听器接收信号信噪比为0 dB, 由于运动声源辐射噪声LOFAR谱中线谱成分较连续谱成分能量高, 将225, 250和275 Hz线谱声强提高15 dB作为已知线谱声强起伏, 如图 2所示。

可以看出, 简正波之间的相互干涉作用致使运动声源多频点声强随时间发生明暗起伏变化, 产生了类似干涉条纹在特定频率上调制的结果。3条线谱声强干涉起伏随相对时间的变化曲线, 如图 3所示, 图中黑色实线表示频率225 Hz参考线谱声强干涉起伏I₁与经过尺度伸缩变换后的275 Hz线谱声强干涉起伏的3对声强极大值的连线, 它们是波导不变量真值条件下尺度伸缩关系的示意。黑色虚线表示在波导不变量非真值条件下尺度伸缩关系的示意。时间区间[1 500, 2 300]内, 频率225与275 Hz对应I₁与的相关系数随波导不变量的变化曲线如图 4所示。

图 4表明, 一方面, 相关系数在波导不变量真值处达到峰值; 另一方面, 相关系数在某些波导不变量搜索值处存在较高的旁瓣。这主要是由于在波导不变量非真值条件下的尺度伸缩关系指向其他的声强极大值(如图 3黄色线上的红色虚线圈), 此时I₁与具有一定的相似性, 相关系数较高。当存在噪声干扰, 这些旁瓣可能高于真实峰, 成为伪解。因此, 需要对旁瓣进行抑制。

从图 3可以看出, 在波导不变量真值条件下, 黑色实线经过250 Hz线谱声强干涉起伏的极大值(如橙色线的红色实线圈), 而黑色虚线对应非极大值(如橙色线的红色虚线圈), 基于此, 通过分别求取多个频率的参考线谱声强干涉起伏和经过尺度伸缩变换后的线谱声强干涉起伏之间相关系数之和(下文称为多线谱声强干涉起伏相关系数), 可以实现相关系数旁瓣抑制效果。

频率225 Hz和250 Hz的参考线谱声强干涉起伏分别与经过尺度伸缩变换后的275 Hz线谱声强干涉起伏之间的相关系数之和的变化曲线如图 5所示。

从图 5可以看出，相关系数最高的旁瓣由图 4中的0.7降低到了0.4，得到有效抑制。相关系数峰值对应波导不变量为0.949，与之前波导不变量估计结果之间相差0.003。这主要由于不同频率线谱声强干涉起伏结构存在细小差别，致使波导不变量估计结果存在细微差异，认为它们都是真实波导不变量的反映。为不失一般性，给出各个时间区间的波导不变量估计结果，如图 6所示。

其中，时间区间长度保持一定，起始时间从1 500 s变化到1 800 s，步长为20 s。可以看出，各时间区间下的波导不变量估计结果基本一致，以它们的均值0.95作为估计结果。为验证波导不变量结果的准确性，将基于β=0.95得到的时间-频率关系曲线(图 2中红色实线)与干涉条纹进行对比。可见，它与干涉条纹保持一致，表明估计结果的正确性。

3.2 非径向运动情形

假设运动声源在CPA处的时刻为3 800 s，该时刻声源相对水听器的最近距离为800 m。水听器接收信号的LOFAR谱如图 7所示。给出时间区间2 000~3 800 s内频率225, 250和275 Hz线谱声强干涉起伏如图 8a)所示。通过(14)式对原始采样时间的线谱声强进行重采样，声强随新的采样时间的干涉起伏如图 8b)所示。黑色实线表示3条重采样线谱声强干涉起伏的极大值的对应关系，即波导不变量为真值时的尺度伸缩关系。与之前类似，依据此关系通过尺度伸缩变换得到的多线谱声强干涉起伏之间满足正比关系，此时相关系数最大。当距速比和波导不变量的搜索值偏离真实值时，正比关系不满足，相关系数减小。

设距速比的搜索区间为[200, 600]，搜索间隔为1；波导不变量的搜索区间为[0.8, 1.2]，搜索间隔为0.001。根据(14)式和(17)式，利用时间区间[2 700, 3 800]内重采样后的频率225 Hz和250 Hz参考线谱声强干涉起伏，分别与尺度伸缩变换后的重采样275 Hz线谱声强干涉起伏之间的相关系数如图 9所示。

从图 9可以看出，相关系数主瓣集中分布在真实距速比和波导不变量附近，其中相关系数峰值对应波导不变量和距速比分别为0.948和397，相较3.1节估计结果0.95仅存在微小偏差。

图 10给出了不同时间区间的波导不变量和距速比估计结果。其中，时间区间从[2 300, 3 500]到[2 600, 3 800]，步长为20 s。蓝色星号表示各时间区间的参数估计结果。红色方框表示通过K-means聚类算法求得各区间参数估计结果分布的质心，对应距速比和波导不变量估计值分别为398.8和0.947，后者与3.1节估计结果0.95几乎一致，表明所提方法能有效实现波导不变量估计。

4 实测数据验证

数据选自SWellEx-96水声试验^[15]，由于本文假设声源在线谱频率处辐射能量是稳定的，故采用试验S5单频信号数据对波导不变量估计方法有效性进行验证。垂直线列阵(vertical line array, VLA)位置处水深约为213 m，海底有23.5 m的泥沙层，覆盖在800 m厚的岩石层上，声速剖面选自第5个站点CTD数据，如图 11所示。声源深度为9 m，接收器为VLA的第10号水听器(以下简称为水听器)，深度大约为150 m。

声源沿着200 m等水深线以2.5 m/s速度自南向北运动，行进时长75 min。根据安装在声源和VLA的GPS系统，声源-接收器之间的实际距离随时间变化如图 12所示。声源相对水听器的最近通过时间约为3 540 s，最近通过距离约为900 m。图 12中的模拟距离表示运动声源与CPA位置之间的距离随时间的变化。蓝色括号表示0~30 min时间区间，此区间内声源相较水听器距离较远，模拟距离与实际距离基本一致，最近通过距离的影响可以忽略，可以将声源看作沿着水听器径向运动；绿色括号表示45~59 min时间区间，声源随着时间增加逐渐接近水听器，最近通过距离的影响不可忽略，模拟距离与实际距离之间的差异明显。下面将使用上述2段时间区间内声强干涉起伏对第3节中2种运动模型估计波导不变量方法的有效性进行验证。

水听器采样频率是1 500 Hz，约以1.36 s为间隔(采样点数为2 048)，重叠率为50%。图 13给出了水听器记录时间区间0~75 min内接收信号LOFAR谱。由于声源-水听器相对运动产生的多普勒效应，选取声源相邻频率窗中最大声强作为该时刻声强值。

利用15~27 min时间区间内频率127和145 Hz参考线谱声强干涉起伏，通过(10)式得到尺度伸缩变换后的频率163 Hz线谱声强干涉起伏，多线谱声强干涉起伏相关系数随波导不变量搜索值的变化曲线如图 14所示。

考虑不同时间区间下的波导不变量估计结果。时间区间长度保持不变，起始时间从10~20 min，步长为0.4 min。不同时间区间的多线谱声强干涉起伏相关系数峰值对应波导不变量估计结果如图 15所示。可以看出，不同时间区间波导不变量估计结果基本稳定在1.11左右，相较图 6估计结果存在起伏，可能是由于试验海域地形不平坦、海深变化对波导不变量产生影响。

时间区间45~57 min内的频率127和145 Hz重采样参考线谱声强干涉起伏与尺度伸缩变换后的频率163 Hz重采样线谱声强干涉起伏之间的相关数随波导不变量变化曲线不变量估计结果强干涉起伏计算得到相关系数系数如图 16所示。其中，波导不变量估计结果为1.126，与径向运动模型下的估计结果1.11基本保持一致。保持时间区间长度不变，起始时间从45~47 min，步长为0.2 min，各时间区间的波导不变量和距速比估计结果，如图 17所示。蓝色星号表示各时间区间的参数估计结果。红色方框表示各区间参数估计结果分布的质心。对应距速比和波导不变量估计值分别为312.3和1.15。表明所提方法能有效实现波导不变量估计。

5 结论

本文研究了低信噪比下利用线谱声源提取波导不变量的问题，根据参考线谱声强干涉起伏与尺度伸缩变换后的线谱声强干涉起伏存在的正比关系包含波导不变量这一事实，提出了一种利用时间域线谱声强干涉起伏相关系数估计波导不变量方法。通过理论分析、仿真实验和实测数据处理得到了如下结论：①对于声源相对水听器径向和非径向2种运动模型，通过对原始采样时间的线谱声强干涉起伏在新的采样时间下进行重采样，参考线谱声强干涉起伏与经过尺度伸缩变化的线谱声强干涉起伏满足正比关系。无需声源相对水听器径向运动，也不要求已知声源运动速度；②相关系数旁瓣由线谱声强干涉起伏准周期结构导致，通过采取多线谱处理一定程度上降低了相关系数旁瓣；③相比宽带连续谱，机会声源线谱声强干涉起伏可以在更低信噪比下实现波导不变量估计。