基于SVD-TLS算法的非预测欧拉反褶积_王海青.pdf
第41卷第5期 2013年IO月 煤田地质与勘探Vol. 41 No.5 Oct. 2013 COALGEOLOOY 2.中南大学信息物理工程学院,湖南长沙410083 摘要针对地下有多个异常源时,单一预测构造指数难于表征多个异常源.采用非预测欧拉反精 积以避免可能错误确定构造指数使得欧拉解过度发散的问题;同时针对欧拉反精积超定方程组的 条件数很大,致使欧拉反精积解集中良解占优率低等解的非唯一性和解的不稳定性等局限性,采 用奇异值分解总体最小二乘法(SVD-TLS算法),以降低由于奇异值分析不当造成计算欧拉解非唯 一性和解的不稳定性的问题,并利用SVD-TLS的截断误差构造闽值函数对解集进行过滤.数值结 果表明了算法的有效性和可靠性. 关键词奇异值分解总体最小二乘法;非预测;欧拉反柑积;多源 中图分类号P631文献标识码ADOI 10.3969/j.issn.1001-1986.2013.05.017 Euler deconvolution based on SVD-TLS WANG Haiqing1, QING Jianke2 I. Geological Exploration Technologies Institute of Anhui Province, Bengbu 233000, China; 2. School of Info-physics and Geomatics Engineering, Central South University, Changsha 410083, China Abstract It is difficult to characterize complicated anomaly sources with the single prescribed structure index and to determine structural index when these types of causative sources are unknown, and solutions of Euler deconvo- Iution with divergent位endcaused by error set or wrongly estimated structural index. Pathological equations of Euler deconvolution results in good solution. In this paper, we use SVD-TLS algorithm to reduce instability and non-uniqueness of Euler solution due to improper singul缸valueanalysis, then use truncation error of SVD-TLS to construct threshold function to filter the erroneous solutions. The proposed is verified by one numerical example and results show that the approach is effective. Key words SVD-TLS; unprescribed; Euler deconvolution; multiple-source 对于重力异常导数的快速反演/成像算法,位场 的成像技术有张量不变量(1-3]、图像自动边界识别与增 强(4-5)、位场偏移(6-7]、井中成像[町、概率密度成像(9-11]、 解析信号[也13]等。但大多数方法仅提供位场场源的指 引信息,并没有给出完整的位场物性分布。其在联合 反演上可以起到局部约束作用。 欧拉反槽积作为一种快速高效、适应性强的位场 解释方法,其欧拉超定方程组广义逆条件数极大,解 的扰动也大(14),致使欧拉反榴积方法存在解的非唯 一性和解的不稳定性等局限性(15]。通常构造指数依 据人工经验和地质先验信息选取(16-17],一方面可能带 人很大的人为干扰因素,另一方面也不容易一次性选 准构造指数,还可能会导致反演结果的不准确(15]。 实际上,异常源多为复合地质体,形态复杂,这就为 经验选取构造指数造成了困难一一地质体的复杂性 导致的构造指数的不确定性(14-18]。 收稿日期2012-08-25 相对于非预测构造指数的欧拉反槽积,采用预测 构造指数的计算模式则其解趋于发散,但解集总体优 于前者。如何在众多的深度点中分辨和确定有效深度 点一直是困扰人们的一个难题。传统的提取正确解的 策略基于FitzGeratdl川,但其不能估计整个欧拉解集 合的优良性或者确定哪些相类似的解标示同一个异 常源(20]。姚长利等针对欧拉反演中的问题提出了水 平梯度滤波(21]、基于主体异常距离准则的有效性统 计筛选和构造指数有效性筛选等改进措施(16,22)。 Gerovska2003)采用了欧拉超定方程组的标准差作为 评价标准,以过滤欧拉反榴积解集中的谬解,但其解 集中的大部分解源自条件数极大的欧拉超定方程,致 使解集中良解占有率低(4 针对欧拉反榴积超定方程组的逆条件数很大,欧 拉反槽积解集中良解占优率低等解的非唯一性和解 的不稳定性等局限性,从计算欧拉反榴积方法源头的 作者简介王海青(1966一),男,河南南阳人,高级工程师,从事地球物理研究. ChaoXing 第5期王海青等基于SVD-TLS算法的非预测欧拉反禧积 81 奇异值分析出发,对非预测欧拉反槽积以避免可能错 误确定构造指数使得欧拉解过度发散的问题,采用奇 异值分解总体最小二乘法(Singularvalue decomposi- tion of the total least squares ,简写SVD-TLS, 以降低由于奇异值分析不当造成计算欧拉解非唯一 性和解的不稳定性。最后通过数值例子进一步表明了 方法的有效性。 1 欧拉反槽积反演 重力位场满足Euler齐次方程,则三维欧拉方程 可表示为[23-24] if aJ aJ x-x0)一+(Y-Yo)一+(z-z0)一=-NJ1 ax句F 式中f fx,y,z)为位场函数;(x,y,z)为测量点位 置;(Xo,Yo,zo)为待求解位场异常源位置;N为位场 异常随场源深度衰减变化“陡缓”的量度,特定的地质 构造有特定的衰减率,即构造指数,可作为未知参数 或作为一个预设参数,采用一定步长在0~3范围内枚 举的构造指数。由于枚举结构指数将使得解集呈线性 增大,不利于过滤谬解,因此将构造指数N作为待 求参数,故式(1)重写为 aJ司f与fu.r司f司ffJf 1m x一+y一+z一+Xo一+Yo一+句一+NJJ2 &钞缸’ax句yt泣 式中B为消除区域背景场的影响,而引人一个代表 区域背景异常的参数,当基于剩余密度计算时,B数 量级很小,故此B略去。3个坐标方向的梯度值 司flθx,y,z可以利用空间域或波数域的一般位场变换 计算出来。利用nxn的滑窗在研究区域滑动,将区域 中的n2个测点代入相应的方程组,利用如奇异值分 解法或其它数值解法求取异常空间(x,y,z)及构造指 数No 2 奇异值分解总体最小二乘法(SVD-TLS 由式(2)组成的方程组由如下形式 Axb 3 I aJ aJ aJ l矿aJ 式中A为I------f I; b为Xoγ+均γ+ |缸dv z[l,x]T o 对增广矩阵B进行奇异值分解 BU“f.VT 其中U为西矩阵;Z是半正定对角矩阵,“I. I \ V,, v., diag(σ11,,σii,,σ,0,0,,OJ; VI .. .“ I 飞PPI l V21 V 町IeRnxk V22εRnxk; k为矩阵A的秩。 在最小二乘法中,构造的自相关矩阵的有效秩p 可由如下公式确定 σtt=σtt Iσ11, 1三k运h6 选择一个接近于零的正数做阔值,并把等式左边所取 的小于此阔值的最大正整数k取为矩阵的有效秩p. 其中h为所得到的奇异值的个数。 才=[vk,j, vkl,j,,vkp,jr 其中才是西矩阵Y的第j列,则有公式 Sptn艺Pσ2vi川l x;s-iI,11s-l,1, il,,p 在GEROVSKA2003)的基础上,对式(5)引人一 个问值吨以判定欧拉反榴积解的稳定性[14,26-27] σr=εI z0 其中ε为式(5)的截断误差[28];Zo为欧拉解预测 深度。 3 模型试验 由于简单模型有解析解,且其构造指数容易确 定,故实验效果良好。利用解析解计算简单规则模型, 并给正演数据添加3的噪声,列出了几个简单地质 体模型实验结果。其中模型的测高为om,测区为X 一2500025 000 m, Y一2500025 000 m;测网大小 为250mx250 m。 3.1 单个正方体 为验证算法的有效性,设定模型的立方体大小为 1 000 mxl 000 mxl 000 m,质心(0,o, 2 500),剩余 密度为2.27g/ cm3。 图1基于SVD-TLS截断误差估计的欧拉反榴积 解集完全处于正演模型正中央,这与理论解析解完全 一致,证明算法是有效、可行的[29]。 3.2 两个立方体的组合模型 立方体大小分别为1000 mx 1 000 mx 1 000 m和 800 mx800 mx800 m,质心分别(0,0, 2 500)和 -5000, -5 000, 1 200),剩余密度分别为2.27g/cm3 和1.50g/cm3。 ChaoXing 0, 2 500)、(-5创泊,一5创胁,1200)、(0,-4 500, 650; 剩余密度分别为2.27g/cm.3、1.50g/cm.3、2.00g/cm3。 从图5能够看出,在未对欧拉反槽积解进行过滤的 前提下,表层的欧拉反槽积解受噪声干扰比较明显, 尤其是在测网的四周,但其影响深度有限,对数据解 释无影响;浅部的异常源受到附近深部异常源的影响, 其欧拉深度有下延的趋势。相对于2个异常源,3个 异常源的情况使得这种下延更为明显,且倾角更大, 这也说明3个异常源情况比2个异常源更为复杂。为 此,利用该密度估计和SVD-TLS阑值概率密度曲线 的相关分析,确定采用阔值1.1-2.5作为过滤标准。利 用该算法计算得到的欧拉反榴积解虽然有一定的倾角 存在干扰,但图6的结果表明,经SVD-TLS阑值过 第41卷 3 000 I 000 煤田地质与勘探 500 2000 82 是2500 慰 图1基于SVD-TLS的欧拉反榴积解云图 Cloud point of Euler deconvolution based on SVD-TLS 们喇川咀 .。”1430/1.弘3 14302860/1.53 3.0 去 运2.0 剿 1.0 2.5 Fig. 1 从图2能够看出,在未对欧拉反榴积解进行过滤 的前提下,表层的欧拉反槽积解受噪声干扰比较明 显,但其影响深度有限,对数据解释无影响;浅部的 异常源受到附近深部异常源的影响,其欧拉深度有下 延的趋势,且有一定的倾角。这对布置地质钻孔有一 定的误导。为此利用核密度估计和SVD-TLS阑值概 率密度曲线的相关分析(图3C29l,确定采用阑值 1.5-2.3作为过滤标准。图4表明欧拉反槽积解云图 恰处于正演模型中。 3.3 3个立方体的组合模型 立方体大小分别为1侧mxl刷刷刷m,则mx 800 mx800 m, 500 mx500 mx500 m;质心分别(0, 4 图4欧拉反裙积解云图(过滤后) Cloud point of Euler deconvolutionafter Filtered 02 ♂、J . .,, .,、‘ 、‘.” 嘟 。 去1.0 运2.0 剧3.0 斗 Fig. 4 0-1497/01.5 .。~1497/1.5~3 14972995/01.5 φ1497-2995/1.5-3 一5 。」‘0 .... 0 500 lα)() 毒1500 慰12000 2500 3000 -3 2 图2欧拉反裙积解云图(未过滤) Cloud point of Euler deconvolution before Filtered 0-1497 /0-1.5分别为欧拉深度/构造指数 Fig. 2 图5欧拉反精积解云图(未过滤) Cloud point of Euler deconvolution before Filtered 。~l饵”'/0~1.5 ;他如1阴1.5-3 19993998/01.5 19993998/1.5-3 A斟U 。 g己之隧剿 Fig. 5 1.4 1.2 倒0.8 衡 革0.6 0.4 I 1.5 2 结构指数(N 图3结构指数及截断SYD孔S阑值函数的概率密度估计 Fig. 3 Probability density estimation of structure index and truncated SVD-TLS threshold function 2 3 3 2.5 0.5 0.2 图6欧拉反精积解云图(过滤后) Cloud point of Euler deconvolutionbefore Filtered Fig. 6 ChaoXing 第5期王海青等基于SVD-TLS算法的非预测欧拉反槽积 83 滤,其对异常源的指示是有效的。 5结论 不难看出,该算法在浅部和测区边界的欧拉反槽 积解受数据噪声干扰影响,但并不影响欧拉反槽积解 的解释。在多个异常源存在时,浅部异常源受深部异 常源的影响致使欧拉反榴积解有下延和稍有向外倾 覆的趋势;在异常源较浅时,此种情况更为明显。由 于采用常规的欧拉反槽积求取方法难于确定构造指 数,且单一的构造指数难于表征地下复杂的地质构 造,而SVD-TLS阑值概率密度曲线和构造指数的概 率密度曲线相关分析的相关域,取峰值邻域作为过滤 欧拉反榴积的阑值,其结果显示对单一及多个异常源 共存的情况,可确定多异常源的空间位置信息。 参考文献 [I] SAAD A H. Understanding gravity g回ients-atutorial[.月.The Leading Edge, 2006, 258 942-949. [2] MATARAGIOJ, KIELEYJ. Application of full tensor gradient inv四antsin detection of intrusion-hosted sulphide minerali7Jltion 泣句“cationsfor deposition mechanisms[]]. Mining Geoscience, 2009, 271 95-98. [3] PEDERSEN L B, RASMUSSEN T M. The gradient能nsorof po旬出alfield anomalies some in『plicati创IS侃侃侃collec挝.onand data processing ofmaps[J]. Geophys阳,1”。,55121558-15“ - [4] ZHANG Lili, HAO Tianyao, WU Jiansheng, et al. Application of image enhancement techniques to阴阳出lfield缸”[巧.Appli“ Geophysics, 2005, 23 145- 152. [5] CUMANI A. Edge detection in multiψectra1 images[]]. Graphical Models and Image阶ocessing,1991, 531 40-51. [6] LAMMERMANN T, BADER BL, MONK.LEYS J, et al. Rapid leukocyte migration by integrin-ind句览ndentflowing and 呵ueezing[J].Natur宫,2008,4537191 51- 55. [7] LnJ Zhixin, LU H血,LuH.Res回rchon correlation imaging of gra叫tyg回ientbased on ffi[巧.Procedia Earth and Planetary Sci en,印,2011,其11203-209. [8] ROBBINS S L. R田xaminationof也ues used as constants in calculating rock density命。mborehole gra叫tyda饱问.Geophy- sics, 1981, 462 208-210. [9] NlCOLASD, FEDELER, MAMAN’KO, et al. Tomographic- probability description of solitons in bose-eins能incondensates[]]. The European Physical Journal B-Condensed Mat阳andComplex Systems, 2003, 363 385-390. [IO] MAURIELLO P, PATELLA D. Resistivity tensor probab凶ty tomography[]]. Progr部sIn Electromagnetics Res咄咄,2008, 84 129- 146. [ll] GUO Lianh时,SHILei, MENG Xiaohong. 30 correlation imaging 。fmagnetic total field anomaly and its ver幅“lgradient阴.Journal of Geophysics and Engineering, 20118 287-293. [12] ANSARI AH, ALA岛IDARK.R创uctionto也epole of magnetic anomalies using analytic signal[]]. World Applied Sciences Journal, 2009, 74 405- 409. [ 13] HSU S K, SIBUET J C, SHYU C T. High-resolution detection of geologic boun也riesfrom po忧ntial-fieldanomalies; an enhanced analy世csignal technique[]]. Geophysics, 1996, 612 373-386. [14] GEROVSKAD,ARAUZO-BRAVO M J. Automatic Interpretation 。fmagnetic data based on euler deconvolution with unprescribed structural index[]]. Compu馆”&Geosciences, 2003 , 298 949-960. [15]鲁宝亮,范美宁,张原庆.欧拉反裙积中构造指数的讨算与优 化选取[巧.地球物理学进展,20093 1027- 1031. [16]姚长利,管志宁,吴其斌,等.欧拉反演方法分析及实用技术 改进[巧.物探与化探,2004,282 150- 155. [17]张季生.位场自动反演技术的研究现状及意义[巧.地球学报, 2创始,276 609-612. [18] PAASCHE H, EBERLE D. Automated compilation of pseudo阳 H也ologymaps叠。mgeophysical也usets a comparison of Gustafson-Kessel and fuzzy c-means cluster algori也ms[月- ExplorationGeopysics, 2012, 424 275-285 [19] FITZGERALD D, REID A, MCINERNEY P. New discrimin- ation techniques for Euler deconvolution[]]. Computers姐d Geosciences, 2004, 30句461-469. [20] UGALDE H , MORRIS W A. Cius忧rAn且lysisof Euler d悦。,nvolutionsolutions new filtering techniques and geologic 归keD础nnin血on[J].Geophysics, 2010, 753 L6J- L70. [21]范美宁.欧拉反裙积方法的研究及应用[DJ.长春吉林大学, 2创6. [22]王家林.对我国石油重磁勘探发展的几点思考[几勘探地球物 理迸展,2006,292 82- 86. [23] REID A B, ALLSOP J M, GRANSER H , et创.Magnetic 姐“rpretationin three dimensions using Euler d民onvolution[J]. Geophysics, 19佣,55180-91. [24] TH。如fPSOND T. Euldph a new technique for m冰山gcomputer- assisted d句也estimatesfrom magnetic data[J]. Geophysics, 1982, 471 31- 37. [25]范美宁,孙运生,田庆君.关于欧拉反裙积方法计算中的一点 改进[巧.物探化探计算技术,2005,272 171-174. [26] STAVREV P Y. Euler deconvolution using differential similarity 位ansationsof gra叽tyor magnetic anomalies[]]. Geophysical Prospecting, 1997, 452 207-246. [27] STAVREV P, GEROVSKA D, ARAUZO-BRAVO M J. Depth and shape estimates企omsimultaneous inversion of magnetic fields and th町gradientcomponents using di岱rentialsimilarity transs[]]. Geophysical Prospecting, 2009, 574 707-717. [28] HANSEN PC. Regul缸izationtools a matlab package for analysis and solution of discrete ID-posed problems[]]. Numerical algori也ms,1994, 61 1- 35. [29]曹书锦,朱自强,鲁光银.基于自适应模糊聚类分析的重力张 量欧拉反槽军fl_.[月.中南大学学报(自然科学版),2012, 433 1033-1039. ChaoXing