2018, Vol. 61
2. 中国科学院矿产资源研究院重点实验室, 中国科学院地质与地球物理研究所, 北京 100029;
3. 中国科学院地球科学研究院, 北京 100049;
4. 中国科学院大学, 北京 100049;
5. 长安大学地测学院, 西安 710054;
6. 沈阳师范大学数学与系统科学学院, 沈阳 110034
2. The Key Laboratory of Mineral Resources, Institute of Geology and Geophysics. Chinese Academy of Sciences, Beijing 100029, China;
3. Chinese Academy of Sciences University, Beijing 100049, China;
4. Chinese Academy of Sciences University, Beijing 100049, China;
5. School of Geology Engineering and Geomatics, Chang'an University, Xi'an 710054, China;
6. Shenyang Normal University, Shenyang 110034, China
地球物理探测的目的是通过观测地球物理场数据来推断地下地质结构.作为深地探测的重要手段，为了更好地认识地球内部结构与构造、寻找能源、资源和环境监测、灾害预报, 有必要开展地球物理勘探新理论、新方法、新技术以及地球物理反问题研究(杨文采, 2002).
地球物理中的瞬变电磁法(Transient Electromagnetic Method, 简称TEM)属时间域电磁感应方法，它是利用不接地或接地线源向地下发送一次场，在一次场的间歇期间，测量由地质体产生的感应电磁场(二次场)随时间的变化.根据二次场的衰减特征，就可以判定地下不同深度地质体的电性特征及规模大小等(牛之琏, 2007).近年来瞬变电磁法在国内外得到了迅速发展，瞬变电磁偏移成像技术成为目前研究的热点问题(李貅等, 2005).
如何寻找到一个数学上的处理办法，将瞬变电磁场变换成“波场”，即提取出电磁响应中与传播有关的特征，压制或去除电磁波传播过程中与频散、衰减有关的特征，借鉴地震勘探技术提高瞬变电磁法的解释精度，成为重要研究方向.Weidelt(1972)，Kunetz(1972)，以及Levy等(1988)，Lee等(1989)的相关研究成果，都揭示了在层状大地介质中，电磁扩散方程与地震波动方程间存在数学上的等效对应形式，但他们研究问题的着眼点都是将对应地电模型的波场模拟结果变换成时域电磁响应(李貅等, 2005; Xue et al., 2007).其实，研究如何将已知时域电磁场数据转换为虚拟波场数据，并进行拟地震资料解释，这将有利于偏移以及更加复杂的成像技术的应用(Xue, 2011, 2012; Guo et al., 2012; Xue et al., 2013).近年来，为了提高这种计算稳定性, 瞬变电磁波场转换算法不断地发展和改进(戚志鹏等, 2013; 钟华森等, 2016).
奇异值分解法(SVD)是一种矩阵分解(Matrix Decomposition)的方法, 可以作为类似于波场变换病态问题的处理手段，许多学者作出了大量的研究工作.由于奇异值变化范围比较大，直接进行奇异值分解必然引起计算结果的不稳定.由于大奇异值控制计算矩阵的主要信息, 小的奇异值控制计算矩阵的次要信息，在随后研究中, 作为算法的改进, 主要考虑大奇异值, 而忽略小奇异值的影响, 一定程度上提高了计算的稳定性.Hansen(1986)提出了在病态问题中，对于阶梯形的奇异矩阵，截断奇异值正则化的有效性.陈本池等(1999)采用奇异值分解算法实现了瞬变电磁场的波场变换, 并针对影响波场变换的诸因素(阻尼系数α的选取, 时窗宽度, 误差)进行了系统研究和计算.谢正超等(2015)利用TSVD方法解决温度场重建的病态方程.
考虑到小奇异值对计算精度的影响, 对于均匀下降型的奇异值矩阵，王振杰和欧吉坤(2004)借鉴杨文采(1989)思路, 将全部奇异值分成大奇异值和小奇异两部分, 并对这两部分分别进行修正.由于波场转换的强烈不适定性，提出一种新的奇异值计算方案, 奇异值矩阵保持不变，采用岭估计法只对小奇异所对应的波场进行精细计算，从而提高算法的精度.
1 TEM虚拟波场提取方法时域扩散场Zm(r, t)与虚拟波场U(r, τ)之间的关系式为(Lee et al., 1989)
|
(1) |
式中Zm(t)代表瞬变电磁场，U(τ)代表虚拟波场，τ为虚拟时间变量.
(1) 式的离散形式可以表示为如下算子方程：
|
(2) |
式中，Z, A, U分别表示离散后的电磁场数据矩阵、系数矩阵和虚拟波场矩阵.
通过公式(2)，可根据扩散场数据Z，得到虚拟波场数据U.但是这种转换属于典型的非线性不适定问题.下面介绍与之相关的奇异值分解法.
截断奇异值分解法：
方程(2)中的系数矩阵A为Hilbert空间中的一个算子，定义A的奇异值σi(i=1, 2, …, n)为Hermite算子A′A的平方根，且满足σ1≥σ2≥…≥σn，A′为矩阵A的共轭转置矩阵.
构造一个正交归一序列pi，vi，i=1, 2, …, n.满足(谢正超等，2015)：
|
(3) |
其中pi, vi可以看成是系数矩阵A的左奇异向量和右奇异向量, v′ i为vi的转置.
可以直接采用奇异值分解法，根据(3)式和(2)式来计算(1)式中的虚拟波场矩阵U，但是，由于奇异值变化范围比较大且不稳定，直接进行奇异值分解必然引起计算结果的不稳定.为提高计算的稳定性，在随后研究中, 作为算法的改进, 主要考虑大奇异值, 而忽略小奇异值的影响, 采用截断奇异值分解法进行计算(谢正超等，2015)，相应地，系数矩阵可以分成两部分：
|
(4) |
其中，
|
(5) |
式中，k表示截断位置(矩阵A对角线元素的截断阶数).
|
(7) |
|
(6) |
因此，方程(2)可以写为
|
(8) |
于是方程(2)的解U可以表示为
|
(9) |
进一步地，如果分别考虑由大奇异值和小奇异值所对应的波场值，(9)式可写为
|
(10) |
上式右端第一项表示大的奇异值所对应的波场值项，第二项表示小的奇异值所对应的波场值项.
如果近似地将(6)式等于0，即A2=0.截断后的(10)式可以表示成(徐文，2016；Onunwor and Reichel, 2017)：
|
(11) |
虽然采用截断奇异值分解法对(2)式进行反演得到的值趋于稳定，但是由于把那些小的奇异值舍去后，一定程度上损失了反演值的精确性，因此有必要对截断奇异值法进行改进一下，以提高反演精度.
改进截断奇异值法：
为了让整个奇异具矩阵更加紧凑，降低病态程度，杨文采(1989)对大的奇异值趋近于用截断奇异值法进行修正，对小的奇异值趋近于用岭估计进行修正.虽然这样做可以使得奇异值更加紧凑，但是其修改了全部的奇异值，降低了反演数据的准确性，作者根据其思想，从另外一个角度提出了改进方案，即对大的奇异值保持不变，对小的奇异值进行修正.
即保持A1的大小，仅对A2进行处理.
A2对应的奇异值
|
(12) |
令A2=lI，因此，修正后的奇异值可以表示成
|
(13) |
其中，I表示单位矩阵，l为一个常数，也可看作是正则化参数.关于l的选择，本文根据对比反演图像效果，手动选择.因此，公式(4)中的系数矩阵A可近似表示为
|
(14) |
把(14)代入(9)式，可得到
|
(15) |
通过公式(11)—(15)，可以明显看出解U多了一项，弥补了丢失了奇异值信息，增加了解的精确性.
公式(13)以对直接由奇异值分解法得到得反演值进行修正，不但让小的奇异值参与计算，且降低了系数矩阵的病态程度.
表 1列出了不同截断位置的系数矩阵条件数，可以看出：对条件数为无穷的系数矩阵进行修正后，系数矩阵的条件数得到了降低，证明了算法得合理性.
|
表 1 不同截断位置的系数矩阵条件数 Table 1 Number of coefficient matrix bars for different truncation positions |
根据野外实际观测瞬变电磁场数据或者模型仿真数据，分别采用奇异值分解法、截断奇异值分解法和本文提出的改进截断奇异值分解法进行波场反变换，提取出三种方法所对应的虚拟波场值，然后，再利用改进截断奇异值方法计算出相应的虚拟波场值.
具体步骤如下：
步骤1：读取野外实际观测瞬变电磁场数据，形成电磁场数据矩阵Z
步骤2：计算系数矩阵A，并对其进行奇异值分解
步骤3：利用截断奇异值法计算出虚拟波场值Uj，j=1, 2, …, k.选择出k.
|
(16) |
步骤4：根据公式(15)，计算虚拟波长值.
|
(17) |
为验证改进的截断奇异值法虚拟波场的提取效果，采用三种方法(直接奇异值法，截断奇异值法和本文的修正式截断奇异值法)对电性源瞬变电磁仿真数据进行虚拟波场提取，并对比其效果.
图 1为设计的地电模型，为两层水平模型，地层分界面深度为300 m，上层电阻率为100 Ωm, 下层电阻率为200 Ωm，发射源长度500 m，发射电流2 A，接收点位于垂直轴向距发射源中心点1000 m.
|
图 1 电性源瞬变电磁仿真模型 Fig. 1 Grounded-wire source transient electromagnetic method model |
图 2a为垂直磁场衰减曲线，图 2b为采用直接的奇异值分解法进行虚拟波场提取的结果，从图 2b可知，由于系数矩阵的严重病态性，得到的反演结果并不理想，在坐标轴内波动非常剧烈，无法准确判断出异常点.图 2c为截断奇异值分解法计算结果，曲线过于光滑，反演值波峰不够明显，存在丢失信息问题；图 2d为修正截断奇异值分解法计算结果，通过对比可知：修正截断奇异值分解法得到的反演值波峰较为明显.
|
图 2 电性源瞬变电磁仿真结果 (a)衰减曲线; (b)直接奇异值分解法; (c)截断奇异值分解法; (d)修正截断奇异值分解法. Fig. 2 Modeling result (a) Decay curve; (b) Curve from traditional singular value decomposition method; (c) Curve from Truncated Singular Value Method; (d) Curve from Modified Truncated Singular Value Method. |
测区位于山西宁武煤田北部，属于山西高原朔平台地之低山丘陵，全区多为黄土覆盖.测区含煤地层主要为石炭系山西组(C3s)和石炭系山太原组(C3t)地层，以灰白、灰黄色砂岩和页岩为主，夹煤系地层，可采煤层为4#、9#和11#煤层，最大埋深约200 m, 采空区普遍发育.区域内主要含水层有碳酸盐岩岩溶裂隙含水层，碎屑岩裂隙含水层及松散岩层孔隙含水层.碳酸盐岩岩溶裂隙含水层由寒武-奥陶系灰岩组成，碎屑岩裂隙含水层在区内分布广泛，含水层由石炭、二迭系砂岩构造.区内由下至上发育有奥陶系上统上马家沟组、石炭系中统本溪组、上统太原组、二叠系下统山西组、下石盒子组、上统上石盒子组以及新生界第三系静乐组和第四系中、上更新统、全新统地层.
当地下煤体局部被采出后，在岩体内形成一个有一定规模的空间，使周围的应力平稳状态遭受破坏，产生局部的应力集中，采空区顶板在上覆岩层压力的作用下，发生变形、断裂、位移、冒落，形成的冒落带、断裂带、变形弯曲带，在地下水的充填及地表水沿裂缝向采空区渗漏，其电阻率将明显发生变化，形成一个低阻电性体，也与围岩电性形成较明显的差异，这样就为瞬变电磁勘查采空区破坏区提供了前提条件.
在已知采空陷落区(山西省平鲁区某煤矿)进行测试，采空区深度150 m.所布设的测线长度400 m, 测点距离10 m(图 1)300—500号测点位于采空区范围内, 采用仪器为GDP-32Ⅱ，试验装置为定源回线，回线大小为200 m×300 m，工作频率为32 Hz，发射电流8 A.
|
图 3 已知测区测线布置图 Fig. 3 Layout of survey line |
采用修正截断奇异值分解法对观测曲线进行处理计算, 计算结果如图 4所示, 其中, 图 4a表示600号测点(非采空区)虚拟波场转换结果, 图 4b表示400号测点(采空区)虚拟波场转换结果.通过对比图 4(a, b)可知, 未采区测点数据的虚拟波场曲线比采空区虚拟波场曲线复杂, 这与富水采空区形变后的地质结构有关.
|
图 4 修正截断奇异值分解法转换结果 (a) 600号测点(非采空区)转换结果; (b) 400号测点(采空区)转换结果. Fig. 4 Extracted result from modified truncated singular value method (a) Extracted result of 600 survey point; (b) Extracted result of 400 survey point. |
图 5为虚拟波场转换剖面和已知测区视电阻率等值剖面对比图, 图 5(a—c)中红色部分表示正值，绿色部分表示负值.图 5a表示由传统奇异值法计算的虚拟波场剖面, 由于传统奇异值法计算的不完备, 不仅出现很多附加场信息, 同时对采空区的反映不够明显, 图 5b表示采用截断奇异值法计算的虚拟波场剖面, 图像质量已有改善, 但是, 对于采空区与未采区的分界仍不够明显, 图 5c表示采用本文所述修正奇异值法计算结果, 在500号测点附近采空区分界面比较明显, 图 5d表示已知测区实测视电阻率等值线断面图.
|
图 5 虚拟波场转换结果剖面和已知测区视电阻率等值剖面图 (a)传统奇异值法计算结果; (b)截断奇异值法计算结果; (c)修正奇异值法计算结果; (d)实测视电阻率断面图. Fig. 5 Peudo wave-field converted sections and apparent resistivity sectiom (a) Section from traditional singular value decomposition method; (b) Section from Truncated Singular Value Method; (c) Section from Modified Truncated Singular Value Method; (d) Field data section of apparent resistivity. |
瞬变电磁成像在数学上是一个强烈的非线性不适定问题, 观测数据的微小的扰动, 都会使部分参数发生很大的变化.寻求优化的数学解法十分重要.本文主要对瞬变电磁响应虚拟波场提取方法进行了研究.针对传统的奇异值分解法不能确保奇异值矩阵均匀下降问题，本文在截断奇异值法研究基础上，提出采用改进截断奇异值法进行提取，根据野外实际观测瞬变电磁场数据或者模型仿真数据，分别采用奇异值分解法、截断奇异值分解法和本文提出的改进截断奇异值分解法进行波场反变换，提取出三种方法所对应的虚拟波场值.模型计算结果表明：改进截断奇异值法比传统的奇异值分解法得到的转换结果更好，并在某煤矿采空区进行了实验工作，成功分辨出采空区顶界面和分界面.当然, 如果参与计算的截断奇异值矩阵较小时, 所提取出来的虚拟波形比较光滑, 但会丢失有用信息，如果参与计算的截断奇异值矩阵较大时, 所提取出来的虚拟波形会出现强烈起伏的现象, 带来多余的附加信息，所以，选择合适的截断奇异值矩阵阶数需要根据实际情况来确定.
Chen B C, Li J M, Zhou F T. 1999. Wave-field conversion method for transient electromagnetic field. Oil Geophysics Prospecting (in Chinese), 34(5): 539-545. |
Guo W B, Xue G Q, Li X, et al. 2012. Correlation analysis and imaging technique of TEM data. Exploration Geophysics, 43(3): 137-148. DOI:10.1071/EG11034 |
Hansen P C. 1987. The Truncated SVD as a method for regularization. BIT Numerical Mathematics, 27(4): 534-553. DOI:10.1007/BF01937276 |
Kunetz G. 1972. Processing and interpretation of magnetotelluric soundings. Geophysics, 37(6): 1005-1021. DOI:10.1190/1.1440310 |
Lee K H, Liu G, Morrison H F. 1989. A new approach to modeling the electromagnetic response of conductive media. Geophysics, 54(6): 1180-1192. |
Levy S, Oldenburg D, Wang J. 1988. Subsurface imaging using magnetotelluric data. Geophysics, 53(1): 104-117. DOI:10.1190/1.1442393 |
Li X, Xue G Q, Song J P, et al. 2005. An optimize method for transient electromagnetic field-wave field conversion. Chinese Journal of Geophysics (in Chinese), 48(5): 1185-1190. |
Niu Z L. 2007. Time-Domain Electromagnetic Method (in Chinese). Changsha: Central South University Press.
|
Onunwor E, Reichel L. 2017. On the computation of a truncated SVD of a large linear discrete ill-posed problem. Numerical Algorithms, 75(2): 359-380. |
Qi Y F, Yin C C, Wang R, et al. 2015. Multi-transient EM full-time forward modeling and inversion of m-sequences. Chinese Journal of Geophysics (in Chinese), 58(7): 2566-2577. DOI:10.6038/cjg20150731 |
Qi Z P, Li X, Wu Q, et al. 2013. A new algorithm for full-time-domain wave-field transformation based on transient electromagnetic method. Chinese Journal of Geophysics (in Chinese), 56(10): 3581-3595. DOI:10.6038/cjg20131033 |
Weidelt P. 1972. The inverse problem of geomagnetic induction. Geophysical Journal International, 1972, 35(1-3): 379. |
Wang Z J, Ou J K. 2004. A novel singular value correction scheme for ill-posed problems and its application in geodetic survey. Progress in Natural Science, 14(6): 672-676. |
Xu W, Chen Y, You W. 2016. Application of SVD Method in Ill-posed Problem. Bulletin of Surveying and Mapping (in Chinese), (1): 62-63, 83. |
Xue G Q, Yan Y J, Li X. 2007. Pseudo-seismic wavelet transformation of transient electromagnetic response in engineering geology exploration. Geophysical Research Letters, 34(16): L16405. DOI:10.1029/2007GL031116 |
Xue G Q, Yan Y J, Li X. 2011. Control of wave-form dispersion effect and applications in TEM imaging technique for identifying underground objects. Journal of Geophysics and Engineering, 8(2): 195-201. DOI:10.1088/1742-2132/8/2/007 |
Xue G Q, Bai C Y, LI X. 2012. Extracting the virtual reflected wavelet from tem data based on regularizing method. Pure and Applied Geophysics, 169(7): 1269-1282. DOI:10.1007/s00024-011-0392-1 |
Xue G Q, Gelius L J, Li X, et al. 2013. 3D pseudo-seismic imaging of transient electromagnetic data-a feasibility study. Geophysical Prospecting, 61(S1): 561-571. |
Xie Z C, Wang F, Yan J H. 2015. Comparative studies of Tikhonov regularization and truncated singular value decomposition in the three dimensional flame temperature field reconstruction. Acta Physica Sinaca, 64(24): 18-26. |
Yang W C. 2002. Non-linear geophysical inversion methods:review and perspective. Progress in Geophysics (in Chinese), 17(2): 255-261. |
Yang W C. 1989. Geophysical inversion and seismic tomography. Bejing: Geological Publishing House.
|
Zhong H S, Xue G Q, Li X, et al. 2016. Pseudo wavefield extraction in the multi-channel transient electromagnetic (MTEM) method. Chinese Journal of Geophysics (in Chinese), 59(12): 4424-4431. DOI:10.6038/cjg20161204 |
陈本池, 李金铭, 周凤桐. 1999. 瞬变电磁场的波场变换算法. 石油地球物理勘探, 34(5): 539-545. |
韩波, 李莉. 2011. 非线性不适定问题的求解方法及其应用. 北京: 科学出版社.
|
李貅, 薛国强, 宋建平, 等. 2005. 从瞬变电磁场到波场的优化算法. 地球物理学报, 48(5): 1185-1190. DOI:10.3321/j.issn:0001-5733.2005.05.029 |
牛之琏. 2007. 时间域电磁法原理. 长沙: 中南大学出版社.
|
戚志鹏, 李貅, 吴琼, 等. 2013. 从瞬变电磁扩散场到拟地震波场的全时域反变换算法. 地球物理学报, 56(10): 3581-3595. DOI:10.6038/cjg20131033 |
王振杰, 欧吉坤. 2004. 一种新的病态问题奇异值修正方案及其在大地测量中的应用. 自然科学进展, 14(6): 672-676. DOI:10.3321/j.issn:1002-008X.2004.06.012 |
谢正超, 王飞, 严建华, 等. 2015. 炉膛三维温度场重建中Tikhonov正则化和截断奇异值分解算法比较?. 物理学报, 64(24): 18-26. |
徐文, 陈义, 游为. 2016. 奇异值分解法在病态问题中的应用. 测绘通报, (1): 62-63, 83. |
杨文采. 2002. 非线性地球物理反演方法:回顾与展望. 地球物理学进展, 17(2): 255-261. DOI:10.3969/j.issn.1004-2903.2002.02.010 |
杨文采. 1989. 地球物理反演和地震层析成像. 北京: 地质出版社.
|
钟华森, 薛国强, 李貅, 等. 2016. 多道瞬变电磁法(MTEM)虚拟波场提取技术. 地球物理学报, 59(12): 4424-4431. DOI:10.6038/cjg20161204 |