地球物理学报  2012, Vol. 55 Issue (3): 970-980   PDF    
电阻率成像的混合正则化反演算法
韩波1 , 窦以鑫2 , 丁亮3     
1. 哈尔滨工业大学数学系, 哈尔滨 150001;
2. 哈尔滨工业大学航天学院航天工程与力学系, 哈尔滨 150001;
3. 东北林业大学数学系, 哈尔滨 150001
摘要: 利用混合正则化方法对二维分片常值电阻率进行反演研究,该方法结合了全变差正则化方法(total variation regularization)和经典吉洪诺夫正则化方法(classical Tikhonov regularization)的优点. 全变差正则化能够有效地重构分片常值电阻率和识别边界,而经典吉洪诺夫正则化方法则能有效地重构光滑的电阻率,从数值算例中可以看出,将这两种方法有效地结合到一起能够改进重构效果.
关键词: 电阻率成像      经典吉洪诺夫正则化      全变差正则化方法      混和正则化方法     
Electrical resistivity tomography by using a hybrid regularization
HAN Bo1, DOU Yi-Xin2, DING Liang3     
1. Department of Mathematics, Harbin Institute of Technology, Harbin 150001,China;
2. Astronautical Engineering and Mechanics,School of Astronautics,Harbin Institute of Technology, Harbin 150001,China;
3. Department of Mathematics, Northeast Forest University, Harbin 150001,China
Abstract: In this paper, we investigate the reconstruction of a piecewise constant resistivity in 2D by using a hybrid regularization. This hybrid regularization integrates a total variation (TV) regularization with a classical Tikhonov regularization. A TV regularization can reconstruct a piecewise constant resistivity as well as detect jump boundaries, and a classical Tikhonov regularization has the ability to reconstruct a smooth resistivity. Combining a classical Tikhonov regularization with a TV regularization, we can see that the hybrid regularization improves the results of numerical experiments..
Key words: Electrical resistivity tomography      Classical Tikhonov regularization      Total variation regularization      Hybrid regularization     
1 引言

近年来随着高密度电阻率法的快速发展,直流电法在岩土工程、资源、环境等领域的地质勘查中得到广泛的应用[1-4].和传统的电阻率测量模式所不同的是电阻率成像技术能够用于解决比较复杂的地质结构问题,通过在地表或井中的发射源向勘探区域注入直流电,布设在地表或井中的电极系获取关于地下电阻率信息的大量实测数据进行电阻率正反演精确重构地下电阻率图像.

在最近十几年,地下电阻率图像技术取得了飞速的发展,其中的一些成果已经被水文学家、工程师、地质学家成功地应用到地下模型的识别.电阻率的反演问题是从接收到的电位识别地下电阻率的构成.从反演理论的角度来看,这是一个非常具有挑战性的工作.因为电阻率反演成像是非线性、不适定问题[5].反演问题是一个正、反演计算不断迭代优化的过程,相对于反演而言,正演问题的研究相对比较成熟,算法的多样性给了我们更多的选择空间.Dey和Morrison[6]采用了7点差分格式;Spitzer[7]提出了共轭梯度法的有限差分算法,降低了内存的使用并提高了计算速度.强建科和罗延钟[8]提出了直流电阻率有限元法模拟.汤井田和公劲喆[9]进行了有限元-无限元耦合模拟研究.汤井田等[10]利用后验误差估计技术,研究了自适应有限元算法,使得网格剖分更具有灵活性,改善了计算效果.为了提高计算效率,Moucha和Bailey[11]设计了一个有效稳定的多重网格算法.最近,鲁晶津等[12]针对带有跃变电阻率的情况提出了代数多重网格(AMG),并且指出随着网格节点数的增加,代数多重网格算法的高效性更加凸显.近年来,在电阻率反演研究方面具有代表性的成果有,Li和Oldenberg[13]提出了基于Born近似反演,这种方法不计算偏导数矩阵,只在波数域里进行了一维线性反演,通过对一维线性反演结果进行傅氏变换得到近似解.很多学者利用优化算法重构电阻率,例如:吴小平等[14]用梯度迭代技术实现了电阻率反演;宛新林等[15]用最小二乘正交分解法实现电阻率反演;徐海浪和吴小平[16]用神经网络的方法实现了电阻率重构.目前,在简单地质环境下,应用电阻率成像技术在一定程度上能够提供满足工程需要的成像结果.但是,当地下情况较为复杂时,现有的技术却很难提供较好的重构结果.经典的吉洪诺夫正则化方法(带有l2 罚项)的反演策略会将边界过度光滑化,而全变差正则化方法则具有识别边界的能力,但反演结果对初值的依赖较大.本文提出了利用混合正则化方法改进已有的经典吉洪诺夫正则化方法和全变差正则化方法的计算方案,研究了较为复杂地质情况下的电阻率的重构问题.结合数值算例对算法进行了详细的分析,给出了一些具有实际意义的结论.全变差正则化方法是Rudin 等[17]在1992 年首次提出,主要是识别图像边界问题.该方法提出后,很快得到了图像处理领域的重视.最近,Wang 等[18]利用混和正则化方法研究了地震波的反演,Gao 和Zhao[19]在2010年将全变差正则化方法应用到生物荧光成像中(BLT),Burstedde和Ghattas[20]将全变差正则化方法应用到一维全波场反演(full wave for min version)研究,取得了非常好的反演效果.Bachmay和Burger[21]在全变差正则化方法的基础上给出了混和正则化方法的收敛性证明和相应的数值迭代算法.借鉴这些方法和理论,本文引进混和正则化方法进行二维地面电阻率成像研究.众所周知,初值选取对反演成像精确性有重要影响,为了提高反演效果,引进调节参数平衡两个不同的罚项,数值算例表明,这样的算法扩大了初值的选取范围也提高了重构精度.

通过分析较为复杂的二维数值算例,初步掌握了混合正则化方法在非线性反演问题应用时的一些数值特性,指出了算法的一些优势和不足之处,为下一步应用该算法解释三维模型提供了必要的理论依据和算法特征分析.最后给出了该算法解释简单三维地质情况的数值报告.

2 数学模型

描述二维电阻率的数学模型可以归结为泊松方程的稳恒电流源问题:

(1)

其中,σ 为电导率,U为电位,δ 为狄拉克函数,r为源点到边界的向径,Γs 为区域的地面边界,Γ 为区域的地下边界,A是电源空间位置,Ω = [0,Lx]×[0,Lz],LxLz分别表示两个正数.

2.1 正演算法

本文主要考虑的是具有间断的电导率,因此正演采用有限元方法模拟电位空间分布情况.由于采用的激发源是点电源,因此,为了更加准确和快速地模拟电场分布情况,本文采用自适应网格技术.即在电源附近采用加密的网格,而在离电源较远的地方采用粗网格,在足够远的地方,采用更加粗糙的网格.采用这样的网格进行计算主要是有以下两个方面的考虑:(Ⅰ)从物理学角度讲,点电荷的附近电场强度最大,能否有效地扑捉到电源周围的电场,这对于电场分布计算至关重要,因此,在点电源附近采用加密的网格;(Ⅱ)从数学的角度讲,网格的逐渐加密,固然能更加准确地模拟正演电场分布情况,但给反演成像带来了巨大的计算量,不符合实际要求.因此,必须采用自适应网格技术.具体网格参见图 1.网格剖分是基于Matlab工具箱“pdetool"中的程序进行的修改.在二维问题中采用结构化三角单元网格剖分.当区域进行分解后,首先在每个单元中将电位U表示为插值函数的形式:

其中,Ni为插值形函数,ui为电位在节点处的电位值.为了将方程(1)转化为变分形式,选取测试函数空间基底为Ni,也就是说,测试函数空间与基函数空间是一致的.有限元离散形式为

这里,离散电导率σ = (σ1σ2,…,σk),k表示有限元单元数,连续电导率σ 在每个单元上是常值,之所以要在每个单元上是常值,主要是考虑到全变差罚项的计算问题.这样就可以将连续型的方程写成下面的离散型:

(2)

式中,算子A表示离散正演算子,U表示离散电位U= (u1u2,…,uL),L表示有限元节点数,F表示离散电源.本文在求解代数方程组(2)时,采用拟最小余项法(Quasi Minimum Residual, QMR)解决该问题,该方法在Matlab中易于实现.

2.2 反演算法

电导率反演问题是一个严重不适定的问题,即观测值的微小变化都有可能引起电导率的巨大变化.解决这样的问题,一般用经典吉洪诺夫正则化技术进行反演:

(3)

(4)

式中,d表示观测数据;Q是投影算子,将正演计算得到的电位投影到观测点;α 为正则化参数.该方法针对连续形式的电导率能够得到很好的重构结果,但对于具有跳跃的电导率一般都会出现过光滑现象.重构结果是将原来跳跃部分光滑化,而真实地下结构往往是以分片常值分布的,因此不能得到满意的结果.针对分片常值的电导率重构问题,全变差正则化方法是一个自然的选择,具体形式如下:

(5)

(6)

式中,‖Δσl1 可以表示为Sk表示三角形单元一个边长的长度,该表面的全局编号为kL(k),R(k)分别表示与第k个边相邻的两个三角形单元的全局编号.在反演计算过程中,发现该方法具有下面的性质:当已知电导率的分布区域,而重构电导率数值时,能够得到很好的重构效果;当不知分布区域,但是初值选取接近真值时,该算法也能得到较好的重构效果;但是,如果既不已知分布区域,初值选取又远离真值时,该算法就不能得到满意的重构效果.分析造成这种现象的原因为经典吉洪诺夫正则化对于连续型电导率重构有效,而全变差正则化对于分片常值重构有效.但是,当电导率初值远离真实电导率并且不同电导率所在区域也是未知的时候,分片常值电导率应理解即具有间断的特性,同时在每个区域上的常值电导率又具有连续特性.因此,当电导率的大小和所在区域都未知的时候,应该结合两个正则化方法的优点重构电导率.鉴于此,本文采用下面的混合正则化方法重构电导率[18]:

(7)

(8)

式中,αγ 表示正则化参数,按照反演理论,这两个参数的选取应该与噪音水平有关,起到平衡数据拟合项和罚项之间的作用,而两者的比值起到了调节两个不同罚项的作用.在数值算例中将(7)式写成如下形式:

(9)

(10)

式中参数γ 是用来调节全变差罚项‖Δσl1l2罚项‖σl22 在重构中的作用,该参数的选取将在数值模拟中介绍.本文采用梯度法求解罚函数(7)式的最小值,即:

(11)

参数τ 表示梯度法的步长,采用线性搜索法,确保不等式J(σk+1)< J(σk)成立.为了计算J′(σ),将J(σ)分解为两部分:J1(σ)=1/2‖QU-dl22J2(σ)=α(γΔσl1+(1-γ)‖σl22).对于第一部分的导数可以表示为:J1(σ)= Q* (∂U/∂σ) *(QU-d),该导数可以通过对偶方法来求得,并不会带来计算困难,计算量大约是正演计算的两倍;为了计算J2(σ)的导数,首先分析电导率连续的情况.因为对于连续电导率σ,‖ΔσL1 的导数在Δσ =0处是不存在的,为了克服这样的困难,采用了下面的近似表达式:代替,因此第二部分的导数可以表示为

需要说明的是,采用这样的近似会产生一个计算的困难,当参数ε选取过小的时候,会使计算格式不稳定,当参数ε选取过大的时候,会使重构的结果过于光滑化.但是在数值算例中可以参考图像处理中关于该参数的选取方法,该参数建议取为10-8 ~10-6之间.针对电导率离散的情况,σl1 范数可以表示为:Sk= (0,…,0,Sk,0,…,-Sk,0,…,0),非零元素Sk,-Sk分别位于第L(k),R(k)个位置,L(k),R(k)分别表示与第k个边长相邻的两个三角形单元(二维)的全局编号,Sk表示边长.注意到离散的散度算子Δh•和离散的梯度算子Δh之间有如下关系式: Δh•=-ΔhT,上角T表示转置,因此J2(σ)=α(-γSTE-1Sσ+(1-γ)σ),E=diag(ηk),ηk= ,矩阵S的第k行是由行向量Sk组成,ST,E-1 分别表示SE的转置矩阵和逆矩阵.

2.3 数值模拟

假设初值和真值选取比较接近时,如果已知电导率σ 是分片常值函数,那么参数γ 应该接近于1,突出全变差罚项的作用;如果已知电导率σ 是连续函数,那么参数γ 应该接近于0,突出l2 罚项的作用.在电导率σ 是分片常值函数情况下,如果初值远离真值时,既要兼顾连续性也要考虑间断性,因此,本文提出下面的反演策略:

第一步:首先考虑电导率σ 的连续性,将罚函数(9)式中的正则化参数分别设置为α=10-6 (对应于观测数据不含噪音),γ = 0(此时的罚函数就是经典吉洪诺夫罚函数),利用优化算法计算罚函数(9)式的最小值,得到反演结果.一般来说,这个反演结果在电导率σ 间断处虽然比较光滑,但是还是有一定的跳跃现象存在;

第二步:再将上面得到的结果作为初值,进行下一步的迭代,根据第一步重构的结果,选择参数γ,目的就是有效地平衡电导率σ 的间断性和连续性,最后通过计算罚函数(9)式的最小值重构电导率.本文采用的数值算法可以简单的描述为下面的流程图形式,见图 2.流程图中的N 表示的是γ∈ [0, 1]的剖分数量.在本文的数值算例中,采用的剖分是均匀剖分,剖分数量N=10.

图 2 算法流程图 Fig. 2 Flow chart of the proposed algorithm
2.4 停止准则

反演问题中停止准则直接影响到重构效果.针对外层迭代(即参数γ)采用的停止准则分为含有噪音和无噪音两种情况:(I)含有噪音:停止准则选为‖QU-dobsl2τδnδn是观测数据的噪音水平,即:‖dobs-dtruel2 =δndobs, dtrue 分别表示仪器观测数据和客观真实数据.在本算例中建议1.0≤τ≤1.1,采用的是τ=1.05;(II)无噪音:这时噪音水平δn=0,选取‖QU-dobsl2 中最小值所对应的电导率作为最后结果.在内层迭代中(即固定参数γ),停止准则包含两个部分:第一部分,如果满足停止准则立即停止迭代,此时的重构的电导率即为所求,第二部分是为确保算法能够在有限步停止,加上另外一个停止准则:

这里j表示内层迭代步数.

2.5 有关罚函数(9)式的实现细节

通过大量的数值算例发现,全变差罚项与l2 罚项之间由于存在尺度上的差异,因此需要将罚函数(9)式进行适当修正.在第n+1 次外层迭代过程中,利用离散形式的电导率σγnrec计算σγn+1rec 时(这里σγnrec表示第n次外层迭代的重构电导率,γn表示第n次外层迭代所对应的γ 值),罚函数(9)式的具体修正形式为

(12)

(13)

这样修改的目的是避免由于 或者存在,使得难于通过参数γ 的分割有效地平衡两个不同的罚项.修改后的罚函数(12)式使得两个不同罚项在同一尺度下,因此,通过简单的等分参数γ 就可以平衡两者之间的作用.当γ =0时,优化问题(12)、(13)式改写为标准的经典吉洪诺夫正则化方法:

(14)

(15)

3 数值模拟

在数值算例中,采用的激发源是符合实际的点电源.参数αγ 的选取对于反演十分重要,可以说,这两个参数选取本身就是一个重要的研究课题,已经超出了本文的研究内容,在数值算例中,给出了针对具体数值算例的参数选取策略和物理解释,而并没有提出严格的数学理论解决这两个参数选取的问题.数值算例表明,尽管这两个参数选取很重要,但是,在实际问题中参数的选取具有鲁棒性,这个性质给处理实际数据带来了很大的方便.正则化参数α是平衡数据拟合部分和罚项部分之间的权重.一般来说,如果选取的过大,会造成过稳定的现象发生,使得后面的罚项严重影响数据拟合部分的计算.通过数值算例发现针对于电导率反演,参数α 选取一般为10-4~10-6之间.当然,该值并不是固定不变的,在处理实际地质资料的时候,参数α 的选取可以通过实验室的模拟结果等手段进行选取.根据经验来看,该值的选取应该针对不同问题进行具体问题具体分析,这样才能更好地处理实际数据.在本文的数值算例中,考虑了参数α 的选取同噪音水平有直接关系,因此,观测数据没有噪音情况下,参数α=10-6,在含有噪音的情况下,参数α=10-4.

在数值算例部分,将对二维地表电阻率成像进行研究.二维正演算法的Matlab 代码是在Adler和Lionheart的EIT(electrical impedance to mography)基础上进行的改进.本文主要考虑的是较为复杂的二维地表电阻率成像问题,数值算例包含观测数据无噪音和含有噪音两部分.分为两部分的目的在于:通过在观测数据无噪音情况下电阻率的重构,可以分析混和正则化方法、全变差正则化方法同经典吉洪诺夫正则化方法之间的区别和特点,这些特点在有噪音的时候,很难说清是算法本身所特有的还是由于噪音引起的;而通过研究噪音数据,可以说明算法的稳定性和处理实际数据的可行性,另外通过这些比较还能分析正则化参数α 的选取方法,以及验证停止准则的有效性.

3.1 二维浅地表电阻率成像

首先针对观测数据无噪音,三个异常体位于浅地层的情况进行研究.考虑到实际地下介质分布复杂,高低电阻率并存的情况比比皆是,因此,考虑的地下电阻率分布情况如图 1 所示.区域Ω =[-12, 12]×[-15, 0].在本文中,表示空间点的坐标单位为m.背景介质电阻率是125Ωm, 而蓝色所围成区域的异常体电阻率是25Ωm, 红色所围成区域的异常体电阻率是250 Ωm.地表布置21 个电极,则单极-单极视电阻率数据有378个.图 1是真实电阻率分布示意图.地下含有两个低电阻率异常体和一个高电阻率异常体.为了保证正演的精确性,采用细网格进行正演计算.远离电极区域的三角单元直角边长度为1,在距离电极较近区域的直角边为0.5,而在点电极周围的直角边是0.25.为了提高反演的惟一性和稳定性,将反演参数限定在[-8, 8]× [-11.5,-0.5]范围之内,并采用粗网格对反演参数进行剖分.图 3展示的是利用经典吉洪诺夫正则化方法(对应于γ=0)重构的结果(初始电阻率选为背景介质电阻率125Ωm).可以看出,重构的效果并不是很理想,尽管提供了大体的位置信息,但是与真实模型仍有较大的差别,特别是判断异常体深度方面,提供了错误信息.尽管如此,重构的地下电阻率能够捕捉到真实异常体存在的区域,并且在该区域的重构电阻率与真实电阻率十分接近.将该结果作为下一步的初值是一个自然的思路.图 4 表示的是全变差正则化方法重构结果.从图 3图 4中的双箭头线段可以看出全变差能够削弱重构结果的连续性,增强间断性.从图 4中还可以看出全变差正则化方法所具有的一个显著特点-阶梯效应(stair casing effect)[22],该效应主要体现在当反演的函数具有连续部分的时候,或者说是分片光滑函数的时候,就会出现阶梯状的重构结果.阶梯效应是全变差正则化方法所特有的现象,理论分析可以参见文献[22].从图 4 中可以看出,虽然全变差正则化方法能够增强重构结果的间断性,但是并没有对经典吉洪诺夫正则化方法有明显改进,可能的原因是:在非线性反演问题中,全变差正则化方法的收敛域小于经典吉洪诺夫正则化方法的收敛域[20].比较图 3图 4后,还得到如下结论:两个正则化方法能够有效地判断异常体横向位置,但是,对于识别异常体纵向位置显得无能为力,表现为不能准确定位异常体的下边界.

图 1 真实电阻率分布图 Fig. 1 The distribution of the true resistivity
图 3 经典吉洪诺夫正则化重构电阻率(无噪音) Fig. 3 Reconstructed resistivity using the classical Tikhonov regularization (noise free)
图 4 全变差正则化重构电阻率(无噪音) Fig. 4 The reconstructed resistivity by using the TV regularization (noise free)

图 5表示的是γ=0.5时,利用混和正则化方法重构的结果.初步分析混和正则化方法之所以能够得到最好的结果,除了全变差罚项的作用以外,l2罚项还能起到扩充收敛域的作用.另外参数γ 的灵活选择也有效地平衡两个不同罚项.与前两个独立的正则化方法相比,混和正则化方法有两个突出的优点:第一,能够减少全变差正则化方法的阶梯效应,这也是调节全变差罚项和l2 罚项之间的效果,需要注意的是,阶梯效应依然存在,并没有完全消失;第二,提高异常体纵向位置的识别能力,可以从图 5中看出,异常体的下边界能够有效地识别,这说明,测量数据含有足够的识别信息,因此重构异常体需要的是有效的数值算法提取这些信息.从上面的分析还可以看出,在识别浅地表异常体分布情况时,混合正则化方法在一定程度上改进了经典吉洪诺夫正则化方法和全变差正则化方法.

图 5 混合正则化重构电阻率γ=0.5(无噪音) Fig. 5 The reconstructed resistivity by using the hybrid regularization γ=0.5 (noise free)

针对不同的参数γ 给出定量的分析,从中发现混和正则化方法在电阻率成像问题中所具有的数值特性.为此在表 1中给出不同γ 所对应的数据参差(即:‖QU-dobsl2).从数据参差的走势可以看出,随着参数γ 的增加,l2 罚项的作用会减弱,数据参差也会随之下降,说明全变差罚项的作用起到了重要作用;当γ 超过了0.5 后,随着参数γ 的增加,l2 罚项的作用会继续减弱,但是,数据参差却出现了不减反增的情况,说明两个不同罚项之间正远离最佳的平衡状态.产生这种现象的可能原因是:在本算例中针对参数γ的剖分过大,在γ∈[0.0,0.5]范围内,l2 罚项起到了主导的作用,使得算法的收敛范围较大;但是,当γ∈ [0.5,1.0]范围内,全变差罚项起到了主导的作用,使得算法的收敛范围较小,当γ =0.5时重构的结果σγ=0.5 不能为γ=0.6的重构提供较好的初值时,从而造成这种数值现象的发生.事实上,通过数值实验发现,当在此处均匀加密后,数据参差继续减小,但是,经过计算一段时间后数据参差还是继续增加,从而还可以得到一个重要的特征:随着γ 趋于1,使得混合正则化算法的全变差罚项作用增强,从而算法的收敛范围逐渐减小,因此,建议自适应剖分参数γ (关于自适应剖分参数γ 的研究已经超出了本文的研究范围).

表 1 混合正则化方法对应不同参数7的数据参差(无噪音) Table 1 The data residual corresponding to different / using the hybrid regularization (noise free)

为了测试算法的稳定性,在接收数据上附加1%噪音.由于反演电阻率从数学的角度讲是一个严重不适定的问题,因此如果噪音强度过高,则不会达到满意的重构效果.在应用该算法处理实际数据时,如果已知噪音强度很大,建议先进行去噪处理.数值模拟的物理参数与上面相同.图 6展示的是利用经典吉洪诺夫正则化方法(对应于γ=0)重构的结果(初始电阻率选为背景介质电阻率125Ωm).考虑到正则化参数的选取应与噪音有关,因此正则化参数α 选为10-4.与数据无噪音的情况类似,重构的效果并不是很理想,尽管提供了大体的位置信息,但是与真实模型仍有较大的差别,特别是判断异常体深度方面,提供了错误信息.由于数据存在噪音,经典吉洪诺夫正则化方法在重构异常体深度方面几乎失去了作用,这点可以通过比较图 3图 6看出,虽然经典吉洪诺夫正则化方法在无噪音情况下(见图 3),探测异常体深度的能力十分有限,但是,两端的低电阻率异常体(蓝色区域)随着深度增加而接近背景值(白色区域),中间部分重构的高电阻率异常体(红色区域)随着深度增加,也越来越接近背景值(白色区域).这种逐渐过度的特性正是经典吉洪诺夫正则化方法所具有的数值特性,即光滑作用;从图 6中可以看出,噪音对于重构效果有严重影响,首先,对于异常体的深度判定几乎不可能,因为,从上到下几乎在数值上没有明显变化,虽然重构结果仍含有边界光滑的现象,但是,噪音对于重构效果的影响已远远超过了这种光滑现象.从图 36的分析可以看出,噪音对于经典吉洪诺夫正则化方法重构异常体有严重影响.当然,正则化参数α 的适当选取能够减少噪音的影响,本文仅仅考虑了α=10-2,10-3,10-4,10-5,10-6这5 个参数,通过数值模拟发现α=10-4是最好的选择,本文是在这种意义下讨论的问题.

图 6 经典吉洪诺夫正则化重构电阻率(噪音数据) Fig. 6 Reconstructed resistivity using the classical Tikhonov regularization (noisy data)

图 7是数据含有噪音情况下的全变差正则化方法的重构结果.通过比较图 4图 7可以得到下面结论:就非线性反演问题而言,全变差正则化方法对噪音很敏感,1%噪声使得重构的结果与无噪音结果相比大相径庭,但是,经典吉洪诺夫正则化方法对于噪音相比要稳定,针对有噪音和无噪音的情况,重构结果保持很多相似性.因此,调节两个罚项,期望能够改善重构结果.

图 7 全变差正则化重构电阻率(噪音数据) Fig. 7 The reconstructed resistivity by using the TV regularization (noisy data)

一种直观的改进方法就是以经典吉洪诺夫正则化方法为起点,通过逐步调节两个罚项的作用来改进重构效果.本数值试验的噪音水平为‖dobs -dtruel2 = 1.358479 (本文采用的范数是‖dobs -dtruel2 = ,整数I为观测数据的数量,标量ditrue是向量dtrue 的第i个元素,一些文献评价噪音水平采用下面的形式:‖dobs -dtrue‖std =ΣIi= ,这时的误差水平为0.0036).当γ=0.2的时候,重构结果满足停止准则.图 8表示的是γ =0.2时的混和正则化方法的重构结果,与数据无噪音的结果(图 5)相比,两者的重构结果很接近,这就说明算法对于噪声具有鲁棒性,这也是本文研究混和正则化方法的目的,因为具有稳定性质的算法对于处理实际数据至关重要.值得注意的是,尽管重构结果基本一致,但是,重构结果所对应的γ 值是不同的,这说明l2 罚项能够有效地抑制噪音带来的不稳定现象,当有噪音的时候,应该增加l2罚项的作用.表 2 给出不同γ 所对应的数据参差.经典吉洪诺夫正则化方法(γ=0)所对应的数据参差是1.8639,随着γ 的增加数据参差随之下降,从数据拟合角度来看,混和正则化方法也改进了经典吉洪诺夫正则化方法.

图 8 混合正则化重构电阻率γ=0.2(噪音数据) Fig. 8 The reconstructed resistivity by using the hybrid regularization /=0.2 (noisy data)
表 2 混合正则化方法对应不同参数y的数据参差(含有噪音) Table 2 The data residual corresponding to different / using the hybrid regularization (with noise)
3.2 二维深地表电阻率成像

本节主要研究两个异常体分别位于浅地层和深地层的情况,观测数据含有0.5%噪声.研究这种情况主要目的是初步掌握这种情况下两个不同罚项之间如何平衡的问题,并且说明本文算法的局限性和存在的问题.图 9是真实电阻率分布示意图.在反演过程中,正则化参数α 选为10-5.图 10是经典吉洪诺夫正则化方法重构的结果,从中可以看出,重构的结果基本可以包含真实的区域,而且数值上很接近.同3.1节相类似,仍然不能对异常体的深度信息进行有效恢复.图 11 是混和正则化方法最好的恢复结果,这时对应的γ 为0.8,这表明全变差罚项的权重大于l2 罚项的权重.比较图 10图 11可以看出:混合正则化方法有能力将高电阻率的异常体(位于深底层的红色区域)的横向边界有效地识别出来,但是,纵向边界的识别仍然不够理想.总的说来,针对这种情况,混和正则化方法改进了经典吉洪诺夫正则化方法,但是改进的程度是有限的,因此,深入研究混和正则化方法是十分有必要的.自适应选取γ 也许是重要的研究内容.另外,是否应该根据空间位置平衡混和正则化的两个不同罚项(即:参数γ 是否与空间坐标有关)也是研究的内容之一.本文采用的是局部收敛反演算法,研究具有较好全局收敛特性的信赖域算法(TrustRegion Method)也是重要的问题[23-25].这些问题的研究也许是混和正则化方法应用到三维问题的关键所在.

图 9 真实电阻率分布示意图 Fig. 9 The distribution of the true resistivity
图 10 经典吉洪诺夫正则化重构电阻率(噪音数据) Fig. 10 Reconstructed resistivity using the classical Tikhonov regularization(noisy data)
图 11 混合正则化重构电阻率γ=0.8 (噪音数据) Fig. 11 The reconstructed resistivity by using the hybrid regularization /=0.8 (noisy data)
4 结论和展望

本文利用混和正则化方法研究了二维电阻率成像问题,数值算例表明,该方法能够在一定程度上改善经典吉洪诺夫正则化方法和全变差正则化方法的重构效果.为了平衡两个不同的罚项,引进了调节参数γ.并且针对不同情况分析了混和正则化方法的特点,为实际应用提供了必要的数值分析结果.通过文中的算法分析,可以看出,混和正则化方法将两个单独的正则化方法有机的结合到一起,通过有效地调节各自的作用,使得两者都能够发挥各自的数值特点.另外,正演应该利用自适应网格技术,这样才能快速准确地模拟电场分布情况.尽管该方法已经改善了重构效果,但是,仍然存在一些缺点和问题.在探测异常体分别位于浅地层和深地层时,混和正则化方法的改进十分有限,这说明对于该方法在复杂情况下的数值特性研究十分重要.本文仅仅是初步讨论了混和正则化方法的一些简单的数值特性,因此对于调节参数γ 的剖分采用的是非常简单的等距剖分,从数值实验中可以看出,这样的剖分并不是最优的,因此,如何自适应的调节该参数也是重要的研究内容,我们推测,参数γ 的研究很有可能是混合算法数值特性分析的关键(当然正则化参数α的作用也十分重要,但是,α 是所有带有罚项正则化方法所具有的,并不是混和正则化方法所特有的).通过研究一些三维问题发现,混和正则化方法应用到三维问题存在一些难点.首先计算量的问题,三维问题的计算量是很大的,如何通过有效的数值计算手段改进混和正则化方法,提高计算效率,是很重要的研究问题,近年来,很多学者研究了提高计算效率的优化算法[26-28],在三维电阻率成像中如何将这些优化算法有机的与混和正则化算法融合到一起,从而提高计算效率也是重要的研究内容.另外,混和正则化方法作为一个重构技术,并不能完全解决三维电阻率重构这样一个严重不适定问题.通过一些简单的数值算例(本文不包含这些算例)发现,当利用地-地装置进行三维电阻率重构时,由于测量数据的信息量和需要重构的信息量相差悬殊,因此混和正则化技术并没有明显改善经典吉洪诺夫正则化的重构结果,这是由问题本质所决定.而利用井-地装置时却能改进重构效果,但是,这些仅仅是通过一些简单的数值算例看到的现象,三维情况下的混和正则化方法的特性较为复杂,很多因素影响参数γ的取值,例如:异常体的相对位置,异常体的大小,异常体的深度,异常体的电阻率值以及噪音强度等,因此,如何将混和正则化方法应用到三维问题也是一个有意义的研究方向.通过重构电阻率的研究表明,混和正则化思想很有可能推广到波动方程反演等一系列非线性地质参数反演问题中,应用到这些领域将遇到更加复杂的问题和数值特性.从数值计算角度来看,在求解‖ΔσL1 导数时采用了 近似,该方法的稳定性和准确性受参数ε 影响较大,利用内点法(interior-pointmethod)求解优化问题(10,11)可以避免参数ε的选取问题.内点法的实现和混和正则化方法的分辨率将是今后重要的研究内容,该研究方向将需要更多的数学理论工具和算法分析技术.本文采用的优化算法是局部收敛的算法,考虑到问题的复杂性,特别是研究三维情况的时候,如何将一些全局优化算法(例如:信赖域算法[23-25],同伦延拓算法[29-31])结合起来对于该问题的研究也许会发挥重要的作用.

致谢

作者特别感谢Adler, Lionheart 和Polydorides[3233]所提供的针对EIT 的二维正演Matlab代码,提高了本论文的进展速度.衷心感谢两位匿名审稿人仔细阅读了本文初稿,对算法提出了建设性的意见,纠正了一些表述不当之处,提高了本文写作水平.

参考文献
[1] Slater L, Binley A M, Daily W, et al. Cross-hole electrical imaging of a controlled saline tracer injection. Journal Applied Geophysics , 2000, 44(2-3): 85-102. DOI:10.1016/S0926-9851(00)00002-1
[2] Slater L, Binley A M, Versteeg R, et al. A 3D ERT study of solute transport in a large experimental tank. Journal of Applied Geophysics , 2002, 49(4): 211-229. DOI:10.1016/S0926-9851(02)00124-6
[3] Kemna A, Kulessa B, Vereeckena H. Imaging and characterisation of subsurface solute transport using electrical resistivity tomography (ERT) and equivalent transport models. Journal of Hydrology , 2002, 267(3-4): 125-146. DOI:10.1016/S0022-1694(02)00145-2
[4] Labrecque D J, Yang X J. Difference inversion of ERT data: a fast inversion method for 3-D in situ monitoring. Journal of Environmental and Engineering Geophysics , 2001, 6(2): 83-89. DOI:10.4133/JEEG6.2.83
[5] Pidlisecky A, Knight R, Haber E. Cone-based electrical resistivity tomography. Geophysics , 2006, 71(4): G157-G167. DOI:10.1190/1.2213205
[6] Dey A, Morrison H F. Resistivity modeling for arbitrarily shaped three-dimensional structures. Geophysics , 1979, 44(4): 753-780. DOI:10.1190/1.1440975
[7] Spitzer K. A 3-D finite-difference algorithm for DC resistivity modelling using conjugate gradient methods. Geophysical Journal International , 1995, 123(3): 903-914. DOI:10.1111/gji.1995.123.issue-3
[8] 强建科, 罗延钟. 三维地形直流电阻率有限元法模拟. 地球物理学报 , 2007, 50(5): 1606–1613. Qiang J K, Luo Y Z. The resistivity FEM numerical modeling on 3-D undulating topography. Chinese J. Geophys. (in Chinese) , 2007, 50(5): 1606-1613.
[9] 汤井田, 公劲喆. 三维直流电阻率有限元-无限元耦合数值模拟. 地球物理学报 , 2010, 53(3): 717–728. Tang J T, Gong J Z. 3D DC resistivity forward modeling by finite-infinite element coupling method. Chinese J. Geophys. (in Chinese) , 2010, 53(3): 717-728.
[10] 汤井田, 王飞燕, 任政勇. 基于非结构化网格的2.5-D直流电阻率自适应有限元数值模拟. 地球物理学报 , 2010, 53(3): 708–716. Tang J T, Wang F Y, Ren Z Y. 2.5-D DC resistivity modeling by adaptive finite-element method with unstructured triangulation. Chinese J. Geophys. (in Chinese) , 2010, 53(3): 708-716.
[11] Moucha R, Bailey R C. An accurate and robust multigrid algorithm for 2D forward resistivity modelling. Geophysical Prospecting , 2004, 52(3): 197-212. DOI:10.1111/gpr.2004.52.issue-3
[12] 鲁晶津, 吴小平, SpitzerK. 直流电阻率三维正演的代数多重网格方法. 地球物理学报 , 2010, 53(3): 700–707. Lu J J, Wu X P, Spitzer K. Algebraic multigrid method for 3D DC resistivity modelling. Chinese J. Geophys. (in Chinese) , 2010, 53(3): 700-707.
[13] Li Y G, Oldenberg D W. Approximate inverse mappings in DC resistivity problems. Geophysical Journal International , 1992, 109(2): 343-362. DOI:10.1111/gji.1992.109.issue-2
[14] 吴小平, 徐果明. 利用共轭梯度法的电阻率三维反演研究. 地球物理学报 , 2000, 43(3): 420–427. Wu X P, Xu G M. Study on 3-D resistivity inversion using conjugate gradient method. Chinese J. Geophys. (in Chinese) , 2000, 43(3): 420-427.
[15] 宛新林, 席道瑛, 高尔根, 等. 用改进的光滑约束最小二乘正交分解法实现电阻率三维反演. 地球物理学报 , 2005, 48(2): 439–444. Wan X L, Xi D Y, Gao E G, et al. 3D resistivity inversion by the least-squares QR factorization method under improved smoothness constraint condition. Chinese J. Geophys. (in Chinese) , 2005, 48(2): 439-444.
[16] 徐海浪, 吴小平. 电阻率二维神经网络反演. 地球物理学报 , 2006, 49(2): 584–589. Xu H L, Wu X P. 2D resistivity inversion using the neural network method. Chinese J. Geophys. (in Chinese) , 2006, 49(2): 584-589.
[17] Rudin L I, Osher S, Fatemi E. Nonlinear total variation based noise removal algorithms. Physica D: Nonlinear Phenomena , 1992, 60(1-4): 259-268. DOI:10.1016/0167-2789(92)90242-F
[18] Wang Y F, Cui Y, Yang C C. Hybrid regularization methods for seismic reflectivity inversion. Int. J. Geomath. , 2011, 2(1): 87-112. DOI:10.1007/s13137-011-0014-1
[19] Gao H, Zhao H K. Multilevel bioluminescence tomography based on radiative transfer equation Part 2: total variation and l1 data fidelity. Optics Express , 2010, 18(3): 2894-2912. DOI:10.1364/OE.18.002894
[20] Burstedde C, Ghattas O. Algorithmic strategies for full waveform inversion: 1D experiments. Geophysics , 2009, 74(6): WCC37-WCC46. DOI:10.1190/1.3237116
[21] Bachmayr M, Burger M. Iterative total variation schemes for nonlinear inverse problems. Inverse Problems , 2009, 25(10): 105004. DOI:10.1088/0266-5611/25/10/105004
[22] Strong D, Chan T. Edge-preserving and scale-dependent properties of total variation regularization. Inverse Problems , 2003, 19(6): S165-S187. DOI:10.1088/0266-5611/19/6/059
[23] Wang Y F, Cao J J, Yang C C. Recovery of seismic wavefields based on compressive sensing by an l1-norm constrained trust region method and the piecewise random subsampling. Geophysical Journal International , 2011, 187(1): 1999-213.
[24] Wen Z W, Wang Y F. A new trust region algorithm for image restoration. Science in China Series A , 2005, 48(2): 169-184. DOI:10.1360/03ys0178
[25] Wang Y F, Yuan Y X. A trust region method for solving distributed parameter identification problem. Journal of Computational Mathematics , 2003, 21(6): 759-772.
[26] Haber E. A parallel method for large scale time domain electromagnetic inverse problems. Applied Numerical Mathematics , 2008, 58(4): 422-434. DOI:10.1016/j.apnum.2007.01.017
[27] Kwon K, Yazici B, Guven M. Two-level domain decomposition methods for diffuse optical tomography. Inverse Problems , 2006, 22(5): 1533-1559. DOI:10.1088/0266-5611/22/5/002
[28] 李振华, 王彦飞, 杨长春. 正则化偏移成像的全局优化快速算法. 地球物理学报 , 2011, 54(3): 828–834. Li Z H, Wang Y F, Yang C C. A fast global optimization algorithm for regularized migration imaging. Chinese J. Geophys. (in Chinese) , 2011, 54(3): 828-834.
[29] 傅红笋, 韩波. 二维波动方程速度的正则化-同伦-测井约束反演. 地球物理学报 , 2005, 48(6): 1441–1448. Fu H S, Han B. A regularization homotopy method for the inverse problem of 2D wave equation and well log constraint inversion. Chinese J. Geophys. (in Chinese) , 2005, 48(6): 1441-1448.
[30] 陈勇, 韩波. 时间推移地震反演的连续模型与算法. 地球物理学报 , 2006, 49(4): 1164–1168. Chen Y, Han B. Continuous models and algorithms for time-lapse seismic inversion. Chinese J. Geophys. (in Chinese) , 2006, 49(4): 1164-1168.
[31] 窦以鑫, 韩波. 山体表面重构数值反演的同伦算法. 地球物理学报 , 2011, 54(7): 1893–1899. Dou Y X, Han B. Numerical inversion of reconstruction of mountain surface by homotopy method. Chinese J. Geophys. (in Chinese) , 2011, 54(7): 1893-1899.
[32] Adler A, Lionheart W R B. Uses and abuses of EIDORS: an extensible software base for EIT. Physiological Measurement , 2006, 27(5): S25-S42. DOI:10.1088/0967-3334/27/5/S03
[33] Polydorides N, Lionheart W R B. A Matlab toolkit for three-dimensional electrical impedance tomography: a contribution to the Electrical Impedance and Diffuse Optical Reconstruction Software project. Measurement Science and Technology , 2002, 13(12): 1871-1883. DOI:10.1088/0957-0233/13/12/310