地球物理学报  2012, Vol. 55 Issue (12): 4058-4068   PDF    
基于褶积模型的地球物理反演模型增强
左博新 , 胡祥云 , 韩波     
中国地质大学(武汉)地球物理与空间信息学院, 武汉 430074
摘要: 地球物理数据在采集和处理过程中, 由于存在噪声、模型误差、以及数据离散化误差等系统误差, 导致了异常体边界模糊和模型分辨率降低等一些不可避免的不良系统退化效应的产生.本文提出了一种新的地球物理反演模型增强方法, 通过消除反演估计模型中的系统误差, 压制模型中的不良系统退化效应, 增强反演模型的分辨率.文章从理论上分析了数据中存在的系统误差对模型求解的影响, 提出了一个新的系统误差褶积退化模型, 并根据该模型提出了一种基于混合范数总变分正则化的盲反褶积模型增强算法.最后, 文章通过1D线性反演增强试验和2D大地电磁反演增强试验, 验证了所提出的地球物理系统退化模型的正确性, 以及盲反褶积增强算法的有效性.试验结果表明, 方法可以有效地提高反演参数模型的分辨率.
关键词: 系统误差      点扩展函数      反褶积      增强      分辨率     
The geophysical model enhancement based on the convolution model
ZUO Bo-Xin, HU Xiang-Yun, HAN Bo     
Institute of Geophysics and Geomatics, China University of Geosciences (Wuhan), Wuhan 430074, China
Abstract: Systematic errors, such as noise, errors in forward and inversion models, and discretization error are inevitably present in acquisition and processing processes of geophysical data.They induce systematic degradation, such as blurring and resolution decreasing, to the geophysical models.This article proposes a new enhancement method for geophysical inversion models.Through reducing the systematic error, the method removes the blurring degradations and increases the model resolution.We analyze the influence of systematic errors on the model estimating process, and a new systematic error convolution degradation model is proposed.A mix-norm TV regularization deconvolution algorithm is also proposed to enhance the inversion model according to the proposed convolution model.Finally, we use the 1D linear inversion and 2D magnetotelluric inversion enhancement experiments to prove the correctness of the theoretical basis and the effectiveness of the proposed blind deconvolution enhancement algorithm.The results show that the proposed algorithms can effectively enhance the resolution of inversion models..
Key words: Systematic error      Point spread function      Deconvolution      Enhancement      Resolution     
1 引言

地球物理参数估计模型是对实际地质结构的客观反映.但在获取参数估计模型的各步骤中,会不可避免的存在各种系统误差.例如,在地球物理正演过程中,受正演模型精度和模型离散化误差的影响,会导致一个褶积模型的系统退化过程[1].类似的系统退化过程也存在于模型估计的其它步骤中.可以认为在真实模型与估计模型之间存在有一系列的级联系统退化效应[2],其对估计模型的影响由数据采集、正演计算和反演计算等步骤中存在的系统误差决定[3],这些系统误差导致了模型分辨率的降低.因此,通过消除估计模型中系统误差,可以有效减小系统误差对于估计模型的影响,由此增强参数估计模型的分辨率.

系统误差主要包括:噪声、测量误差、离散化误差、正反演模型误差、仪器退化误差等[3-4].因此,地球物理反演问题是一个病态问题,大多数情况下模型估计问题存在多解性和解不确定性.真实模型和反演估计模型间通常存在较大的差异[5].研究表明,虽然使用正则化方法可以减小反演问题解的病态性,但解的不确定性问题始终存在[6].本文根据地球物理反演理论,分析了系统误差对于估计模型分辨率的影响,根据推导的系统误差退化数学模型得出:即使在理想状况下得到的无偏参数估计模型,同样受到系统误差的影响而存在退化效应.

Backus[7-10]提出了地球物理模型分辨率的概念,提出了由点扩展函数PSF(PointSpreadFunction)来评价反演估计参数模型的理论.通过PSF 的数值分布分析可以评价各种地球物理方法对真实数据模型的分辨能力.在PSF 定义的基础上,Jackson[11]和Menke[12]定义了参数分辨率矩阵MRM(ParameterResolutionMatrix).Daily[13]等人根据实际的直流电阻率方法普查应用,将PSF 和MRM 对反演模型的分辨率评价结果,与实际普查情况进行了对比验证分析.结果显示:在MRM 矩阵的对角线上,较大的值对应了较差的分辨率,较小的值对应较好的分辨率.对比结果与PSF和MRM 的分辨率评价理论推导相符.

Alumbaugh[2]等人根据PSF 研究了在2D 和3D 电阻率反演模型中不同空间位置的空间分辨率分布,较好地解释了在实际的电阻率反演中反演模型分辨率的差异.Miller[14]等人分析了由Oldenburg[15]提出的RDI(RegionofDataInfluence)与PSF在数学上的关系,并进行了线性和非线性试验的验证,结果表明:PSF 的分辨率评价结果能够很好地符合RDI准则对于数据影响区域(RegionofData)的预测.Oldenborger[3]通过PSF有效地评价了3DERT(ElectricalResistivityTomography)电阻率成像的反演结果.并且分析认为PSF 是叠加在真实数据模型上的冲击响应.Alumbaugh[2]提出:PSF是一个连接在真实模型和估计模型之间的退化滤波器,并且设想根据PSF 从估计模型中提取出更多的真实数据模型相关信息.

根据以往对于MRM 矩阵的研究表明,MRM的数值分布受反演正则化项、噪声水平、正反演模型等因素的影响.Oldenborger[3]提出MRM 可以用来观测系统误差所导致的系统退化效应.本文进一步分析了MRM 矩阵与系统误差之间的关系,并建立了系统误差与MRM 的关系模型.

反褶积算法被广泛地应用于地球物理参数模型的计算和校正[16-18].在类似的相关领域,反褶积方法也有广泛的应用,例如:图像与视频的增强[19]、天文数据的增强[20]、遥感影像的增强[21],以及生物医学图像的增强[22]等.此外,盲反褶积技术还被成功地应用在层析成像数据的增强中[23-24].层析成像模型计算过程也存在由于系统误差所导致的系统退化效应,并受不同的物理仪器和不同的模型重建算法的影响.反褶积方法可以很好地消除层析成像方法中系统误差所带来的影响.本文根据相关领域对于系统误差、系统退化效应的研究,以及地球物理模型反演理论,提出了地球物理反演参数褶积退化模型,以及地球物理反演参数模型的反褶积增强算法.方法可以有效地消除估计模型中存在的系统误差,提高地球物理反演模型的分辨率,为后期资料的解释提供更为清晰和直观的参考.

2 模型分辨率矩阵与退化点扩展函数 2.1 模型分辨率矩阵

系统误差存在于观测数据中,本文忽略观测数据采集的具体过程和方式,直接讨论含有系统误差的采集数据的处理.将非线性模型估计分析进行线性化近似,定义如下[3]:d=F(m),其中:F(·)表示正演算子,d表示观测数据向量,m是模型向量,假设反演为病态欠定问题,则构建正则化最小二乘目标函数如下:

(1)

其中μ 为正则化参数,Wd 为观测数据的加权矩阵,Wm 为参考模型mref 的加权矩阵.线性化式(1)并将其进行Taylor展开,并忽略高阶项ο(‖Δm2).则迭代更新估计模型为:

(2)

其中S为敏感度矩阵,H为Hess 矩阵,H=[STWdTWdS+μWmTWm].为了建立估计模型和真实模型的关系,将观测模型和真实模型$\hat{m}$的关系表达为:d=F($\hat{m}$)+eS+eTeS 是测量中的随机误差项,eT 是离散化误差和计算模型误差引起的系统误差项,在非线性情况下,归并误差项eS+eT=e[25-27],得到

(3)

其中:mb 代表无偏估计参数模型,分辨率矩阵R可定义为:

(4)

由于式(3)中叠加在真实模型$\hat{m}$上的RRS-1e的影响,导致在无偏估计模型mb 中出现了退化效应.在地球物理反演过程中,受噪声和正则化项等因素的影响,通常R的数值分布并不完全规则,而且,在多维情况下,R是一个相对规模较大的矩阵,直接计算R所需的空间复杂度和实践复杂度过大,硬件无法满足需求.

因此,对R的数值分布进行了进一步的讨论,以简化退化数学模型.

2.2 点扩展函数

将模型分辨率矩阵R按行和列,划分为行向量al和列向量pk,则模型分辨率矩阵可以表示为M个向量:

(5)

列向量pk被称为PSF,PSF 数值分布受观测噪声、反演中加入的先验信息等因素影响,是加载在真实模型上的类delta扰动函数.如果PSF 有较宽和较大的旁瓣存在,估计模型的分辨率则相应较低[28].

3 点扩展函数褶积近似模型

理想的情况下,R是一个单位矩阵,此时模型的分辨率达到最大值,是对真实模型的无误差反映.但在实际应用中,受系统误差因素,例如离散化误差、正演模型误差、以及受压制噪声稳定解而使用的正则化项的影响,通常估计模型的分辨率矩阵的主对角线存在一定宽度的旁瓣,因此,无偏估计模型无法达到最大分辨率的理想状态.

为从理论上推导系统退化模型,以及分析提出的近似褶积模型的误差,首先将矩阵R扩展为Teoplitz矩阵,将式(3)、式(4)中模型分辨率矩阵与真实模型的乘积退化过程,近似为点扩展函数褶积退化过程.

式(5)中,行向量al称为核函数,在空间域上,核函数al是一个加权平均函数.核函数al与真实模型$\hat{m}$的乘积等于无偏估计模型参数mb.因此,拓展后的分辨率矩阵R与真实模型$\hat{m}$(l=1,…,M)的乘积过程,在空间域上是一个加权平均过程,反映了由系统误差而引起的系统退化效应.

根据式(5)定义,MRM 矩阵中的PSF 存在类似于核函数al的系统退化作用.在分辨率矩阵中,PSF的数值通常沿对角线对称分布,假定PSF 支持域大小为k,则

(6)

式中,pij的数值以矩阵R的对角线为中心对称分布.假设:矩阵主对角线附近,在列方向上不为零的元素个数为N个.将矩阵R在顶部向上延拓N/2行,矩阵R底部向下延拓N/2行,则矩阵大小变为(M+N)×M.将延拓前矩阵R第1 行到第N/2行,第N/2列到第N列之间正方形区域的元素值,对应赋值给拓展后矩阵的第1行到第N/2行,第1列到第N/2列之间正方形区域的元素.同理,对矩阵底部下延拓N/2 行,赋值过程类似.具体的上延拓、赋值过程如图 1所示.

图 1 MRM 的Teoplitz矩阵拓展近似 Fig. 1 Teoplitz matrix expansion approximation of MRP

图 1表示了模型分辨率矩阵的向上拓展、赋值,近似Teoplitz矩阵的过程.黑色小方格代表矩阵中数值不为零的元素,白色小方格代表数值为零的元素.其中黑色虚线组成的方框代表向上拓展后增加的部分,白色虚线表示未拓展矩阵的主对角线,两个红色正方形包围的区域,表示对应的赋值区域.以图中单个元素示例:将A 位置的元素值赋值给A′位置的元素,将B 位置的元素值赋值给B′位置的元素.以此类推下拓展、赋值过程,就完成了整个模型矩阵的拓展和赋值.

拓展后的矩阵Rexp是一个标准Teoplitz矩阵.根据相关的褶积定理,Teoplitz矩阵与一个向量的乘积,等于一个褶积点扩展函数与向量的卷积.因此,Rexp与真实模型向量$\hat{m}$的乘积,等于褶积点扩展函数与真实模型向量$\hat{m}$的褶积.褶积后的矩阵大小为(M+N,1).为使褶积数据长度与mb 相符,截取数据中心区域范围M长度的数据.

由于N$\ll $M,近似Teoplitz矩阵的误差和褶积有限长度所带来的误差可以忽略不计,得到近似关系(其中*表示部分褶积):

(7)

其中PSFe 表示褶积点扩展函数,PSFe=pk.根据式(3)和式(7),可以得到:

(8)

以上的推导将R近似为一个褶积矩阵PSFe(褶积点扩展函数),褶积点扩展函数的规模远小于R的.从计算效率和计算精度上考虑,褶积退化数学模型更为适合模型增强计算.在以上的推导过程中,假定退化点扩展函数为空不变的,这样假设主要依据为:首先,在以往的退化增强的研究中[19],通常假设点扩展函数是空不变的,这种假设所带来的误差和增强效果能满足实际需求,算法的计算效率也较高,空不变假设有可能简化了实际的物理模型,但空不变模型通常可以近似代表实际的物理模型.此外,这种假设具有一定的普遍性.虽然空变的点扩展函数较为复杂且更为准确,对于某种具体地球物理的应用,如果可以估计出其空变的点扩展函数,采用其进行增强,相比空不变点扩展函数,其增强的精度更高,但是对于更广泛意义上的实际应用,不存在普遍意义上的较为稳定的空变点扩展函数.因此,在更广泛的适用范围上,空不变的假设具有一定的普遍性.以上讨论是以一维模型近似为例,二维模型的Teoplitz矩阵近似与一维近似过程相似,二维模型的分辨率矩阵相对一维规模较大,数值分布相似,同样可以近似为一个Teoplitz矩阵,根据Teoplitz矩阵估计得到二维的褶积点扩展函数,具体过程与一维情况相同,本文不再单独论述.

此外,二维数据试验表明,通过试验直接求取的分辨率矩阵通常是不规则的,难以从数值上严格地近似为Teoplitz矩阵.由于误差的影响,直接从分辨率矩阵近似褶积点扩展函数,会导致褶积点扩展函数存在较大的误差,而褶积点扩展函数误差会在反褶积模型增强过程中被放大,严重影响增强效果.此外,由于还存在其它的叠加系统褶积退化效应,其无法通过点扩展函数来估计.因此,从试验结果可以看出,直接从分辨率矩阵中估计出总体的系统退化点扩展函数是不可行的.

另一个误差的不确定性影响是在系统误差中包含的随机误差成分,以类似于噪声叠加的方式存在,且通常相对数值较小[3].因此,本文将RS-1e近似为随机噪声项n.

除根据分辨率矩阵推导出的退化褶积点扩展函数外,一些其它的地球物理过程,例如数据采集、器件退化等,也会产生一系列类似的级联的褶积退化.可以将所有的褶积退化表示为系统级联退化过程,如式(9)所示.

(9)

根据式(9),地球物理无偏估计模型mb,可以表示为级联的卷积退化所形成的整体褶积退化函数PSF褶积真实模型而形成.因此,可根据褶积退化过程,对地球物理无偏估计模型进行反褶积.从无偏估计模型中恢复出更多真实模型的信息,减小估计模型的系统误差.

4 反演估计模型的盲反褶积逆滤波增强 4.1 混合范数总变分正则化盲反褶积算法

根据以上的分析,本文提出新的一种混合范数的TV(TotalVariation)正则化的盲反褶积算法,对反演估计模型进行逆滤波增强.

定义总变分正则化模型如下:

(10)

其中,Δ (·)表示后向差分运算,λ 是Lagrange乘子,m表示求解估计模型.R${\hat{m}}$≈PSF*${\hat{m}}$,根据式(9)则可以得到:

(11)

其中·表示点积,FFT(·)表示离散傅里叶变换.则模型(10)可以表达如下:

(12)

Dx表示x方向的梯度模板,Dy表示y方向的梯度模板.

研究表明[18, 29],混合l1l2 范数总变分正则化的逆滤波算法对模型的边缘具有较好的保持性和抗噪性能,要好于l2 范数总变分正则化的逆滤波算法.因此,在总变分正则化项范数的选择上,采用混合l1l2 范数,则总变分项定义如下:

算法为了更精确、更有效地估计总变分正则化项,在原有总变分框架(式(10))的基础上,引入新变量w,v.使用引入变量w近似梯度项‖Δm1v近似梯度项‖Δm2.‖·‖=‖·‖2.建立了一个新的凸的目标函数模型:

(13)

将近似惩罚项加入目标函数中则(13)可以写成:

(14)

其中β 是惩罚系数.式(14)是一个凸问题.明显,在惩罚系数β 趋近于无穷大时,式(14)收敛为式(10).要求其中两项: ‖w- ‖Δm12 与‖v-‖Δm22 的近似误差都足够小,以满足式(14)整体的误差要求.对于式(14)的解所对应的惩罚系数β 越大,项$\int\limits_{\Omega }{\left\.\omega - \right\.\nabla m{{\left\._{1} \right\.}^{2}}dx+}\int\limits_{\Omega }{\left\.\upsilon - \right\.\nabla m{{\left\._{2} \right\.}^{2}}dx}$值越小,式(14)近似式(10)的精度越高.

4.2 反褶积算法迭代求解

采用迭代最小化的方法来求解式(14),由于真实模型${\hat{m}}$和褶积点扩展函数PSF未知,所以需要两个分步骤来交替迭代求解.假设近似估计模型m已知情况下,求解点扩展函数PSF,这个步骤称为PSF步骤.相反,假设近似估计的点扩展函数PSF 已知情况下,求解估计模型m,这个步骤称为m 步骤.由于对总变分项进行了近似,因此在计算估计模型的m步骤中需要估计总变分近似项wi,j,vi,j,所以在m 步骤中嵌套了一个对wi,j,vi,j进行估计的wv步骤.式(14)是一个最小化凸问题.利用总变分方法,最小化式(15)的问题转化为Euler-Lagrange带有Neumann条件方程.对其求偏导,省略推导步骤,得到解的形式如下:

m 步骤的解为:

(15)

其中为垂直和水平方向的差分算子.

PSF步骤的解为:

(16)

其中κ 是一个大于零的小数,表示输入信号的功率谱,Snn表示噪声功率谱.由于噪声估计比较困难,因此κ 不能被精确地确定大小,通常采用一个小值来近似.

wv步骤的解为:

(17)

其中:τ 是为避免除零而设置的一个相对小数,一般设置为[0.5,1]范围内的一个值,其大小不影响求解结果.

算法设计从一个较小的值β开始,逐步地增加β的大小,把当前的β 值估计的近似解作为下一个β值递增前模型未知项估计的初始解.从一个较小β估计值估计的开始,迭代估计未知项的大小,并迭代递增β 值.对于快速求解情况,追求最短的程序运行时间,可设β=βmax,并通过试验选取一个合适的β值,控制模型的细节保存程度和模型中噪声的平滑程度.

5 试验分析

本文根据以上的理论分析和提出的增强算法,分别进行了1D 和2D 的仿真试验,1D 的线性仿真试验主要验证系统退化褶积近似理论的正确性,2D非线性大地电磁反演试验验证提出的TV 正则化逆滤波增强算法的有效性.

5.1 一维线性反演增强

首先,以一维线性反演问题为例:di= (fi,m),其中fi是高斯模糊退化核函数,如图 2b 所示(从100个核函数中等间隔选取10 个核函数显示),采用一个1D 长度为100的合成数据模型,加载5%高斯随机噪声后观测数据、真实模型以及核函数如图 2所示.

图 2 真实模型(实线)、观测数据(虚线)(a)和核函数(b) Fig. 2 (a) True model (solid line) and observed data (dashed line).(b) Kernel function

对观测数据使用均值为0.5的参考模型进行反演,反演方法采用带Tikhonov 正则化的TSVD(TruncatedSingularValueDecomposition)算法[30],反演中的正则化参数的确定使用L-Curve算法.真实模型和反演结果对比如图 3所示.

图 3 真实模型(实线)和反演模型(虚线) Fig. 3 True model (solid line) and inversion model (dashed line)

受噪声的影响,反演中存在病态性问题,本例中,正则化系数α=0.54117,从图 3中可以看出:真实模型的跳变边缘在反演后变得较为平滑,估计模型的分辨率较真实模型有所降低,反演模型中出现了由于噪声等因素引起的数值扰动效应.

根据已知的核函数建构模型分辨率矩阵如图 4所示.根据提出的模型分辨率矩阵PSF 褶积近似理论,对反演估计的模型分辨率矩阵进行PSF 褶积近似,根据式(10)的推导,从图 4模型分辨率矩阵中提取褶积退化PSF,如图 5所示.

图 4 模型分辨率矩阵 Fig. 4 Model resolution matrix
图 5 褶积退化点扩展函数 Fig. 5 Convolution degradation point spread function

为分析矩阵近似过程对模型精度的影响,将得到的分辨率矩阵拓展为Teoplitz矩阵,从而得到褶积点扩展函数PSF.将分辨率矩阵与真实模型${\hat{m}}$的乘积,与近似褶积PSF和真实模型${\hat{m}}$的褶积结果对比,如图 6所示.

图 6 fi${\hat{m}}$(黑色虚线)与PSF*${\hat{m}}$(灰色线)的对比(表示褶积) Fig. 6 Comparison of fi${\hat{m}}$(black dashed line) and PSF * ${\hat{m}}$ (gray line) ( * means convolution)

图 6中,核函数与真实模型乘积的结果(fi${\hat{m}}$),与点扩展函数和真实模型褶积PSF*${\hat{m}}$的数值曲线基本重合.但根据第3节的理论推导可知,由于褶积近似过程和褶积有限长度的影响,因此,fi${\hat{m}}$与PSF*${\hat{m}}$ 之间存在误差距离,将两者相减,结果如图 7所示.

图 7 fi${\hat{m}}$与PSF*${\hat{m}}$的误差 Fig. 7 Difference between fi${\hat{m}}$ and PSF * ${\hat{m}}$

fi${\hat{m}}$与PSF*${\hat{m}}$误差绝对值的最大值为0.0011,与fi${\hat{m}}$的最大值0.7520相比,误差为0.15%,可以认为褶积近似误差较小,在反演增强中可以忽略.

针对本例图 3所示的反演模型,根据提出的TV正则化增强算法对反演结果进行增强,TV 算法使用的初始PSF如图 5所示,正则化因子:α=7.3564×10-4,增强结果如图 8所示.

图 8 真实模型(黑色线),反演模型(虚线),反演增强模型(灰色线) Fig. 8 True model (black line) ,inversion model (dashed line) and enhanced inversion model (gray line)

图 8中,反演估计模型经过TV 算法增强后,真实模型的陡峭的边界细节得到了恢复,其更接近于真实模型,同时反演估计模型中零值附近的数据扰动在增强后也得到了很好的抑制.估计参数模型增强后与真实模型的平均距离误差减小了58.50%,同时分辨率得到了提高.

5.2 二维非线性反演增强

以2D MTOccam 反演为例,真实模型采用合成模型如图 10所示:在地表下5~10km 深度、水平方向-5到5km范围内有一个10Ωm的低阻体,周围围岩为100Ωm 的高阻体.频段范围从0.5Hz~0.2×10-2 Hz,观测数据加入了10%的高斯随机噪声,采用不均匀网格划分方法,在异常体附近和靠近地面位置网格划分较密,参考模型电阻率值为100Ωm.

根据反演过程中噪声的偏差,设置反演正则化因子值为0.32.经过10次反演迭代后,得到的估计模型如图 10所示.反演后周围高阻围岩的电阻率整体降低,而低阻体的阻值整体升高.并且,低阻体的范围扩大、边界模糊,接近一个类椭球体.

图 10 反演估计模型 Fig. 10 Inversion estimated model

图 11,深度5km 到10km,宽度-5km 到5km的低阻异常体,其反演估计模型受各种噪声和误差的影响,数值分布由真实模型的矩形分布变为类高斯状分布.

图 11 异常体反演后的数据分布 Fig. 11 Value distribution of anomalous areas after inversion

反演估计模型最低电阻率值为36.981Ωm,大于真实模型的数值10Ωm.异常区域内平均电阻率值为57.981 Ωm,相比真实模型的平均电阻率值10Ωm误差较大.估计模型在数值分布连续区域,如高阻异常体精度较高,受系统退化误差的影响较小,而在估计模型数值变化区域,如高低阻交界处,存在较为严重的系统退化效应.

本例反演的模型分辨率矩阵是一个较大的稀疏矩阵,由于噪声和网格划分、以及反演算法误差的影响,估计出的反演模型分辨率矩阵数值分布形式较为复杂,难以从中近似估计出严格的2D PSF 褶积阵.对图 10所示的反演估计模型,采用本文提出的盲反褶积TV 正则化增强算法进行增强,增强的反演模型如图 12所示.

图 12 增强的反演模型 Fig. 12 Enhanced inversion model

估计模型经过反褶积算法增强后,低阻异常区域的边界变得更为清晰,异常区域的形状更接近真实模型,从视觉上更利于数据的解释和分析,如图 13所示.

图 13 异常体区域增强后的数据分布 Fig. 13 Value distribution of anomalous areas afterenhancement

对比图 9图 10以及图 12增强后的模型,异常区域内电阻率值最小值为12.983Ωm,相比反演估计模型36.981Ωm,更接近真实模型电阻率值10Ωm.同时,增强后异常区域内的平均电阻率值为15.484Ωm,相比反演模型平均值57.981Ωm,误差减小了88.57%.

图 9 合成模型 Fig. 9 Synthetic model

总体上,反演估计模型的数值范围为(36.981,117.119)Ωm,反演增强估计模型的数值范围为(12.983,114.968)Ωm,更接近真实模型范围(10,100)Ωm,增强后模型数值精度得到提高.

以7.5km 深度的各水平位置模型的数值分布为例,对比真实模型、估计模型和反演增强模型,如图 14所示.

图 14 7.5km 深度水平位置模型对比 真实模型(黑线),反演模型(蓝色虚线),反演增强模型(红线) Fig. 14 Comparison of the horizontal position models at 7.5 km depth True model (black line) ,inversion model (blue dashed line),enhancement model (red line)

增强后模型中的异常体的边界相对增强前更为清晰,模型的视觉分辨率得到了提高.真实模型与反演估计模型电阻率的距离加和均值为14.0997Ωm,真实模型与增强模型电阻率的距离加和均值为12.7991Ωm.由此可得,在7.5km 深处各水平位置的反演估计模型,经过增强后,误差减少9.22%.

对比水平方向0km 处,真实模型、估计模型和反演增强模型在各个深度的数值分布,如图 15所示.

图 15 水平方向0km 处各深度模型对比 真实模型(黑线),反演模型(虚线),反演增强模型(灰线) Fig. 15 Comparison of the depth model at 0 km horizon position True model (black line),inversion model (dashed line),enhancement model (gray line)

水平方向0km 处各深度的估计模型与真实模型电阻率的距离均值为13.9074Ωm,真实模型与增强模型电阻率的距离均值为10.7118Ωm.水平方向0km 处各深度反演估计模型,经过增强后,误差减少22.98%.

图 14图 15 中,在高阻围岩部分,由于反演误差引起的模型数值扰动,在经过反褶积增强后,得到了很好的平滑,增强结果与真实数据的变化趋势更为相似.

总体上,真实模型与反演模型电阻率的距离加和均值为10.2412Ωm,真实模型与增强反演模型电阻率的距离加和均值为9.8809Ωm.经过反演增强后,反演模型整体上误差减少3.5%.由于区域内大量数据属于连续的高阻体,反演后整体的系统退化误差并不大,系统退化主要出现在高阻和低阻交界的异常区域.由于低阻异常模型区域所占整体模型区域的比重较小,所以整体上精度提高了3.5%,提高精度幅度较小.但对于异常体区域的模型,增强后模型的精度提高88.57%.所以,方法对于估计模型中存在较多系统退化效应的模型数值跳变区域,可以更为显著地减少系统退化效应,提高模型的分辨率和数值精度.

6 结论

文章根据地球物理反演理论提出了一个新的反演估计参数的近似退化褶积模型,以及基于该模型的盲反褶积增强算法.通过1D 和2D 的反演增强试验,文章验证了褶积近似PSF 理论的正确性,以及提出的反褶积增强算法的有效性.试验表明,反演模型在增强处理后误差整体减小,特别是对于数值跳变的异常区域,可以较好地增强异常体边界,减少系统误差所带来的系统退化效应,有效地增加模型分辨率和数值精度.

参考文献
[1] Ganse A A. A Geophysical Inverse Theory Primer. Princeton: Princeton University Press, 2008 .
[2] Alumbaugh D L, Newman G A. Image appraisal for 2-D and 3-D electromagnetic inversion. Geophysics, , 2000, 65(5): 1455-1467. DOI:10.1190/1.1444834
[3] Oldenborger G A, Routh P S. The point-spread function measure of resolution for the 3-D electrical resistivity experiment. Geophys.J.Int. , 2009, 176(2): 405-414.
[4] 于波, 翟国君, 刘雁春, 等. 噪声对磁场向下延拓迭代法的计算误差影响分析. 地球物理学报 , 2009, 52(8): 2182–2188. Yu B, Zhai G J, Liu Y C, et al. Analysis of noise effect on the calculation error of downward continuation with iteration method. Chinese J.Geophys. (in Chinese) , 2009, 52(8): 2182-2188.
[5] Tikhonov A N, Arsenin V I A. Solution of Ill-Posed Problems. Washington: W.H.Winston and Sons, 1977 .
[6] ŽdanovM S. Geophysical Inverse Theory and Regularization Problems. New York: Elsevier Science, 2002 .
[7] Backus G E. Inference from inadequate and inaccurate data, I. Proceedings of the National Academy of Sciences , 1970a, 65(1): 1-7. DOI:10.1073/pnas.65.1.1
[8] Backus G E. Inference from inadequate and inaccurate data, II. Proceedings of the National Academy of Sciences , 1970b, 65(2): 281-287. DOI:10.1073/pnas.65.2.281
[9] Backus G E. Inference from inadequate and inaccurate data, III. Proceedings of the National Academy of Sciences , 1970c, 67(1): 282-289. DOI:10.1073/pnas.67.1.282
[10] Backus G E, Gilbert J F. Numerical applications of a formalism for geophysical inverse problems. Geophys.J.Roy.Astr.Soc. , 1967, 13(1-3): 247-276.
[11] Jackson D D. Interpretation of inaccurate, insufficient and inconsistent data. Geophys.J.Roy.Astr.Soc. , 1972, 28(2): 97-109.
[12] Menke W. Geophysical Data Analysis: Discrete Inverse Theory. Orlando: Academic Press Inc, 1984 .
[13] Daily W, Ramirez A. Electrical resistance tomography during in-situ trichloroethylene remediation at the Savannah River Site. Journal of Applied Geophysics , 1995, 33(4): 239-249. DOI:10.1016/0926-9851(95)90044-6
[14] Miller C R, Routh P S. Resolution analysis of geophysical images: Comparison between point spread function and region of data influence measures. Geophysical Prospecting , 2007, 55(6): 835-852. DOI:10.1111/gpr.2007.55.issue-6
[15] Oldenburg D, Li Y. Estimating depth of investigation in dc resistivity and IP surveys. Geophysics , 1999, 64(2): 403-416. DOI:10.1190/1.1444545
[16] Yu J H, Hu J X, Schuster G T, et al. Prestack migration deconvolution. Geophysics , 2006, 71(2): S53-S62. DOI:10.1190/1.2187783
[17] 张季生, 高锐, 李秋生, 等. 欧拉反褶积与解析信号相结合的位场反演方法. 地球物理学报 , 2011, 54(6): 1634–1641. Li G S, Gao R, Li Q S, et al. A combined Euler and analytic signal method for an inversion calculation of potential data. Chinese J.Geophys. (in Chinese) , 2011, 54(6): 1634-1641.
[18] 薛国强, 李貅, 戚志鹏, 等. 瞬变电磁拟地震子波宽度压缩研究. 地球物理学报 , 2011, 54(5): 1384–1390. Xue G Q, Li X, Qi Z P. Study of sharpening the TEM pseudo-seismic wave-form. Chinese J.Geophys. (in Chinese) , 2011, 54(5): 1384-1390.
[19] Almeida M S C, Almeida L B.Blind deblurring of natural images.//Proceedings of IEEE International Conference on Acoustics, Speech and Signal (ICASSP).Las Vegas, NV: IEEE, 2008: 1261-1264.
[20] Bertero M, Boccacci P. Image restoration methods for the large binocular telescope (LBT). Astron.Astrophys.Suppl. , 2004, 147: 323-333.
[21] Müller J P. Digital Image Processing in Remote Sensing. New York: Taylor & Francis, 1988 .
[22] Taxt T. Three-dimensional blind deconvolution of ultrasound images. IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control , 2001, 18(4): 867-871.
[23] Villain N, Goussard Y, Idier J, et al. Three-dimensional edge-preserving image enhancement for computed tomography. IEEE Transactions on Medical Imaging , 2003, 22(10): 1275-1288. DOI:10.1109/TMI.2003.817767
[24] Jefferies S, Schulze K, Matson C, et al. Blind deconvolution in optical diffusion tomography. Opt.Exp. , 2002, 10(1): 46-53.
[25] Bronstein A M, Bronstein M M, Zibulevsky M, et al.Quasi-maximum likelihood blind deconvolution of images acquired through scattering media.// Proceedings of EEE International Symposium on Biomedical Imaging: Nano to Macro (ISBI).Arlington, VA: IEEE, 2004: 352-355.
[26] Oldenborger G A, Routh P S, Knoll M D. Sensitivity of electrical resistivity tomography data to electrode position errors. Geophys.J.Int. , 2005, 163(1): 1-9.
[27] van Wijk K, Scales J A, Navidi W, et al. Data and model uncertainty estimation for linear inversion. Geophys.J.Int. , 2002, 149(3): 625-632.
[28] Friedel S. Resolution, stability and efficiency of resistivity tomography estimated from a generalized inverse approach. Geophys.J.Int. , 2003, 153(2): 305-316.
[29] Fu H Y, Ng M K, Nikolova M, et al. Efficient minimization methods of mixed l2-l1 and l1-l1 norms for image restoration. SIAM J.Sci.Comput. , 2006, 27(6): 1881-1902.
[30] Hansen P C. Truncated singular value decomposition solutions to discrete ill-posed problems with ill-determined numerical rank. SIAM J.Sci.and Stat.Comput. , 1990, 11(3): 503-518.