地震走时和航空电磁联合反演三维速度和电阻率
王月 张捷
摘要：地震走时数据可以反演近地表的速度，但不能反演出隐蔽层和低速层。航空电磁数据可以反演近地表的高电阻率和低电阻率，但是对垂直方向的分辨率低。联合三维地震走时数据和航空电磁数据反演了近地表的速度和电阻率结构，并采用棋盘模型测试和实际数据应用展示了联合反演出的速度模型优于单独的走时反演出的速度模型。结果表明：该算法可以应用到处理数据量大的勘探地震成像中，能提供优化的速度结构。
关键词：联合反演；航空电磁；地震走时；近地表
中图分类号：P313.23，P315.2 文献标识码：A 文章编号：1000-0666（2018）01-0022-10
0 引言
地球物理反演是指利用接收到的信号推测地下介质的物理性质及结构分布，不同信号能够反演出不同的物理参数信息，为地球物理解释提供证据，有助于寻找矿产、油藏等资源。近地表地质结构复杂，分布着沼泽、沙漠、戈壁、风化层等松软结构，而且存在山地等起伏地形，给地下介质反演带来很大困难。因此，利用有限的数据资料反演出近地表的地质结构，有利于地球物理深部成像。
联合反演在地球物理学的发展始于由Vozoff和Jupp（1975）提出的利用直流电磁测深数据和大地电磁测深数据联合反演一维层状电阻率模型，采用的反演算法为迭代二阶马奎特阻尼最小二乘法。随后，联合不同数据的反演逐渐发展，Gomez-Tre-vino和Edwards（1983）联合可控源电磁数据和直流电测深数据利用奇异值分解法（SVD）反演了地层的电阻率和厚度，Roth和Holliger（1998）利用高分辨率的瑞雷波和导波联合反演的P波速度和S波速度，Bosch和McGaughey（2001）在先验条件下利用重力数据和电磁数据联合反演二维密度结构和磁导率分布，Tryggvason等（2002）利用P波走时和S波走时信息联合反演纵横波速度和地震定位，Pieta和Bala（2008）联合直流电测深数据和电磁数据应用蒙特卡罗全局优化算法实现了并行计算反演视电阻率等。
地震数据和非地震数据之间的联合反演，也得到广泛发展。周辉等（1994）利用广义线性反演方法及非线性反演的预条件最速下降法开展了一维地震一大地电磁测深资料反演方法研究，得出顺序反演效果优于地震、电磁单独反演效果，而同时反演的效果最优。王西文（1997）采用剥离法对重力、地震资料联合反演目的层密度值，进而预测油气藏。Hering等（1995）利用一维直流电测深数据和地震面波数据联合反演得到近地表的电阻率模型和面波速度模型。Misiek等（1997）利用电法数据和勒夫波数据联合反演了Borsod地区近地表的速度和电阻率结构。Meju和Gallardo（2003）验证了近地表地区电阻率和纵波速度存在结构一致的现象，提出了交叉梯度算法，实现了二维直流电测深数据和地震走时数据的联合反演（Gallardo，Meju，2003，2004），以及二维大地电磁和地震走时数据的联合反演（Gallardo，Meju，2007）。交叉梯度在联合反演中得到极大地推广（Linde et al，2006；Fregoso，Gallardo，2009；Ogunbo，Zhang，2014）。
本文利用地震走时数据和航空电磁数据，应用交叉梯度对模型的结构限制，联合反演了近地表的速度和电阻率模型。地震走时正演是三维模拟，航空电磁正演是一维模拟，通过三维Tik-honov正则化（Tikhonov，Arsenin，1977），联合反演出三维速度模型和三维电阻率模型，再对合成棋盘模型进行了測试，并应用到一组实际数据中，检验反演结果。
1 研究方法
1.1 地震走时正演模拟
地震波走时正演模拟的目的是计算波从源到接收器的时间以及射线路径，主要有基于射线的方法和基于网格的方法。基于射线的方法简单高效，能解决各向同性和各向异性介质，模拟的到时相对准确，比如打靶法（Anderson，Kak，1982；Wil-lianMSon，1990）和弯曲法（Julian，Gubbins，1977；Graeber，Asch，1999）。基于网格的方法种类繁多，应用原理各不相同，可以计算模型中每个网格点的到时，能够稳定处理各向同性介质，并且能同时高效计算出到时和射线路径，如求解程函方程的有限差分法（Vidale，1988，1990）和应用费马原理的最短路径法（Nakanishi，Yamaguchi，1986；Moser，1991）。
求解程函方程和最短路径的方法可以提供初至波到时，而打靶法和弯曲法只能得到一个不能确定是初至波的走时，因此，本文采用了求解程函方程的方法正演模拟地震走时和射线路径。射线追踪的程函方程为（Vidale，1988，1990）：式中：s为慢度；t为走时。此方法是采用有限差分向外扩展的方式，利用已知点的到时寻找到达下一个网格点的用时最短的路径。同时，可以根据最终走时的梯度从接收点回溯到源，得到射线路径。这种方法得到的地震走时准确而且计算效率高，计算速度与定义的网格点个数成正比。由于这种立方体向外扩展计算波前的方法不总是遵循波在介质中传播的方向，存在波前不连续的问题，比如，初至波可能在传回立方体之前需要对立方体外的区域再进行扩展计算。但是随着后续的发展，Vidale提出的方法在地震波走时计算中得到广泛的应用（VanTrier，Symes，1991；Qin et al，1992；Hole，Zelt，1995）。
1.2 航空电磁正演模拟
航空电磁包括频率域系统和时间域系统，时间域电磁场通常可提供比频率域地球更深部的信息（Holladay et al，1997），而频率域电磁场能得到高分辨率的浅层信息，因此，在近地表研究中，本文采用频率域电磁场数据反演浅层的电阻率结构。
时空域电磁场的正演计算根据麦克斯韦定律：式中：D=εE；B=μH；D为位移电流；E为电场；H为磁场；B为磁通量；σ为介质的电导率；ε为介质的电容率；μ为介质的磁导率；Jm为磁流；Je为电流；r和t为空间坐标和时间坐标。
横断电场是由垂直磁偶极子源产生的（Ward，Hohmann，1988）。本文在正演模拟中采用水平圆环发射器，利用横断电场计算瞬时电压。一维电磁场正演模拟计算效率高，得到了广泛发展（Seng-piel，2010；Beard，2000；Siemon，2009）。由于近地表区域结构复杂，横向变化大，一维电阻率结构不能很好地反映地下结构，因此需要二维或者三维反演。三维电磁数据可以提供丰富的地下结构信息，但计算量大且耗时（Cox et al，2010；Haber etal，2004）。Auken和Christiansen（2004）首次提出用一维正演模拟，加入水平方向约束反演出伪二维电阻率模型。Schultz和Ruppel（2005）通过多组一维数据，在反演目标函数中引入二维正则化项，同时反演出二维电阻率模型，并应用到实际数据中。为了提高计算效率，本文根据一维频率域航空电磁数据的正演模拟理论，提出了应用Tikhonov正则化矩阵加强相邻网格点之间的约束，实现反演三维电阻率的伪三维算法。
1.3 联合反演方法
地震走时层析成像可以解决近地表地面起伏大的问题，得到地下速度结构，但不能反演出隐藏的低速层。电磁数据反演可以得到高电阻率区和低电阻率区的分布，但反演结果的垂直分辨率较低。如果将地震走时数据和电磁数据联合起来，同时反演地下的速度和电阻率，能够为地球物理解释提供更加丰富的信息，降低结果的非唯一性。
本文联合三维地震走时计算和一维航空电磁模拟，通过三维Tikhonov正则化的应用，在三维交叉梯度的结构约束作用下，同时反演出近地表三维速度模型和电阻率结构。联合反演的目标函数为：式中：mr为电阻率模型；ωr表示电磁数据的权重；dr表示观测到的数据；Gr表示正演模拟的数据；L为Tikhonov正则化矩阵；τr为正则化项的平滑系数；MS为慢度模型；Gs为正演模拟的走时数据；ds为观测到的走时数据；ωs为数据的权重；τs为模型的平滑系数；t是三维二参量的交叉梯度；τ为交叉梯度的权重。交叉梯度因子根据下列公式计算（Gallardo，Meju，2007）：
t（x，y，z）=▽log（mr（x，y，z））×▽MS（x，y，z）=txi+tyj+tzk（4）为了得到迭代反演的模型更新量，对公式（3）进行泰勒展开，得到：式中：β为加在联合反演上平衡mr和MS模型更新量的参数，定义为：式中，mws和mwr分别代表模型的权重，由于电阻率的对数比慢度近似高3个量级，因此，模型权重比值乘以系数1000。三维模型反演的参数多，计算数据量大，因此利用共扼梯度算法求解公式（5），得到模型更新量。
2 棋盘模型测试
为了对比联合反演和单独的走时反演、电磁反演的结果，本文对棋盘模型进行测试，真实模型如图1所示。图1a～d展示了不同深度速度模型的4个切面，可以看出速度模型内部共有9个异常块体，5个高速体和4个低速体交错放置。图1e～h展示了不同深度的电阻率模型的4个切面，同样可以看出，在电阻率模型内部，高速体位置对应高电阻率异常，低速体位置对应低电阻率异常。速度模型的观测系统，炮点在地下10m，检波器在地表，炮间距为250m，检波器的间距为50m，共有169个炮点，5929个检波器。电阻率模型的观测系统，发射接收装置离地表50m，在X方向和Y方向的间距均是200m，每个发射器伴随着一个接收器，二者之间的水平距离为7.9m。
首先，分别进行单独的地震初至波走时层析成像和伪三维航空电磁反演，将没有异常的背景速度结构和背景电阻率结构作为反演的初始模型，得到三维速度模型和三维电阻率模型（图2）。然后，将单独反演的结果作为联合反演的初始输入模型。进行联合反演之前，我们首先需要讨论在目标函数中应用的权重参数（τr，τs，ωs，β和γ）的选择。根据数据和模型的二阶范数的L曲线，选定速度和电阻率模型的正则化项平滑因子τr和τs分别为0.001和0.5。设定参数ωr=1，ωs=1，β=1，然后测试了不同的交叉梯度的权重下的联合反演。图3为不同γ值反演的地震走时数据不匹配度、频率域航空电磁数据（FrequencyAirborne Electromagnetic，简称FAEM）不匹配度、归一化的交叉梯度值以及目标函数值随迭代次数的变化，交叉梯度在第三次反演迭代时引入。图3c显示随着γ的增大，交叉梯度对反演的作用越来越大，电阻率模型和速度模型在结构上越来越相似，作为补偿，走时数据的不匹配度逐渐增加。在4组测试中，我们选择较合理的γ=3时的反演结果展示在图4中。
从图4a～d中反演的速度模型中可以看出，Z为250m，300m和400m的切面中有明显的低速体，说明联合反演中交叉梯度项起作用，使得电阻率模型的信息传递给了速度模型。Z为450m的切面中低速体不明显，是由于电阻率反演在深度上的分辨率低，且航空磁测的探测深度有限，深部的电阻率反演結果和真实模型有差距，导致传递给速度模型的信息有差异，因此深部的速度模型也没有明显的低速结构。
通过以上合成模型测试，可知地震走时数据和航空电磁数据联合反演的方法能够得到比单独的地震走时反演近地表速度和单独的电磁数据反演近地表电阻率更准确的结果。通过交叉梯度的作用，速度模型和电阻率模型互相传递信息，使同一地区的2个物理特征量保持结构上的一致。
3 实际数据应用
我们将联合地震走时数据和航空电磁数据反演三维速度和电阻率结构的方法应用到野外采集的实际数据中。加拿大阿尔伯塔地区的地表起伏严重，走时数据采集有难度，因此走时数据的观测系统比电磁数据的观测系统较稀疏。地震数据观测系统在X方向覆盖10000m，Y方向覆盖7000m，有3542个炮点，4621个接收器，共接收了6355116条地震记录，炮间距和接收器间距都为50m。电磁数据采集仪器为Fugro Airborne RESOLVE Sys-tem，观测系统共有6933个发射接收装置，每条测线的间距为100m，发射器和接收器水平距离为7.9m，高程随地表变化，平均在地表以上35m的位置，接收到数据包括二次场与一次场的比值的实部和虚部，共有5个频率（0.390kHz，1.787kHz，8.391kHz，40.740kHz，131.630kHz）。
首先根据地震记录在共炮域拾取初至波走时，利用拾取到的走时文件进行初至波走时反演。从图5a～d反演的速度模型4个Y方向切面（4151m，5851m，7351m和8551m处）可以看出，该研究地区地表起伏大，速度分层明显，从地表到400m处大致为一层，400～600m大致为第二层。在浅层分布着一些离散的、小的低速体，速度横向不均匀。利用航空电磁数据进行单独的伪三维反演出电阻率模型，从图5e～h可以看出，在浅层存在明显的低电阻率层，在Y=8551m的切面中，X从60000m到63000m之间凸起地表处有明显的高电阻率异常区。
进一步将单独的反演结果作为联合反演的初始速度模型和初始电阻率模型输入，设置不同的权重因子，比较反演结果。固定wr=1，ωs=1，β=1，测试交叉梯度权重因子γ分别为0.1、0.5、1.0和10.0作用下电磁数据不匹配度、走时数据不匹配度、归一化的交叉梯度和目标函数随迭代次数的变化，如图6所示。由图6a看出，随着交叉梯度权重的增大，电磁数据的不匹配度没有明显的变化，仍然能够降到很小的值。由图6b看出，当γ<1时，走时数据的不匹配度仍然会随着迭代次数的增加而降低，但当γ=10时，在交叉梯度引入的第三次迭代后，走时不匹配度陡增，而且随着反演迭代次数增加而增加，这是由于交叉梯度的权重过大，强烈限制了速度模型和电阻率模型在结构上逐渐一致，而电阻率模型的分辨率较低，传递给速度模型的结构信息不准，所以在走时数据上会出现不匹配度增加的现象。由图6c看出，随着γ的增加，速度模型和电阻率模型在结构上越来越相似，交叉梯度的值越来越小，这也验证了交叉梯度在联合反演中对模型结构上约束的作用。当γ<1时，图6d中目标函数的值仍然随着迭代次数的增加而降低，γ>1时，其值随着迭代次数增加而增加，与走时数据的变化相同。进一步对比γ分别为 0.1和0.5的情况，在第四次迭代之前二者数值几乎重叠，从第五次迭代开始，γ=0.5的曲线整体在γ=0.1曲线的下方，取值更小。也就是说，当γ=1时，交叉梯度的作用变大，联合反演中模型的匹配度开始以牺牲数据的拟合度来实现，速度模型和电阻率模型在结构上的相似度较高。因此，本文认为对这组地震走时数据和航空电磁数据的联合反演，γ=1时反演的速度模型和电阻率模型较为可靠。
联合反演的速度模型和电阻率模型如图7所示。反演的速度模型中，图7a中的浅层部分出现了明显的低速块体，在图7d中低速体的边缘更加清晰。反演的电阻率模型中，图7e中浅层出现横向不均匀现象，与速度模型有相似结构，图7f中在薄的低电阻率结构下分布着不均匀的高电阻率结构，图7g、7h中在X为60000～62000m，凸起的地表处浅层的电阻率值偏低，这是由于在速度结构中此区域为低速区，交叉梯度的作用使得速度结构信息传递给了电阻率结构，使二者在结构上趋于一致。
为了进一步对比单独反演算法和联合反演算法得到的速度结构的差异，我们分别在2个速度模型的基础上对地震数据进行静校正，然后实施叠加处理，如图8所示，箭头指示了叠加剖面变化明显的位置，尤其在图8b和d中，基于联合反演得到的速度结构计算出来的静校正处理得到的叠加剖面图显示了反射层较好的连续性和较少的起伏变化。从对比的结果可以看出，整体上联合反演得到的速度结构产生的叠加剖面要优于单独的走时反演得到的速度结构产生的叠加剖面。因此，联合航空电磁数据和地震走时数据反演近地表的电阻率和速度结构可以提供更准确的速度值来计算静校正，有利于我们进行后续的地震数据处理和地下深部成像。
4 结论与讨论
本文提出了利用地震初至波走时数据和航空电磁数据联合反演近地表三维速度和电阻率结构的方法，该方法应用了三维地震走时反演，一维电磁正演，结合三维Tikhonov正则化，并利用三维交叉梯度算子对结构的约束，实现了伪三维反演，得到近地表的三维速度结构和电阻率结构。这种伪三维反演算法大大提高了计算效率，节省了存储空间。电磁反演能对高电阻率区和低电阻率区成像，初至波走时反演对低速区成像有困难，而电磁数据和走时数据可以互相补充，联合反演同一区域的速度和电阻率结构，提高速度成像质量。通过对理论模型测试和实际数据应用，验证了这种联合反演方法的可行性和结果的可靠性。
本文仔细讨论了不同的权重因子在联合反演中的作用。交叉梯度的权重过小，联合反演的模型之间传递的结构信息不足以影响联合反演结果；交叉梯度的权重过大，数据的不匹配度增加，反演结果出现假象。因此，在联合反演中，需要根据走时数据、电磁数据、交叉梯度以及目标函数随反演迭代次数变化的曲线判断权重因子的合理取值。在实际数据的应用中，本文对比了联合反演的速度结果和单独的走时反演的速度结果，然后分别基于这2种速度结构對数据做了静校正和叠加处理。叠加剖面显示，联合反演的速度结构优于单独反演结果。因此，电磁数据和走时数据的联合反演有助于得到更接近地下真实结构的模型。
航空电磁数据具有采集方便、数据量大的特点，能够反映地下电阻率的分布。航空电磁数据和地震走时数据联合反演近地表的电阻率结构和速度结构，为地球物理解释提供了更丰富的信息，也降低了反演的非唯一性问题。非地震数据和地震数据的结合有效地扬长避短，在勘探地球物理成像领域具有广泛的应用前景。因此，在未来联合更多的数据，比如重力数据、水文监测数据、波形数据等，对研究地球物理成像有极大的推动作用。
参考文献：
王西文.1997.用重力、地震资料联合反演直接预测油气藏的方法[J].石油地球物理勘探，32（2）；221-228.
周辉，杨宝俊，刘财.1994.地震—大地电磁测深资料综合反演[J].地球物理学报，37（增刊1），471-485.
Auken E，Christiansen A V.2004.Layered and laterally constrained 2Dinversion of resistivity data[J].Geophysics，69（3）：752-761.
Anderson A H，Kak A C.1982.Digital ray tracing in two-dimensional re-fractive fields[J].Journal of the A coustical Society of America，72（5）：1593-1606.
Beard L P.2000.Comparison of methods for estimating earth resistivityfrom airborne electromagnetic measurements[J].Journal of AppliedGeophysics，45（4）：239-259.
Bosch M，McGaughey J.2001.Joint inversion of gravity and magnetic dataunder lithologic constraints[J].The Leading Edge，20（8）：877-881.
Cox L H，Wilson G A，Zhdanov M S.2010.3D inversion of airborne elec-tromagnetic data using a moving footprint[J].Exploration Geophys-ics，41（4）：250-259.
Fregoso E，Gallardo L A.2009.Cross-gradients joint 3D inversion withapplications to gravity and magnetic data[J].Geophysics，74（4）：L31-L42.
Gallardo L A，Meju M A.2003.Characterization of heterogeneous near-surface materials by joint 2D inversion of DC resistivity and seismicdata[J].Geophysical Research Letters，30（13）：1658.
Gallardo L A，Meju M A.2004.Joint two-dimensional DC resistivity andseismic travel time inversion with cross-gradients constraints[J].Journal of Geophysical Research，109（B3），B03311，1-11.
Gallardo L A，Meju M A.2007.Joint 2D cross-gradients imaging of mag-netotelluric and seismic travel-time data for structural and litholog-ical classification[J].Geophysical Journal International，169（3）：1261-1272.
Graeber F M，Asch G.1999.Three-dimensional models of P wave veloci-ty and P-to-S velocity ratio in the southern central Andes by sim-ultaneous inversion of local earthquake data[J].Journal of Geophys-ical Research，104（B9）：20237-20256.
Gomez-Trevino E，Edwards R N.1983.Electromagnetic soundings in thesedimentary basin of southern Ontario-A case history[J].Geophys-ics，48（3）：311-330.
Haber E，Ascher U，Oldenburg D.2004.Inversion of 3D electromagneticdata in frequency and time domain using an inexact all-at-onceapproach[J].Geophysics，69（5）：1216-1228.
Hering A，Misiek R，Gyulai A，et al.1995.A joint inversion algorithm toprocess geoelectric and surface wave seismic data.Part Ⅰ：Basis ideas[J].Geophysical Prospecting，43（2）：135-156.
Hole J A，Zelt B C.1995.3-D finite-difference reflection travel times[J].Geophysical Journal International，121（2）：427-434.
Holladay J S，Lo B，Prinsenberg S J.1997.Bird orientation effects inquantitative airborne electromagnetic interpretation of pack ice thick-ness sounding[C].Oceans Conference，Marine Technology Society/Institute of Electrical and Electronics Engineers，Proceedings，2：1114-1119.
Julian B R，Gubbins D.1977.Three-dimensional seismic ray tracing[J].Journal of Geophysics，43：95-113.
Linde N，Binley A，Tryggvason A，et al.2006.Improved hydrogeophysicalcharacterization using joint inversion of cross-hole electrical resist-ance and ground-penetrating radar traveltime data[J].Water Re-sources Research，42（12）：W12404.
Meju M A，Gallardo L A.2003.Evidence for correlation of electrical resis-tivity and seismic velocity in heterogeneous near-surface materials[J].Geophysical Research Letters，30（7）：1373.
Misiek R，Liebig A，Gyulai A，et al.1997.A joint inversion algorithm toprocess geoelectric and surface wave seismic data.Part Ⅱ：applica-tions[J].Geophysical Prospecting，45（1）：65-85.
Moser T J.1991.Shortest path calculation of seismic rays[J].Geophys-ics，56（1）：59-67.
Nakanishi I，Yamaguchi K.1986.A numerical experiment on nonlinearimage reconstruction from first-arrival times for two-dimensionalisland arc structure[J].Journal of Physics of the Earth，34（2）：195-201.
Ogunbo J N，Zhang J.2014.Joint seismic traveltime and TEM inversionfor near surface imaging[C].SEG：84th Annual International Meet-ing，2104-2108.
Pieta A，Bala J.2008.Joint inversion of direct current（DC）and electro-magnetic（EM）measurements in parallel computing environment[C].Barcelona；14th EAGE European Meeting of Environmentaland Engineering Geophysics.
Qin F，Luo Y，Olsen K B，et al.1992.Finite-difference solution of theeikonal equation along expanding wavefronts[J].Geophysics，57（3）.478-487.
Roth M，Holliger K.1998.Joint inversion of Rayleigh and guided waves inhigh-resolution seismic data using a genetic algorithm[C].SEG：67th Annual International Meeting，1570-1573.
Schultz G，Ruppel C.2005.Inversion of inductive electromagnetic data inhighly conductive terrains[J].Geophysics，70（1）：G16-G28.
Sengpiel K-P.2010.Approximate inversion of airborne EM data from amulti-layered ground[J].Geophysical Prospecting，36（4）：446-459.
Siemon B.2009.Levelling of helicopter-borne frequency-domain elec-tromagnetic data[J].Journal of Applied Geophysics，67（3）：206-218.
Tikhonov A N，Arsenin V Y.1977.Solutions of Ill-Posed Problems[M].New York：Wiley.
Tryggvason A，Rognvaldsson S Th，Flovenz G.2002.Three-dimensionalimaging of the P-and S-wave velocity structure and earthquake lo-cations beneath Southwest Iceland[J].Geophysical Journal Interna-tional，151（3）：848-866.
Van Trier J，Symes W W.1991.Unwind finite-difference calculation oftraveltimes[J].Geophysics，56（6）：812-821.
Vidae J E.1990.Finite-difference calculations of traveltimes in three di-mensions[J].Geophysics，55（5）：521-516.
Vidale J E.1988.Finite-difference calculations of traveltimes[J].Bulle-tin of the Seismological America，78（6）：2062-2076.
Vozoff K，Jupp D L B.1975.Joint inversion of geophysical data[J].Geo-physical Journal International，42（3）：977-991.
Ward S H，Hohmann G W.1988.Electromagnetic theory for geophysicalapplications[M].SEG，131-311.
William P R.1990.Tomographic inversion in reflection seismology[J].Geophysical Journal International，100（2）：255-274.
摘要：地震走时数据可以反演近地表的速度，但不能反演出隐蔽层和低速层。航空电磁数据可以反演近地表的高电阻率和低电阻率，但是对垂直方向的分辨率低。联合三维地震走时数据和航空电磁数据反演了近地表的速度和电阻率结构，并采用棋盘模型测试和实际数据应用展示了联合反演出的速度模型优于单独的走时反演出的速度模型。结果表明：该算法可以应用到处理数据量大的勘探地震成像中，能提供优化的速度结构。
关键词：联合反演；航空电磁；地震走时；近地表
中图分类号：P313.23，P315.2 文献标识码：A 文章编号：1000-0666（2018）01-0022-10
0 引言
地球物理反演是指利用接收到的信号推测地下介质的物理性质及结构分布，不同信号能够反演出不同的物理参数信息，为地球物理解释提供证据，有助于寻找矿产、油藏等资源。近地表地质结构复杂，分布着沼泽、沙漠、戈壁、风化层等松软结构，而且存在山地等起伏地形，给地下介质反演带来很大困难。因此，利用有限的数据资料反演出近地表的地质结构，有利于地球物理深部成像。
联合反演在地球物理学的发展始于由Vozoff和Jupp（1975）提出的利用直流电磁测深数据和大地电磁测深数据联合反演一维层状电阻率模型，采用的反演算法为迭代二阶马奎特阻尼最小二乘法。随后，联合不同数据的反演逐渐发展，Gomez-Tre-vino和Edwards（1983）联合可控源电磁数据和直流电测深数据利用奇异值分解法（SVD）反演了地层的电阻率和厚度，Roth和Holliger（1998）利用高分辨率的瑞雷波和导波联合反演的P波速度和S波速度，Bosch和McGaughey（2001）在先验条件下利用重力数据和电磁数据联合反演二维密度结构和磁导率分布，Tryggvason等（2002）利用P波走时和S波走时信息联合反演纵横波速度和地震定位，Pieta和Bala（2008）联合直流电测深数据和电磁数据应用蒙特卡罗全局优化算法实现了并行计算反演视电阻率等。
地震数据和非地震数据之间的联合反演，也得到广泛发展。周辉等（1994）利用广义线性反演方法及非线性反演的预条件最速下降法开展了一维地震一大地电磁测深资料反演方法研究，得出顺序反演效果优于地震、电磁单独反演效果，而同时反演的效果最优。王西文（1997）采用剥离法对重力、地震资料联合反演目的层密度值，进而预测油气藏。Hering等（1995）利用一维直流电测深数据和地震面波数据联合反演得到近地表的电阻率模型和面波速度模型。Misiek等（1997）利用电法数据和勒夫波数据联合反演了Borsod地区近地表的速度和电阻率结构。Meju和Gallardo（2003）验证了近地表地区电阻率和纵波速度存在结构一致的现象，提出了交叉梯度算法，实现了二维直流电测深数据和地震走时数据的联合反演（Gallardo，Meju，2003，2004），以及二维大地电磁和地震走时数据的联合反演（Gallardo，Meju，2007）。交叉梯度在联合反演中得到极大地推广（Linde et al，2006；Fregoso，Gallardo，2009；Ogunbo，Zhang，2014）。
本文利用地震走时数据和航空电磁数据，应用交叉梯度对模型的结构限制，联合反演了近地表的速度和电阻率模型。地震走时正演是三维模拟，航空电磁正演是一维模拟，通过三维Tik-honov正则化（Tikhonov，Arsenin，1977），联合反演出三维速度模型和三维电阻率模型，再对合成棋盘模型进行了測试，并应用到一组实际数据中，检验反演结果。
1 研究方法
1.1 地震走时正演模拟
地震波走时正演模拟的目的是计算波从源到接收器的时间以及射线路径，主要有基于射线的方法和基于网格的方法。基于射线的方法简单高效，能解决各向同性和各向异性介质，模拟的到时相对准确，比如打靶法（Anderson，Kak，1982；Wil-lianMSon，1990）和弯曲法（Julian，Gubbins，1977；Graeber，Asch，1999）。基于网格的方法种类繁多，应用原理各不相同，可以计算模型中每个网格点的到时，能够稳定处理各向同性介质，并且能同时高效计算出到时和射线路径，如求解程函方程的有限差分法（Vidale，1988，1990）和应用费马原理的最短路径法（Nakanishi，Yamaguchi，1986；Moser，1991）。
求解程函方程和最短路径的方法可以提供初至波到时，而打靶法和弯曲法只能得到一个不能确定是初至波的走时，因此，本文采用了求解程函方程的方法正演模拟地震走时和射线路径。射线追踪的程函方程为（Vidale，1988，1990）：式中：s为慢度；t为走时。此方法是采用有限差分向外扩展的方式，利用已知点的到时寻找到达下一个网格点的用时最短的路径。同时，可以根据最终走时的梯度从接收点回溯到源，得到射线路径。这种方法得到的地震走时准确而且计算效率高，计算速度与定义的网格点个数成正比。由于这种立方体向外扩展计算波前的方法不总是遵循波在介质中传播的方向，存在波前不连续的问题，比如，初至波可能在传回立方体之前需要对立方体外的区域再进行扩展计算。但是随着后续的发展，Vidale提出的方法在地震波走时计算中得到广泛的应用（VanTrier，Symes，1991；Qin et al，1992；Hole，Zelt，1995）。
1.2 航空电磁正演模拟
航空电磁包括频率域系统和时间域系统，时间域电磁场通常可提供比频率域地球更深部的信息（Holladay et al，1997），而频率域电磁场能得到高分辨率的浅层信息，因此，在近地表研究中，本文采用频率域电磁场数据反演浅层的电阻率结构。
时空域电磁场的正演计算根据麦克斯韦定律：式中：D=εE；B=μH；D为位移电流；E为电场；H为磁场；B为磁通量；σ为介质的电导率；ε为介质的电容率；μ为介质的磁导率；Jm为磁流；Je为电流；r和t为空间坐标和时间坐标。
横断电场是由垂直磁偶极子源产生的（Ward，Hohmann，1988）。本文在正演模拟中采用水平圆环发射器，利用横断电场计算瞬时电压。一维电磁场正演模拟计算效率高，得到了广泛发展（Seng-piel，2010；Beard，2000；Siemon，2009）。由于近地表区域结构复杂，横向变化大，一维电阻率结构不能很好地反映地下结构，因此需要二维或者三维反演。三维电磁数据可以提供丰富的地下结构信息，但计算量大且耗时（Cox et al，2010；Haber etal，2004）。Auken和Christiansen（2004）首次提出用一维正演模拟，加入水平方向约束反演出伪二维电阻率模型。Schultz和Ruppel（2005）通过多组一维数据，在反演目标函数中引入二维正则化项，同时反演出二维电阻率模型，并应用到实际数据中。为了提高计算效率，本文根据一维频率域航空电磁数据的正演模拟理论，提出了应用Tikhonov正则化矩阵加强相邻网格点之间的约束，实现反演三维电阻率的伪三维算法。
1.3 联合反演方法
地震走时层析成像可以解决近地表地面起伏大的问题，得到地下速度结构，但不能反演出隐藏的低速层。电磁数据反演可以得到高电阻率区和低电阻率区的分布，但反演结果的垂直分辨率较低。如果将地震走时数据和电磁数据联合起来，同时反演地下的速度和电阻率，能够为地球物理解释提供更加丰富的信息，降低结果的非唯一性。
本文联合三维地震走时计算和一维航空电磁模拟，通过三维Tikhonov正则化的应用，在三维交叉梯度的结构约束作用下，同时反演出近地表三维速度模型和电阻率结构。联合反演的目标函数为：式中：mr为电阻率模型；ωr表示电磁数据的权重；dr表示观测到的数据；Gr表示正演模拟的数据；L为Tikhonov正则化矩阵；τr为正则化项的平滑系数；MS为慢度模型；Gs为正演模拟的走时数据；ds为观测到的走时数据；ωs为数据的权重；τs为模型的平滑系数；t是三维二参量的交叉梯度；τ为交叉梯度的权重。交叉梯度因子根据下列公式计算（Gallardo，Meju，2007）：
t（x，y，z）=▽log（mr（x，y，z））×▽MS（x，y，z）=txi+tyj+tzk（4）为了得到迭代反演的模型更新量，对公式（3）进行泰勒展开，得到：式中：β为加在联合反演上平衡mr和MS模型更新量的参数，定义为：式中，mws和mwr分别代表模型的权重，由于电阻率的对数比慢度近似高3个量级，因此，模型权重比值乘以系数1000。三维模型反演的参数多，计算数据量大，因此利用共扼梯度算法求解公式（5），得到模型更新量。
2 棋盘模型测试
为了对比联合反演和单独的走时反演、电磁反演的结果，本文对棋盘模型进行测试，真实模型如图1所示。图1a～d展示了不同深度速度模型的4个切面，可以看出速度模型内部共有9个异常块体，5个高速体和4个低速体交错放置。图1e～h展示了不同深度的电阻率模型的4个切面，同样可以看出，在电阻率模型内部，高速体位置对应高电阻率异常，低速体位置对应低电阻率异常。速度模型的观测系统，炮点在地下10m，检波器在地表，炮间距为250m，检波器的间距为50m，共有169个炮点，5929个检波器。电阻率模型的观测系统，发射接收装置离地表50m，在X方向和Y方向的间距均是200m，每个发射器伴随着一个接收器，二者之间的水平距离为7.9m。
首先，分别进行单独的地震初至波走时层析成像和伪三维航空电磁反演，将没有异常的背景速度结构和背景电阻率结构作为反演的初始模型，得到三维速度模型和三维电阻率模型（图2）。然后，将单独反演的结果作为联合反演的初始输入模型。进行联合反演之前，我们首先需要讨论在目标函数中应用的权重参数（τr，τs，ωs，β和γ）的选择。根据数据和模型的二阶范数的L曲线，选定速度和电阻率模型的正则化项平滑因子τr和τs分别为0.001和0.5。设定参数ωr=1，ωs=1，β=1，然后测试了不同的交叉梯度的权重下的联合反演。图3为不同γ值反演的地震走时数据不匹配度、频率域航空电磁数据（FrequencyAirborne Electromagnetic，简称FAEM）不匹配度、归一化的交叉梯度值以及目标函数值随迭代次数的变化，交叉梯度在第三次反演迭代时引入。图3c显示随着γ的增大，交叉梯度对反演的作用越来越大，电阻率模型和速度模型在结构上越来越相似，作为补偿，走时数据的不匹配度逐渐增加。在4组测试中，我们选择较合理的γ=3时的反演结果展示在图4中。
从图4a～d中反演的速度模型中可以看出，Z为250m，300m和400m的切面中有明显的低速体，说明联合反演中交叉梯度项起作用，使得电阻率模型的信息传递给了速度模型。Z为450m的切面中低速体不明显，是由于电阻率反演在深度上的分辨率低，且航空磁测的探测深度有限，深部的电阻率反演結果和真实模型有差距，导致传递给速度模型的信息有差异，因此深部的速度模型也没有明显的低速结构。
通过以上合成模型测试，可知地震走时数据和航空电磁数据联合反演的方法能够得到比单独的地震走时反演近地表速度和单独的电磁数据反演近地表电阻率更准确的结果。通过交叉梯度的作用，速度模型和电阻率模型互相传递信息，使同一地区的2个物理特征量保持结构上的一致。
3 实际数据应用
我们将联合地震走时数据和航空电磁数据反演三维速度和电阻率结构的方法应用到野外采集的实际数据中。加拿大阿尔伯塔地区的地表起伏严重，走时数据采集有难度，因此走时数据的观测系统比电磁数据的观测系统较稀疏。地震数据观测系统在X方向覆盖10000m，Y方向覆盖7000m，有3542个炮点，4621个接收器，共接收了6355116条地震记录，炮间距和接收器间距都为50m。电磁数据采集仪器为Fugro Airborne RESOLVE Sys-tem，观测系统共有6933个发射接收装置，每条测线的间距为100m，发射器和接收器水平距离为7.9m，高程随地表变化，平均在地表以上35m的位置，接收到数据包括二次场与一次场的比值的实部和虚部，共有5个频率（0.390kHz，1.787kHz，8.391kHz，40.740kHz，131.630kHz）。
首先根据地震记录在共炮域拾取初至波走时，利用拾取到的走时文件进行初至波走时反演。从图5a～d反演的速度模型4个Y方向切面（4151m，5851m，7351m和8551m处）可以看出，该研究地区地表起伏大，速度分层明显，从地表到400m处大致为一层，400～600m大致为第二层。在浅层分布着一些离散的、小的低速体，速度横向不均匀。利用航空电磁数据进行单独的伪三维反演出电阻率模型，从图5e～h可以看出，在浅层存在明显的低电阻率层，在Y=8551m的切面中，X从60000m到63000m之间凸起地表处有明显的高电阻率异常区。
进一步将单独的反演结果作为联合反演的初始速度模型和初始电阻率模型输入，设置不同的权重因子，比较反演结果。固定wr=1，ωs=1，β=1，测试交叉梯度权重因子γ分别为0.1、0.5、1.0和10.0作用下电磁数据不匹配度、走时数据不匹配度、归一化的交叉梯度和目标函数随迭代次数的变化，如图6所示。由图6a看出，随着交叉梯度权重的增大，电磁数据的不匹配度没有明显的变化，仍然能够降到很小的值。由图6b看出，当γ<1时，走时数据的不匹配度仍然会随着迭代次数的增加而降低，但当γ=10时，在交叉梯度引入的第三次迭代后，走时不匹配度陡增，而且随着反演迭代次数增加而增加，这是由于交叉梯度的权重过大，强烈限制了速度模型和电阻率模型在结构上逐渐一致，而电阻率模型的分辨率较低，传递给速度模型的结构信息不准，所以在走时数据上会出现不匹配度增加的现象。由图6c看出，随着γ的增加，速度模型和电阻率模型在结构上越来越相似，交叉梯度的值越来越小，这也验证了交叉梯度在联合反演中对模型结构上约束的作用。当γ<1时，图6d中目标函数的值仍然随着迭代次数的增加而降低，γ>1时，其值随着迭代次数增加而增加，与走时数据的变化相同。进一步对比γ分别为 0.1和0.5的情况，在第四次迭代之前二者数值几乎重叠，从第五次迭代开始，γ=0.5的曲线整体在γ=0.1曲线的下方，取值更小。也就是说，当γ=1时，交叉梯度的作用变大，联合反演中模型的匹配度开始以牺牲数据的拟合度来实现，速度模型和电阻率模型在结构上的相似度较高。因此，本文认为对这组地震走时数据和航空电磁数据的联合反演，γ=1时反演的速度模型和电阻率模型较为可靠。
联合反演的速度模型和电阻率模型如图7所示。反演的速度模型中，图7a中的浅层部分出现了明显的低速块体，在图7d中低速体的边缘更加清晰。反演的电阻率模型中，图7e中浅层出现横向不均匀现象，与速度模型有相似结构，图7f中在薄的低电阻率结构下分布着不均匀的高电阻率结构，图7g、7h中在X为60000～62000m，凸起的地表处浅层的电阻率值偏低，这是由于在速度结构中此区域为低速区，交叉梯度的作用使得速度结构信息传递给了电阻率结构，使二者在结构上趋于一致。
为了进一步对比单独反演算法和联合反演算法得到的速度结构的差异，我们分别在2个速度模型的基础上对地震数据进行静校正，然后实施叠加处理，如图8所示，箭头指示了叠加剖面变化明显的位置，尤其在图8b和d中，基于联合反演得到的速度结构计算出来的静校正处理得到的叠加剖面图显示了反射层较好的连续性和较少的起伏变化。从对比的结果可以看出，整体上联合反演得到的速度结构产生的叠加剖面要优于单独的走时反演得到的速度结构产生的叠加剖面。因此，联合航空电磁数据和地震走时数据反演近地表的电阻率和速度结构可以提供更准确的速度值来计算静校正，有利于我们进行后续的地震数据处理和地下深部成像。
4 结论与讨论
本文提出了利用地震初至波走时数据和航空电磁数据联合反演近地表三维速度和电阻率结构的方法，该方法应用了三维地震走时反演，一维电磁正演，结合三维Tikhonov正则化，并利用三维交叉梯度算子对结构的约束，实现了伪三维反演，得到近地表的三维速度结构和电阻率结构。这种伪三维反演算法大大提高了计算效率，节省了存储空间。电磁反演能对高电阻率区和低电阻率区成像，初至波走时反演对低速区成像有困难，而电磁数据和走时数据可以互相补充，联合反演同一区域的速度和电阻率结构，提高速度成像质量。通过对理论模型测试和实际数据应用，验证了这种联合反演方法的可行性和结果的可靠性。
本文仔细讨论了不同的权重因子在联合反演中的作用。交叉梯度的权重过小，联合反演的模型之间传递的结构信息不足以影响联合反演结果；交叉梯度的权重过大，数据的不匹配度增加，反演结果出现假象。因此，在联合反演中，需要根据走时数据、电磁数据、交叉梯度以及目标函数随反演迭代次数变化的曲线判断权重因子的合理取值。在实际数据的应用中，本文对比了联合反演的速度结果和单独的走时反演的速度结果，然后分别基于这2种速度结构對数据做了静校正和叠加处理。叠加剖面显示，联合反演的速度结构优于单独反演结果。因此，电磁数据和走时数据的联合反演有助于得到更接近地下真实结构的模型。
航空电磁数据具有采集方便、数据量大的特点，能够反映地下电阻率的分布。航空电磁数据和地震走时数据联合反演近地表的电阻率结构和速度结构，为地球物理解释提供了更丰富的信息，也降低了反演的非唯一性问题。非地震数据和地震数据的结合有效地扬长避短，在勘探地球物理成像领域具有广泛的应用前景。因此，在未来联合更多的数据，比如重力数据、水文监测数据、波形数据等，对研究地球物理成像有极大的推动作用。
参考文献：
王西文.1997.用重力、地震资料联合反演直接预测油气藏的方法[J].石油地球物理勘探，32（2）；221-228.
周辉，杨宝俊，刘财.1994.地震—大地电磁测深资料综合反演[J].地球物理学报，37（增刊1），471-485.
Auken E，Christiansen A V.2004.Layered and laterally constrained 2Dinversion of resistivity data[J].Geophysics，69（3）：752-761.
Anderson A H，Kak A C.1982.Digital ray tracing in two-dimensional re-fractive fields[J].Journal of the A coustical Society of America，72（5）：1593-1606.
Beard L P.2000.Comparison of methods for estimating earth resistivityfrom airborne electromagnetic measurements[J].Journal of AppliedGeophysics，45（4）：239-259.
Bosch M，McGaughey J.2001.Joint inversion of gravity and magnetic dataunder lithologic constraints[J].The Leading Edge，20（8）：877-881.
Cox L H，Wilson G A，Zhdanov M S.2010.3D inversion of airborne elec-tromagnetic data using a moving footprint[J].Exploration Geophys-ics，41（4）：250-259.
Fregoso E，Gallardo L A.2009.Cross-gradients joint 3D inversion withapplications to gravity and magnetic data[J].Geophysics，74（4）：L31-L42.
Gallardo L A，Meju M A.2003.Characterization of heterogeneous near-surface materials by joint 2D inversion of DC resistivity and seismicdata[J].Geophysical Research Letters，30（13）：1658.
Gallardo L A，Meju M A.2004.Joint two-dimensional DC resistivity andseismic travel time inversion with cross-gradients constraints[J].Journal of Geophysical Research，109（B3），B03311，1-11.
Gallardo L A，Meju M A.2007.Joint 2D cross-gradients imaging of mag-netotelluric and seismic travel-time data for structural and litholog-ical classification[J].Geophysical Journal International，169（3）：1261-1272.
Graeber F M，Asch G.1999.Three-dimensional models of P wave veloci-ty and P-to-S velocity ratio in the southern central Andes by sim-ultaneous inversion of local earthquake data[J].Journal of Geophys-ical Research，104（B9）：20237-20256.
Gomez-Trevino E，Edwards R N.1983.Electromagnetic soundings in thesedimentary basin of southern Ontario-A case history[J].Geophys-ics，48（3）：311-330.
Haber E，Ascher U，Oldenburg D.2004.Inversion of 3D electromagneticdata in frequency and time domain using an inexact all-at-onceapproach[J].Geophysics，69（5）：1216-1228.
Hering A，Misiek R，Gyulai A，et al.1995.A joint inversion algorithm toprocess geoelectric and surface wave seismic data.Part Ⅰ：Basis ideas[J].Geophysical Prospecting，43（2）：135-156.
Hole J A，Zelt B C.1995.3-D finite-difference reflection travel times[J].Geophysical Journal International，121（2）：427-434.
Holladay J S，Lo B，Prinsenberg S J.1997.Bird orientation effects inquantitative airborne electromagnetic interpretation of pack ice thick-ness sounding[C].Oceans Conference，Marine Technology Society/Institute of Electrical and Electronics Engineers，Proceedings，2：1114-1119.
Julian B R，Gubbins D.1977.Three-dimensional seismic ray tracing[J].Journal of Geophysics，43：95-113.
Linde N，Binley A，Tryggvason A，et al.2006.Improved hydrogeophysicalcharacterization using joint inversion of cross-hole electrical resist-ance and ground-penetrating radar traveltime data[J].Water Re-sources Research，42（12）：W12404.
Meju M A，Gallardo L A.2003.Evidence for correlation of electrical resis-tivity and seismic velocity in heterogeneous near-surface materials[J].Geophysical Research Letters，30（7）：1373.
Misiek R，Liebig A，Gyulai A，et al.1997.A joint inversion algorithm toprocess geoelectric and surface wave seismic data.Part Ⅱ：applica-tions[J].Geophysical Prospecting，45（1）：65-85.
Moser T J.1991.Shortest path calculation of seismic rays[J].Geophys-ics，56（1）：59-67.
Nakanishi I，Yamaguchi K.1986.A numerical experiment on nonlinearimage reconstruction from first-arrival times for two-dimensionalisland arc structure[J].Journal of Physics of the Earth，34（2）：195-201.
Ogunbo J N，Zhang J.2014.Joint seismic traveltime and TEM inversionfor near surface imaging[C].SEG：84th Annual International Meet-ing，2104-2108.
Pieta A，Bala J.2008.Joint inversion of direct current（DC）and electro-magnetic（EM）measurements in parallel computing environment[C].Barcelona；14th EAGE European Meeting of Environmentaland Engineering Geophysics.
Qin F，Luo Y，Olsen K B，et al.1992.Finite-difference solution of theeikonal equation along expanding wavefronts[J].Geophysics，57（3）.478-487.
Roth M，Holliger K.1998.Joint inversion of Rayleigh and guided waves inhigh-resolution seismic data using a genetic algorithm[C].SEG：67th Annual International Meeting，1570-1573.
Schultz G，Ruppel C.2005.Inversion of inductive electromagnetic data inhighly conductive terrains[J].Geophysics，70（1）：G16-G28.
Sengpiel K-P.2010.Approximate inversion of airborne EM data from amulti-layered ground[J].Geophysical Prospecting，36（4）：446-459.
Siemon B.2009.Levelling of helicopter-borne frequency-domain elec-tromagnetic data[J].Journal of Applied Geophysics，67（3）：206-218.
Tikhonov A N，Arsenin V Y.1977.Solutions of Ill-Posed Problems[M].New York：Wiley.
Tryggvason A，Rognvaldsson S Th，Flovenz G.2002.Three-dimensionalimaging of the P-and S-wave velocity structure and earthquake lo-cations beneath Southwest Iceland[J].Geophysical Journal Interna-tional，151（3）：848-866.
Van Trier J，Symes W W.1991.Unwind finite-difference calculation oftraveltimes[J].Geophysics，56（6）：812-821.
Vidae J E.1990.Finite-difference calculations of traveltimes in three di-mensions[J].Geophysics，55（5）：521-516.
Vidale J E.1988.Finite-difference calculations of traveltimes[J].Bulle-tin of the Seismological America，78（6）：2062-2076.
Vozoff K，Jupp D L B.1975.Joint inversion of geophysical data[J].Geo-physical Journal International，42（3）：977-991.
Ward S H，Hohmann G W.1988.Electromagnetic theory for geophysicalapplications[M].SEG，131-311.
William P R.1990.Tomographic inversion in reflection seismology[J].Geophysical Journal International，100（2）：255-274.
