基于自适应模糊聚类分析的重力张量欧拉反褶积解
曹书锦，朱自强，鲁光银
(中南大学 信息物理工程学院，湖南 长沙，410083)
摘要：利用将重力全张量数据应用在欧拉反褶积中，规避位场梯度计算的精度问题，引入自适应模糊聚类算法克服聚类数目需要求预先确定、模糊聚类分析局部最优、分类不确定等弱点，并准确的确定多异常源的情况。核密度估计结果表明，张量欧拉反褶积比预设结构参数的欧拉反褶积方法更能表征地下异常类型；反演结果表明，传统欧拉反褶积难于识别在深大型异常源附近的浅部规模相对较小异常源；过滤后的欧拉反褶积解的空间包络基本与初始模型的一致，张量欧拉反褶积在获得多异常源的空间结构信息更具有优势。
关键词：全张量重力梯度；三维欧拉反褶积；自适应模糊聚类分析
中图分类号：P631.1 文献标志码：A 文章编号：1672-7207(2012)03-1033-07
Gravity tensor Euler Deconvolution solutions based on adaptive fuzzy cluster analysis
CAO Shu-jin, ZHU Zi-qiang, LU Guang-yin
(School of Info-physics and Geomatics Engineering, Central South University, Changsha 410083, China)
Abstract: Tensor gravity Euler Deconvolution was used to avoid calculate Potential field’s derivatives. To accurately interpret the case with many causative sources. The adaptive fuzzy clustering analysis to overcome the drawbacks of fuzzy clustering analysis was used, such as (pre-determined number of clusters, local optimization, or classification uncertainty). Kernel density estimation results show that Tensor gravity Euler Deconvolution is a better estimate structural index of causative sources than traditional Euler Deconvolution. The inversion results show that Euler Deconvolution is difficult to identify small and shallow anomaly source，a relatively large-scale and deep anomaly source near it; the Euler solutions after filtered are coincident with the Forward model, which proves that the method of adaptive fuzzy clustering analyzing tensor gravity Euler Deconvolution is better than that analyzed by fuzzy clustering analyzing.
Key words: full gradient tensor gravity; three-dimensional Euler Deconvolution; adaptive fuzzy clustering analysis
欧拉反褶积是一种半自动/自动估算场源位置的反演方法，是针对磁测数据反演场源深度和水平位置而提出的，它以反演时使用参数少，无需知道场源物性(密度或磁性)的先验信息，不需要准确的解释模型，只要事先确定出与场源性质有关的构造指数或通过枚举构造指数。其方法简明、快速，只需利用位场异常及其导数就可以进行反演，快速而有效地圈出异常源的基本轮廓，尤其适合于大面积位场数据的分析和解释。应用欧拉反褶积解决位场问题的前提是需高精度的取得欧拉齐次方程中位场导数计算值，即导数的换算精度问题是影响欧拉反褶积精度的核心问题。为取得高精度的位场导数，Mickus等[1]利用垂向重力梯度经由傅里叶变换计算得到全张量分量，但是其误差最大为3.3×10-9 s-2，远大于重力梯度张量测量精度。
Hinojosa等[2]在傅里叶域导出张量分量间希尔伯特变换上的关系。Mickus等[1]和Wang等[3]利用三次B样条变换在空间域高精度的完成了二维和三维情况下的重力场值和梯度值之间的变换，相对FFT变换仅能满足均匀间隔网格，而三次B样条变换却能适用于三维均匀和非均匀规则网格和二维规则和非规则网格的情况。张凤旭等[4]推导出重力位余弦谱一般表达式以及重力异常各阶导数计算公式，建立了位场余弦谱分析理论。传统重力测量主要依赖于密度的变化，而张量重力比其更敏感于密度差异引起的微小密度梯度变化,重力张量测量精度远比计算梯度的精度高。自2003年其测量精度已远优于同类重力测量精度[5-6]。欧拉张量反褶积远由于传统欧拉反褶积，源于前者利用测量梯度数据，后者利用计算梯度。当异常场中源为多个时(如结构界面或大小深度不一致)，则异常特征可能代表的不仅仅是一个特定的构造指数，即构造指数可能是变化的。同时在实际情况中，地下异常源类型未知，其构造指数很难确定。在欧拉褶积中错误地给定构造指数，使得欧拉解往往具有发散的趋势。为消除欧拉方程解的发散问题，传统的提取正确解的策略基于FitzGerald等[7]，但其不能估计整个欧拉解集合的优良性或者确定哪些相类似的解标示同一个异常源[8]。范美宁[9]突破了以往只用位场异常数据计算的限制，在数据使用方面延展到导数、解析信号、不同高度异常的组合，丰富欧拉反演方法中的数据源，在一定程度上改善了反演效果。姚长利等[10]针对欧拉反演中的问题提出了水平梯度滤波、基于主体异常距离准则的有效性统计筛选和构造指数有效性筛选等改进措施，促进了三维欧拉反演方法的实用化[10-11]。Mikhailov等[12]采用人工智能的方式提取欧拉解。Gerovska[13]在算法Stavrev等[14-15]的基础上采用非聚类和聚类的方式处理欧拉解，并取得了非常好的效果。Ugalde等[8]采用模糊聚类方法来提取欧拉解，但是虽然其避免人为确定构造指数带来的欧拉解得发散的问题，但其过度依赖于人工确定地下异常源的个数。在此，本文作者根据张量数据高精度的特点，利用重力张量数据代替计算位场导数完成重力张量欧拉反褶积，规避了传统欧拉反褶积中对位场导数计算精度高要求的限制。直接求取构造指数，而不采用预设的方式，以减小可能错误确定构造指数使得欧拉解过度发散的问题。引入自适应模糊聚类算法克服了聚类数目要求预先确定、局部最优、分类不确定等弱点，以达到准确地确定了多异常源的目的。最后通过数值例子进一步表明本文方法的有效性。
1 传统欧拉反褶积
三维欧拉方程可表示为[16-17]：
式中：f为位场函数，
；(x, y, z)为测量点位置；(x0, y0, z0)为待求解位场异常源位置；N为位场异常随场源深度衰减变化“陡缓”的量度，特定的地质构造有特定的衰减率，即构造指数，可作为未知参数或作为一个预设参数，采用一定步长在0~3范围内枚举的构造指数。式(1)可写为：
式中：B为消除区域背景场的影响，而引入一个代表区域背景异常的参数；a为失真率，当且仅当N=0是存在[9]。3个坐标方向的梯度值fx,y,z可以利用空间与或波数域的一般位场变换计算出来。但如果直接将观测梯度值用于式(2)则更可取[6, 18]。
2 张量欧拉反褶积
位场张量是位场异常沿坐标系轴向的导数(表明重力张量的优点)。重磁场的一阶导数、二阶导数及延拓后的异常等也同样满足欧拉方程，并且这些重、磁位场数据之间的联系也满足泊松关系[9]。三维张量欧拉反褶积方程[6]由传统的欧拉方程和2个附加的类似方程和组成：
式中：
；
；
为异常背景场；
为重力梯度；
为重力梯度张量。其中，异常背景场很小，将自式~中去除，同时将构造指数N作为待求参数；利用n×n的滑窗在研究区域滑动，并将滑窗中的n2个测点代入相应的方程组，利用奇异值分解法或其他数值解法求取异常空间位置 (x, y, z)及构造指数N。
3 自适应模糊聚类分析
欧拉解的聚类分析就是通过无监督训练将欧拉解集按相似性分类，把相似性大的归为一类，占据特征空间的一个局部区域，而每个局部区域的聚合中心又起着相应类型的代表作用，以指示相应的地质体的空间信息。Ruspini[19]首次将模糊集理论应用到聚类分析中，提出模糊聚类算法。在各类聚类分析中，需要初始化某些聚类参数，诸如希望的聚类数、最小标准差、初始聚类中心、聚类中心间的最小距离等，当聚类数目与实际地质情况不符合时(尤其是先验信息不足，难于确定算法初始化信息时)，由于这些算法是全局搜索优化算法，受样本数据的分布影响较大，很容易陷入局部极值，导致分类错误。因此初聚类中心和样本输入次序对最终结果有很大的影响。常用的试错法不仅工作量巨大，而且不能保证聚类结果的最优性。Gerovska在算法Stavrev等[13-15]的基础上采用非聚类和聚类的方式处理欧拉解，并取得了非常好的效果。Ugalde等[8, 20]采用模糊聚类方法来提取欧拉解，但是虽然其避免人为确定构造指数带来的欧拉解得发散的问题，但其过度依赖于人工确定地下异常源的个数。
特引入一类自适应模糊聚类方法[21-22]，定义如下：数据集
，
，其组成数据集合X，利用聚类分析算将其划分如M个质心为vi，半径为ri的簇，dij为第i与j个质心距离，将各聚类之间的关系表现矩阵U={uij}，uij为第i与j个质心间的关联度，
，m为模糊指数，模糊指数越大，聚类的模糊程度就越大[23]。聚类算法主要是基于目标函数最小化原理，引入如下目标函数(见文献[22])：
自适应模糊聚类方法算法如下[24]。
初始化聚类算法信息，在l=1, 2, …, n内迭代。
(1) 求取每个聚类的质心：
, 
(2) 计算质心半径：
，
(3) 计算簇质心距离
[25]：
,
, 
(4) 更新
[26-27]，对每个特征向量xj，
，i={1, 2, …, c}有
，
(5) 计算簇间相似关系，并确认那两个簇最相 似[28-29]
，
1≤i, j≤M(l-1)
4 应用算例
4.1 算例1
首先考察均匀半空间，检验基于H-自适应重力张量正演(H-FEM)的欧拉反褶积方法的有效性。利用解析解的和H-FEM分别计算简单规则模型(见文献[30])，并给正演数据添加3%的噪声。立方体大小2 000 m×2 000 m×2 000 m，质心(0, 0, 1 500)。球体半径 1 000 m，球心为(0, 0, 1 500)，水平圆柱体半径1 000 m，沿Y轴长2 000 m，轴心位置为(0, 0, 1 000)。剩余密度为2.36 g/cm3, 测高为0 m，测区为X：-10 000~ 10 000 m, Y：-10 000~10 000 m；测网大小为100 m×100 m (蓝色的点为欧拉反褶积解，等高线为Gxy，为显示方便已将等高线移z=2 500的平面上)。
欧拉解基本上落在长方体中心部分，顶部的已标示出长方体的平面上的轮廓(图1(a))；欧拉解完全处于球体的中心部位(图1(b))，且图2曲线2结构指数N的概率密度峰值点的完全处以一点(N=2)，与理论情况完全相符，表明本文正演算法的正确性。
图1 欧拉反褶积反演各简单异常体
Fig.1 Gravity tensor Euler Deconvolution based on Cube(a), Sphere(b) and Horizontal cylinder(c)
图2 欧拉反褶积结构指数的概率密度函数
Fig.2 Probability density function of structure index in Euler Deconvolution
4.2 算例2
模型2，(X0, Y0, Z0)为立方体质心，模型剩余密度均为2.36 g/cm3, 测高为0 m，测区为X：-10 000~ 10 000 m，Y：-10 000~10 000 m；测网大小为100 m× 100 m。5个立方体模型参数如表1所示。
表1 模型参数
Table 1 Model parameters
图3 均匀半空间赋存组合模型张量重力Gxy分量等值线图(含3%高斯随机噪声)
Fig.3 Full tensor gravity gradient contour map buried sphere (Synthetic data contaminated with 3% Gaussian random noise for inversion)
同时由于重力张量分量对模型深度和模型大小较为敏感，模型1和2相对于测网尺寸较小其显示的异常范围较小(图3)，甚至浅部规模较小的异常源模型中没有一个欧拉反褶积解(图4(a))。分析在正演模型内的欧拉反褶积解的个数，表明本文算法能更正确地估计异常体的深度(表2)；规模较大异常体的模型5在总体欧拉解中占较大比率，同时也表明单一的结构指数是难于表征含有多个异常源的情况。若逐步的将噪音权重提高，则发现靠近较大异常源的规模相对较小异常源的将淹没(图4(b), (c)，表2)。
图4 有噪声和无噪声数据算法比较
Fig.4 Comparison optimal algorithm and tensor Euler Deconvolutions
表2 处于正演模型内部欧拉反褶积解个数
Table 2 Number of Euler Deconvolution solutions inside Forward Model
4.3 聚类分析
令m=2，
，
，相似度阈值
，终止标准
，最大迭代次数为500，初始聚类数Nfc=50，Nafc为1~100之间的一个随机数且大于异常源个数，异常源个数可由张量不变量的闭合边界个数确定，如不能确定可选取较大初始值。
图5所示为模糊聚类及自适应模糊聚类聚类的中心点。可以看出：当各个聚类的之间距离的最小可控制距离较小时，在异常源影响较大的区域内，聚类的个数较多。
图5 欧拉反褶积解的聚类分析
Fig.5 Clustering analysis Euler Deconvolution solution
为进一步剔除坏解，针对取得聚类样本再一次进行核密度估计，同时只将核密度估计的最大峰值附近区域作为结构指数N的过滤标准同时利用张量不变量提取到的异常源边界作为辅助提取得到的欧拉反褶积解已能充分反映低下异常源的情况(见图6)。
图6 过滤后的欧拉反褶积解
Fig.6 Solution of Euler Deconvolution after filtered
5 结论
(1) 对比简单规则模型欧拉反褶积解和理论解，及其结构指数密度核估计表明本文方法正确可行。
(2) 单一的结构指数很难表征地下异常源的几何特征，预测单一结构指数常导致欧拉反褶积解过度 发散。
(3) 采用自适应聚类方法克服了聚类数目需要求预先确定，模糊聚类分析局部最优、分类不确定等弱点，并准确地确定了多异常源的情况。
(4) 利用本文算取得欧拉反褶积有效解占优率高，无需对其做更多分析，便可用于地质解释。
(5) 欧拉反褶积解受到测网尺寸和异常体规模深度的影响，规模较大的异常源在整个欧拉反褶积解集合中总体占优，而反演模型相对于测网较小的异常源的可能性较小；若逐步的将噪音权重提高，则发现靠近较大异常源的规模相对较小异常源的将被淹没。
(6) 本文利用边缘检测张量不变量数据提取异常体边界来辅助剔除欧拉反褶积解集中的发散解，是非常高效的位场数据过滤的一种手段。
参考文献：
[1] Mickus K L, Hinojosa J H. The complete gravity gradient tensor derived from the vertical component of gravity: A fourier transform technique[J]. Journal of Applied Geophysics, 2001, 46(3): 159-174.
[2] Hinojosa J H, Mickus K L. Hilbert transform of gravity gradient profiles: Special cases of the general gravity-gradient tensor in the fourier transform domain[J]. Geophysics, 2002, 67: 766-769.
[3] Wang B, Krebes E S, Ravat D. High-precision potential-field and gradient-component transformations and derivative computations using cubic B-splines[J]. Geophysics, 2008, 73(5): I35-I42.
[4] 张凤旭, 孟令顺, 张凤琴, 等. 重力位谱分析及重力异常导数换算新方法——余弦变换[J]. 地球物理学报, 2006, 49(1): 244-248.
ZHANG Feng-xu, MENG Ling-shun, ZHANG Feng-qin, et al. A new method for spectral analysis of the potential field and conversion of derivative of gravity-anomalies: Cosine transform[J]. Chinese J Geophys, 2006, 49(1): 244-248.
[5] Vasilevsky A, Druzhinin A, Evans R, et al. Feasibility of FTG reservoir monitoring[J]. SEG Technical Program Expanded Abstracts, 2003, 22(1): 1450-1453.
[6] Zhang C, Mushayandebvu M F, Reid A B, et al. Euler Deconvolution of gravity tensor gradient data[J]. Geophysics, 2000, 65(2): 512-520.
[7] FitzGerald D, Reid A, McInerney P. New discrimination techniques for Euler Deconvolution[J]. Computers and Geosciences, 2004, 30(5): 461-469.
[8] Ugalde H, Morris W A. Cluster analysis of Euler Deconvolution solutions: New filtering techniques and geologic strike determination[J]. Geophysics, 2010, 75(3): L61-L70.
[9] 范美宁. 欧拉反褶积方法的研究及应用[D]. 长春: 吉林大学地球探测科学与技术学院, 2006: 19-63.
FAN Mei-ning. The study and application of Euler Deconvolution method[D]. Changchun: Jilin University. College of GeoExploration Science and Technology, 2006: 19-63.
[10] 姚长利, 管志宁, 吴其斌, 等. 欧拉反演方法分析及实用技术改进[J]. 物探与化探, 2004, 28(2): 150-155.
YAO Chang-li, GUAN Zhi-ning, WU Qi-bin, et al. An analysis of Euler Deconvolution and its improvement[J]. Geophysical and Geochemical Exploration, 2004, 28(2): 150-155.
[11] 王家林. 对我国石油重磁勘探发展的几点思考[J]. 勘探地球物理进展, 2006, 29(2): 82-86.
WANG Jia-lin. Views on the domestic situation and progress of gravity and magnetic petroleum exploration[J]. Progress in Exploration Geophysics, 2006, 29(2): 82-86.
[12] Mikhailov V, Galdeano A, Diament M, et al. Application of artificial intelligence for Euler solutions clustering[J]. Geophysics, 2003, 64(2): 168-180.
[13] Gerovska D, Araúzo-Bravo M J. Automatic interpretation of magnetic data based on Euler Deconvolution with unprescribed structural index[J]. Computers & Geosciences, 2003, 29(8): 949-960.
[14] Stavrev P. Euler Deconvolution using differential similarity transformations of gravity or magnetic anomalies[J]. Geophysical Prospecting, 1997, 45(2): 207-246.
[15] Stavrev P, Gerovska D, Arauzo-Bravo M J. Depth and shape estimates from simultaneous inversion of magnetic fields and their gradient components using differential similarity transforms[J]. Geophysical Prospecting, 2009, 57(4): 707-717.
[16] Reid A B, Allsop J M, Granser H, et al. Magnetic interpretation in three dimensions using Euler Deconvolution[J]. Geophysics, 1990, 55(1): 80-91.
[17] Thompson D T. Euldph: A new technique for making computer-assisted depth estimates from magnetic data[J]. Geophysics, 1982, 47(1): 31-37.
[18] 许家姝. 漠河盆地重磁场特征与沉积构造研究[D]. 长春: 吉林大学地球探测科学与技术学院, 2007: 38-41.
XU Jia-shu. Study on the features of gravi-magnetic field and the sedimentary structure in Mohe Basin[D]. Changchun: Jinlin University. College of GeoExploration Science and Technology, 2007: 38-41.
[19] Ruspini Z H. A new approach to clustering[J]. Information and Control, 1969, 15(1): 22-32.
[20] Balasko B, Abonyi J, Feil B. Fuzzy clustering and data analysis toolbox[M]. University of Veszprem, 2002: 8-12.
[21] Pham D L, Prince J L. An Adaptive fuzzy C-means algorithm for image segmentation in the presence of intensity inhomogeneities[J]. Pattern Recognition Letters, 1999, 20(1): 57-68.
[22] Xie X L, Beni G. A validity measure for fuzzy clustering[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1991, 13(8): 841-847.
[23] Pal N R. On cluster validity for the fuzzy C-means model[J]. IEEE Transactions on Fuzzy Systems, 1995, 3(3): 370-379.
[24] Kaymak U, Setnes M. Fuzzy clustering with volume prototypes and adaptive cluster merging[J]. IEEE Transactions on Fuzzy Systems, 2002, 10(6): 705-712.
[25] 张敏, 于剑. 基于划分的模糊聚类算法[J]. 软件学报, 2004, 15(6): 858-869.
ZHANG Min, YU Jian. Fuzzy partitional clustering algorithems[J]. Journal of Softwave, 2004, 15(6): 858-869.
[26] Martino F, Sessa S. Implementation of the extended fuzzy C-means algorithm in geographic information systems[J]. Journal of Uncertain Systems, 2009, 3(4): 298-306.
[27] Martino F, Loia V, Sessa S. Extended fuzzy C-means clustering algorithm for hotspot events in spatial analysis[J]. International Journal of Hybrid Intelligent Systems, 2008, 5(1): 31-44.
[28] Dubois D, Prade H. A unifying view of comparison indices in a fuzzy set-theoretic framework[C]//Recent Development in Fuzzy Set and Possibility Theory. New York: Pergamon Press, 1982: 3-13.
[29] Kosko B, Burgess J C. Neural networks and fuzzy systems[J]. The Journal of the Acoustical Society of America, 1998, 103(6): 3131.
[30] 曹书锦, 朱自强, 鲁光银, 等. 基于单元细分 H-自适应有限元全张量重力梯度正演[J]. 地球物理学进展, 2010, 25(3): 1015-1023.
CAO Shu-jin, ZHU Zi-qiang, LU Guang-yin, et al. Forward modelling of full gravity gradient tensors based H-Adaptive mesh refinement[J]. Progress in Geophysics, 2010, 25(3): 1015-1023.
[31] 鲁宝亮, 范美宁, 张原庆. 欧拉反褶积中构造指数的计算与优化选取[J]. 地球物理学进展, 2009, 24(3): 1027-1031.
LU Bao-liang, FAN Mei-ning, ZHANG Yuan-qing. The calculation and optimization of structure index in Euler Deconvolution[J]. Progress in Geophys, 2009, 24(3): 1027-1031.
[32] 马国庆, 杜晓娟, 李丽, 等. 利用斜角总水平导数和欧拉反褶积法研究四川盆地线性构造[J]. 地球物理学进展, 2011, 26(3): 916-921.
MA Guo-qing, DU Xiao-juan, LI Li, et al. Utilizing total horizantal derivate of tittangle and Euler Deconvolution to study the linear structure of the Sichuan Basin[J]. Progress in Geophys, 2001, 26(3): 916-921.
(编辑 陈爱华)
收稿日期：2011-05-25；修回日期：2011-07-22
基金项目：国家高科技研究发展计划(“863”计划)项目(2007AA06Z102)；国家自然科学基金资助项目(41174061)
通信作者：鲁光银(1976-)，男，湖北宜昌人，教授，从事地球物理正反演研究；电话：13975894898；E-mail: luguangyin2006@126.com
