地球物理学报  2016, Vol. 59 Issue (11): 4310-4322   PDF    
基于交叉梯度交替结构约束的二维地震走时与全通道直流电阻率联合反演
高级1,2 , 张海江1     
1. 中国科学技术大学地球和空间科学学院 地震与地球内部物理实验室, 合肥 230026;
2. 安徽万泰地球物理技术有限公司, 合肥 230026
摘要: 在利用不同的地球物理勘探方法对地下复杂介质成像时,因观测系统的非完备性及数据本身对某些岩石物性的不敏感性,单独成像的结果存在较大的不确定性和不一致性.对于地震体波走时成像与直流电阻率成像,均面临着成像阴影区问题.对于地震走时成像,地震射线对低速区域覆盖较差形成阴影区,造成低速区域分辨率降低.对于电阻率成像,电场线在高阻区域分布较少,造成高阻区域分辨率较低.为了提高地下介质成像的精度,Gallado和Meju(2003)提出了基于交叉梯度结构约束的联合地球物理成像方法.在要求不同的物性模型拟合各自对应的数据同时,模型之间的结构要求一致,即交叉梯度趋于零.为了更有效地实现基于交叉梯度的结构约束,我们提出了一种新的交替结构约束的联合反演流程,即交替反演不同的数据而且在反演一种数据时要求对应的模型与另一个模型结构一致.新的算法能够更容易地把单独的反演系统耦合在一起,而且也更容易建立结构约束和数据拟合之间的平衡.基于新的联合反演流程,我们测试了基于交叉梯度结构约束的二维跨孔地震走时和直流电阻率联合成像.合成数据测试表明,我们提出的交替结构约束流程能够很好地实现基于交叉梯度结构约束的联合成像.与单独成像结果相比,地震走时和全通道电阻率联合成像更可靠地确定了速度和电阻率异常.
关键词: 体波走时成像      电阻率成像      结构相似性      交叉梯度      联合反演     
Two-dimensional joint inversion of seismic velocity and electrical resistivity using seismic travel times and full channel electrical measurements based on alternating cross-gradient structural constraint
GAO Ji1,2, ZHANG Hai-Jiang1     
1. Laboratory of Seismology and Physics of Earth's Interior, School of Earth and Space Sciences, University of Science and Technology of China, Hefei 230026, China;
2. Anhui Wantai Geophysical Technology Co. Ltd., Hefei 230026, China
Abstract: When imaging complex near surface structures using different geophysical exploration methods, because of the incompleteness of the observation system and insensitivity of geophysical data to some geophysical properties of media, there exist relatively large uncertainties and significant inconsistency in imaging results from separate inversions. For seismic travel time tomography and electrical resistivity tomography, both methods are subject to the issue of imaging shadow zones. For seismic travel time tomography, seismic rays tend to travel through high velocity zones, causing poor ray path coverage and thus low resolution of low velocity zones. For electrical resistivity tomography, the distribution of electric field lines is relatively sparse in high resistivity zones, leading to lower resolution. In order to reliably image underground structure, Gallardo and Meju (2003) proposed a joint geophysical inversion method based on cross-gradient structure constraint, which requires separate models fitting the data and at the same time structurally consistent by having approximately zero cross gradients. To more effectively realize the cross-gradient based structure constraint, we propose a new joint inversion scheme that alternatively inverts one type of data but structurally constrained by another model. With this new joint inversion scheme, we test two-dimensional cross-hole seismic travel time tomography and full channel resistivity tomography. Synthetic results show that the new scheme can efficiently realize the cross-gradient based joint inversion. Compared to separate inversions, the joint inversion can more reliably determine seismic velocity and electrical resistivity anomalies..
Key words: Seismic travel time tomography      Electrical resistivity tomography      Structure similarity      Cross-gradient      Joint inversion     
1 引言

不同的地球物理勘探方法对地下介质的成像各有其优缺点,而且都存在某种程度的不确定性和多解性(胡祥云等,2006).引起地球物理勘探不确定性和多解性的原因包括:有限的观测数据、由于环境干扰和仪器精度存在的观测误差、地球物理正演方程和真实物理过程存在的偏差、不同地质模型观测的地球物理场相同情况的存在、反演算法存在的偏差等因素.因此为了更可靠地刻画地下介质,多种地球物理数据之间的联合反演解释就成为一种有效的途径(Gallardo and Meju,20032004; Gallardo et al.,2012; Abubakar et al.,2012).

地下介质具有不同的物理性质,如密度、电阻率、极化率和磁性.一类的联合反演是基于多属性参数的耦合即建立地下介质不同物理性质之间的经验关系(如Berge et al.,2000; 于鹏等,2006).例如因为速度和密度之间存在较好的经验关系,Maceira和Ammon(2009)发展了面波和重力的联合成像方法.但是对于某些地球物理属性,例如地震速度和电阻率,很难建立起二者之间普遍适用的经验关系.在这种情况下,发展了基于结构一致性的另一类联合反演方法(如,Lines and Lines,1988; Haber and Oldenburg,1997; Musil et al.,2003; Gallardo-Delgado et al.,2003; Hu et al.,2009).上述两类联合反演方法均存在一定的问题(王俊,2014):第一类方法需要预先知道不同物性参数基于岩石物理学的联系关系才能实现联合反演,而多数情况下两种属性之间的联系难以确定;对于第二类方法,如何解决模型结构的融合一直是制约该方法发展的重大障碍.Gallardo和Meju(20032004)提出了基于模型交叉梯度的结构约束实现直流电法与地震走时的联合反演,取得了较好的效果.该方法利用梯度场来刻画模型的结构,如果两个模型的梯度一致即结构也一致.在反演的过程中,一方面要求各个模型拟合数据,同时通过改变模型让二者之间的交叉梯度趋于零实现结构的一致.该联合反演方法无需明确模型之间的物性关系,目前得到了地球物理学界广泛的关注.例如,Abubakar等(2012)进一步结合岩石物理性质及结构相似性实现了大地电磁与面波的联合反演,Bennington等(2015)发展了基于归一化交叉梯度约束的三维天然地震和二维大地电磁的联合反演算法,李兆祥等(2015)通过交叉梯度约束实现了电阻率与极化率的联合反演,并取得了较好的效果.

对于不同地球物理数据的反演,模型的离散化和数据对模型的敏感度是不一样的,因此把对应于不同数据的反演系统与交叉梯度结构约束完全耦合成一个大的联合反演系统是一个复杂的具有挑战性的工程(Gallardo and Meju,2007).为了降低实现基于交叉梯度结构约束的难度,Gallardo和Meju(2007)提出了分别进行数据拟合和交叉梯度结构约束的联合反演流程,通过迭代交替进行数据反演和交叉梯度反演实现联合反演.虽然该流程降低了实现联合反演的难度,但是由于数据拟合与结构约束完全分开,因此迭代交替反演时很难平衡数据拟合和结构约束的程度.例如在一次迭代中,利用不同的数据首先对各自对应的模型进行了更新,但再用交叉梯度结构约束对两个模型进行改变使二者结构一致时,可能会导致对数据的拟合在很大程度上变差;另外利用数据反演对模型进行更新之后,二者的结构一致性可能会得到破坏.所以利用Gallardo和Meju(2007)提出的联合反演流程,因为数据拟合和结构约束是分开进行的,因此如何平衡二者对模型的更新是一个挑战.

针对于此问题,本文提出了一种新的基于交叉梯度结构约束的联合反演流程,避免了Gallardo和Meju(2007)提出的数据拟合和结构约束完全分开的缺点,同时还避免了把对应于不同数据的反演系统放在一起的困难.在我们新提出的反演流程中,对一种数据进行反演的时候,同时要求对应的模型与另外一个模型结构一致,两种数据反演交替进行直至收敛.结构的约束还是通过交叉梯度实现,但不同的是只改变一个模型而另一个模型不变使二者的结构一致.本文首先研究了全通道2.5D跨孔直流电阻率层析成像(潘纪顺等,2010;潘克家等,2012; 雷旭友等,2009; Loke and Barker,1996; Loke and Dahlin,2002; Loke,2003),并比较了全通道接收与传统四极法采集反演结果的不同.其次介绍了基于MSFM(Multi-stencils Fast Marching Methods)射线追踪方法的跨孔地震走时反演.然后提出了新的基于交叉梯度交替结构约束的联合直流电阻率和地震走时反演流程.最后基于合成数据,验证了本文提出的联合反演流程的有效性.

2 全通道2.5D跨孔直流电阻率层析成像

全通道直流电阻率法为除去2个供电电极外其余电极均作为采集电极采集数据并用于反演的一种电阻率成像方法.与传统四极采集方式相比,全通道直流电阻率法具有更好的采集方位覆盖与更大的数据采集量.例如若采用32个电极采集,全通道采集得到的电位值数为29270(C322×30=29760),而对应的温纳装置采集数据个数为 155($\underset{n=1}{\overset{10}{\mathop{\sum }}}\,$(32-3n)=155), 即全通道采集数据量为温纳装置的192倍.全通道采集方式包含了常规四极法中的温纳、偶极、二极等装置,同时具有各种常规四极方法的优点.例如具有温纳装置的抗噪能量强、垂直分辨率高和水平分辨率低的优点;具有偶极装置抗噪能力弱、水平分辨率高、垂直分辨率低的优点;具有二极装置抗干扰能力弱、探测深度大的优点(冯锐等,2004).因为全通道采集方式可以对地下结构的方位约束更好并且可获取更为丰富的信息量,因此能获得更可靠的地质模型.

本文在实现电阻率正反演时,首先采用有限差分方法求解2.5D静电场方程(范翠松等,2012; 刘斌等,2012),求取观测位置处的电位值,并基于伴随矩阵(Adjoint method)方法计算非线性灵敏度矩阵(Pidlisecky et al.,2007;吴小平和徐果明,2000; Mcgillivray and Oldenburg,1990),最后利用非线性牛顿共轭梯度反演方法(刘云鹤和殷长春,2013)实现全通道电阻率层析成像(吴小平等,2015).

二维跨孔直流电阻率层析成像装置形式为在两个钻孔内布置供电与接收电极.其供电形式可以采用单极、偶极、三极、四极等装置形式,跑极方式与地表采集相同,也可以采用全通道接收形式.本文在实现直流电阻率与地震波走时联合反演时采用全通道采集方式.

图 1为跨孔电阻率装置及正反演结果图.其中图 1a为电极布置方式,分别在两个钻孔中布置32个供电与采集电极,孔距15 m;图 1b为电阻率模型,其中背景电阻率为100 Ωm,低高阻电阻率异常值分别为10 Ωm和1000 Ωm,异常大小为4 m×4 m;图 1c为两个电极供电时观测区域2D电场分布等值线图,从中可以看出跨孔电阻率电场在空间分布较为均匀,电场线穿过孔间观测区域,具有对孔间异常进行成像的基础;图 1d为电阻率反演结果,从反演结果图中可以看出高阻及低阻异常均具有较好的成像及较少的假异常,但受到电法勘探体积效应的影响反演结果分辨率相对较低.

图 1 二维跨孔全通道电阻率成像 (a)电极分布;(b)真实电阻率模型;(c)二维电位场,V为电位,I为电流;(d)电阻率反演结果. Fig. 1 Two-dimensional cross borehole full channel resistivity tomography (a)Distribution of electrodes;(b)Synthetic resistivity model;(c)2D electric field,V is electric potential and I is electric current; (d)Inverted resistivity model.
3 地震体波初至走时成像

地震体波走时层析成像是利用初至波到时构建地下速度结构的一种重要的方法(Zelt and Smith,1992),该方法可以利用主动源、被动源或同时使用两种震源(Di Stefano and Chiarabba 2002;Ritter et al.,2001).对于体波走时成像,需要追踪震源和检波器之间的射线和计算对应的走时.射线追踪方法主要包括基于程函方程(Vidale,19881990)及基于射线方程(Zelt and Smith,1992; 张风雪等,2010)方法.本文采用不等间距节点法对模型区域网格进行剖分(Sambridge,1990; Sambridge and Rawlinson,2004),走时计算采用有限差分法(Vidale,19881990)求解程函方程的MSFM(Multi-Stencils Fast Marching Methods)方法(Osher and Fedkiw,2003; Popovici and Sethian,2002).反演时为提高反演结果的稳定性,采用二阶吉洪诺夫正则化(Tikhonov,1977)约束(Calvetti et al.,2000),并采用阻尼最小二乘LSQR法迭代求解(Paige and Saunders,1982).

图 2为跨孔地震波走时观测系统及正反演结果图,对应于图 1所示的跨孔电阻率观测系统.其中图 2a为观测系统图,检波器和震源分别放置在两个钻孔中,每个孔中放置16炮30个检波器,炮距2 m、检波距1 m,在一个孔中放炮时另一个孔中30个检波器同时接收;图 2b为速度模型,其中背景速度为2000 m·s-1、异常速度值分别为2500 m·s-1和1500 m·s-1、异常大小为4 m×4 m,其中低速对应低阻而高速对应高阻异常;图 2c为基于MSFM求解的波前走时场等值线图;图 2d为走时反演结果,从反演结果图中可以看出走时成像能确定高低速异常体所在位置,与电阻率成像结果(图 1d)比较,具有较高的纵横向分辨率,但在异常周围产生较多的假异常,产生假异常的主要原因是模型测试采用了符合实际工程施工的跨孔观测系统,即炮点和接收点位于两个孔中,因此观测系统覆盖角度具有一定的局限性,导致反演的速度模型在x方向的分辨率低于z方向分辨率.

图 2 二维跨孔地震走时层析成像 (a)两个孔中震源和检波器的分布;(b)真实速度模型;(c)单炮地震波走时场等值线分布;(d)地震波走时反演速度模型. Fig. 2 Two-dimensional cross borehole seismic travel time tomography (a)Distribution of sources and receivers in two boreholes;(b)Real velocity model; (c)Travel time field using the MSFM ray tracing method;(d)Inversion velocity model.
4 基于交叉梯度约束的联合反演方法

基于交叉梯度结构约束的联合反演由Gallardo和Meju(2003)首次提出,后续被广泛应用于地震数据和大地电磁数据、地震与直流电法数据的联合反演(Gallardo and Meju,200320042007; Gallardo et al.,20052012).国内基于交叉梯度约束的联合反演研究也逐步展开(彭淼,2012;周丽芬,2012王俊,2014).对于三维地震速度模型ms和电阻率模型mr,二者之间的交叉梯度定义为速度梯度和电阻率梯度的叉积:

(1)

在三维模型情况下,垂直于x、y、z平面方向上的交叉梯度分量分别为 txtytz. 其展开形式分别为:txty=-

在二维模型情况下,模型沿y轴方向变化为0,模型之间沿xz方向变化的交叉梯度只存在y方向分量:

(2)

图 3显示了两个不同二维模型如何计算交叉梯度.从图中可以看出,当两个模型梯度方向相同或相反时,即结构一致时交叉梯度值为零,反之不为零.两个模型交叉梯度的这个性质为利用交叉梯度约束反演奠定了基础,即在反演时要求速度与电阻率模型交叉梯度值趋于零,即模型之间具有一致的结构.因此基于交叉梯度结构约束的联合反演要求不同的模型在拟合对应的数据同时它们之间的结构要一致.另外考虑到联合反演系统存在一定的病态性,因此对模型加入一定的平滑约束,据此建立的联合反演目标函数为(Gallardo and Meju,2003Bennington et al.,2015):

(3)

图 3 两个模型交叉梯度的测试 (a)模型1;(b)模型1的梯度场;(c)模型2;(d)模型2的梯度场;(e)模型1与模型2的交叉梯度. Fig. 3 Test of cross gradient for two models (a)Model 1;(b)Gradient field of model 1;(c)Model 2;(d)Gradient field of model 2; (e)Cross gradients between model 1 and model 2.

式中: fs(ms) 为根据地震速度模型 ms 计算的理论地震走时数据, ds 为观测走时数据, μs 为地震走时拟合权重, Cs 为观测走时数据误差的协方差矩阵, Ls 为速度模型平滑因子, λs 为速度模型平滑权重; fr(mr) 为根据电阻模型 mr 计算的理论电位数据, dr 为观测电位数据, μr 为电位数据拟合权重, Cr 为观测电位数据误差的协方差矩阵, Lr 为电阻率模型平滑因子, λr 为电阻率模型平滑权重. t(msmr) 为速度模型与电阻率模型之间的交叉梯度, β 为交叉梯度约束权重.

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

(4)

式中$\frac{\partial t}{\partial {{m}_{s}}}$为交叉梯度对速度的偏导数,$\frac{\partial t}{\partial {{m}_{s}}}$为交叉梯度对电阻率的偏导数, t(ms0mr0) 为参考模型 ms0mr0 之间的交叉梯度.假设地震走时数据的观测误差互不相关,数据点的误差为 (σ1s,σ2s,σ3s,…,σnss), 其中 ns 为走时观测数据个数;同时假设电位数据的观测误差也互不相关,数据点的误差为 (σ1r,σ2r,σ3r,…,σnrr), 其中 nr 为电位数据观测个数.对于式(3)中列出的联合反演目标函数,可以通过线性化并通过阻尼最小二乘反演矩阵系统进行求解(Aster et al.,2005):

(5)

式中: Gs 为地震走时数据对地震速度模型 ms 的灵敏度矩阵, Ws=diag(1/σ1s,1/σ2s,…,1/σnss),Gr 为电位数据对电阻率模型 mr 的灵敏度矩阵, Wr=diag(1/σ1r,1/σ2r,…,1/σnrr),Δms 为速度模型的更新量, Δmr 为电阻率模型的更新量, Δds 为观测走时数据与计算走时数据之间的残差, Δdr 为观测电位数据与计算电位数据之间的残差, εs 为速度模型阻尼权重, εr 为电阻率模型阻尼权重, I 为单位矩阵, t0 为当前速度模型和电阻率模型的交叉梯度.

Bennington等(2015)进行三维地震走时和二维大地电磁联合反演时,就是基于公式(5)的反演矩阵系统,把地震走时反演和大地电磁反演以及交叉梯度结构约束真正融合在一起.该种方式虽然在数学上具有严格性,但在实际编程实现时具有挑战性,因为把复杂的单个反演系统融合成为一个大的反演系统容易出现问题.Gallardo和Meju(2007)提出了一种解决方案,避免了融合不同反演系统存在的挑战性,如图 4所示的反演流程图.例如对于地震走时和电阻率的联合反演,在Gallardo和Meju(2007)提出的反演流程中,首先基于独立的地震走时成像和电阻率成像系统分别利用地震走时数据和电位数据进行反演得到新的地震速度模型和电阻率模型,然后利用交叉梯度结构约束对两个模型进行更新得到新的模型,接着再重复上述的过程直至拟合数据和二个模型结构一致.这种联合反演方式把结构约束与数据拟合孤立地进行,虽然可以直接利用现有的反演系统进行联合反演,但是造成结构约束权重难以选择,即结构约束太强会导致数据拟合发散,若结构约束太弱会导致其失去作用,难以选取合适的交叉梯度结构约束权重因子.

图 4 Galldado和Meju(2007)提出的基于交叉梯度结构约束的联合反演流程图 Fig. 4 Flow chart of joint inversion based on cross-gradient structure constraint proposed by Gallardo and Meju(2007)

本文为克服此缺陷,提出一种新的交叉梯度结构约束与数据拟合进行联合反演的方式,避免了交叉梯度结构约束及两种数据拟合之间权重选择困难的问题.式(5)所示的联合反演矩阵系统可以分开为两个交替进行的子反演系统:

(6)

(7)

公式(6)所列出的反演系统表示的是基于电阻率模型结构约束的地震走时成像,在这里交叉梯度结构约束只用来改变地震速度模型的结构使其与电阻率模型结构一致,同时还要拟合地震走时数据,二者之间的平衡通过权重系数 βs 来控制.同样,式(7)所列出的反演系统表示的是基于地震速度模型结构约束的电阻率成像,在这里交叉梯度结构约束只改变电阻率结构使其与地震速度模型结构一致,二者之间的平衡通过权重系数 βr 来控制.与公式(5)给出的大的联合反演矩阵系统相比较,我们给出的基于某个模型结构对另外一个模型结构进行约束的反演算法避免了不同反演系统反演参数和敏感矩阵数量级不一致的情况.与Gallardo和Meju(2007)的方法相比较,基于式(6)和(7)的反演同时考虑了结构约束和数据拟合之间的一个平衡.具体实现时只需要把$\frac{\partial t}{\partial {{m}_{s}}}、\frac{\partial t}{\partial {{m}_{r}}}$通过βsβr 调整为和数据灵敏度矩阵相同数量级即可.

图 5给出了本文提出的基于交叉梯度结构约束的地震走时和电阻率联合反演流程:

图 5 本文提出的基于交叉梯度交替结构约束的联合反演流程图 Fig. 5 Flow chart of joint inversion with alternating cross-gradient based structure constraint proposed in this study

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

(2) 基于公式(6)和(7)分别进行基于交叉梯度结构约束的地震走时和电阻率成像.

(3) 更新速度与电阻率模型.

(4) 计算两个模型的交叉梯度以及地震走时数据和电位数据的残差.

(5) 判断数据拟合和交叉梯度值是否满足迭代终止条件,若满足直接输出结果,若不满足跳回步骤(2).一般迭代终止条件包括,迭代次数、每次迭代残差下降的梯度、残差剩余量占初始残差百分比、模型更新量等.本文迭代终止及收敛条件为电阻率或地震走时残差下降量低于某一百分比时迭代终止(如迭代残差下降量小于1%).

5 合成数据测试

我们利用图 1图 2所建立的地震速度与电阻率模型以及对应的观测系统,对本文提出的基于交叉梯度的交替结构约束联合反演算法进行了测试.对于公式(6)和(7)中出现的平滑系数λ,阻尼因子ε和数据相对权重β可以采用L-curve方式进行选择(Hansen and O′Leary,1993; Aster et al.,2005).例如对于基于电阻率模型结构约束的地震走时速度成像,通过L-curve分析,得到对应于曲线拐点的最优权重因子为4×106,即在数据拟合和结构一致性方面达到了平衡,如图 6所示.

图 6 交叉梯度与拟合误差L-curve曲线 Fig. 6 Cross-gradient and fitting error L-curve

图 7显示了地震走时及电阻率联合反演得到的速度模型和电阻率模型.与单独反演得到的速度模型比较,联合反演得到的速度模型很大程度上消除了低速和高速异常周围出现的由于拖尾效应引起的假异常.对于电阻率模型,联合反演得出的异常中心位置更靠近真实模型,尤其是低阻异常,而且反演结果数值更接近真实模型数值.

图 7 单独反演和联合反演得出的速度和电阻率比较 (a)真实速度模型;(b)真实电阻率模型;(c)单独走时反演得到的速度模型;(d)单独电阻率反演得到的电阻率模型;(e)速度和电阻率联合反演得到的速度模型;(f)速度和电阻率联合反演得到的电阻率模型. Fig. 7 Comparison of seismic velocity and electrical resistivity from separate and joint inversions (a)True seismic velocity model;(b)True electrical resistivity model;(c)Seismic velocity model from separate travel time inversion;(d)Resistivity model from separate resistivity tomography;(e)Seismic velocity model from joint inversion;(f)Resistivity model from joint inversion.

图 8为地震走时及电阻率单独反演及联合反演走时残差和电位残差随着迭代次数的变化曲线.对于联合反演和单独反演来说,电位残差随着迭代次数变化的趋势基本一致,但最终联合反演得出的电阻率模型能更好地拟合数据,也对应着恢复更好的电阻率模型.相比较而言,联合反演得到的速度模型对地震走时数据的拟合比单独反演变差,最终模型对应的走时残差变大.从速度模型(图 7c)可以看出,由于地震走时反演系统的病态性,单独地震走时反演过度地拟合了数据,导致在真正的异常周围出现很多的假象,但利用联合反演,速度模型的结构得到了电阻率模型的约束,因此速度模型不能过于“自由”的变化.对于本理论模型测试,可以看出电阻率结构起了相对主导的作用.

图 8 单独反演与联合反演地震走时和电位数据残差随着迭代次数变化的曲线 (a)地震走时数据;(b)电位数据. Fig. 8 Variation curves of residuals for seismic travel times and electrical potentials along with iteration times for separate and joint inversions (a)Seismic travel time data;(b)Electric potential data.

图 9比较了单独反演与联合反演速度和电阻率模型的交叉梯度.从图中可以看出,单独反演模型的交叉梯度值分布范围为-5.0~+5.0,联合反演后模型交叉梯度值分布范围为-0.5~+0.5,即联合反演后模型之间交叉梯度值下降到约为单独反演的10%(图 9a9b),同时从整个模型空间看,联合反演之后,交叉梯度值都较小,而且数值分布较均匀(图 9c9d),说明利用交叉梯度的结构约束使得模型之间具有了更好的结构一致性,符合地质模型结构唯一的特性,因此基于结构一致性的联合反演提高了反演结果的合理性.

图 9 单独与联合反演速度和电阻率模型交叉梯度的比较 (a)单独反演模型交叉梯度在每个深度上在水平方向的变化;(b)与(a)类似但是为联合反演情况; (c)单独反演模型交叉梯度在模型空间的分布;(d)联合反演模型交叉梯度在模型空间的分布. Fig. 9 Comparison of cross-gradient variations for seismic velocity and resistivity models from separate and joint inversions (a)Cross-gradient variations along the horizontal direction at each depth for models from separate inversions;(b)Same as(a)but for joint inversion; (c)Distribution of cross-gradients in model space for models from separate inversions;(d)Same as(c)but for joint inversion.

在某些地质条件下,地震速度和电阻率二者之间可能存在结构不完全一致的情况.例如当岩层空隙局部充水时,对岩层电阻率有较大影响,造成对应区域的电阻率降低,但可能对地震波传播速度的影响较小,造成两种模型具有不同异常结构.为此我们也测试了算法在地震速度和电阻率具有不一致结构情况下的联合反演效果,同时为了检验本文算法在观测数据存在噪声时反演的稳定性,我们在走时观测数据与电位观测数据中分别加入了3%的随机误差.

图 10为地震速度模型与电阻率模型具有不同结构情况下的联合反演结果.地震速度和电阻率模型中的高速和高阻异常具有结构一致性,但是低速和低阻异常的结构不一致.从图中可以看出,与单独反演得到的速度模型比较,联合反演得出的速度模型很大程度上消除了高速异常周围的假异常,同时速度反演结果在对应低电阻率异常位置没有产生假异常.其原因为当一个模型的梯度为零,不论另一个模型的结构如何,交叉梯度也为零.在低电阻率异常区域,对应的速度模型变化非常小接近常数,梯度变化接近于零.因此在进行交叉梯度约束反演时,对应的交叉梯度值亦趋近于零.所以低电阻率异常结构不会对速度结构起到约束的作用而导致假异常的出现.同样速度模型中的低速异常结构也没有在电阻率模型对应区域带来假电阻异常.

图 10 速度模型与电阻率模型结构不一致情况下单独反演和联合反演结果比较 (a)真实速度模型;(b)真实电阻率模型;(c)单独走时反演得到的速度模型;(d)单独电阻率反演得到的电阻率模型; (e)联合反演得到的速度模型;(f)联合反演得到的电阻率模型. Fig. 10 Comparison of seismic velocity and electrical resistivity models from separate and joint inversions when two models have different structures (a)True seismic velocity model;(b)True electrical resistivity model;(c)Seismic velocity model from separate travel time inversion;(d)Resistivity model from separate resistivity tomography;(e)Seismic velocity model from joint inversion;(f)Resistivity model from joint inversion.

在速度模型与电阻率模型具有相同结构的区域,速度模型与电阻率模型之间得到较好的相互结构约束.高速异常得到高电阻率异常的x方向上的结构约束,使高速异常区域具有更好的形态恢复.对于电阻率模型来说,高电阻率异常区域得到速度模型中高速异常在z方向结构上的约束,因此联合反演提高了电阻率模型高电阻率异常在z方向上的分辨率.

6 讨论和结论

多物理属性的联合反演,目标是利用地下介质的多种岩石物理性质来减小模型的多解性(Vozoff and Jupp,1975).不同的地球物理方法对不同岩石物理性质具有不同的分辨率,如电阻率法对低阻异常区敏感、速度成像对高速异常区敏感,同时高阻体一般为高速,低速体一般为低阻,因此电阻率成像与速度成像对介质敏感性是互补的.另外不同地球物理数据采集时受不同来源噪声的影响是不同的,如电阻率法对电性干扰较为敏感,但地震数据受环境噪声影响较大.基于上面两点,电阻率成像和地震走时成像的联合反演系统可以更好地刻画地下介质.

针对完全耦合在一起的联合反演系统存在多属性数据数量级不同以及不同反演系统耦合在一起存在的技术实现问题,我们提出基于交叉梯度交替结构约束的联合反演算法.新的反演方法避免了不同属性数据直接耦合在一起的挑战而且还可以兼顾数据拟合和结构约束之间的平衡.基于一个二维的地震速度和电阻率模型,我们利用新提出的联合反演流程进行了测试,表明电阻率成像与速度成像是一种较为有效互补的成像方式.通过基于交叉梯度约束的结构一致性联合成像,较单独反演能够更好地恢复速度和电阻率模型.合成模型测试还表明,因为电阻率反演可以得到较稳定平滑模型,因此有助于消除地震走时成像中出现的次生干扰异常.同时我们还测试了地震速度和电阻率模型结构不一致的情况,结果表明基于交叉梯度结构约束的联合反演只对两个模型具有一致结构的区域产生结构约束,而没有在结构不一致的地方产生假的异常.这主要是因为在一个模型的梯度为零的情况下不论另外一个模型的异常结构形态如何变化,两个模型的交叉梯度也是零.

在联合反演迭代过程中,要适当控制单一模型的收敛速度,即如果一个模型收敛速度较快使异常周围模型梯度为零,两个模型的交叉梯度在异常周围也为零,此时交叉梯度无法起到对另一模型的结构进行结构约束的作用.另外因为地震走时和电阻率成像具有不同的收敛性和分辨率,可以采用在联合反演迭代第一阶段选择以地震速度模型结构约束为主,提高电阻率反演的分辨率,随着迭代可以逐渐增大电阻率模型结构约束的权重,提高地震走时反演的稳定性.

参考文献
Aster R C, Borchers B, Thurber C H. Parameter Estimation and Inverse Problems. San Diego: Academic Press, 2005 .
Bennington N L, Zhang H J, Thurber C H, et al. 2015. Joint inversion of seismic and magnetotelluric data in the Parkfield region of California using the normalized cross-gradient constraint. Pure and Applied Geophysics , 172 (5) : 1033-1052. DOI:10.1007/s00024-014-1002-9
Berge P A, Berryman J G, Bertete-Aguirre H, et al. 2000. Joint inversion of geophysical data for site characterization and restoration monitoring. LLNL Rep. UCRL-ID-128343, Proj, 55411.
Calvetti D, Morigi S, Reichel L, et al. 2000. Tikhonov regularization and the L-curve for large discrete ill-posed problems. Journal of Computational and Applied Mathematics , 123 (1-2) : 423-446. DOI:10.1016/S0377-0427(00)00414-3
Di Stefano R, Chiarabba C. 2002. Active source tomography at Mt. Vesuvius:Constraints for the magmatic system. Journal of Geophysical Research:Solid Earth , 107 (B11) . DOI:10.1029/2001JB000792
Fan C S, Li T L, Yan J Y. 2012. Research and application experiment on 2. 5D SIP inversion. Chinese J. Geophys. (in Chinese) , 55 (12) : 4044-4050. DOI:10.6038/j.issn.0001-5733.2012.12.016
Feng R, Li Z M, Li Z W, et al. 2004. Resistivity tomography technology. Earthquake Research in China (in Chinese) , 20 (1) : 13-30.
Gallardo L A, Meju M A. 2003. Characterization of heterogeneous near-surface materials by joint 2D inversion of DC resistivity and seismic data. Geophys. Res. Lett. , 30 (13) : 1658.
Gallardo L A, Meju M A. 2004. Joint two-dimensional DC resistivity and seismic travel time inversion with cross-gradients constraints. Journal of Geophysical Research:Solid Earth , 109 (B3) . DOI:10.1029/2003JB002716
Gallardo L A, Pérez-Flores M A, Gómez-Treviño E. 2005. Refinement of three-dimensional multilayer models of basins and crustal environments by inversion of gravity and magnetic data. Tectonophysics , 397 (1-2) : 37-54. DOI:10.1016/j.tecto.2004.10.010
Gallardo L A, Meju M A. 2007. Joint two-dimensional cross-gradient imaging of magnetotelluric and seismic traveltime data for structural and lithological classification. Geophysical Journal International , 169 (3) : 1261-1272. DOI:10.1111/gji.2007.169.issue-3
Gallardo L A, Fontes S L, Meju M A, et al. 2012. Robust geophysical integration through structure-coupled joint inversion and multispectral fusion of seismic reflection, magnetotelluric, magnetic, and gravity images:Example from Santos Basin, offshore Brazil. Geophysics , 77 (5) : B237-B251. DOI:10.1190/geo2011-0394.1
Gallardo-Delgado L A, Pérez-Flores M A, Gómez-Treviño E. 2003. A versatile algorithm for joint 3D inversion of gravity and magnetic data. Geophysics , 68 (3) : 949-959. DOI:10.1190/1.1581067
Gao G Z, Abubakar A, Habashy T M. 2012. Joint petrophysical inversion of electromagnetic and full-waveform seismic data. Geophysics , 77 (3) : WA3-WA18. DOI:10.1190/geo2011-0157.1
Haber E, Oldenburg D. 1997. Joint inversion:a structural approach. Inverse Problems , 13 (1) : 63-77. DOI:10.1088/0266-5611/13/1/006
Hansen P C, O'Leary D P. 1993. The use of the L-curve in the regularization of discrete ill-posed problems. SIAM Journal on Scientific Computing , 14 (6) : 1487-1503. DOI:10.1137/0914086
Hu W Y, Abubakar A, Habashy T M. 2009. Joint electromagnetic and seismic inversion using structural constraints. Geophysics , 74 (6) : R99-R109. DOI:10.1190/1.3246586
Hu X Y, Yang D K, Liu S H, et al. 2006. The developing trends of environmental and engineering geophysics. Progress in Geophysics (in Chinese) , 21 (2) : 598-604.
Jupp D L B, Vozoff K. 1975. Stable iterative methods for the inversion of geophysical data. Geophysical Journal International , 42 (3) : 957-976.
Lei X Y, Li Z W, Zhe J P. 2009. Applications and research of the high resolution resistivity method in exploration of caves, mined regions and Karst region. Progress in Geophysics (in Chinese) , 24 (1) : 340-347.
Li Z X, Tan H D, Fu S S, et al. 2015. Two-dimensional synchronous inversion of TDIP with cross-gradient constraint. Chinese J. Geophys. (in Chinese) , 58 (12) : 4718-4726. DOI:10.6038/cjg20151232
Liu B, Li S C, Li S C, et al. 2012. 3D electrical resistivity inversion with least-squares method based on inequality constraint and its computation efficiency optimization. Chinese J. Geophys. (in Chinese) , 55 (1) : 260-268. DOI:10.6038/j.issn.0001-5733.2012.01.025
Liu Y H, Yin C C. 2013. 3D inversion for frequency-domain HEM data. Chinese J. Geophys. (in Chinese) , 56 (12) : 4278-4287. DOI:10.6038/cjg20131230
Loke M H, Barker R D. 1996. Rapid least-squares inversion of apparent resistivity pseudosections by a quasi-Newton method. Geophysical Prospecting , 44 (1) : 131-152. DOI:10.1111/gpr.1996.44.issue-1
Loke M H, Dahlin T. 2002. A comparison of the Gauss-Newton and quasi-Newton methods in resistivity imaging inversion. Journal of Applied Geophysics , 49 (3) : 149-162. DOI:10.1016/S0926-9851(01)00106-9
Maceira M, Ammon C J. 2009. Joint inversion of surface wave velocity and gravity observations and its application to central Asian basins shear velocity structure. Journal of Geophysical Research:Solid Earth (1978-2012) , 114 (B2) . DOI:10.1029/2007JB005157
McGillivray P R, Oldenburg D W. 1990. Methods for calculating Fréchet derivatives and sensitivities for the non-linear inverse problem:a comparative study. Geophysical Prospecting , 38 (5) : 499-524. DOI:10.1111/gpr.1990.38.issue-5
Musil M, Maurer H R, Green A G. 2003. Discrete tomography and joint inversion for loosely connected or unconnected physical properties:Application to crosshole seismic and georadar data sets. Geophysical Journal International , 153 (2) : 389-402. DOI:10.1046/j.1365-246X.2003.01887.x
Osher S, Fedkiw R. 2003. Level Set Methods and Dynamic Implicit Surfaces. New York:Springer. http://cn.bing.com/academic/profile?id=2093834886&encoded=0&v=paper_preview&mkt=zh-cn
Paige C C, Saunders M A. 1982. LSQR:An algorithm for sparse linear equations and sparse least squares. ACM Transactions on Mathematical Software (TOMS) , 8 (1) : 43-71. DOI:10.1145/355984.355989
Pan J S, Ge W Z, Zhe J P. 2010. High-density electrical resistivity imaging technique of ground/borehole-to-surface/inter-well. Journal of North China Institute of Water Conservancy and Hydroelectric Power (in Chinese) , 31 (2) : 74-78.
Pan K J, Tang J T, Hu H L, et al. 2012. Extrapolation cascadic multigrid method for 2. 5D direct current resistivity modeling. Chinese J. Geophys. (in Chinese) , 55 (8) : 2769-2778. DOI:10.6038/j.issn.0001-5733.2012.08.028
Peng M. 2012. Joint inversion of magnetotelluric and teleseismic data[Ph. D. thesis] (in Chinese). Beijing:China University of Geosciences (Beijing).
Pidlisecky A, Haber E, Knight R. 2007. RESINVM3D:A 3D resistivity inversion package. Geophysics , 72 (2) : H1-H10. DOI:10.1190/1.2402499
Popovici A M, Sethian J A. 2002. 3-D imaging using higher order fast marching traveltimes. Geophysics , 67 (2) : 604-609. DOI:10.1190/1.1468621
Rawlinson N, Sambridge M. 2003. Seismic traveltime tomography of the crust and lithosphere. Advances in Geophysics , 46 : 81-199. DOI:10.1016/S0065-2687(03)46002-0
Ritter J R R, Jordan M, Christensen U R, et al. 2001. A mantle plume below the Eifel volcanic fields, Germany. Earth and Planetary Science Letters , 186 (1) : 7-14. DOI:10.1016/S0012-821X(01)00226-6
Sambridge M S. 1990. Non-linear arrival time inversion:constraining velocity anomalies by seeking smooth models in 3-D. Geophysical Journal International , 102 (3) : 653-677. DOI:10.1111/gji.1990.102.issue-3
Tikhonov A N, Arsenin V Y. 1977. Solutions of Ill-Posed Problems. Washington, DC:Vh Winston.
Treitel S, Lines L. 1988. Geophysical examples of inversion (with a grain of salt). The Leading Edge , 7 (11) : 32-35. DOI:10.1190/1.1439464
Vidale J. 1988. Finite-difference calculation of travel times. Bulletin of the Seismological Society of America , 78 (6) : 2062-2076.
Vidale J E. 1990. Finite-difference calculation of traveltimes in three dimensions. Geophysics , 55 (5) : 521-526. DOI:10.1190/1.1442863
Wang J. 2014. Research of the theory of joint inversion of gravity and magnetic data based on cross-gradient function[Master's thesis] (in Chinese). Beijing:China University of Geosciences (Beijing).
Wu X P, Xu G M. 2000. Study on 3-D resistivity inversion using conjugate gradient method. Chinese J. Geophys. (in Chinese) , 43 (3) : 420-427.
Wu X P, Liu Y, Wang W. 2015. 3D resistivity inversion incorporating topography based on unstructured meshes. Chinese J. Geophys. (in Chinese) , 58 (8) : 2706-2717. DOI:10.6038/cjg20150808
Yu P, Wang J L, Wu J S, et al. 2006. Review and discussions on geophysical joint inversion. Progress in Exploration Geophysics (in Chinese) , 29 (2) : 87-93.
Zelt C A, Smith R B. 1992. Seismic travel time inversion for 2-D crustal velocity structure. Geophysical Journal International , 108 (1) : 16-34. DOI:10.1111/gji.1992.108.issue-1
Zhang F X, Wu Q J, Li Y H, et a1. 2010. Application of FMM ray tracing to forward and inverse problems of seismology. Progress in Geophysics (in Chinese) , 25 (4) : 1197-1205. DOI:10.3969/j.issn.1004-2903.2010.04.008
Zhou L F. 2012. Two dimensional joint inversion of MT and seismic data[Master's thesis] (in Chinese). Beijing:China University of Geosciences (Beijing).
范翠松, 李桐林, 严加永. 2012. 2.5维复电阻率反演及其应用试验. 地球物理学报 , 55 (12) : 4044–4050. DOI:10.6038/j.issn.0001-5733.2012.12.016
冯锐, 李智明, 李志武, 等. 2004. 电阻率层析成像技术. 中国地震 , 20 (1) : 13–30.
胡祥云, 杨迪坤, 刘少华, 等. 2006. 环境与工程地球物理的发展趋势. 地球物理学进展 , 21 (2) : 598–604.
雷旭友, 李正文, 折京平. 2009. 超高密度电阻率法在土洞、煤窑采空区和岩溶勘探中应用研究. 地球物理学进展 , 24 (1) : 340–347.
李兆祥, 谭捍东, 付少帅, 等. 2015. 基于交叉梯度约束的时间域激发极化法二维同步反演. 地球物理学报 , 58 (12) : 4718–4726. DOI:10.6038/cjg20151232
刘斌, 李术才, 李树忱, 等. 2012. 基于不等式约束的最小二乘法三维电阻率反演及其算法优化. 地球物理学报 , 55 (1) : 260–268. DOI:10.6038/j.issn.0001-5733.2012.01.025
刘云鹤, 殷长春. 2013. 三维频率域航空电磁反演研究. 地球物理学报 , 56 (12) : 4278–4287. DOI:10.6038/cjg20131230
潘纪顺, 葛为中, 折京平. 2010. 地面/井地/井间超高密度电阻率成像技术. 华北水利水电学院学报 , 31 (2) : 74–78.
潘克家, 汤井田, 胡宏伶, 等. 2012. 直流电阻率法2. 5维正演的外 推瀑布式多重网格法. 地球物理学报 , 55 (8) : 2769–2778. DOI:10.6038/j.issn.0001-5733.2012.08.028
彭淼. 2012. 大地电磁与天然地震数据联合反演研究[博士论文]. 北京:中国地质大学(北京).
王俊. 2014. 基于交叉梯度重磁联合反演方法技术研究[硕士论文]. 北京:中国地质大学 (北京). http://cdmd.cnki.com.cn/article/cdmd-11415-1014238597.htm
吴小平, 徐果明. 2000. 利用共轭梯度法的电阻率三维反演研究. 地球物理学报 , 43 (3) : 420–427.
吴小平, 刘洋, 王威. 2015. 基于非结构网格的电阻率三维带地形反演. 地球物理学报 , 58 (8) : 2706–2717. DOI:10.6038/cjg20150808
于鹏, 王家林, 吴健生, 等. 2006. 地球物理联合反演的研究现状和分析. 勘探地球物理进展 , 29 (2) : 87–93.
张风雪, 吴庆举, 李永华, 等. 2010. FMM射线追踪方法在地震学正演和反演中的应用. 地球物理学进展 , 25 (4) : 1197–1205. DOI:10.3969/j.issn.1004-2903.2010.04.008
周丽芬. 2012. 大地电磁与地震数据二维联合反演研究[硕士论文]. 北京:中国地质大学 (北京). http://cdmd.cnki.com.cn/article/cdmd-11415-1012364756.htm