0 引言
海洋声学研究表明，海水声速是温度、盐度和压强(深度)等物理化学参数的复杂函数，还受洋流和海洋生物等因素影响，具有三维和随机时变的特点。然而，无论在技术上还是在工作量上，利用观测手段研究海水声速的三维随机变化都具有很大的挑战性，只是在近年来才被提到有关研究日程上来。另一方面，利用理论手段研究海水声速的三维时空变化更具有挑战性，基本上没有研究结果问世。因此，对海水声速的观测和理论研究一般都是停留在给定的研究点上，所得到的研究结果是海水声速随深度的变化，称为海水声速剖面[1-2]。
数学上，海水声速剖面主要是以经验或回归公式的形式给出，只有著名物理海洋学家Walter H. Munk[3]1974年发表的公式是同时利用海水声速沿深度方向的梯度与海水沿深度方向的温度、盐度和压力梯度之间的经验关系(回归公式)和Brunt-Väisälä频率导出的。
此外，与经验公式以盐度、温度以及压强或深度为变量不同，Munk公式给出海水声速随深度的变化。鉴于在海洋反射地震观测中一般并不采集海水温度和盐度等海洋环境参数，所以只有Munk公式可以被直接用于解决海洋反射地震中与水体速度有关的问题，例如海水速度分布对地震波传播和偏移成像的影响研究等[4]。
在温度、盐度和压强等参数中，温度对海水声速剖面的影响最大[1, 3]。研究表明，海水中的温度分布具有明显的垂直分层性，最典型的是三层结构，即由上至下为混合层、过渡层和恒温层结构。其中：混合层下边界的深度位于100~200 m，其温度主要受太阳辐射影响，基本不变；过渡层位于混合层之下，其温度一般具有阶跃性变化，也称为温跃层；恒温层位于过渡层之下，其温度保持不变。
Munk公式是海水温度三层结构中海水声速的典型数学表达，适合于世界上的大部分海域[1, 3]。
在浅海和中深海中，声速一般随深度呈负梯度递减。在深海中，声速将首先随深度下降，直至一个极小值；然后再以一定梯度上升，直至海底，形成声波波导，称为深海声波波导、深海声道或水下声道。类似地，在浅海区中也会出现声波波导现象，称为浅海声道或表面声道[1-3]。Munk公式用指数函数描述深海声波波导的形态，既没有考虑混合层中声波速度随深度的变化，也没有考虑深海水体中随机水团对深海水体声速的影响。因此，Munk公式所给出的深海水体声速剖面是一种在理想状态下的声速剖面。
在常规海洋地震勘探数据处理中，一般将海水声速取为常数，即1 500或1 450 m/s[2]。在浅水区(0~300 m)，由于海水声速变化幅度很小，所以做上述选取不会引入明显的误差。然而，在深水区采用这样简单的选取将会引起明显的误差。事实上，深海声道的存在对地震波的射线路径、射线走时和射线振幅均会产生不同程度的影响，进而使偏移成像结果发生变化[4-6]。此外，海水声速随时间的变化会使海上反射地震和时移地震的处理结果偏离正确的结果，所以在实际的资料处理工作中必须考虑海水速度时空变化对观测数据的影响[7]。
近十几年来，海水速度变化对深水地震观测数据和偏移成像的影响逐渐得到学术界和工业界的重视，涌现了一批研究结果[4-8]。就深海水体速度建模来讲，一般是通过速度分析和层析来实现的[9-10]。在公开的文献中，没有利用Munk公式进行深海水体速度建模的报道。然而，对水层速度进行速度分析的前提条件是存在明显的水层反射[9]，而利用层析对海水层进行速度建模一般假设速度在深度方向上保持不变[10]。
本文的目的是建立利用Munk公式和海底反射走时反演深海水体速度的基本理论。与常规速度反演不同，深海水体声速剖面的大致趋势是已知的，而未知的只是速度沿深度变化的精确形态。具体到Munk声速剖面，速度随深度的变化规律完全由其所含的4个参数(简称为Munk参数)决定。因此，在利用Munk公式进行水体速度反演时只要对其所含的4个参数进行反演即可，没必要对整个水体进行离散并在每个深度点上进行反演。因此，利用Munk公式实现深海水体速度反演将会大幅度减少深海水体速度建模的计算量，这对实际应用具有重要的意义。
作为水下地层的盖层，水体的速度分布会直接影响水下地层中的地震照明度和偏移像场的位置，也会极大地影响复杂海域速度建模的精度和效率。由于海底反射的运动学及动力学特点仅与海水层和海底的物理性质有关，所以在利用海底反射进行速度反演时可以不考虑水下地层的几何形态和物理性质，也可以对海底的声学边界条件进行简化处理，从而降低反演的难度和工作量。如果利用水下地层反射进行水体速度反演，则必须对水层速度和地层速度同时进行调整和更新。这显然要比只利用海底反射复杂得多。况且，海底反射和海底位置是比较容易识别和提取的信息，不会存在无法克服的困难。
在下文中，将首先对Munk公式进行变分分析，以得到Munk参数，即扰动参数、声道轴处的速度、声道轴深度和波导宽度等对海水声速影响程度的显式表达。其次，将在海底深度为已知的条件下分别建立面向海底反射走时的射线走时反演和波动方程走时反演理论。具体地讲，是建立理论走时和观测走时之差的Fréchet导数和海森(Hessian)矩阵公式，以及建立对Munk参数的更新公式。最后，给出本文的结论并对本文所涉及的假设条件和一些与实际应用有关的问题进行简要讨论。
1 Munk公式及其变分本节的目的是建立Munk声速剖面相对于其4个参数的变分和变化率公式。这些公式构成了利用Munk公式反演深海水体速度理论的基础，具有重要的意义。
根据定义，一个物理量的变分(用δ表示)与其增量(用Δ表示)等价。因此，为了求出速度增量，既可以采用变分法，也可以采用微扰法和微分法。然而，如果在不同的方法中采用不同的近似，则所得到的结果会不尽相同。
数学上，Munk公式的外观形式[3]为
(1) 式中：v(z)为深海水体速度；z为海水深度；v1为声道轴处的速度；ε为扰动参数，ε=BγA/2(B为声道宽度，一般可设其等于声道轴深度z1；γA为绝热速度梯度，γA=0.0114/km[3])；η为复合参数，η=2(z-z1)/B，物理上，η代表相对于z1的归一化距离，无量纲。
不难看出，Munk公式给出的深海水体速度由两部分组成：一部分为常数，即声道轴处的速度v1；另一部分为校正项。如果将v1、ε和η视为参数，则校正项是v1和ε的双线性函数，是η的非线性函数。此外，校正项是深度的非线性函数。
1.1 变分法对式(1)作相对于ε、v1和η的变分，得到:
(2) 将e-η展开为关于η=0的Taylor级数并取其前三项，可将(2)式简化为
(3) 根据定义，η=2(z-z1)/B，因此，
(4) 将式(4)代入式(3)，得到v(z)关于v1、ε、η和B的变分为
(5) 根据式(5)：如果只有v1或ε发生变化，则δv(z)呈线性变化；如果v1和ε同时发生变化，则δv(z)呈双线性变化；如果只有z1发生变化，则δv(z)呈非线性变化。同样的结论对于波导宽度B也成立。
根据式(5)，可得到速度相对于Munk参数的变化率公式为：
(6) 
(7) 
(8) 
(9) 为得到上述变分在模型参数空间中的一点α0=(ε=ε0, v1=v1, 0, z1=z1, 0, B=B0)的值，只需要在式(6)—(9)中用η0(η0=2(z-z1, 0)/B0)替换η即可。
如果B=z1，则
(5') 此时，式(6)和式(7)保持不变，而式(8)和式(9)合并为δv/δz1，即
(10) 式中：
(11) 
(12) 至于v(z)相对于v1、ε和η的2阶及高阶变分，可以利用式(5)式进行计算，此处不再进行进一步的讨论。
1.2 微扰法设B=z1、ε=ε0+Δε、v1=v1, 0+Δv1以及z1=z1, 0+Δz1，构建速度差函数：
式中：v=v(z; ε0+Δε, z1, 0+Δz1, v1, 0+Δv1)；v0=v(z; ε0, z1, 0, v1, 0)。
对指数函数e-η进行关于η=η0的Taylor级数展开并舍去高次项，得到
(13) 由此得出：
(14) 
(15) 
(16) 将式(14)—(16)中的指数函数展开为Taylor级数并取到3次项，得到
(14') 
(15') 
(16') 利用式(11)和式(12)可将式(14')—(16')表示为z、ε0、z1, 0和v1, 0的函数。
对于2阶及高阶变分，可以采用类似的方法进行计算，具体结果在此从略。
1.3 Taylor级数展开法令Δv=v(ε, v1, z1)-v(ε0, v1, 0, z1, 0)。如果对式(1)在ε=ε0、v1=v1, 0以及z1=z1, 0处进行Taylor展开并仅取其线性项，则得到
(17) 式中：α0=(ε0, v1, 0, z1, 0)；Δε=ε-ε0；Δz1=z1-z1, 0；Δv1=v1-v1, 0。
(18) 
(19) 
(20) 或者
(18') 
(19') 
(20') 对于2阶及高阶变分，可以利用式(13)进行计算，具体结果在此从略。
1.4 3种方法比较表 1汇总了3种方法在3阶近似下的结果。可以看出，3种方法所得到的结果完全相同。值得指出的是，如果在变分法中不先把函数e-η展开成级数，则可得到与式(14)—(16)或式(18)—(20)完全相同的结果。
射线走时反演利用射线理论反演深海水体的声速剖面。为此，首先要计算出海底反射的射线走时相对于速度的变分，然后再利用复合函数求导的链式法则建立射线走时关于Munk参数的Féchet导数、目标函数和海森矩阵，最后给出Munk参数的梯度和共轭梯度法更新公式。如果海底是水平的，则可以得到免射线追踪的反演公式。反之，必须在已完成射线追踪的基础上讨论问题。考虑到射线追踪是个纯粹的计算问题且已经有很多研究结果，所以在此不对其进行讨论。
2.1 射线走时及其变分首先考虑崎岖或倾斜海底条件下的一般公式。令x=(x, z)为地下任意点的坐标且u(x)=1/v(x)为慢度。根据文献[11-12]，射线走时公式为
(21) 式中：σ为沿射线的参数，dσ=ds/u=vds；ds为射线弧长微元；σ0为σ在射线起点处的值。令u0(x)和δu(x)分别为背景慢度和一远远小于u(x)的慢度扰动，则u(x)=u0(x)+δu(x)。如果射线端点保持不变，则根据Fermat原理，射线路径不受速度扰动的影响。此时，
(22) 为了求出慢度扰动δu(x)与速度扰动δv(x)的关系，考虑
(23) 式中，v0(x)是背景速度且|δv(x)|≪v0(x)。根据定义，有
(24) 将式(24)代入到式(22)，得到
(25) 这是δt关于δv的1阶近似。可以证明，δt关于δv的2阶近似为
(26) 公式(25)和(26)给出了走时扰动与速度扰动之间的普遍关系。如果在这两个公式中用Δv代替δv，可得到在δv为有限时的走时时差Δt=t-t0(t0为速度为v0时的射线走时)。
如果速度只随深度变化，即v=v(z)，则可得到不显式依赖于射线路径的走时公式。事实上，若令θ为所考虑射线的出射角(射线在源点处与垂直轴的夹角)，p=(sin θ)/v为射线参数(注意：p沿给定射线保持不变)，则对于一个给定的p(对于一条给定的射线或对于一个给定的θ)，透射射线的走时t(p)和射线终点距源点的水平距离X(p)分别为[11-12]：
(27) 
(28) 式中：zP为射线终点的深度，其最大值为射线转折点的深度，由公式c(zP)=1/p确定。
对于水平海底，与给定p相对应的反射走时和射线在海面上的出射点距源点的水平距离分别为2t(p)和2X(p)。对于倾斜和复杂海底，要根据反射定律对入射射线和反射射线分别进行计算。如果存在射线转折点，则同样需要对t(p)和X(p)进行分段计算。
现在考虑t(p)和X(p)相对于速度的变分。利用直接变分法，有
(29) 
(30) 令
在经过简单的计算后，得到时差函数Δt(p)关于Δv的2阶近似为
(31) 对于X(p)，令ΔX(p)=X(p)-X0(p)(X0(p)是当速度为v0时射线终点相对于射线始点的水平距离)。将v=v0+Δv代入到式(28)，得到X0(p)关于Δv的2阶近似为
(32) 公式(29)给出了射线走时相对于速度的变分。利用链式法则δt(m)=(δvt)(δmv)(m为Munk参数矢量，δv和δm分别表示对v和对m的变分)，可将式(29)变换为射线走时关于m=(ε, v1, z1, B)T=(m1, m2, m3, m4)T的变分。
2.2 Fréchet导数为求出射线走时关于速度的Fréchet导数，引入走时积分算子KL和KZ以及其1阶和2阶Fréchet导数(KL)'、(KL)″、(KZ)'以及(KZ)″。其中，
式中：L表示沿射线路径积分；Z表示沿垂直坐标积分。此外，KL只作用在慢度u(x)上；KZ只作用在式(27)中的被积函数上。
类似地，引入水平距离积分算子KX，以及其1阶和2阶Fréchet导数(KX)'和(KX)″。其中，
KX只作用在式(28)中的被积函数上。注意：虽然KX和KZ在形式上完全一样，但是其作用的被积函数不同，因此有必要分开定义。
根据定义和公式(26)，在用Δt替换δt以及用Δv替换δv后得到：
(33) 
(34) 类似地，利用公式(29)和(30)，得到:
(35) 
(36) 
(37) 
(38) 为得到射线走时相对于Munk参数矢量m=(ε, v1, z1, B)T的Fréchet导数，可将公式(33)和(35)分别代入链式法则(∂t(m)/∂m)=(∂t(m)/∂v)(∂v/∂m)。由于m为矢量，所以在利用链式法则后所得到的是矢量而不是标量。
2.3 目标函数和海森矩阵在Munk公式(式(1))中，令B=z1，且令m=(m1, m2, m3)=(ε, v1, z1)T为模型空间中的Munk参数矢量。定义反射走时的理论值(计算值)tcal与观测值tobs之间的差为Δt，Δt=tcal-tobs。为求出深海水体的速度剖面，计算函数
(39) 的极小值。式中：Δd=(Δt1, Δt2, …, ΔtN)T是定义在数据空间中的列向量；‖·‖为Hilbert空间中的L2范数；N为参加反演的地震道道数。
现在求J(m)的梯度。对于∇J(m)的第j个分量，有
令
(40) 则
(41) 式中，yT=FTΔd为列向量，矩阵F称为Fréchet导数矩阵，其矩阵元素为理论走时tcali的Fréchet导数。
现在考虑J(m)的2阶导数。根据定义，
写成矩阵形式，有
(42) 式中，H=FTF+(∂FT/∂m)Δd为海森矩阵。一般来讲，为了避免求∂FT/∂m，可取H≈FTF。
将公式(39)中的J(m)在m=m0处展开并令Δm=m-m0，得到
(43) 式中，yT和H均取其在m=m0处的值。
在实际计算中，利用公式(40)和(42)计算F和H中的矩阵元。
2.4 最优解与梯度类更新公式根据定义，目标函数J(m)的极小值给出射线走时反演的最优解Δm。为了求出这个最优解，对式(41)求梯度并令其为0，得到
因此，
(44) 如果利用梯度法实现式(44)的迭代求解，则可利用Pica等[13]文章中的公式(A-9)和公式(A-10)。此时，m的更新公式为
(45) 式中：αk为更新步长；φk为搜索方向，通过公式J(mk+1)=[tcal(mk)-tobs+αkFkφk]T定义。一般来讲，可以取φk为最速下降方向γk，即取∇J(m)的负方向。在实际反演中，也可利用文献[13]中的公式(A-9)，即取mk+1=mk+χφk(χ为一个趋近于零的小参数，一般取max{χ·φk}≤max{mk}/100)。
对于更新步长αk，根据文献[13]中的公式(A-10)，有
(46) 式中，Fk是在当前模型中算出的Fréchet导数矩阵。
如果利用共轭梯度法对m进行更新，则其搜索方向pk和mk的更新公式[14]为：
(47) 
(48) 式中：∇J由公式(41)给出；αk和βk按下列公式进行计算[14]：
(49) 
(50) 波动方程走时反演(wave equation traveltime inversion, WTI)利用波动方程计算已知速度模型中的点源辐射场(理论波场)，并利用互相关提取计算走时和观测走时之间的走时差[15]。一旦得到了时差，就可以利用目标函数
进行最优化计算。式中：mWTI是模型空间中的矢量，由地下离散网格节点上的速度定义；ΔdWTI是时差，与公式(39)中Δd的定义一样，唯一不同的是ΔdWTI是利用数值模拟场与观测场之间的互相关求得的，而公式(39)中的Δd是利用射线走时计算的；s表示对源点求和(对接收点的求和已隐含在式(39)中)。
本节的核心内容是对文献[15]中的有关公式进行变换，以得到数值模拟场相对于Munk参数的Fréchet导数。
3.1 基本公式在波动方程走时反演中，利用数值模拟数据和观测数据之间的互相关函数f(xr, τ; xs)计算目标函数JWTI(mWTI)中的ΔdWTI，其定义[15]为
(51) 式中：xs为源点；xr为接收点；τ是观测场和数值模拟场在时间轴上的相对位移；p(xr, t+τ; xs)obs是经过相对时移τ的观测场；p(xr, t; xs)cal是数值模拟场；A(xr; xs)obs是观测场p(xr, t; xs)obs中最大的振幅，起归一化的作用。f(xr, τ; xs)极大值的位置给出p(xr, t+τ; xs)obs和p(xr, t+τ; xs)cal之间的时差Δτ，即
(52) 式中：T是p(xr, t; xs)obs和p(xr, t; xs)cal之间时差的最大值。因此，Δτ必须满足下列方程:
(53) 式中，
在Luo等的WTI论文[15]中，用梯度法对速度进行更新，即vk+1=vk+αkξk。式中：αk为步长；ξk是目标函数JWTI(v)沿其最陡下降方向的值。
(54) 
(55) 式中：求和号下的s和r分别表示对源点和接收点求和；*表示关于时间的褶积；g(x, t; x', t')是标量波动方程的格林函数；因子E由下列积分给出。
(56) 在利用Munk速度公式的波动方程走时反演中，需要把文献[15]中的速度更新公式vk+1=vk+αkξk变换为对Munk参数矢量m=(ε, v1, z1)T的更新公式。为此，要首先求出Δτ关于m的Fréchet导数。利用链式法则和公式(55)有
(57) 式中，∂v/∂m可以根据第1节中的有关结果进行计算。
利用式(57)，可将公式(40)写为下列形式：
(58) 将公式(58)代入到第2节的有关公式中，即可得到针对m的梯度和共轭梯度更新公式。
3.2 数值模拟波场计算为得到数值模拟波场p(xr, t; xs)cal，理论上可以采用任何一种数值模拟方法。然而，海底反射是来自于水下地层顶部的反射，只与海底和水体的声速和密度差有关，与水下地层的物理性质无关。另外，从上文可知，在WTI中，对数值模拟场的走时要求很高，而对其振幅的精度要求不高。因此，可以采用几何射线理论等快速数值模拟方法进行数值模拟场计算[11]。如果海底形态复杂，则可用高斯波束法取代几何射线法[11]。为了能够准确地模拟海底散射，可以采用Kirchhoff法或散射积分法，例如Born近似或广义Born近似等[16-19]。除了这些在地球物理中常用的方法之外，还可采用计算海洋学中的一些常用方法，例如简正模式展开法等[1, 20]。
4 走时反演的问题与对策前面三节分别建立了利用Munk公式和海底反射走时构建深海水体声速剖面的一般性理论，包括Fréchet导数计算、射线走时反演和波动方程走时反演，没有涉及如何利用海洋声学中的有关研究结果制定走时反演策略和走时反演方案等问题，也没有讨论在利用Munk公式和海底反射走时进行深海水体速度建模时将会遇到的问题和挑战，以及解决问题和迎接挑战的方法与对策。本节的目的是要处理这些在前文中没有涉及和没有讨论的问题。为此目的，将利用海洋声学中的有关研究结果，在前面三节所建立的一般性理论基础之上，讨论利用Munk公式和海底反射走时构建深海水体声速剖面的具体算法和策略。考虑到篇幅的限制和问题的复杂性，计划专门成文讨论与数值实现有关的一些具体问题。因此，在本节中只论及策略而不讨论算例。
对于射线走时反演，两个具有挑战性的问题是海底定位和海底反射走时拾取。本文假设海底深度和形状为已知，主要是根据我国近海海域有比较详细的水深资料这一事实。如果在勘探海域没有公开的水深资料，则可以考虑同时反演水体速度和海底，包括海底的深度和形状。对于这一方案，将在后续的研究中给予考虑。
对于海底反射走时拾取，理论上可以通过人工拾取、自动拾取以及智能拾取来实现。考虑到在文献中已经有数种自动拾取反射走时的方法，例如互相关法和瞬时信号法等[15, 21-23]，所以在此不再对海底反射走时拾取问题进行进一步的讨论。
射线走时反演中存在的另一个问题是反射点的精确定位问题。由于假设海底位置和形状为已知，所以可以采用由下而上的透射射线族追踪方案，即从海底向海面进行波前构建或有限差分走时计算。对于前者，可以同时获得射线路径和射线走时。对于后者，必须在得到走时后利用反向梯度法计算射线路径。如果采用波前构建法计算，则在得到射线路径后必须对射线进行挑选，即从射线族中选出满足反射定律的反射射线对。在射线分布密度足够大的条件下，可对波前构建的结果进行射线插值，以使透射射线的终点与源点或接收点重合。或者，可采用反向梯度插值法直接从源点和接收点进行射线插值。对于基于有限差分走时的反向梯度法，可以直接从所考虑的源点和接收点进行反向梯度计算，直接得到反射射线[16-17]。
此外，还可以利用射线路径和射线走时的精确或高精度近似公式代替常规的射线追踪和有限差分计算所给出的数值结果。之所以这样考虑，是因为深海水体在深度上的不均匀性相对于水下地层的速度不均匀性来讲很弱，所能引起的走时差并不大，有可能会淹没在数值计算所引起的误差里。考虑到篇幅的限制，将专门成文讨论如何建立海底反射走时的精确或高精度解析近似公式，故而在此不再赘述。
因此，利用Munk公式和海底反射走时进行深海水体速度反演的第一个策略就是采用精确或高精度近似公式计算射线路径和射线走时，以减小数值计算误差对走时反演的影响。或者，采用由下至上的射线追踪方案，即将点源置于海底，计算由海底向海面的射线路径和射线走时，并利用插值技术确定海面处满足Snell定律的炮-检对位置。
与射线走时反演对计算机资源要求不高不同，波动方程走时反演对计算机资源有很高的要求。因此，为了能走向实用，必须要设法提高数值模拟的实现速度和降低有关算法对内存和输入输出的要求。鉴于本文所引用的波动方程走时反演理论对振幅的精度要求不高，所以可选用成熟的非网格快速正演方法，例如射线法、高斯波束法、WKBJ(Wentzel-Kramers-Brillouin-Jeffreys)法、Kirchhoff法、屏近似法和Born近似法等近似方法。在具有复杂海底的海域，为避免在计算Green函数时出现焦散，可利用高斯波束叠加替代几何射线Green函数以及采用广义Born近似等。
波动方程走时反演的另一个问题是海底边界条件的选取。由于波动方程走时反演对振幅信息的要求不高，所以可以将水下地层视为一个具有一定厚度的均匀层并采用海洋声学常用的边界条件(即绝对软边界条件、绝对硬边界条件、混合边界条件以及密度间断边界条件，外加辐射条件[24-25])，而不采用常用的、考虑实际地层物理模型的衔接条件。这样做的好处是可以不考虑水下地层的物理性质，只对海洋水体进行离散(如果采用网格法实现的话)。
因此，利用Munk公式和海底反射走时进行深海水体速度反演的第二个策略就是选用简单、高效的数值模拟方法和最简单的海底边界条件，把波动方程数值模拟的计算区域限定为水体而不用考虑模型的水下部分。
在本文中，将v1、ε、z1、B这4个参数统称为Munk参数。由于一般可以取B与z1近似相等，所以实际上Munk公式中只有3个参数，即v1、ε和z1。对于世界上的绝大多数海域，声道轴的深度和变化规律是已知的(纬度越高，声道轴深度越浅，越靠近海面。例如：我国南海深水区的声道轴深度在1 000 m附近；大西洋中部的声道轴深度位于1 100~1 400 m；地中海、黑海和日本海以及温带太平洋中，声道轴深度位于100~300 m[24-25])。这暗示着，真正需要反演的Munk参数只有2个。因此，在利用了海洋声学中的有关研究结果之后，原来的四参数反演问题即可转化为双参数反演问题。
进一步，对v1和ε的反演可以采用交替方式实现。例如，可先固定v1而反演ε，然后再利用所得到的ε反演v1。反之亦然。
因此，利用Munk公式和海底反射走时进行深海水体速度反演的第三个策略就是要利用已知的深海声道深度数据和交替反演策略来减少反演参数的数量，先从4个变成2个(从四参数反演问题变为双参数反演问题)，然后再从2个变成1个(从双参数反演问题化为单参数反演问题)。
近年来，深海水体声速随方位和时间的变化逐渐得到关注[24-25]，这对应着反射地震中的三维走时反演问题。由于Munk公式描述的是海水声速随深度的变化，所以在前文中所建立的2个基于海底反射走时的速度反演理论均为一维理论，没有考虑水体速度在横向和方位上的变化。对于这个缺陷，可在实际应用中采用逐点反演的方式加以克服。事实上，根据海洋声学的有关研究，深海水体声速在横向上的变化是很缓慢的。因此，在一个CDP(共深度点)数据体(道集)内，水体速度在横向上一般不会出现显著的变化。因此，如果利用本文所提出的理论进行一个CDP接一个CDP的逐点反演，即可实现深海水体的三维速度建模。
因此，利用Munk公式和海底反射走时进行深海水体速度反演的第四个策略就是由点到面，利用一维声速剖面公式实现深海水体的三维速度建模。对于海水声速随时间的变化，可以利用不同时间采集的海底反射数据进行反演，从而不构成实质性的问题和困难。
对于实现射线走时反演理论的最优化算法，前文只讨论了梯度类方法。这主要是因为所选为基础之一的波动方程走时反演方法是利用梯度类方法实现的。理论上，可以选用任意一种最优化方法计算目标函数(式(39))的极小点。然而，考虑到在实际实施过程中只需要讨论单参数反演问题，所以应该尽量降低实现算法的复杂程度。
因此，利用Munk公式和海底反射走时进行深海水体速度反演的第五个策略就是利用最简单的算法实现走时反演中的最优化计算。
5 结论与展望本文建立了利用Munk公式和海底反射走时构建深海水体速度剖面的基本理论，包括射线走时反演理论和波动方程走时反演理论。对于射线走时反演，主要是建立射线走时关于Munk参数(扰动参数ε、声道轴速度v1、声道轴深度z1和声道宽度B)的Fréchet导数和海森矩阵；对于波动方程走时反演，主要是对文献[15]中的有关公式进行变换，以使其能从对速度的更新变换为对Munk参数的更新。由于Munk速度剖面随深度的变化关系完全由4个Munk参数控制，所以在利用本文所建立的理论进行深海水体速度建模时不必对水体进行离散，而且其计算量也与水深无关。事实上无论水体的厚度是多少(无论水深是多少)，本文所提理论的反演参数最多为4个，从而可以大大地提高反演的计算速度。这与常规的速度反演完全不同。
在计算射线走时的Fréchet导数时，分别采用了直接变分法、微扰法和Taylor级数展开法。结果表明，这3种方法不仅在理论上完全等价，而且在结果上也完全相同。
为实现利用Munk公式和海底走时对海水声速剖面的反演，提出了5个策略：①利用精确或高精度近似公式计算射线路径和射线走时，或者采用由下至上的射线追踪方案，精确计算海面处满足Snell定律的炮-检对位置；②选用简单、高效的数值模拟方法和最简单的海底边界条件，把波动方程数值模拟的计算区域限定为水体；③利用已知深海声道深度和交替反演策略将四参数反演问题简化为单参数反演问题；④由点到面，利用一维声速剖面公式实现深海水体的三维速度建模；⑤用最简单的算法实现走时反演中的最优化计算。
与常规水下地层速度反演已具有数十年的研究和应用历史不同，深海水体速度建模在近十年内才逐步得到学术界和工业界的重视。因此，作为利用Munk公式进行水体建模的第一步，肯定无法覆盖射线走时和波动方程走时反演的全部理论问题。换句话说，在将来肯定还必须要进行很多理论和计算方面的研究工作才能使利用Munk公式进行深海水体速度建模的理论走向实用。
| [1] |
Jensen F B, Kuperman W A, Poter M B, et al. Computational Ocean Acoustics[M]. New York: Springer, 2011.
|
| [2] |
Jones E J W. Marine Geophysics[M]. London: University College London, 1999.
|
| [3] |
Munk W H. Sound Channel in an Exponentially Stratified Ocean, with Application to SOFAR[J]. Journal of Acoustical Society of America, 1974, 55(2): 220-226. |
| [4] |
韩复兴, 孙建国, 王坤. 深海声道对波场传播的影响[J]. 石油地球物理勘探, 2014, 49(3): 444-450. Han Fuxing, Sun Jianguo, Wang Kun. Deep Sea Channel Influence on Wave Field Energy Propagation[J]. Oil Geophysical Prospecting, 2014, 49(3): 444-450. |
| [5] |
Han F, Sun J, Wang K. The Influence of Sea Water Velocity Variation on Seismic Traveltimes, Ray Paths, and Amplitude[J]. Applied Geophysics, 2012, 9(3): 319-325. DOI:10.1007/s11770-012-0344-2 |
| [6] |
姬莉莉, 林缅. 中尺度涡影响下目的层地震成像的畸变与影响因素分析[J]. 地球物理学报, 2013, 56(1): 195-203. Ji Lili, Lin Mian. The Distortion of Seismic Imaging in the Presence of Mesoscale Eddies and Analysis of Influencing Factors[J]. Chinese Journal of Geophysics, 2013, 56(1): 195-203. |
| [7] |
Bertrand A, MacBeth C. Seawater Velocity Variations and Real-Time Reservoir Monitoring[J]. The Leading Edge, 2003, 22(4): 351-355. DOI:10.1190/1.1572089 |
| [8] |
韩复兴, 孙建国, 王坤, 等. 海水速度差异对中深层偏移成像的影响分析[J]. 地球物理学报, 2015, 58(9): 3439-3447. Han Fuxing, Sun Jianguo, Wang Kun, et al. Analysis of the Influence of Sea Velocity Difference on Migration Imaging in Middle and Deep Layers[J]. Chinese Journal of Geophysics, 2015, 58(9): 3439-3447. |
| [9] |
鲁统祥, 张文祥, 邓勇, 等. 针对海洋水层反射的处理方法研究[J]. 地球物理学进展, 2018, 33(3): 1247-1254. Lu Tongxiang, Zhang Wenxiang, Deng Yong, et al. Method for Processing the Seismic Reflection of Sea Water[J]. Progress in Geophysics, 2018, 33(3): 1247-1254. |
| [10] |
Ritter G L da S. Water Velocity Estimation Using Inversion Methods[J]. Geophysics, 2010, 75(1): U1-U8. DOI:10.1190/1.3280232 |
| [11] |
Červen' V. Seismic Ray Theory[M]. Cambridge: Cambridge University Press, 2001.
|
| [12] |
Goldin S V. Seismic Traveltime Inversion[M]. Tulsa: SEG, 1989.
|
| [13] |
Pica A, Diet J P, Tarantola A. Nonlinear Inversion of Seismic Reflection Data in a Laterally Invariant Medium[J]. Geophysics, 1990, 55(3): 284-292. DOI:10.1190/1.1442836 |
| [14] |
张贤达. 矩阵分析与应用[M]. 北京: 清华大学出版社, 2004. Zhang Xianda. Matrix Analysis and Applications[M]. Beijing: Qinghua University Press, 2004. |
| [15] |
Luo Y, Schuster G. Wave-Equation Traveltime Inversion[J]. Geophysics, 1991, 56(5): 645-653. DOI:10.1190/1.1443081 |
| [16] |
孙建国, 苗贺. 基于Chebyshev走时逼近的三维多次反射射线计算[J]. 吉林大学学报(地球科学版), 2018, 48(3): 890-899. Sun Jianguo, Miao He. Computation of Three Dimensional Multi-Reflection Rays Based on Traveltimes Numerical Approximation Using Chebyshev Polynomials[J]. Journal of Jilin University (Earth Science Edition), 2018, 48(3): 890-899. |
| [17] |
孙建国, 李懿龙, 孙章庆, 等. 基于模型参数化的地震波走时与射线路径计算[J]. 吉林大学学报(地球科学版), 2018, 48(2): 343-349. Sun Jianguo, Li Yilong, Sun Zhangqing, et al. Computation of Seismic Traveltimes and Ray Path Based on Model Parameterization[J]. Journal of Jilin University (Earth Science Edition), 2018, 48(2): 343-349. |
| [18] |
孙建国. 高频渐近散射理论及其在地球物理场数值模拟与反演成像中的应用:研究历史与研究现状概述以及若干新进展[J]. 吉林大学学报(地球科学版), 2016, 46(4): 1231-1259. Sun Jianguo. High-Frequency Asymptotic Scattering Theories and Their Applications in Numerical Modeling and Imaging of Geophysical Fields:An Overview of the Research History and the State-of-the-Art, and Some New Developments[J]. Journal of Jilin University (Earth Science Edition), 2016, 46(4): 1231-1259. |
| [19] |
孙建国. 声波散射数值模拟的两种新方案[J]. 吉林大学学报(地球科学版), 2006, 36(5): 863-868. Sun Jianguo. Two New Schemes for Numerical Modeling of Acoustic Scattering[J]. Journal of Jilin University (Geoscience Edition), 2006, 36(5): 863-868. |
| [20] |
刘巍, 王勇献, 张理论. 数值海洋声学[M]. 北京: 科学出版社, 2019. Liu Wei, Wang Yongxian, Zhang Lilun. Numerical Ocean Acoustics[M]. Beijing: Science Press, 2019. |
| [21] |
Sun B, Alkhalifah T. Adaptive Traveltime Inversion[J]. Geophysics, 2019, 84(4): U13-U29. DOI:10.1190/geo2018-0595.1 |
| [22] |
Luo Y, Ma Y, Wu Y, et al. Full-Traveltime Inversion[J]. Geophysics, 2016, 81(5): R261-R274. DOI:10.1190/geo2015-0353.1 |
| [23] |
Saragiotis C, Alkhalifah T, Fomel S. Automatic Traveltime Picking Using Instantaneous Traveltime[J]. Geophysics, 2013, 78(2): T53-T58. DOI:10.1190/geo2012-0026.1 |
| [24] |
刘伯胜, 黄益旺, 陈文剑, 等. 水声学原理[M]. 3版. 北京: 科学出版社, 2019. Liu Bosheng, Huang Yiwang, Chen Wenjian, et al. Principles of Water Acoustics[M]. 3rd Ed. Beijing: Science Press, 2019. |
| [25] |
汪德昭, 尚尔昌. 水声学[M]. 2版. 北京: 科学出版社, 2017. Wang Dezhao, Shang Erchang. Water Acoustics[M]. 2nd Ed. Beijing: Science Press, 2017. |