基于广义交叉验证的电测深曲线一维全自动迭代反演

施昕祎，刘海飞，柳建新，李 星，汪强强

(1.中南大学 地球探测与信息物理学院，长沙 410083；2.有色金属与地质灾害探查湖南省重点实验室,长沙 410083)

0 引言

Hadamard[1]在研究偏微分方程解的稳定性时提出对于给定的偏微分方程与定解条件，定解问题中解的存在性、唯一性和稳定性是理论上的三个基本问题，并将这三个基本问题统称为定解问题的适定性。之后前苏联学者Tikhonov[2]与其团队对不满足上述适定条件的定解问题(即反问题)的处理方法进行了深入研究，提出了沿用至今的Tikhonov正则化方法。随后国外学者Engl[3]和Hansen[4]在正则化方法的理论和算法方面做了相关研究，给出了可用的线性和非线性不适定问题正则化方法的结果并提出奇异值分解算法作为不适定线性最小二乘问题的正则化方法，如Miller正则化和Tikhonov两大类正则化方法，后者在理论上已经发展成为较完善的理论体系。我国在反问题方面的研究起于1980年，国内学者开始了数学领域关于反问题的正则化处理方法研究。栾文贵[5]首先将正则化方法应用于地球物理反演中，通过对正则化后的泛函取极小值实现反问题的正则化最小二乘处理，当采用Tikhonov正则化时，正则化因子的选取会对反演结果造成重要影响;阮百尧[6]实现了激电测深一维最优反演。当初始模型层数的选取不适当时会导致反演结果较差，且初始模型的层厚和层阻也需要手动输入。当线性反演方程组病态程度高时，用一般解方程组方法不能取得稳定解，特别是在求解电磁法二维、三维反演问题时，解病态方程组是不可避免的。因此求解大型、超定、病态的方程组成为反演研究的首要问题。数值计算学家Forsythe教授在用奇异值分解的正则化算法解决病态问题方面做出了较好应用，基于奇异值分解算法来选择正则化因子，奇异值算法的关键在于适当的修改奇异值(即施加正则化)来增强求解稳定性[7-11]。戴亦军[12]为了提高野外大地电磁测深数据的处理效率和初步解释的精度,提出了大地电磁测深数据的一维正则化反演进行拟二维反演解释方法,用Bostick反演的深度来控制层参数,采用阻尼高斯-牛顿算法进行反演计算,并将Bostick反演结果作为反演计算的初始模型，得到的反演结果更加稳定;朱肖雄[13]采用最小二乘正则化反演迭代,实现自然电场场源的二维反演,得到了收敛稳定的满意结果;张彭[14]以长导线电性源半航空电磁法为例,将自适应正则化反演算法应用于半航空时间域电磁数据反演中,给出最平缓模型约束条件下的半航空时域电磁数据自适应正则化反演算法，对典型地电模型理论正演加噪数据进行并行反演计算，减少计算时间提高效率。笔者在前人的基础上基于奇异值分解算法(Singular value decomposition algorithm，SVD)，利用奇异值分解法先将病态线性方程组进行分解，然后采用正则化技术削弱小奇异值的影响，可以使求得的解更好地逼近真实解，能够达到有效求解病态方程组的目的。综合奇异值分解法与广义交叉验证自适应正则化技术求解线性反演方程，实现根据极距数自动构建反演初始模型并自适应迭代反演计算，通过典型地电模型验证了算法的准确性和稳定性。

1 电测深完全自动迭代反演

1.1 基本原理

由于在建立的地电模型中，选择层数自动为观测到的极距数，进而模型参数多于实测数据，这会导致线性方程组欠定，以至于无法对其进行求解。在野外实际数据处理中，超定、欠定问题是常见的，因此构建目标函数方程的韧性要充分估计到。为得到较优的解估计并增强反演的稳定性，通常在模型空间引入稳定化泛函(即对模型参数施加先验约束)，施加光滑模型约束条件在目标函数中。目标函数φ被施加了局部光滑约束条件为式(1)。

φ=‖Δd-AΔm‖2+λ‖CΔm‖2

(1)

(2)

将式(1)两端对ΔmT求导，并且使它等于零，得到方程：

(ATA+λCTC)Δm=ATΔd

(3)

式(3)的最小二乘解等价形式为式(4)。

(4)

使用最小二乘的正则化方法解超定方程组(4)，求得模型修正量Δm:

m(k+1)=m(k)+Δm(k)

(5)

为便于进行偏导数矩阵计算的参变量的无量纲化处理此处不使用取对数处理。

施加光滑约束的线性反演方程构建完毕之后，将在方程组的求解过程中引入奇异值分解法与广义交叉验证自适应正则化方法。

1.2 计算正则化因子的广义交叉验证法

在求解构建好的线性反演方程时应尽量选择稳定性好抗病态能力强的解法。奇异值分解算法具有抗病态能力强的优点且无需判定矩阵A是否对称正定，是反演中解病态线性方程组的最有效的方法之一，在地球物理反演中得了较好地应用。如何修改奇异值是使用关键，为了达到增强方程组解稳定性的目的，需要修改奇异值。只有奇异值得到恰当修改，才能使奇异值分解算法达到降低方程病态程度的要求。

对于线性方程组Ax=b,经过奇异值分解算法分解处理后形式可变为A=USVT(U=[u1,u2,…,un]，V=[v1,v2,…,vn]，S为递减对角阵)，解向量正则化表达式为式(6)。

(6)

其中：fi为滤波函数。广义交叉验证法就是 Golub[15]由统计学中提出的一种对SVD法的正则化手段，它是Tikhonov滤波正则化压制矩阵奇异性方法的其中之一。

正则化方法削弱矩阵的奇异性，是通过采用零阶Tikhonov正则化方式来实现的，则式(6)中滤波函数fi可表示为为式(7)。

(7)

通过改变λ的值可以使滤波系数fi从“0”到“1”光滑变化。将奇异值分解后存在的较小的奇异值压制掉，需选择合理的正则化因子λ(即用于修改小奇异值的正则化因子λ)加入矩阵对角线上(即采用零阶Tikhonov正则化方式)。笔者采用广义交叉验证法来给定最优正则化参数λ。

Golub[15]提出了广义交叉验证法(generalized cross-validation,GCV)，该方法用于估计正则化参数。其思想来源于统计学，当在数据项b中去掉一个分量bi后，则由此产生的新模型解也能较好的预测b中被去掉的那一个分量bi。基于这一思想。Golub给出了关于参数λ的广义交叉验证函数为式(8)。

(8)

其中：tr为方阵对角线各元素的和，即矩阵的迹。使GCV(λ)函数达到极小时λ的值就是“最优”的正则化参数。从式(7)中可以看出，分子作为正则解的残差是比较容易计算的。而方阵的迹采用直接方法进行计算不仅困难而且计算量比较大，但根据矩阵A的奇异值分解形式，可以较方便的得出迹估计如下:

(9)

通过给定一系列的λ值，“最优”的正则化参数就是使GCV(λ)函数取极小值时对应的参数λ值。选定奇异值分解法与广义交叉验证正则化方法后，要选择合适的初始模型来进行反演才能进一步提升反演效果。

1.3 自动构建反演初始模型

为了达到反演过程中无需手动输入初始模型参数的目的，实现了一种根据反演数据自动化构建的多层等效模型构建方法。在多层等效模型中，层数就是电测深曲线的极距数,层厚则是相邻电极距的差,层阻为对应极距的视电阻率。由于电测深浅部分辨率高、深部分辨率低的特点，所以采用这种方式创建模型。层厚有极距差表示恰好与电测深曲线的纵向分辨规律相符合，而层阻由对应极距的视电阻率表示，可以确保等效模型的电性参数与电测深曲线在深度方向上具有相同的变化规律。除此之外，由于电测深曲线所记录的视深度(AB/2)总是大于所反映地层的真实深度，故在反演之前，初始层厚要乘以一个系数c(c<1,如0.8)，以改善初始模型从而达到加快收敛速度的目的。而且只有确保电测深曲线的首、尾支完整，才能自动构建多层等效模型，否则很大程度上会使分层与真实情况不符，导致反演结果不准确。考虑到实际电测深方法的分辨特点，在反演过程中以极距数自动构建初始模型可以避免任意设置初始模型层数带来的人为不确定性并且可在一定程度上改善反演质量使其尽可能接近真实地质条件。

1.4 计算偏导数矩阵

完成反演初始模型的构建后，还需解决参数反演问题中的一个难点，即计算观测数据关于模型参数的偏导数矩阵。阮百尧[6]用差分方法计算偏导数矩阵，通过无量纲化偏导数矩阵、解向量以及方程右端项中的元素，解决了视电阻率数据和电阻率模型参数变化范围大且参变量量纲不同的问题,前人曾用它们的对数值来进行计算,但是这并没有解决厚度模型参数与电阻率模型参数量纲不同问题,因此在反演过程中,出现了视电阻率测深曲线形态虽然拟合很好,但预测模型与真实模型参数数值相差很大的情况，故使用差分法来实现电阻率测深一维最优化反演。此方法在初始模型层参数不是很准确的情况下，均能既快又稳的收敛到最优解。利用差分法计算偏导数矩阵的过程如下：

取模型改正量Δmj=0.1mi，则视电阻率对模型参数得偏导数可近似为:

ρci(m1,…,mj,…,mNM)]/0.1mj

(10)

通过对式(10)右端第一项中Δm求偏导，使其等于零，经整理得:

(11)

两端同时除以ρci(m1,…,mj,…,mNM),使其两端的向量元素变为便于处理的无量纲化令：

(12)

(13)

Bi=ρai/ρci-1

(14)

图1 反演流程图Fig.1 Inversion flowchart

式(11)可以写成矩阵形式:

近10 a来，许多国家尤其是我国进行了大量的探地雷达(Ground Penetratin Rador，简称GPR)应用于铁路线路翻浆冒泥病害的相关检测与试验[1-5]，GPR勘探方法开始在铁路线路可见和隐伏的翻浆冒泥病害的探测与圈定应用中被广泛地接受[1，6]，已经成为铁路线路翻浆冒泥病害检测、圈定的重要手段。

Ax=b

(15)

上述方程组中,左端系数矩阵A中各系数Aij与未知向量x中各变量xi以及右端向量b中各系数Bi都无量纲,从而解决了参变量量纲不同的问题。至此便完成了基于广义交叉验证的电测深曲线一维全自动迭代反演相关基本理论的阐述，下面将对反演过程做简要说明。

1.5 反演的计算流程

由于反演过程无需输入初始模型参数，使得整个反演过程实现了全自动化，其反演流程如下：

1)将同一断面的激电测深曲线保存在文件中作为反演输入文件。格式为两列，其中第一列第一行为总极距数(AB/2)，后续第二行第一列为从小到大递增的极距列表，第二列则为与相应极距对应的视电阻率数值。

2)读取数据文件，程序根据读取的极距与视电阻率进行深度修正系数自适应选择，构建相应反演初始模型，由GCV法进行用于每次迭代中修改奇异值算法的正则化因子自适应选取，计算偏导数矩阵得到模型修改量后用新模型正演并计算拟合差，当RMS小于5%并随着迭代次数不在变化时终止迭代。

3)全自动迭代反演计算结束，将反演结果保存为三个dat类型文件，广义交叉验证函数数据、实测和模拟数据拟合信息、最终模型与反演结果，利用Surfer软件绘制图件。

为了验证算法的准确性，将采用典型的地电模型K型和四层模型与实测数据进行试算。

2 模型与算例

2.1 算例一

选用K型地电模型正演理论数据进行反演计算，其真实的模型参数为见表1。

表1 K型地电模型理论参数Tab.1 Theoretical parameters of K-type geo-electric model

给出一个测深点上各极距反演视电阻率与理论视电阻率数据见表2。

表2 视电阻率数据及反演所得模型的理论视电阻率Tab.2 Comparison of apparent resistivity data observed with inverted

图2(a)中GCV(λ)函数曲线随着迭代次数的增加而呈现出缓慢下降的趋势。在达到第六次迭代时程序自动终止反演迭代，最后一次迭代的GCV(λ)函数曲线首支趋于平坦，程序从函数曲线的后段向下突变处自适应选出最小GCV(λ)函数值对应横轴上的正则化因子λ(0.001附近范围)；从图2(a)、图2(b)与表2，经过6次反演迭代，反演所得模型的理论视电阻率与实际视电阻率的相对误差为0.35%；反演所得模型参数为：第一层电阻率为50 Ω·m，厚度为20.5 m；第二层电阻率为98.4 Ω·m，厚度为10.8 m；第三层电阻率为50 Ω·m，与真实模型吻合较好。各次迭代视电阻率拟合相对误差为157.68%、126.16%、66%、11.34%、7.41%、0.36%，可见收敛速度较快，反演时间15 s。

图2 K型模型电测深一维自动迭代反演Fig.2 One-dimensional automatic iterative inversion result of K-model electric sounding(a) 广义交叉验证曲线;(b) 实测和模拟数据的拟合差曲线;(c)电测深一维反演结果图

2.2 算例二

选用四层地电模型正演理论数据进行反演计算，以验证模拟的准确性与有效性。其真实的模型参数为见表3。

表3 四层地电模型理论参数Tab.3 Four-layer geo-electric model theoretical parameters

四层模型正演理论数据经过正则化自动迭代反演得图3(a)广义交叉验证曲线、图3(b)实测与模拟数据拟合差曲线和图3(c)电测深一维反演结果图。

图3 四层模型电测深一维自动迭代反演Fig.3 One-dimensional automatic iterative inversion result of four-layer model electrical sounding(a)广义交叉验证曲线;(b)实测和模拟数据的拟合差曲线;(c)电测深一维反演结果图

给出一个测深点上各极距反演视电阻率与理论视电阻率数据见表4。

表4 视电阻率数据及反演所得模型的理论视电阻率Tab.4 Comparison of apparent resistivity data observed with inverted

从图3(a)可以看出GCV(λ)函数曲线整体在第六次迭代时较前五次迭代数值更小，GCV(λ)曲线前半部分平滑，后半部分开始缓慢下降，而后尾部缓慢增大上升。可见第六次迭代GCV(λ)函数值最小时对应的正则化因子λ取值在0.01附近，算法据此自适应选定最优正则化因子。从图3(a)，图3(b)中可以看出，经过6次反演迭代，反演所得模型的理论视电阻率与实际视电阻率的相对误差为0.42%；反演所得模型参数为：第一层电阻率为50 Ω·m，厚度为9.8 m；第二层电阻率为8.9 Ω·m，厚度为22.8 m；第三层电阻率为96.8 Ω·m，厚度为8.3 m；第四层电阻率为494 Ω·m与真实模型吻合较好。各次迭代视电阻率拟合相对误差为168.24%、146.44%、48.15%、7.64%、2.59%、0.42%，可见收敛速度较快，反演时间30 s。

2.3 实例

靶区位于野马泉附近，离矿区约9 km。工作区位于东准噶尔盆地东缘，气候干燥炎热，降水稀少。区内地表水主要接受大气降水补给，其次为构造裂隙水。物探工作主要查明断层产状和断层含水性，为后续矿区供水方案提供依据。区域内出露地层主要为志留系、泥盆系。志留系地层分布在工作区以南的广大地区，由志留系库布苏群下亚群(Skpa)组成。

选取工作区中某测线上的两点2530号点,2540号点,2520号点进行反演算法的试算。从图4(a)、图4(b)、图4(c)三个实测点处理后的广义交叉验证曲线可以看出经迭代后可较明显的识别出函数值最小时对应的正则化因子；而在图4(d)、图4(e)、图4(f)拟合差曲线图中可见三个测点的拟合误差均在允许范围内；在图4(g)、图4(h)、图4(i)三测点的电测深一维反演结果图中也可看出大致分层状况，进而为激电测深二维反演提供依据。在前期对资料进行初步推断解释的基础上，在Y3线2522号点布设钻孔ZK4-1，0 m～36 m以下219 mm套管对地表水进行封闭，36 m以下的井径为180 mm。0 m～25 m为第四系松散沉积物，25 m～70 m凝灰质砂岩，70 m以下为花岗岩。在深度130 m处打到构造破碎带，含水量较大，由于水压过大，导致冲击钻探方法难于施工，并于135 m深度处终孔，未打穿断层破碎带，终孔的孔径为110 mm。采用32 m3/h水泵进行抽水试验，经计算孔出水量为46.3 m3/h(1 211.2 m3/d)，不能满足选场1 500 m3/d的水量要求。故在2526号点增加钻孔ZK4-2,距ZK4-1钻孔4 m，孔径为315 mm，在深度165 m～180 m范围内见断层破碎带。采用32 m3/h水泵进行抽水试验，水位降深为156.5 m，经计算ZK4-2孔出水量为16.05 m3/h(385.2 m3/d)，在2526处钻井进行抽水试验时，2522钻井的水位变化很小，说明两钻孔之间可能存在泥质充填物堵塞。综合考虑，将工区作为矿山今后用水的主要取水区，应加以保护。

图4 2530号点、2540号点、2520号点实测数据电测深一维自动迭代反演Fig.4 2530 point,2540 point，2520 point electrical sounding one-dimensional automatic iterative inversion(a)2520号测点广义交叉验证曲线;(b)2530号测点广义交叉验证曲线;(c)2540号测点广义交叉验证曲线;(d)2520号测点拟合差曲线;(e)2530号测点拟合差曲线;(f)2540号测点广义交叉验证曲线;(g)2520号测点电测深一维反演结果图;(h)2530号测点电测深一维反演结果图;(i)2540号测点电测深一维反演结果图

3 结论

1)提出了一种电测深一维自适应正则化反演方法，在整个反演过程中根据广义交叉验证方法自适应选取正则化因子，压制了反演的奇异性并使反演过程更加收敛稳定。

2)依据电测深方法浅部分辨率高、深部分辨率较低的特点，实现依据电测深曲线极距数与对应的视电阻率自动创建多层等效模型，使所构建等效模型的电性参数与电测深曲线尽可能与实际地质情况贴近并在深度方向上具有相似的变化规律。

3)由于反演过程中无需手动输入初始模型参数，使得整个反演过程实现了全自动化。通过对模型模拟数据与野外实测数据进行反演试算，结果表明本文提出的正则化方法反演方法可行有效且可以对实际生产提供指导。