引言

在勘探地震学中，由于近地表（0—500 m）分布的沙漠、沼泽、戈壁和风化层等松软结构的影响，地震波在传播过程中，出现明显的衰减，给地球物理成像带来了很大干扰。传统的全波形反演方法忽略了衰减对波形的影响，但通过品质因子Q则可以判断地下介质对波形的衰减能力，Q值越高，衰减能力越弱。如果已知具体的Q值分布，则可提高速度反演的精度以及地球物理解释的可靠性（Pramanik et al，2000 ；Wang，2002；Zhu et al，2014 ）。除此之外，Q值也可以作为地下天然气存储的指示依据（Odebeatu et al，2006 ）。

由波动方程可知，地震波在地球介质中传播时，速度和衰减是耦合的（Kjartansson，1979；Aki，Richards，1980；Zhu，Carcione，2014），因此，利用波形同时反演速度和衰减难度较大。目前已经发展了较多的利用地震数据反演同时得到速度和衰减结构的方法，例如，Stewart（1983）利用频率域上行波和下行波的波形同时反演出速度和地下介质品质因子Q值。Q值模型的反演需要以较准确的速度数据为基础，才能从多重外在因素所导致的衰减中独立出固有Q值的作用，故Kamei和Pratt （2008）率先提出利用波形首先反演出纵波速度结构，然后进行速度和Q值联合反演的方法。Smithyman等（2009）则利用声波波形反演了二维P波速度及其衰减。Zhu和Harris （2015）利用井间地震走时数据联合反演纵波速度、横波速度和Q值。

地球物理学方面的联合反演始于Vozoff和Jupp （1975），其利用直流电磁测深数据和大地电磁测深数据联合反演一维层状电阻率模型。自此之后，联合反演在地震数据领域得到了广泛的发展。杨文采和焦富光（1987）利用阻尼最小二乘法联合反射波走时和均方根速度反演出速度结构；张树林等（1993）利用井间地震和逆垂直地震剖面（vertical seismic profiling，简写为VSP）联合层析反演三维速度结构；Tryggvason等（2002）利用P波和S波走时信息联合反演纵横波速度，并用其进行地震定位；Gallardo和Meju （2003）提出利用交叉梯度联合反演的方法，并得到了广泛应用（Gallardo，Meju，2004，2007；Tryggvason，Linde，2006；Colombo et al，2008 ）。

通常我们会利用一维弹性波正演模拟数据反演一维速度结构，利用二维弹性波正演模拟数据反演二维结构，或者利用三维弹性波正演模拟数据反演三维结构（Zhou et al，1995 ；Pratt，Shipp，1999）。然而由于一维弹性波波形反演计算效率高，但由于近地表结构复杂，一维模型很难准确地描述实际的地质结构；二维或三维弹性波波形反演能够提供更加精确的地下结构，但计算耗时多，而且对计算机的存储有较高的要求，因此，如何利用一维弹性波模拟反演得到地下二维结构是值得探讨的。Auken和Christiansen （2004）提出伪二维反演电阻率的方法，即利用一维正演模拟，通过增加水平约束条件来反演二维结构。目前地球物理联合反演依然面临着不同数据权重选择的难题（Moorkamp et al，2007 ；Colombo et al，2016 ）。为了提高运算效率，降低计算所用存储空间，本文拟提出一种利用一维弹性波正演，通过二维Tikhonov正则化矩阵（Tikhonov，Arsenin，1977）对模型参数的约束，来反演二维结构的方法。该方法具体步骤为：首先，分析纵波速度、横波速度、Q_P和Q_S等4个参量联合反演的可行性，并设计应用交叉梯度联合反演的目标函数；其次，进行理论合成模型测试，讨论不同的权重因子对反演结果的影响；最后，对某地区的地震波数据进行处理，联合反演近地表的速度和Q值。

1.   弹性波正演

地震波是能量的传播，通常包括体波和面波两种。体波是指脉冲短且能够传到地下深部的波，包含P波和S波，遵循由地下不同区域不同的弹性模量和密度引起的反射、折射定律；面波通常指在地表传播的波，包括瑞雷波和勒夫波，其频率越低，能够到达的地下深度越大，而且面波的振幅随深度的增加衰减得很快。本文主要研究地震波中的P波和瑞雷波。

由于地球介质的黏弹性或内在摩擦，地震波能量在传播过程中逐渐衰减（Shearer，1999）。衰减的强度用1/Q来表示，即：

式中，E表示峰值应变能，−ΔE表示一个周期内能量的衰减，因此Q值越大，衰减越小。Knopoff （1964）认为高频地震波的Q值与频率有关，而低频地震波的Q值受频率影响则较小。

在地震波传播过程中，Q值通过

与速度结合，影响弹性波传播（Aki，Richards，1980）。式中，v（ω）表示频率为ω时的速度值，v₀表示参考频率为f₀时的速度值，故非线性波形正演计算可以简化为

式中：d为计算出的波形数据；m为模型参数向量，包含v_P，v_S，Q_P和Q_S；G表示正演计算函数。随后，利用有限差分公式分别求取弹性波对v_P，v_S，Q_P和Q_S的偏微分，并用波形将敏感度体现出来，其表达式为

式中，d（m_i，t）表示t时刻的地震记录，m_i表示第i个模型的参数，δm_i表示对m_i的扰动。用扰动后模型正演得出的波场减去扰动前模型正演得出的波场，再除以扰动量，便可得到波场对扰动模型的偏微分。

2.   弹性波反演

弹性波的初至波含有v_P和Q_P信息，面波含有v_S和Q_S的信息（Wang，Zhang，2017）。勘探地球物理中，初至波包含详细的地下结构信息，而面波对浅层结构信息更敏感。因此，综合初至波信息和面波信息，可以同时反演v_P，v_S，Q_P和Q_S。尽管目前的研究尚未发现Q_P与Q_S之间存在解析关系，但由于在地质学分析中，二者在结构上存在一定的联系，对于干燥或部分饱和的岩石，Q_P和Q_S的值的变化趋势相同（Winkler，Nur，1979），所以，我们考虑在联合反演过程中增加结构约束，使Q_P与Q_S的模型在结构上保持一致，并与速度模型相近。因此，所使用交叉梯度的目标函数为

This page contains the following errors:

Below is a rendering of the page up to the first error.

对式（6）添加扰动，可以得到

根据式（7），可得到反演过程中模型更新量的公式，即

式中，J_P和J_S分别为G_P和G_S的雅克比矩阵，可以通过式（4）计算出来。B为交叉梯度对v_P，Q_P，v_S和Q_S的偏微分矩阵，则

目前已发展了多种算法来反演式（8），例如模拟退火法、遗传算法、奇异值分解法、共轭梯度法等。利用全局最优化算法或半全局最优化算法，可以得到方程的最优解，但计算量大，所需计算机的存储空间大，在处理少量模型参数的情况下有优势，例如一维模型的反演；但在实际勘探地震学中，大量的波形数据以及模型参数需要解决，如果使用模拟退火法、遗传算法等会极大地增加计算时间（Beaty et al，2002 ）；而共轭梯度法在求解方程解时，仅需要动态保存雅克比矩阵和反演的模型参数，在处理勘探地震数据时更加有效实用（Hestenes，Stiefel，1952）。因此本文采用共轭梯度法求解式（8）得到模型的更新量。最后，通过迭代反演的模型正演计算的数据与观测到的数据求取差值，当误差达到预设的条件时，迭代停止；或者，当迭代次数达到预设迭代次数N时，反演结束。从目标函数的表达式式（5）可以看出，波形数据的权重以及交叉梯度的权重直接对反演产生作用，因此，我们将在后面的理论模型测试中详细讨论不同权重因子对反演结果的影响。

3.   理论模型测试

建立一个层状模型，如图1所示。v_P变化范围为1 900—2 200 m/s，v_S根据v_P/v_S＝2.0得到，故其变化范围为950—1 100 m/s，Q_P的取值范围为30—60，Q_S的取值范围为20—60。模型具有明显的分层：在180 m以下为均匀半空间；第二层在水平方向350—550 m之间和850—1 050 m之间设置一凸起的异常体和一凹陷的异常体。采用震源时间函数为雷克子波，中心频率为8 Hz，震源深度设置在10 m，接收器位于地表，偏移距为800 m，且记录了2.048 s的地震波序列，采样率为0.004 s。为了使联合反演过程中数据对速度和衰减的响应能够合理地体现出来，设置速度和衰减的参考量（v_P0，v_S0，Q_P0和Q_S0）来平衡波形的敏感度。通过测试得到一组比较合理的参考量为v_P0＝100 m/s；v_S0＝100 m/s；Q_P0＝1；Q_S0＝10。初至波和面波对v_P/v_P0，v_S/v_S0，Q_P/Q_P0以及Q_S/Q_S0的敏感度，如图2所示。

图  2  弹性波对参数v_P/v_P0 （a），v_S/v_S0 （b），Q_P/Q_P0 （c） 和Q_S/Q_S0 （d）的敏感度

Figure  2.  Sensitivity of elastic waveform to parameters v_P/v_P0 （a），v_S/v_S0 （b），Q_P/Q_P0 （c） and Q_S/Q_S0 （d）

下载: 全尺寸图片 幻灯片

图  1  真实理论测试模型

（a） v_P模型；（b） v_S模型；（c） Q_P模型；（d） Q_S模型

Figure  1.  True models for synthetic tests

（a） v_P model；（b） v_S model；（c） Q_P model；（d） Q_S model

下载: 全尺寸图片 幻灯片

建立如图3所示的联合反演初始模型。通常波形反演对初始速度模型的要求较高，因此，我们根据观测数据的初至波到时建立初始P波速度模型，然后根据纵横波速度比为2，得到初始S波速度模型。由于初始模型正演模拟得到的波形与真实波形之间的相位差不能过大，否则波形反演会陷入周期跳跃问题，不可得到准确的反演结果，又因波形对纵横波衰减的敏感度较小，联合反演对Q_P和Q_S初始模型的要求则较低，因此我们建立如图3c和图3d所示的各向均匀的衰减模型作为初始Q_P和Q_S模型。

图  3  联合反演的初始输入模型

（a） v_P模型；（b） v_S模型；（c） Q_P模型；（d） Q_S模型

Figure  3.  Initial models for joint inversion

（a） v_P model；（b） v_S model；（c） Q_P model；（d） Q_S model

下载: 全尺寸图片 幻灯片

随后测试不同权重因子对反演结果的影响。首先，确定初至波和面波的波形数据权重为1，v_P，v_S，Q_P和Q_S模型的光滑因子均为0.1，测试不同交叉梯度权重α_t对反演结果的影响。当α_t分别为0，0.1，1，3和10的情况下，波形不匹配度随反演迭代次数的变化如图4所示。可以看出，5组参数下的联合反演均收敛，而且随着权重因子的增大，交叉梯度对模型的约束从反演结果来看越来越强，在联合反演中的作用越来越大，致使数据的不匹配度可收敛至很小。

图  4  不同交叉梯度权重下初至波（a）和面波（b）的波形不匹配度随迭代次数的变化

Figure  4.  Misfit variation of initial arrival data （a） and surface wave data （b） with iteration number by different cross-gradient weightings

下载: 全尺寸图片 幻灯片

α_t＝0表示联合反演中不加入交叉梯度算子的限制，经过20次迭代反演得到的速度模型及其衰减模型如图5所示。横波速度结构的异常体反演清晰（图5a）；而纵波速度结构异常体的轮廓模糊，出现假构造，在垂直方向的分辨率很低，例如深度约为500 m时第二层的凸起异常区延伸至地表和第三层（图5b）；Q_P模型与真实模型的相差较大（图5c）；Q_S模型第二层的起伏变化明显，但第三层与第二层的分界模糊（图5d）。由于未引入交叉梯度算子，横波信息的反演结果比纵波信息的反演结果更接近真实模型，所以我们结合图4中波形不匹配度的变化认为，α_t＝10的反演结果要优于α_t＝0时的结果。对比α_t＝10时反演的速度及衰减模型（图6）与无交叉梯度的反演结果（图5）可以看到：v_P在图6a中第二层有起伏变化；v_S在图6b中的分层较图5b中更加清晰，第二层的起伏更加明显；Q_P和Q_S在图6c和d也显示出第二层具有明显的起伏变化且在深度为100 m处形成分界线。此后在进行真实波形与初始模型的波形对比以及真实波形与经过20次迭代后反演所得波形的对比（图7），可见初至波与面波拟合得较好。

图  5  交叉梯度权重α_t= 0时理论数据的联合反演结果

（a） v_P模型；（b） v_S模型；（c） Q_P模型；（d） Q_S模型

Figure  5.  The synthetic data inversion results with cross-gradient weighting α_t= 0

（a） v_P model；（b） v_S model；（c） Q_P model；（d） Q_S model

下载: 全尺寸图片 幻灯片

图  7  理论模型测试的初至波和面波的真实数据（黑线）与模拟数据（红线）的对比

图（a）中红线表示基于初始模型正演计算的数据，图（b）中红线表示基于交叉梯度为α_t=10反演模型正演计算的数据，下同

Figure  7.  The comparison of true data （black line） and calculated data （red line） for early arrivals and surface waves in the synthetic tests

The red line in Fig.（a） represents the forward calculation date associ-ated with the initial model，the red line in Fig.（b） represents the forward calaulation data resulted from inversion model with cross-gradient weighting α_t=10，the same below

下载: 全尺寸图片 幻灯片

图  6  交叉梯度权重α_t=10时理论数据的联合反演结果

（a） v_P模型；（b） v_S模型；（c） Q_P模型；（d） Q_S模型

Figure  6.  The synthetic data inversion results with cross-gradient weighting α_t=10

（a） v_P model；（b） v_S model；（c） Q_P model；（d） Q_S model

下载: 全尺寸图片 幻灯片

由前文得到的反演结果可知，交叉梯度因子能够为反演模型提供结构上的约束，使反演的4个模型在结构上趋于一致。联合反演的模型中，有关S波信息的横波速度和衰减结构与真实模型的拟合度较高，而有关P波信息的纵波速度和衰减结构与真实模型拟合度较低。因此，下文将简单讨论不同波形数据的权重因子对联合反演结果的影响。

确定交叉梯度的权重为10，S波的权重为1，分别设置ω_P为0.2，0.5，1.0和5.0等4种不同权重值进行联合反演，计算反演结果与真实模型之间的误差，结果如图8所示。由图可见：随着P波数据权重的增加，反演的P波速度结构与真实结构的差异逐渐变小，但当ω_P＝5.0时，差异在第5次迭代时最小，而随着迭代次数的增加而变大（图8a）；S波速度的反演结果随着ω_P的改变无明显变化，与真实结构之间的差异降至较低水平（图8b）；当ω_P＝1.0时，反演得到的Q_P模型与真实模型之间的差异达到最小（图8c）；经过20次迭代，ω_P＝0.5时的Q_S模型与真实模型之间差异最小，但ω_P＝1.0的结果在可接受的范围内（图8d）。综合对比v_P，v_S，Q_P和Q_S与真实模型之间的差异，可知ω_P＝1.0时的反演结果较好。

图  8  ω_P取不同值时反演的模型与真实模型之间的差异

（a） v_P模型；（b） v_S模型；（c） Q_P模型；（d） Q_S模型

Figure  8.  Differences between the true model and the inverted model with different ω_p

（a） v_P model；（b） v_S model；（c） Q_P model；（d） Q_S model

下载: 全尺寸图片 幻灯片

由图1—8的理论模型测试结果可知，联合初至波和面波反演速度与衰减参数的方法收敛速度较快，而且伪二维反演算法能够处理含有微弱起伏变化的近似水平层状结构。通过测试不同的交叉梯度权重对反演结果的影响（图4—8）可得：随着α_t的增加，反演的4个模型在结构上越来越相似，有效地弥补了波形数据对Q_P和Q_S敏感度低的缺陷。尽管得到的衰减模型较为光滑，且Q_P值和Q_S值与真实模型存在差异，但对于近地表地区，分辨率较低的Q_P和Q_S结构也足够为我们提供分析近地表的地质构造信息。通过对不同数据权重因子的测试（图8）表明：当减小P波数据的权重时，S波数据占主导地位，S波的相关参量v_S和Q_S，均能反演出较为准确的结构，而P波的相关参量v_P和Q_P的模型更新量受交叉梯度的限制作用变大，趋向于与v_S和Q_S结构上的一致，但是忽略P波数据的拟合度，容易导致反演的v_P和Q_P模型不合理；当增大P波数据的权重时，P波数据在反演中占主导地位，而S波数据的拟合被压制，因此，在图8a和8c中，ω_P＝0.5时分别在第5次和第7次迭代横型的误差中降至最小值。然而，由于初至波对v_P和Q_P的分辨率较低，反演的结构与真实结构存在差异，因此交叉梯度的作用引起v_S和Q_S结构随着不准确的v_P和Q_S结构变化，最终导致反演的模型拟合度不高。理论模型测试结果表明，根据反演的数据不匹配度曲线和模型的差异曲线选择最优权重因子，有利于得到合理可靠的反演结果。

4.   实际数据应用

将伪二维联合反演近地表的速度和衰减结构的方法应用到实际数据中。首先，根据到时文件拾取数据的初至到时，利用射线追踪走时反演方法得到一个初始纵波速度模型；其次，对地震波数据进行滤波、截取波形等预处理，滤波范围是2—15 Hz。本文测试了几种频段，包括2—10 Hz，2—15 Hz，2—20 Hz以及2—30 Hz，其中利用2—15 Hz频段滤波的波形反演得到的模型较合理；最后，选取偏移距为500 m、间隔为80 m的地震记录，炮点埋深为41 m，检波器置于地表，假定纵横波速度比为1.732，根据走时反演得到的纵波速度，得到一个横波速度模型。由于近地表浅层存在风化层，沉积层土壤松软，浅层的速度值很低，不完全遵循假定的纵横波速比，因此，我们根据面波到时，不断修改横波速度进行正演计算，直到正演模拟的面波到时与观测数据的面波到时基本吻合。由于地震波对Q_P和Q_S敏感度较低，我们建立均匀的Q_P和Q_S模型作为初始模型，如图9所示。

图  9  实际测试数据的初始输入模型

（a） v_P模型；（b） v_S模型；（c） Q_P模型；（d） Q_S模型

Figure  9.  Initial models for real data tests

（a） v_P model；（b） v_S model；（c） Q_P model；（d） Q_S model

下载: 全尺寸图片 幻灯片

联合反演中不同交叉梯度权重下，数据不匹配度随迭代次数的变化如图10所示。可见：随着α_t的增加，数据的不匹配度变大，但基本随着迭代次数的增加而降低；当α_t小于10时，数据的不匹配度基本保持在相同水平；当α_t为10时，数据的匹配度和交叉梯度的权重基本平衡了；当α_t为100时，数据的不匹配度陡然增高。α_t为0和10时联合反演的速度及其衰减模型，如图11和图12所示，对比可知，图12反演的Q_P和Q_S模型与速度在结构上有很大相似度。基于初始输入模型计算的波形与图12中真实波形的对比以及基于联合反演结果正演计算的波形与真实波形的对比如图13所示。结果显示，初至波拟合得较好，而面波由于实际采集的数据质量较差，最终拟合效果较初至波差。

图  11  交叉梯度权重α_t= 0时实际数据的联合反演结果

（a） v_P模型；（b） v_S模型；（c） Q_P模型；（d） Q_S模型

Figure  11.  The real data inversion results with cross-gradient weighting α_t= 0

（a） v_P model；（b） v_S model；（c） Q_P model；（d） Q_S model

下载: 全尺寸图片 幻灯片

图  12  交叉梯度权重α_t=10时实际数据的联合反演结果

（a） v_P模型；（b） v_S模型；（c） Q_P模型；（d） Q_S模型

Figure  12.  The real data inversion results with cross-gradient weighting α_t=10

（a） v_P model；（b） v_S model；（c） Q_P model；（d） Q_S model

下载: 全尺寸图片 幻灯片

图  13  实际数据测试的初至波和面波的真实数据（黑线）与基于初始模型（a）和反演模型（b）模拟数据（红线）的对比

Figure  13.  The comparison of true data （black line） and calculated data （red line） for early arrivals and surface waves associated with initial models （a） and inversion models （b）

下载: 全尺寸图片 幻灯片

图  10  不同交叉梯度权重下初至波（a）和面波（b）的不匹配度随迭代次数的变化

Figure  10.  Misfit variation of initial arrival data （a） and surface wave data （b） with iteration number by different cross-gradient weightings

下载: 全尺寸图片 幻灯片

研究区域地表风化严重，沉积层地质结构疏松，因此浅层介质对地震波吸收较多，尤其是地表至60 m深度的地层，Q值较小。由于面波的信噪比较低，因此，在反演过程中可利用交叉梯度的作用，得到与纵波速度和Q_P结构类似的横波速度和Q_S模型。通过对不同交叉梯度权重因子的讨论（图10），结果显示：如果交叉梯度权重过小，则4个参量之间的结构限制作用变小，联合反演结果与不引入交叉梯度的反演结果类似；如果交叉梯度的权重过大，则联合反演被交叉梯度因子主导从而忽略数据对反演的作用。因此，在联合反演中，要通过对比不同交叉梯度权重作用下的反演结果，综合分析波形的拟合度，以求得到最优的权重取值。

5.   讨论与结论

联合初至波和面波反演近地表的速度及其衰减，研究结果表明，Q值的影响在波形反演中不可被忽视，尤其是地质结构复杂的近地表区域，介质的Q值越小，对波形的衰减作用越大。利用弹性波的初至波和面波联合反演v_P，v_S，Q_P和Q_S等4个参量，不仅可以为地震解释提供丰富的结构信息，而且可以降低地球物理解释的非唯一性。本文在反演的目标函数中，引入交叉梯度算子，利用其对v_P，v_S，Q_P和Q_S这4个参量在结构上的约束，实现了弹性波全波形反演近地表的速度及其衰减结构。

在理论模型测试中，本文讨论了目标函数中的不同数据权重因子和交叉梯度权重因子对反演结果的影响。如果交叉梯度权重因子过大，反演的4个参数虽然在结构上具有较高的一致性，但是不能体现出地震波数据的影响，致使反演的结果可靠性降低；如果交叉梯度权重因子过小，反演的结果对波形数据的质量依赖性较强，而且波形对衰减的敏感度远低于对速度的敏感度，造成在反演迭代过程中只更新速度结构而衰减结构无明显的变化，反演结果也不可靠。因此，适当提高交叉梯度的权重，而又不忽略地震波数据在联合反演中的主导地位，导致联合反演得到的4个参量更接近真实模型。实际数据处理过程的难点在于对数据和交叉梯度权重因子的选择。由于研究区域近地表风化严重，浅层的衰减强，弹性波数据，尤其面波数据，质量较差，因此，在反演中适当降低面波数据的权重，同时增加交叉梯度的权重，使反演的横波速度和Q_S倾向于与纵波速度和Q_P的结构一致，有利于得到可靠的结果。实际操作中，需要通过大量的测试，依据反演结果与真实数据之间的拟合度，选择合理的权重值。

本文提出的伪二维反演是基于一维共中心点假设，利用二维正则化矩阵增加反演模型网格点之间的联系，实现二维结构反演。对于第4节中实际数据的处理，采用15个核心并行计算，反演迭代一次约需6分钟，20次迭代需2个小时，而真正的二维弹性波波形反演迭代一次的时间远大于1个小时。这种算法可以处理近似水平层状结构或者有微小起伏变化的结构。当地震波在一维介质中传播时，地表接收到的波形可以用来反演震源和接收器之间的中心点下方的结构，而如果速度或者Q模型在水平方向上有明显的变化，那么这种利用一维正演反演二维结构的伪二维算法可能会失败。然而，这种算法相比于一维反演结果能够提供更加丰富的地下结构水平方向的变化，相比于真正的二维弹性波反演，提高了计算效率，节省了计算时间。因此，这种算法在处理数据量大的勘探地球物理反演问题方面有很大的发展前景。