Warping变换是近年来广泛应用于水声领域的信号处理方法，它适用于浅海波导中宽带脉冲信号的简正波分离^[1]。Warping变换是一种域转换的计算方法，是一个酉变换。通过Warping变换可以将时域的宽带信号转换为Warping域中的单频信号，在对转换后的信号进行时频分析可使各阶简正波在时频平面上分离，再利用频域滤波分别提取各阶简正波的信息后进行Warping逆变换即可实现单阶信号的还原^[2]。Warping变换最初由Richard G. Baraniuk^[3]提出，与时频分析相结合并应用于信号处理。之后，Warping变换作为一种信号处理工具被广泛应用于水下被动测距、海洋环境参数反演等水声研究领域^[4-8]。

Warping变换公式的推导基于理想波导的简正波理论，是一个鲁棒性很强的变换^[9]。因此，Warping变换在水声领域的大部分应用都基于理想波导的Warping算子^[10-12]，这在一定程度上影响了Warping变换的处理效果。Warping变换在本质上是对接收信号的相位进行处理，非理想波导和理想波导之间在相位上的差别主要是非理想波导的相位考虑了来自于波导中海底参数的影响。在使用Warping变换时，为了消除或者减弱相位中来自海底参数的影响，Zhou等^[13]对接收信号进行了自相关处理，因为不同阶简正波之间的相互关系，自相关处理后，信号的相位只保留了与环境无关的频率不变量部分。Wang^[14]对接收时域信号的能量谱进行了Warping变换处理，接收信号能量谱分为相同阶简正波相干部分和不同阶简正波的相干部分，前者相位为零，而后者相位中特征频率在数值上等于理想波导中相干的两号简正波的截止频率差，只与海水中平均声速和海深有关, 与海底参数无关。以上两种方法虽然可以消除Warping变换在使用过程中来自海底参数等环境因素的影响，但却不能直接对接收信号进行处理，无法对接收信号的单阶简正波进行分离提取。因此，针对非理想波导，国内外学者对Warping变换公式进行了修正。G. Le Touzé^[15]引进了Pekeris波导下的Warping算子，实现了在时频域上对各阶模态的分离，该算子基于波导中关于群速度和相速度以及海水中声速的近似关系(v_pv_g=v₀²), 该算子虽然适用于理想波导，但应用于pekeris波导却没有相应的理论依据。Niu^[16]依据射线理论，利用海底反射系数和本征射线循环距离，推导出了适用于pekeris波导的更加精确的频散公式，对Warping变换公式进行了修正。此后，又针对不同海洋环境展开研究，推导出了适用于非理想波导的Warping变换的修正因子^[17]。但文献[17]中得到的修正公式需要已知详细的海底参数，包括海底声速、海底密度等，这在实际测量中难度较大，而海底相移参数P是一个包含海底参数信息却相对容易获得的参数^[18]。

针对海底参数获取难度大的问题，本文在文献[18]的基础上，利用海底相移参数P和海底环境参数之间的关系，提出了基于海底相移参数P的具有均匀半无限海底的非理想波导的Warping变换修正公式。

1 Warping变换的修正公式 1.1 Warping变换的基本理论

Warping变换适用于浅海波导且能够对单水听器的接收信号实现准确有效的分离^[10]。其基本理论如下^[9-10]。

根据简正波理论，海洋理想波导中水听器接收到的声压信号的时域表达式可以表示为

$ p\left( t \right) = \sum\limits_{n = 1}^N {{B_n}\left( t \right){{\rm{e}}^{{\rm{j}}2{\rm{ \mathsf{ π} }}{f_{cn}}\sqrt {{t^2} - t_r^2} }}} $

(1)

式中：$ {t_r} = {t_r} = {r_0}/c, {r_0}$为声源传播距离，c为波导中的平均声速，B_n(t)为第n阶简正波的瞬时幅度，f_cn是第n阶简正波的截止频率。

Warping算子h(t)表示为

$ h\left( t \right) = \sqrt {{t^2} + t_r^2} $

(2)

将Warping算子及式(1)代入Warping变换的计算公式，可得到Warping变换的结果：

$ {p_h}\left( t \right) = {W_h}p\left( t \right) = \sqrt {\left| {h'\left( t \right)} \right|} \sum\limits_{n = 1}^N {{B_n}\left( {h\left( t \right)} \right){{\rm{e}}^{{\rm{j}}2{\rm{ \mathsf{ π} }}{f_{cn}}}}} $

(3)

$\sqrt {|h\prime \left( t \right)|} $保证了变换前后信号的能量不变。Warping变换后的各阶简正波变成了以该阶简正波的截止频率为信号频率的单频信号，信号转换成了若干先后到达的单频信号的组合。

上述推导基于理想波导，虽然Warping变换具有很强的稳健性，对于大部分浅海波导均适用，但对于非理想波导，来自环境的影响会在一定程度上影响Warping变换的处理效果。因此，本文在已有的研究基础上利用海底相移参数P，给出了适用于具有均匀半无限海底的非理想波导的Warping变换修正公式。

1.2 Warping变换的修正

首先，将海水声速剖面的表达式为

$ c\left( z \right) = \bar c\left[ {1 - a\left( z \right)} \right] $

(4)

式中：$ \bar c = \frac{1}{H}\int_0^H {c\left( z \right){\rm{d}}z} $为平均声速，H为海深，a(z)为随深度变化的量，满足$ \int_0^D {a\left( z \right){\rm{d}}z = 0} $。此时，海水的波数为

$ k\left( z \right) = \omega /c\left( z \right) \approx \bar k\left[ {1 + a\left( z \right)} \right] $

(5)

其中$\bar k = \omega /\bar c $。

由波束位移的射线模型(beam displacement ray mode，BDRM)理论^[19]，循环距离的表达式为

$ {S_n} = 2\int_0^H {\frac{{{\mu _n}}}{{\sqrt {{k^2}\left( z \right) - \mu _n^2} }}{\rm{d}}z} $

(6)

将水平波数k(z)在$ {\bar k}$处进行泰勒展开，保留一阶项，并考虑到$\int_0^H {a\left( z \right){\rm{d}}z = 0} $，忽略高阶小项后，式(6)可化为

$ {S_n} \approx 2\int_0^H {\left[ {\frac{{{\mu _n}}}{{\sqrt {{{\tilde k}^2} - \mu _n^2} }} - \frac{{{\mu _n}\bar k\left( {k\left( z \right) - \bar k} \right)}}{{\sqrt {{{\left( {{{\bar k}^2} - \mu _n^2} \right)}^3}} }}} \right]{\rm{d}}z} $

(7)

同理，经过一个循环距离的时间可近似为

$ {T_n} = 2\int_0^H {\frac{{k\left( z \right)}}{{\sqrt {{k^2}\left( z \right) - \mu _n^2} }}{\rm{d}}z} \approx \frac{{2H\bar k}}{{\bar c\sqrt {{{\bar k}^2} - \mu _n^2} }} $

(8)

简正波群速的表达式可写为

$ v_n^g = \frac{{{S_n} + {\delta _{n1}} + {\delta _{n2}}}}{{{T_n} + {\tau _{n1}} + {\tau _{n1}}}} $

(9)

式中：δ_n1、δ_n2分别代表上、下界面处的波数位移，τ_n1、τ_n1为上、下界面处的时延，当波导为自由海面时，有δ_n1=0、τ_n1=0，且$ {\delta _{n2}} \ll {S_n}, {\tau _{n2}} \ll {T_n}$，所以群速度为

$ v_n^g = \frac{{{S_n} + {\delta _{n2}}}}{{{T_n} + {\tau _{n1}}}} \approx \frac{{{S_n}}}{{{T_n}}}\left[ {1 + \frac{{{\delta _{n2}}}}{{{S_n}}} - \frac{{{\tau _{n1}}}}{{{T_n}}}} \right] = \frac{{{S_n}\left( {1 + \varepsilon } \right)}}{{{T_n}}} $

(10)

式中：$\varepsilon = \frac{{{\delta _n}}}{{{S_n}}} - \frac{{{\tau _n}}}{{{T_n}}}, {\delta _n}、{\tau _n} $分别为相应的波束位移和时间延迟。群速度可表示为

$ v_n^g = \frac{{{r_0}}}{t} $

(11)

则将式(7)~(8)、(10)代入式(11)，可得到水平波数关于时间t的近似表达式：

$ {\mu _n} \approx \frac{{\bar k{r_0}}}{{\bar ct'}} $

(12)

其中t′=t(1+ε)。

根据射线理论，海面海底反射(surface-reflected bottom reflected, SBRB)类简正波的本征方程考虑小入射角的情况，在WKB近似下，该类波导的频散方程可以写为^[18]

$ 2\int_0^H {\sqrt {{k^2}\left( z \right) - u_n^2} {\rm{d}}z} + {\phi _b} - {\rm{ \mathsf{ π} }} = 2n{\rm{ \mathsf{ π} }} $

(13)

式中：u_n为第n阶简正波的水平波数，$ {\phi _b}$为海底反射相移。在$ \theta \ll 1$ (即小掠射角)的情况下，有如下近似：

$ {\phi _b} \approx - {\rm{ \mathsf{ π} }} + P\theta $

(14)

式中P为海底反射相移参数。将(14)式代入式(13)可得

$ \begin{array}{*{20}{c}} {2\int_0^H {\sqrt {{k^2}\left( z \right) - \xi _n^2} {\rm{d}}z} + P{\theta _n} = 2n{\rm{ \mathsf{ π} }}}\\ {n = 1,2,3, \cdots } \end{array} $

(15)

θ_n为第n阶简正波的掠射角，因为考虑的是小掠射角的情况，所以有θ_n≈sinθ_n，则式(15)可以化为

$ 2\int_0^H {\sqrt {{k^2}\left( z \right) - u_n^2} {\rm{d}}z} + P\sin \theta = 2n{\rm{ \mathsf{ π} }} $

(16)

引入有效深度的概念，则式(16)可以化为

$ 2\int_0^{H + \Delta H} {\sqrt {{k^2}\left( z \right) - u_n^2} {\rm{d}}z} = 2n{\rm{ \mathsf{ π} }} $

(17)

式中：ΔH为有效深度，在不考虑剪切波的情况下，已知有效深度和海底相移参数之间的关系如下

$ \Delta H = \frac{P}{{2k}} $

(18)

将式(17)的左侧按照式(5)进行泰勒级数展开并保留第一项后，可以得到为^[20]

$ \left( {H + \Delta H} \right)\sqrt {{{\bar k}^2} - u_n^2} = n{\rm{ \mathsf{ π} }} $

(19)

由式(12)可知

$ {u_n} \approx \frac{{\bar k{r_0}}}{{\bar ct'}} \approx \frac{{\bar k{r_0}}}{{\bar ct}} $

(20)

将式(18)、(20)代入到式(19)中，则有^[21-22]

$ \left( {H + \frac{P}{{2\bar k}}} \right)\bar k\sqrt {1 - \left( {\frac{{{r_0}}}{{\bar ct}}} \right)} = n{\rm{ \mathsf{ π} }} $

(21)

则第n阶简正波的瞬时频率表达式为

$ f\left( t \right) = \frac{{n\bar ct}}{{2H\xi \left( t \right)}} - \frac{{P\bar c}}{{4{\rm{ \mathsf{ π} }}H}} $

(22)

其中$\xi \left( t \right) = \sqrt {{t^2} - \frac{{r_0^2}}{{{{\bar c}^2}}}} $。

已知在t′≈t时，瞬时频率的表达式为^[6]

$ f\left( t \right) = \frac{{\bar c}}{{4H}}\left[ {\frac{{\left[ {\left( {2n - 1} \right){\rm{ \mathsf{ π} }} + \frac{{2{\phi _b}}}{{\rm{ \mathsf{ π} }}}} \right]t}}{{\xi \left( t \right)}}} \right] $

(23)

所以，式(22)等于式(23)，故：

$ \frac{{\bar c}}{{4H}}\left[ {\frac{{\left[ {\left( {2n - 1} \right){\rm{ \mathsf{ π} }} + \frac{2}{{\rm{ \mathsf{ π} }}}{\phi _b}} \right]t}}{{\xi \left( t \right)}}} \right] = \frac{{n\bar ct}}{{2H\xi \left( t \right)}} - \frac{{P\bar c}}{{4{\rm{ \mathsf{ π} }}H}} $

(24)

得到海底相移的表达式

$ \begin{array}{l} {\phi _b} = \frac{{\rm{ \mathsf{ π} }}}{2} - \frac{P}{{2t}}\xi \left( t \right)\\ {{\phi '}_b} = - \frac{{P \cdot {{\left( {\frac{{{r_0}}}{{\bar c}}} \right)}^2}}}{{2{t^2}\xi \left( t \right)}} \end{array} $

(25)

文献[17]对Warping变换的修正因子为

$ \gamma \left( t \right) = \frac{{\bar c}}{{4H}}\frac{2}{{\rm{ \mathsf{ π} }}}{{\phi '}_b}\xi \left( t \right) $

(26)

将式(25)代入式(26)得到：

$ \gamma \left( t \right) = - \frac{{\bar c}}{{4H}}\frac{2}{{\rm{ \mathsf{ π} }}} \cdot \frac{{P \cdot {{\left( {\frac{{{r_0}}}{{\bar c}}} \right)}^2}}}{{2{t^2}\xi \left( t \right)}}\xi \left( t \right) $

(27)

结合式(2)、(27)可以得到修正后的瞬时频率的表达式为

$ f\left( t \right) = \frac{{n{\rm{ \mathsf{ π} }}\bar ct}}{{2{\rm{ \mathsf{ π} }}H\xi \left( t \right)}} - \frac{{P\bar c}}{{4{\rm{ \mathsf{ π} }}H}} - \frac{{\bar c}}{{4H}}\frac{2}{{\rm{ \mathsf{ π} }}} \cdot \frac{{P \cdot {{\left( {\frac{{{r_0}}}{{\bar c}}} \right)}^2}}}{{2{t^2}\xi \left( t \right)}}\xi \left( t \right) $

(28)

利用式(28)求得瞬时相位为

$ {\psi _n}\left( t \right) = 2{\rm{ \mathsf{ π} }}\int\limits_{{t_r}}^t {f\left( t \right){\rm{d}}t} = 2{\rm{ \mathsf{ π} }}\left[ {{{\tilde f}_{cn}}\xi \left( t \right) + \tilde \chi \left( t \right)} \right] $

(29)

其中

$ \begin{array}{l} {{\tilde f}_{cn}} = \frac{{n\bar c}}{{2H}}\\ \tilde \chi \left( t \right) = - \frac{{P\bar c}}{{4{\rm{ \mathsf{ π} }}H}}\left( {t - {t_r}} \right) + \frac{{\bar cP \cdot t_r^2}}{{4H{\rm{ \mathsf{ π} }}}}\left( {\frac{1}{t} - \frac{1}{{{t_r}}}} \right) \end{array} $

(30)

则接收时域声压信号可以写为

$ p\left( t \right) = \sum\limits_{n = 1}^N {{B_n}\left( t \right){{\rm{e}}^{{\rm{j}}2{\rm{ \mathsf{ π} }}\left[ {{{\tilde f}_{cn}}\xi \left( t \right) + \tilde \chi \left( t \right)} \right]}}} $

(31)

式中：与简正波阶数有关的是第一项${{\tilde f}_{cn}}\xi \left( t \right) $，使得变换之后每号简正波的频率对应于单频频率$ {{\tilde f}_{cn}}$，由于第二项$\tilde \chi (t) $与简正波号数无关，故可以通过指数项的形式进行补偿。修正因子为

$ \left( {{M_q}x} \right)\left( t \right) = x\left( t \right){{\rm{e}}^{{\rm{j}}2{\rm{ \mathsf{ π} }}q\left[ t \right]}} $

(32)

其中$ q\left( t \right) = - \tilde \chi (t)$。

则修正后的非理想波导的Warping函数的形式：

$ \left( {{W_h}{M_q}p} \right)\left( t \right) = {\left| {h'\left( t \right)} \right|^{1/2}}p\left[ {h\left( t \right)} \right]{{\rm{e}}^{{\rm{j}}2{\rm{ \mathsf{ π} }}q\left[ {h\left( t \right)} \right]}} $

(33)

式(33)只适用于具有液态半无限空间的海底环境，同时未考虑海底吸收，研究包含海底吸收的修正公式需要更加详细的海洋环境参数，这将是接下来的研究工作之一。

2 修正公式准确性的仿真验证

为了说明修正公式的正确性，本文对修正公式利用计算软件kraken进行了理论仿真，仿真含海洋环境为pekeris波导和含跃层的水体环境。两种仿真环境下的收发距离r=10 km，取前五阶简正波进行处理。分别利用式(3)和式(33)对接收信号进行处理，pekeris波导下的修正结果

2.1 pekeris波导

pekeris波导海深H=100 m, 海水中声速c₀=1 500 m/s，海底声速为c_b=1 600 m/s，海底密度为ρ_b=1.5 g/cm³。仿真结果图 1所示。

2.2 含跃层波导

波导海深H=100 m, 海水中声速如图 2(a)所示，在海深20~60 m处存在声速跃层，海底衰减为α=0.02 dB/λ，其他环境参数与pekeris波导中一致。仿真结果如图 2所示。

由图 1和图 2的仿真结果可以看出，修正后的Warping变换和传统的Warping变换相比，各阶简正波在时频域上的分离更加明显。对接收信号进行Warping变换处理后，信号由宽频信号变为各阶简正波的单频信号，由图 1(e)和图 2(e)可以看出，经过修正后，各阶简正波所对应的特征频率为${{\tilde f}_{cl}} = \frac{{l\bar c}}{{2H}} $，与未修正的结果相比更准确。由图 1和图 2相比可知，式(33)的修正效果对pekeris波导更好。

2.3 修正因子中参量对修正结果的影响

由式(32)可知，修正因子中包含三个参量，即海深、收发距离以及海底相移参数P。在修正过程中，这三个参量的变换会影响修正结果，为了说明这三个参数对修正结果的影响，本文进行了如下验证：在同一海洋环境下的同一接收信号，分别改变上述三个参量，使其在一定范围内变化，观察经过式(33)处理后对应的频谱。

由图 3可以看出，修正因子中的三个参量在±10%范围内波动时，虽然会引起各阶简正波特征频率的变化，但对Warping变换后各阶简正波的分离效果影响不大，具体表现为各阶简正波频谱的左右移动。而三个参量的变换对特征频率的影响主要有以下几个方面：1)对收发距离r以及海深H而言，特征频率会随着他们值的增大而变小，相反，当海底相移参数变大时，特征频率的值则会变大。2)海深和相移参数的变化对特征频率的影响较小，而收发距离的变化对Warping变换后特征频率的影响则相对较大，所以修正因子对收发距离的变化较为敏感，因此，在利用式(33)对接收信号进行Warping变换修正时，如果对各阶简正波的特征频率的准确度要求较高，则需要使用较为准确的收发距离。

3 实验数据处理

为了验证本文提出方法在实际应用的有效性，对实测实验数据进行处理。实验数据为黄海某海域的海上试验获得的爆炸声信号。实验环境声速剖面为CTD测量获得，如图 4(a)所示，海深约为27.6 m，信号带宽取10~200 Hz，海底为液态半无限海底，海底声速为1 700 m/s，海底密度为ρ_b=1.5 dB/λ，利用式(3)和式(33)对接收信号进行了处理。

由实验结果可以看出，经过修正后，接收信号在时频域上有两阶明显的简正波信号，与未修正时相比(图 4 (c))，分离效果变好，由宽带信号变成对应其特征频率的两阶单频信号。同时，两种修正方法吻合较好，证明了两种修正方法的准确性。

4 结论

1) 在不考虑海底衰减的情况下，修正因子与海深H、海水中的平均声速c、收发距离r以及海底相移参数P四个参量有关。经过修正，Warping变换后的信号频谱表现为一系列特征频率为$\frac{{l\bar c}}{{2H}} $的单频信号，但对于存在海底衰减的海底，该结论依然具有一定的适用性。

2) 与已有的Warping变换修正公式相比，本文提出的修正公式不需要已知详细的海底参数，计算方便。

3) 仿真结果表明，修正因子中参量存在误差时会影响修正结果，特征频率会发生相应改变，但对Warping变换后各阶简正波的分离效果影响不大，具体表现为各阶简正波频谱的左右移动。通过仿真和实验数据对本文提出的方法进行了验证，修正结果良好，说明了本方法的有效性。研究海底具有衰减时的修正公式并利用修正结果对海洋环境参数进行反演将是接下来需要研究的重要内容。