1 引言

浅层地震勘探与直流电阻率法是解决浅地表工程地质问题最常用的地球物理勘探方法.地震走时成像在气田(Korenaga et al., 1997)、油藏(瞿辰等, 2010)、二氧化碳存储(Morency et al., 2011; Zhang et al., 2012; Carcione et al., 2012)、岩溶勘察(段成龙等, 2013; Pepe et al., 2015)、含水层结构探测(Ehosioke and Fechner, 2014)、孤石探测(郭云峰等, 2015)、矿产资源勘探(Menu et al., 2015)等领域得到较多的应用.电阻率成像技术在隧道超前预报(李术才等, 2014)、孤石探测(王俊超等, 2012; 李术才等, 2015)、岩溶勘察(张文俊等, 2014)、岩体完整性评价(Bellmunt et al., 2012; 曹权等, 2013)、大坝渗漏探测(Greve et al., 2012)、障碍物探测(Leontarakis and Apostolopoulos, 2013)、二氧化碳地下储存监测(Schmidt-Hattenberger et al., 2014)等领域也有广泛的应用.

采用地震走时法或电阻率法进行单独成像，因观测覆盖角度、数据采样率、人文干扰、成像算法中存在问题导致单一地球物理方法存在多解性，无法对地下结构进行精确成像(胡祥云等, 2006).例如，地震走时数据对高速异常体敏感，因而对其成像效果较好，一般用于解决存在速度异常的构造问题；直流电阻率数据对低阻异常体敏感，因而对低阻体具有较好的分辨能力，一般用于解决地层的富水性问题.然而这两种方法均存在较为明显的缺点，即地震走时成像对低速异常体、电阻率成像对高阻异常体存在成像阴影问题.

为了解决单一反演方法存在的多解性问题，很多学者对多数据联合反演地下介质多属性的方法进行了较为深入的研究.Lines (1988)采用岩石物理属性之间的关系建立两种勘探方法数据之间的联合反演，Haber和Oldenburg (1997)采用两种物理模型之间的结构相似性联合反演，Bosch (1999)运用两种属性之间的统计学关系进行联合反演，这三种联合反演方式可归纳为两类：(1) 建立两种或多种模型之间属性关系的联合反演；(2) 采用两种模型的结构相似作为约束条件的联合反演.当属性之间的对应关系较为复杂或没有对应关系、关系较弱时，采用结构相似进行约束的联合反演更为合理(Haber and Oldenburg, 1997; Gallardo and Meju, 2004).Gallardo和Meju (2004)首次提出以模型之间交叉梯度为结构一致约束方式的二维直流电阻率与地震折射数据的联合反演，极大促进了基于交叉梯度结构相似约束联合反演方法的发展(Gallardo and Meju, 2004; Gallardo, 2007; Fregoso and Gallardo, 2009; Moorkamp et al., 2011; 彭淼, 2012; 周丽芬, 2012; Rovetta et al., 2013; Ogunbo and Zhang, 2014; Sánchez and Delgado, 2015; Bennington et al., 2015; 李兆祥等, 2015; Wang et al., 2015; 高级和张海江, 2016).

因不同物性模型参数及灵敏度矩阵元素的数值变化范围大，可能存在量级上的巨大差异，给联合反演中求解大型灵敏度矩阵造成极大困难.Bennington等(2015)采用对基于交叉梯度约束建立的大型联合反演系统正则化的方法，实现了三维天然地震和二维大地电磁的联合反演，一定程度解决了联合反演中大型矩阵求解不稳定性问题，但在实现过程中需要显式计算地震走时与大地电磁的灵敏度矩阵，需要较大的存储空间并对单独反演程序进行修改，造成该方法实现较为困难.Um等(2014)提出了分别进行数据拟合和交叉梯度结构约束的联合反演流程，通过交替迭代进行数据约束和交叉梯度结构约束实现联合反演.虽然该流程降低了实现联合反演的难度，但是由于数据拟合与结构约束完全分开，因此迭代交替反演时较难平衡数据拟合和结构约束的权重.

针对联合反演矩阵存储空间巨大、求解不稳定及难以实现的问题，本文提出一种新的基于模型更新量及交叉梯度结构约束结合的联合反演流程，通过把单独反演得到的模型更新量加入到交叉梯度控制的结构约束之中，解决了Um等(2014)把数据拟合和交叉梯度结构约束分开反演权重难以选择的问题，同时实现数据拟合与结构约束的平衡.在具体实现过程中，把交叉梯度约束从联合反演大矩阵中分离，解决了Bennington等(2015)需要较大存储空间的问题，而且无需大幅度修改单一方法数据反演程序，为多方法数据联合反演地下介质的多属性提供了一个新的思路.

在实现基于模型更新量与交叉梯度约束的联合反演算法的基础上，本文采用了三维跨孔地震走时与直流电阻率数据，测试了反演流程的稳定性，验证了新的反演策略的有效性.相比于单独反演结果，联合反演结果减少了地震走时成像的假异常，提高了直流电阻率成像的分辨率.另外，在单独反演中可以根据L曲线等选择最优反演参数，而联合反演只需要考虑结构约束与数据拟合程度的权重大小，有效降低了联合反演权重选择的困难.

2 交叉梯度结构约束方法

对于三维速度模型m_s和电阻率模型m_r，二者之间的交叉梯度定义为速度梯度和电阻率梯度的叉积：

(1)

在三维模型情况下，垂直于Y-Z、X-Z、X-Y平面即平行于X、Y、Z方向上的交叉梯度分量分别为t_x、t_y、t_z.其展开形式分别为：

在二维模型情况下，因模型沿Y方向变化为0，模型之间沿X与Z方向变化的交叉梯度只存在Y方向分量：

(2)

图 1显示了三维速度与电阻率模型交叉梯度性质，速度与电阻率模型大小均为64 m×48 m×32 m.图 1a为速度模型，该模型为均匀背景中包含一个高速与一个低速异常，背景速度为2 km·s^-1，高低速异常体的速度分别为3 km·s^-1和1 km·s^-1，大小为10 m×10 m×10 m.因异常体与背景速度具有差异性，所以其周围均具有梯度值，梯度如图中箭头方向所示.图 1b为电阻率模型，电阻率数值分布范围为100 Ωm至260 Ωm，电阻率数值只沿Z方向变化，电阻率变化规律为从Z=0 m的100 Ωm至Z=32 m的260 Ωm，变化梯度为5 Ωm·m^-1.由于电阻率只在垂向上(Z)变化，所以电阻率模型的梯度只在Z方向上有非零值.速度模型与电阻率模型之间的交叉梯度如图 1c所示，从该图可以看出交叉梯度仅在X与Y方向上有值，即两个模型梯度方向相同或相反时交叉梯度值为零，反之不为零.该性质为利用交叉梯度约束反演奠定了基础，即在反演时要求速度与电阻率模型交叉梯度值趋于零，模型之间具有一致的结构.

3 基于模型更新量约束的联合反演方法

联合反演方法主要包括建立反演目标函数、选择反演迭代方法、联合反演具体流程、迭代终止条件等，下面对这几部分分别予以阐述.

3.1 单独反演目标函数的建立

地震走时成像目标函数定义为(Aster et al., 2011)：

(3)

式中：f (m_s)为根据地震速度模型m_s计算的理论地震走时数据，d_s为观测走时数据，L为速度模型平滑因子，α为速度模型平滑权重，I为单位矩阵，β为模型阻尼权重.目标函数中采用L曲线调整权重因子α、β的大小从而达到数据拟合与模型平滑之间的平衡.在实现地震走时成像时，我们采用基于射线理论的层析成像方法.首先对模型区域网格进行剖分(Sambridge, 1990)，采用有限差分法求解程函方程(Vidale, 1988; Vidale, 1990)和MSFM (Multi-Stencils Fast Marching Methods)方法计算走时场(Osher and Fedkiw, 2006).为了提高反演结果的稳定性，采用二阶范数的吉洪诺夫正则化(Tikhonov and Arsenin, 1977)约束(Calvetti, 2000)，并采用阻尼最小二乘LSQR法迭代求解(Paige and Saunders, 1982).

电阻率成像目标函数φ(m)可分为两个部分，即数据拟合和模型正则化，

(4)

式中：f (m_r)为根据电阻率模型m_r计算的理论电位数据，d_r为观测电位数据，为参考模型权重，m_r^ref为电阻率参考模型，W_d为观测数据权系数矩阵、W_m为参考模型权系数矩阵，为数据拟合和参考模型之间的平衡权重因子.(4) 式中第二项为引入的正则化约束条件(Tikhonov and Arsenin, 1977).在直流电阻率成像中，首先采用有限差分法求解模型网格节点电位，然后采用伴随矩阵方法求解一阶和二阶偏导数即Jacobian和Hessian矩阵(Mcgillivray and Oldenburg, 1990; 吴小平和徐果明, 2000; Pidlisecky et al., 2007)，最后采用基于共轭梯度的高斯牛顿法迭代求解电阻率模型更新量.

3.2 联合反演算法

对于交叉梯度函数t (m_s, m_r)，利用一阶泰勒展开的方法进行线性化，可写成矩阵形式为：

(5)

式中为交叉梯度对速度的偏导数，为交叉梯度对电阻率的偏导数，t (m_s⁰, m_r⁰)为参考模型m_s⁰和m_r⁰之间的交叉梯度.

对于速度和电阻率模型，可以对二者进行改变使两个模型之间的交叉梯度为零，即最小化交叉梯度值，这样两个模型的结构就变得一致.在公式(5) 中，如果模型之间的交叉梯度t (m_s, m_r)为零，针对不同的交叉梯度分量可以写成下面的矩阵形式：

(6)

式中Δm_s、Δm_r分别为速度与电阻率模型更新量, 分别为交叉梯度在X、Y、Z三个方向上的分量对速度模型的偏导数，分别为交叉梯度在X、Y、Z三个方向上的分量对电阻率模型的偏导数，t_x⁰、t_y⁰、t_z⁰分别为上一次迭代结束时速度模型与电阻率模型之间的交叉梯度在X、Y、Z三个方向上的分量.

基于交叉梯度结构约束的联合反演系统，是通过把对应于数据拟合的单个反演系统和对应于结构约束的单个反演系统(公式(6))结合而建立起来的(Gallardo and Meju, 2004).但是因为不同数据所对应的灵敏度矩阵可能存在量级上的差异，而且不同数据的反演程序相差很大，因此很难把不同数据的反演系统以及交叉梯度反演系统融合起来实现联合反演.为了解决这个问题，Um等(2014)提出了一种新的联合反演策略，即首先基于不同数据的单独反演系统对模型进行更新对数据进行拟合，然后利用交叉梯度结构约束反演系统进一步更新模型使得两个模型的结构一致.这两个步骤交替进行直至数据拟合和交叉梯度最小化.这种联合反演方式把结构约束与数据拟合分开交替地进行，虽然可以直接利用现有的反演系统进行联合反演，但是在该流程实际实现的过程当中，很难平衡数据的拟合和结构一致的约束.为了解决这个问题，高级和张海江(2016)提出了一种新的基于交叉梯度交替结构约束的联合反演流程，即在进行地震走时反演时，用电阻率模型同时约束速度模型的结构，在直流电阻率反演时，同时用速度结构约束电阻率结构，这两个步骤交替进行直至收敛.这种交替结构约束的联合反演策略很容易实现结构约束和数据拟合之间的平衡.该反演策略要求把单独反演系统的灵敏度矩阵提取出来，和交叉梯度约束对应的灵敏度矩阵组合在一起，但是在某些情况下如三维直流电阻率反演很难显式提取出灵敏度矩阵.

所以为了解决上面的问题，基于Um等(2014)的联合反演策略，我们提出一种新的联合反演策略来实现数据拟合和结构一致性的平衡.把单独地震走时反演与直流电阻率反演得到的速度更新量Δm_s⁰与电阻率更新量Δm_r⁰与只基于交叉梯度结构约束的反演系统(6) 结合起来，即把Δm_s⁰、Δm_r⁰加入到基于交叉梯度的反演系统中，用于约束基于交叉梯度最小化的线性反演系统，得到新的速度更新量Δm_s和新的电阻率更新量Δm_r，实现数据拟合和结构一致性的平衡.在三维模型情况下，新的联合反演系统可以写成式(7)：

(7)

式中：α_s为速度更新量权重，α_r为电阻率更新量权重，β_t为结构一致性约束权重.在二维模型情况下，式(7) 可以简化为式(8) 的形式：

(8)

在公式(7) 和(8) 中，通过调整α_s，α_r和β_t三个参数的相对大小，就可以达到数据拟合和结构一致约束之间的平衡.

3.3 联合反演流程

图 2为本文提出的基于模型更新量及交叉梯度

结构约束的地震走时和电阻率联合反演流程图，包括以下几个主要步骤：

(1) 输入初始速度模型与电阻率模型，并计算二者之间的交叉梯度；

(2) 根据公式(3) 和(4) 定义的目标函数分别进行地震走时和电阻率单独成像，求解出速度模型更新量与电阻率模型更新量；

(3) 根据交叉梯度结构约束、速度和电阻率模型更新量，基于公式(7) 或者(8) 所表示的联合反演矩阵，求解得到新的速度模型与电阻率模型更新量；

(4) 更新速度与电阻率模型，并计算更新后的速度模型与电阻率模型之间的交叉梯度以及地震走时数据和电位数据的残差；

(5) 判断数据拟合和交叉梯度值是否满足迭代终止条件，若满足直接输出结果，若不满足返回步骤(2).

通常采用的迭代终止条件包括：迭代次数、每次迭代残差下降的梯度、残差剩余量占初始残差百分比、模型更新量大小等.本文采用的迭代终止即假设迭代过程已经收敛的条件为电阻率或地震走时残差下降量低于某一百分比时迭代终止(例如迭代残差下降量小于1%).

由于不同反演方法得到的模型更新量数量级具有较大差异，因此在联合反演中，不同方法之间模型更新量的权重选择极为重要.地震走时成像中参与计算的模型参数为慢度，其数量级为10^-3，模型更新量数据数量级为10^-4；电阻率成像中参与计算的参数为电导率的自然对数，其数量级为10¹，更新量数据量级为10⁰.基于两种模型更新量数量级的差异，在设置地震慢速与电阻率模型更新量之间平衡的权重参数时，需要使α_s比α_r高4个数量级.本文设置的权重参数α_s为10⁰，α_r为10^-4，这样就改善了地震走时与电阻率数据联合反演时两种参数的数量级及变化范围不同的问题，为联合反演权重选择提供了参考依据.β_t从10^-10至10¹⁰以间隔10倍关系分别计算地震走时、电阻率观测数据拟合残差与交叉梯度值，根据L曲线选择合适的结构约束权重β_t(Aster et al., 2011).本文对于前5次迭代将β_t设置为10²，第6次至第20次迭代将β_t设置为10³.

综合起来，两种不同成像方法的联合反演参数选择分为两步：

(1) 首先根据两种模型参数的更新量数量级，确定两种模型之间合适的权重参数使两种模型均起到结构约束的作用.

(2) 根据L曲线等权重选择方法确定结构约束与数据拟合之间的平衡.

4 合成数据测试

为了验证本文提出的基于模型更新量的联合反演流程，我们利用三维跨孔地震走时与三维跨孔直流电阻率联合反演进行数值试验.由于本文提出的联合反演方法不需要修改单个反演程序，因此，可以非常便捷地实现两者的联合反演.

图 3为地震走时与直流电阻率跨孔观测装置布置示意图.在该算例中，我们设计了符合实际工程施工条件的4个钻孔式的跨孔观测装置.假设模型地表一个角为(0, 0) 点，则4个钻孔的坐标分别为：Bore1(3, 3)、Bore2(14, 3)、Bore3(14, 16)、Bore4(3, 16)，坐标的单位为m，孔深均为12 m.

在地震走时跨孔装置中，每个孔中设置11个激发点、11个接收点，激发点和接收点间距均为1 m；地表布设56个接收点，间距为1 m×1 m.数据采集时，每个孔中的激发对应的接收点为除去激发点所在钻孔之外的所有接收点.

在电阻率跨孔装置中，每个孔中放置10个电极，电极距为1 m；地表布设56个电极，间距为1 m×1 m.数据采集时，地表和孔中所有的电极任取两个作为供电电极，其余电极全部采集观测电位.

我们首先通过均匀速度模型2000 m·s^-1和均匀电阻率模型100 Ωm来验证三维地震走时场与三维电位场的数值计算精度.比较数值计算结果与解析结果发现，两者误差均达到0.01%的误差水平.由于实际反演迭代终止条件一般选择数据拟合为1%的拟合水平，因此，本文采用的地震走时与直流电阻率正演算法能够满足反演精度要求.图 4为均匀背景中存在一个高速异常体的理论速度模型与对应的三维地震走时场，其中背景速度为2000 m·s^-1，高速异常体速度为3000 m·s^-1，异常体大小为3 m×3 m×2 m；图 5为均匀背景电阻率中存在一个高阻异常体的理论电阻率模型与对应的三维电位场，其中背景电阻率为100 Ωm，高阻体电阻率为1000 Ωm，异常体大小为3 m×3 m×2 m.

基于图 4、图 5中所示的理论速度模型和理论电阻率模型，通过正演计算得到了理论地震走时和电位数据，然后均加入了3%的高斯噪声构建了用于模型测试的数据.地震走时和电阻率成像所用的初始速度模型和初始电阻率模型均为均匀模型，初始速度为2000 m·s^-1，初始电阻率为100 Ωm.单独地震走时、单独电阻率成像以及联合反演得到的速度模型和电阻率模型的三维立体图、等深度切片及沿Y方向的侧向切片分别显示在图 6、图 7和图 8.从图中可以看出，地震走时成像与电阻率成像均能较好地分辨出异常的水平中心位置及形态；但在垂向上，地震走时成像具有较高的分辨率，而电阻率成像得到的异常体向地表延伸，无法较好地刻画异常体的纵向形态.相比较来说，通过联合反演提高了电阻率成像对异常体的纵向分辨率，同时由于电阻率成像产生干扰异常较少，可以通过交叉梯度的约束减少地震走时成像的拖尾异常.

图 7和图 8显示的联合反演结果表明联合反演比单独反演具有较大的优势.在该测试中，两个模型中对应的异常为正比例关系，即高速对应高阻.为了进一步研究基于交叉梯度结构约束联合反演的适用范围，本文设计了异常结构属性呈反比例对应关系的速度和电阻率模型，其中速度模型与图 4中速度模型一致，而电阻率模型为在均匀背景电阻率(100 Ωm)中存在一个低阻异常体，低阻异常体的电阻率值为10 Ωm，大小为3 m×3 m×2 m.图 9和图 10为反比例模型对应的单独反演和联合反演得到的速度模型和电阻率模型在深度上的切片和沿Y方向的侧向切片.从图中可以看出，与正比例异常模型一样，联合反演对于反比例异常模型也能得到较好的反演结果.

在实际地质工程勘探时，若勘探区内有钻孔时，在钻孔中测井或取样可以测得垂向岩层的速度和电阻率值.通过孔中已知的速度和电阻率值，可以构建一维层状参考速度和电阻率模型.在进行联合反演时，可以利用钻孔中已知的垂向上的界面信息来约束三维速度和电阻率模型的反演.式(7) 中阻尼权重系数α_s和α_r分别替换为diag(∂m_s^priori/∂z)和diag(∂m_r^priori/∂z)得到式(9)，其中m_s、m_r分别为上一次迭代得到的速度模型和电阻率模型.m_s^priori和m_r^priori分别根据已知速度和电阻率建立的层状参考速度和电阻率模型，∂m_s^priori/∂z和∂m_r^priori/∂z分别为先验速度模型和电阻率模型在垂向上的梯度，因为先验模型为层状模型，速度和电阻率只在模型垂向上变化，所以∂m_s^priori/∂z和∂m_r^priori/∂z计算得到的垂向梯度就保留了已知参考模型的上下边界信息，然后把∂m_s^priori/∂z和∂m_r^priori/∂z分别对角化并与单位矩阵I相乘，在通过式(9) 求解新的速度更新量Δm_s和电阻率更新量Δm_r时，就保留了上一次反演得到的速度更新量Δm_s⁰和电阻率更新量Δm_r⁰中的边界信息，如式(9) 所示：

(9)

式中：diag(∂m_s^priori/∂z)为先验速度模型m_s^priori纵向梯度对角化，diag(∂m_r^priori/∂z)为先验电阻率模型m_r^priori纵向梯度对角化.

图 11和图 12分别为加入垂向层状速度与电阻率先验模型约束的三维跨孔地震走时与直流三维电阻率联合反演模型的等深度切片与侧向剖面图.该算例采用的观测装置与图 3所示的一致，在均匀速度(2000 m·s^-1)、电阻率(100 Ωm)模型中设置了两个异常体，分别为高速高阻体(2500 m·s^-1, 1000 Ωm)和低速低阻体(1500 m·s^-1, 10 Ωm).异常体大小为3 m×3 m×3 m，距离地表 3 m.速度和电阻率反演初始模型分别为2000 m·s^-1和100 Ωm.从加入先验信息的联合反演结果可以看出，相比较于地震走时数据与电位数据单独反演的结果，联合反演得到的速度与电阻率模型均具有更好的形态及更准确的异常体中心成像，与真实模型相比速度异常与电阻率异常均得到了较好的数值恢复.

5 讨论和结论

地球物理联合反演一般基于多方法之间的互相约束实现，包括基于岩石物理属性联系与结构一致两种约束方式的联合反演.本文在实现三维地震走时与直流电阻率联合反演时采用了基于交叉梯度结构约束的方法.但是传统的基于交叉梯度结构约束的联合反演可能存在如下问题：

(1) 构建联合反演矩阵时，反演系数矩阵较为庞大，需要较大的计算及存储成本.

(2) 不同地球物理方法对岩石属性灵敏度不同，造成不同方法的反演系数矩阵数量级差距较大，致使联合反演系数矩阵的求解存在不稳定性.

(3) 为了构建联合反演矩阵，需要对单一反演程序进行修改，给联合反演程序实现带来巨大困难.

针对上述联合反演存在的问题，本文提出了基于单个数据单独反演得到的模型更新量与交叉梯度结构约束相结合的联合反演策略，实现数据拟合和结构一致之间的拟合.该反演策略具有以下优势：

(1) 不需要对单个数据的反演程序进行改动，只需要给出其迭代一次的模型更新量即可，极大地简化了联合反演的实现难度.采用本文提出的新的联合反演流程，可以方便地实现多属性联合反演.例如实现直流电测深与多道面波的联合反演，可提高直流电测深的纵向分辨率及多道面波的横向分辨率.

(2) 在单个数据反演中可以有效地选择数据拟合权重与模型平滑权重，而在联合反演中只需要平衡结构约束与数据拟合即可，解决了传统的联合反演系统中权重难以选择的问题(Bennington et al., 2015; Zhou et al., 2015).

(3) 与Um等(2014)的反演策略相比较，本文提出的联合反演策略可以很容易地实现结构一致和数据拟合之间的平衡.

通过三维跨孔地震走时与直流电阻率联合反演对本文提出的联合反演流程进行了测试，结果证明了该联合反演策略的正确性与有效性.通过电阻率结构约束有效减少了速度反演结果中的拖尾等干扰异常，通过速度结构的约束有效提高了电阻率反演结果的分辨率，进而总体提高了地震走时与电阻率的成像效果.

致谢

感谢评审专家对论文提出的宝贵修改意见.