CN104122598A - 一种从物探重力异常中提取断层异常的结构函数法 - Google Patents

Abstract

本发明公开了一种从物探重力异常中提取断层异常的方法，属于地球科学中的勘探地球物理学。现有的分离区域异常与局部异常的方法(趋势分析方法、频率域滤波方法、滑动平均法、上延等)存在边界效应明显和没有考虑区域异常的空间结构特性的不足，导致分离提取出来的区域异常精度不高，为了能有一个精度高的分离断层区域异常与局部异常的方法，本发明给出了结构函数法。在已知重力观测数据的空间分布几何特征、协方差函数和区域异常形式下，它用求最优、线性、无偏内插估计量的方法求出断层区域异常。趋势分析方法只是结构函数法的一个特例。本发明适合分离提取具有线性结构的断层重力异常。实验证明本发明是有效性的并且在实施方面是可行的。

Description

一种从物探重力异常中提取断层异常的结构函数法

一.技术领域：本发明涉及到一种从物探重力异常中提取断层异常的方法，按中华人民共和国国家标准“学科分类与代码”，本发明属于地球科学中的“勘探地球物理学(编码：1702065)”，即，应用地球物理(物探，地球物理勘探)。 

二.背景技术

断层在地质学研究中的重要性不言而喻，例如，活动的深大断层带一般是地震多发的地震带(如环太平洋地震带和四川521地震带)。勘探地球物理学(简称物探)中的重力异常一般都能反映这种深大断层带，表现为重力异常梯度带。这种断层重力异常具有明显的空间结构性，即线性的，具有很长的走向(例如几十公理甚至上百公理，但在垂直方向上也就是几公理甚至不到1公理)，在其走向方向上物探重力数据变化缓慢且是相关的，但在其垂直方向上数据变化剧烈且相关性很小。大断层重力异常大多数情况下表现为线性的区域异常，小断层重力异常在个别场合下表现为线性的局部异常。这样分离提取断层重力异常就归结为如何从物探重力异常图中分离提取区域异常与局部异常的问题。分离提取区域异常与局部异常，在勘探地球物理学(简称物探)重力资料的分析处理中占有十分重要的地位，因为这种分离是很多处理方法的基础。多年来人们一直在研究这个问题，直到最近还有人在拟文发表这方面的研究成果。然而，这个问题直到现在也没有得到较好的解决。分离区域异常与局部异常目前常用的方法是趋势分析方法，但是，用趋势分析方法求得的区域异常值在边界附近与真实区域异常值差别较大，即边界效应明显。此外，在一个工区用一个大的统一的趋势面代表复杂的区域异常是不太合理的，因为那样往往会引起假异常，或丢掉一些有用异常。其它一些划分区域异常与局部异常的方法，如频率域的滤波方法，空间域的滑动平均法、上延等，对划分区域异常与局部异常有一定的效果，但也有其缺点。这些方法共同的一个缺点是没有考虑物探重力异常的空间结构特性(例如，线性的，或等圆状的，矩形状的等)。以断层物探重力异常为例，空间结构特性是指断层物探重力异常有非常明显的走向，是线性的，在其走向方向上物探重力数据变化缓慢且是相关的，但在其垂直方向上数据变化剧烈且相关性很小。由于没有考虑物探重力异常的空间结构特性，这就导致应用它们(常规方法)求得的区域异常值与真实区域异常值有时差别较大。为了能有一个精度高且效果好的分离区域异常与局部异常的方法，本文给出一个方法----结构函数法。结构函数法考虑了物探重力异常的空间结构特性并对边界效应给出了解决方法，使得分离出的区域异常与局部异常的精度更高、更能反映真实情况。 

三.发明内容

0发明目的(阐明本发明所要解决的问题，从而归纳出本发明的目的)：如背景技术中所讲，现有的分离物探重力区域异常与局部异常的方法存在没有考虑重力区域异常的空间结构特性和边界效应明显的不足，为了能有一个精度高且效果好的分离重力区域异常与局部异常的方法，本文给出一个方法----结构函数法。结构函数法考虑了物探重力异常的空间结构特性并对边界效应给出了解决方法，使得分离出的区域异常与局部异常的精度更高、更能反映真实情况。下面是本发明----结构函数法的工作原理、技术方案和优点。 

1工作原理 

1.1基本思想 

概括的讲，结构函数法就是一种估值方法。更具体地说，结构函数法就是在考虑了求区域异常所用物探重力观测数据的形状、大小及其与待估区域异常值的点相互之间的空间分布位置等几何特征，以及重力观测数据的空间结构信息后，为了达到线性、无偏和最小估计方差的估计，而对每个求区域异常所用重力观测数据分别赋予一定的权系数，最后进行加权平均来估计出待估区域异常点处的重力区域异常值。从重力观测数据中减去所求区域异常就可得到局部异常。例如，已知平面上一个不大的有限邻域D(x₀)内S₀，S_PS₁，S₂，S₃，S₄六处的物探重力观测数据分别为z₀，z_P，z₁，z₂，z₃，z₄，要估计出S₀和S_P点处的重力区域异常与局部异常。对具有某种结构的物探重力观测数据来说，S₀点待估的重力区域异常值为局部异常为而对S_P点，则有一组不同的系数λ₁，λ₂，λ₃，λ₄与其对应，使得区域异常局部异常为诸系数λ₁，λ₂，λ₃，λ₄一般是利用求条件极值的拉格朗日乘数法求得的。 

从数学角度抽象地说，结构函数法就是一种在已知物探重力观测数据z(x)的协方差函数(具体求法下面将说明)和z(x)的区域异常形式下，求最优、线性、无偏内插估计量(Best linear unbiased Estimator，简写为BLUE)的方法。无偏性(泛性)是指对某个量反复观测后，该量观测值的平均值与各观测值的差的数学期望值为零。 

结构函数法要解决的问题是如何求出估计区域异常值公式m^*＝λ₁z₁+λ₂z₂+λ₃z₃+…+λ_nz_n中的诸系数λ₁，λ₂，λ₃，…λ₄。为此，需要先说明几个相关的概念。 

1.2相关的几个概念 

1.2.1物探重力观测数据的两重性 

物探重力观测数据z(x，y)观测前可看作是一个随机场，观测后可看作是一个普通的二元实值函数(或空间点函数)。有时把它理解为随机场，有时把它理解为普通的空间点函数，需视具体场合而定。 

物探重力观测数据z(x，y)即有结构性又有随机性，这是因为：一方面，当空间一点(x，y)固定后，z(x，y)就是一个随机变量，这就体现了其随机性；另一方面，在空间两个不同点(x_p，y_p)及(x_p+Δx，y_p+Δy)处的重力观测数据具有某种程度的相关性，这种相关性反映了重力观测数据的某种连续性和关联性(例如，在断层重力异常的走向方向上数据是相关的)，也就体现了其结构性的一面。 

由于重力观测数据具有这种随机性和结构性，经典概率统计方法(只考虑随机性)是不够用的，需要引入一个基本工具-----结构函数，来较好地研究重力观测数据z(x，y)。 

1.2.2结构函数 

先考虑简单的一维情况，假设空间点x只在一维的x轴上变化，把变量z(x)在x、x+h点处的值之差的方差之半定义为z(x)在x轴方向上的结构函数，记为γ(x，x+h)，简记为γ(x，h)，即： 

$\mathrm{\γ}(x,h)=\frac{1}{2}\mathrm{Var}[z\left(x\right)-z(x+h)]=\frac{1}{2}E{[z\left(x\right)-z(x+h)]}^2-\frac{1}{2}{\left\{E\right[z\left(x\right)]-E[z(x+h)\left]\right\}}^2....\left(1\right)$

其中Var[z(x)]为z(x)的方差，E[z(x)]为z(x)的数学期望。 

从式(1)可见，结构函数γ(x，h)一般是依赖于x和h两个自变量的。如果又知道γ(x，h)与x取值无关，只依赖于h(滞后、或间隔、或步长)，则可把结构函数γ(x，h)写为γ(h)，此时以h为横坐标，以γ(h)值为纵坐标做出的图形就叫结构图。 

把上述在一维x轴上的点x和间隔h，改为在二维平面上的空间点x_P和间隔h_p(Δx，Δy)，则可类似地得到二维平面上的结构函数的定义。为了简化叙述，在不致发生混淆处，把空间点x_P仍记为x，间隔h_p(Δx，Δy)仍记为h。 

1.2.3结构函数的计算 

一维结构函数γ^*(h)的计算公式(二维情况下的结构函数的求法类似)如下。 

可以把x轴上相隔为h的N(h)对点x_i和x_i+h(i＝1，2，…，N(h))处的N(h)对观测值z(x_i)和z(x_i+h)(i＝1，2，…，N(h))看成是z(x)和z(x+h)的N(h)对实现。于是可使用z(x)的算术平均值代替z(x)的数学期望的方法来计算γ^*(h)： 

${\mathrm{\γ}}^{*}\left(h\right)=\frac{1}{2N\left(h\right)}\underset{i=1}{\overset{N\left(h\right)}{\mathrm{\Σ}}}{[z\left(x_i\right)-z(x_i+h)]}^2-\frac{1}{2}{\left\{\frac{1}{N\left(h\right)}\right[z\left(x_i\right)-z(x_i+h)\left]\right\}}^2......\left(2\right)$

这样，对不同的滞后h，可算出相应的γ^*(h)值来。经计算后，在h-γ^*(h)的座标图上把诸点(h，γ^*(h))标出，再将相邻各点用线段连接起来，就得到实验结构函数图。对此图用曲线拟合的方法(具体拟合时，有许多技巧)可得到结构函数γ(h)与h的具体函数关系。 

1.2.4基本假设 

(1)协方差函数 

若z(x)是个随机场，在空间两点x和x+h处的两个随机变量z(x)和z(x+h)的二阶中心混合矩： 

Cov{z(x)，z(x+h)}＝C(x，x+h)＝E[z(x)×z(x+h)]-E[z(x)]×E[z(x+h)]称为随机场z(x)的自协方差函数，简称协方差函数。 

(2)基本假设 

假设z(x)具有一个非平稳的数学期望(即，区域异常)和一个平稳的协方差函数，即： 

$\left\{\begin{array}{c}E\left[z\left(x\right)\right]=m\left(x\right)\\ E\left[z\left(x\right)z(x+h)\right]=m\left(x\right)m(x+h)+C(x,x+h)\end{array}\right.---\left(3\right)$

其中m(x)就是物探重力观测数据z(x)的区域异常，而协方差函数C(x，x+h)存在且平稳(即只依赖于h而与x无关)，可记为C(h)。 

(3)协方差函数与结构函数的关系 

①协方差函数与结构函数之间有一个重要关系式： 

C(x，x+h)＝C(0)-γ(x，x+h)  (4) 

根据协方差函数C(x，x+h)与结构函数γ(x，x+h)的定义可以证明式(4)，证明如下， 

因为假设式(3)中E[z(x)z(x+h)]＝m(x)m(x+h)+C(x，x+h)，而且C(x，x+h)只依赖于h而与x无关，所以，当h＝0时，E[z(x)]²＝m²(x)+C(0)： 

在上式中，令x＝x+h，则E[z(x+h)]²＝m²(x+h)+C(0)。 

根据结构函数γ(x，x+h)的定义有： 

2γ(x，x+h)＝D[z(x)-z(x+h)]＝E[z(x)-z(x+h)]²-E[z(x)-z(x+h)]E[z(x)-z(x+h)]＝E[z(x)]²+E[z(x+h)]²-2E[z(x)×z(x+h)]-[m(x)-m(x+h)]²＝C(0)+m²(x)+C(0)+m²(x+h)-2[C(x，x+h)+m(x)m(x+h)]-[m(x)-m(x+h)]²＝2C(0)-2C(x，x+h) 

所以，C(x，x+h)＝C(0)-γ(x，x+h) 

如果又知道γ(x，x+h)与x取何值无关，只依赖于h(滞后、或间隔、或步长)，则可把结构函数γ(x，x+h)记为γ(h)，这样则有： 

C(h)＝C(0)-γ(h)  (5) 

②另外，协方差函数还有一个重要的性质：当|h|→∞时，C(h)→0，或写成C(∞)＝0。实际上，C(h)是反映变量z(x)和z(x+h)之间相关程度的，当距离h变得太大时，这种相关就消失了，即h→∞时，C(h)→0，记作C(∞)＝0。 

式(4)、式(5)的重要性在于通过γ(h)就可求得C(h)：C(0)＝γ(∞)，而结构函数γ(h)是可以通过对实际物探重力观测数据作计算并拟合求得(上面已介绍)。 

1.3结构函数法方程组 

1.3.1区域异常和局部异常 

设物探重力观测数据z(x)的数学期望E[z(x)]＝m(x)，这里把m(x)称为区域异常，也即一般意义上的区域异常。 

对于非平稳(或说有区域异常)的物探重力观测数据z(x)，假设它们可分解为区域异常和局部异常两部分，即： 

z(x)＝m(x)+R(x) 

其中，m(x)＝E[z(x)]为在点x处z(x)的数学期望，而R(x)＝z(x)-m(x)称为局部异常，它的数学期望值为零，因为： 

E[R(x)]＝E[z(x)-m(x)]＝E[z(x)]-m(x)＝0 

区域异常m(x)通常采用多项式形式： 

一维时m(x)＝a₀+a₁x+a₂x²+…， 

二维时m(x，y)＝a₀+a₁x+a₂y+a₃x²+a₄xy+a₅y²÷… 

区域异常还有次数的区别，有一次(线性)区域异常，二次区域异常及更高次区域异常(视多项式的次数而定)。由于在实际工作中对z₀点的区域异常的估计，只是根据其周围一个不大的有限邻域D(x₀)内的数据来进行的，为了表示D(x₀)内的区域异常，一般只需要用一次或二次多项式就足够了，也就是说只需要用一次(线性)区域异常和二次区域异常就足够了。至于三次区域异常，由于其曲面形状较复杂，在实际工作中，在一个有限邻域内很难遇到，况且我们研究区域异常的目的总是要了解物探重力观测数据z(x)在空间的变化总趋势，用一、二次区域异常，由于其形状较规则，有利于达到此目的，而三次或三次以上区域异常则不行。 

即使在一次和二次区域异常中，还是要尽可能多的采用线性区域异常，因为线性区域异常要简单得多。那么，何时可采用线性区域异常，何时须采用二次区域异常呢？这要根据以下两方面的情况来确定： 

①要看样品数据的密集程度。当我们估计某一点的区域异常时，需要使用该点周围一定范围内的数据，这就要确定一个适当大小的“窗口”，如果数据很密集，原“窗口”内数据很多，则可适当缩小“窗口”。由于窗口小，所以窗口内的区域异常就可近似地看作为线性区域异常：反之，如果原“窗口”内数据不足，就应将“窗口”放大，以便有足够的数据用于估计，这时，这种较大范围的区域异常就可能要用二次的了。 

②要看数据递增或递减的速率。当物探重力观测数据沿某方向上升或下降得很快时，则要应用二次区域异常：反之，则用线性区域异常。由于在小范围内，物探重力观测数据变化缓慢，适合于用一次多项式拟合。 

1.3.2用结构函数法求区域异常和局部异常时的假设条件及问题的提出 

用结构函数法求区域异常和局部异常时的假设条件可概括为如下四点： 

①   $\left\{\begin{array}{c}E\left[z\left(x\right)\right]=m\left(x\right)\\ E\left[z\left(x\right)z(x+h)\right]=m\left(x\right)m(x+h)+C(x,x+h)\end{array}\right.$

②假设z(x)可分解为区域异常m(x)和局部异常R(x)之和：z(x)＝m(x)+R(x)。 

③假设区域异常m(x)在人们工作尺度下变化缓慢，且可以表示为k+1个x的单项式函数f_l(x)(l＝0，1，2，…，k)(假设f_l(x)是已知的)的线性组合： 

$m\left(x\right)=\underset{l=0}{\overset{k}{\mathrm{\Σ}}}{a}_l{f}_l\left(x\right)$

其中x是邻域v中任一点，a_l为未知的系数。 

④假设事先已知道协方差函数C(x，x+h)。 

C(x，x+h)＝C(0)-γ(x，x+h)，通过γ(x，x+h)求C(x，x+h)，而γ(x，x+h)可通过对z(x)作实验结构图来求得(在上面结构函数的计算部分中已讲过)。 

问题的提出：在上述假设下，我们要求出： 

区域异常m(x)的最优线性估计量，或进一步说，要求出区域异常的各个未知系数a₀，a₁，a₂，…a_k的最优线性估计量。 

1.3.3结构函数法方程组 

在上述假设下，又设已知在n个信息样品点x_α(α＝1，2，…，n)处的物探重力观测值为z_α(α＝1，2，…，n)，即z(x_α)＝z_α(α＝1，2，…，n)，估计在任一固定点x处的区域异常值m(x)。 

我们设m(x)的估计量m^*(x)是n个已知数据z_α的线性组合，即 

$m^{*}\left(x\right)=\underset{\mathrm{\α}=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ρ}}_{\mathrm{\α}}{z}_{\mathrm{\α}}$

其中ρ_α为权系数。 

这样用结构函数法求区域异常的问题就归结为求权系数ρ_α(α＝1，2，…，n)了。 

本文要这样来选择这些权系数ρ_α，使得对任意a_l(l＝0，1，2，…，k)都能满足无偏性(也称为泛性)和最小方差性(也称最优性)条件。 

①无偏性(泛性)条件：即E[m^*(x)]＝m(x) 

由于一方面， $m\left(x\right)=\underset{l=0}{\overset{k}{\mathrm{\Σ}}}{a}_l{f}_l\left(x\right)$

另一方面又有 

$\begin{array}{c}E\left[m^{*}\left(x\right)\right]=E\left[\underset{\mathrm{\α}=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ρ}}_{\mathrm{\α}}{z}_{\mathrm{\α}}\right]=\underset{\mathrm{\α}=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ρ}}_{\mathrm{\α}}E\left[z_{\mathrm{\α}}\right]=\underset{\mathrm{\α}=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ρ}}_{\mathrm{\α}}E\left[z\left(x_{\mathrm{\α}}\right)\right]=\underset{\mathrm{\α}=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ρ}}_{\mathrm{\α}}m\left(x_{\mathrm{\α}}\right)=\\ \underset{\mathrm{\α}=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ρ}}_{\mathrm{\α}}\left[\mathrm{\α}\underset{l=0}{\overset{k}{\mathrm{\Σ}}}{a}_l{f}_l\left(x_{\mathrm{\α}}\right)\right]==\underset{l=0}{\overset{k}{\mathrm{\Σ}}}{a}_l\left[\underset{\mathrm{\α}=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ρ}}_{\mathrm{\α}}{f}_l\left(x_{\mathrm{\α}}\right)\right]\end{array}$

由于E[m^*(x)]＝m(x)E[m^*(x)]＝m(x)对任一组系数a₀，a₁，a₂，…a_k均成立，故有： 

$\underset{l=0}{\overset{k}{\mathrm{\Σ}}}{a}_l\left[\underset{\mathrm{\α}=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ρ}}_{\mathrm{\α}}{f}_l\left(x_{\mathrm{\α}}\right)\right]=\underset{l=0}{\overset{k}{\mathrm{\Σ}}}{a}_l{f}_l\left(x\right)$

对任一组系数a₀，a₁，a₂，…a_k均成立，因此必定有： 

$\underset{\mathrm{\α}=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ρ}}_{\mathrm{\α}}{f}_l\left(x_{\mathrm{\α}}\right)=f_l\left(x\right),l=\mathrm{0,1},...,k---\left(6\right)$

成立。式(6)共有k+1个式子，它称为无偏性(无泛性)条件组，也就是说，其中条件的个数与基函数f_l(x)的个数相等。 

②最小估计方差(或最优性)条件 

在上述k+1个无偏性条件下，用m^*(x)估计m(x)的估计方差为： 

$\begin{array}{c}{\mathrm{\α}}^2=E{[m^{*}\left(x\right)-m\left(x\right)]}^2=E{\{m^{*}\left(x\right)-E[m^{*}\left(x\right)\left]\right\}}^2=D\left[m^{*}\left(x\right)\right]=D\left[\underset{\mathrm{\α}=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ρ}}_{\mathrm{\α}}{z}_{\mathrm{\α}}\right]\\ E{[\underset{\mathrm{\α}=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ρ}}_{\mathrm{\α}}{z}_{\mathrm{\α}}-\underset{\mathrm{\α}=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ρ}}_{\mathrm{\α}}E\left(z_{\mathrm{\α}}\right)]}^2=E{\left\{\underset{\mathrm{\α}=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ρ}}_{\mathrm{\α}}\right[z_{\mathrm{\α}}-E\left(z_{\mathrm{\α}}\right)\left]\right\}}^2\\ E\left\{\underset{\mathrm{\α}=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ρ}}_{\mathrm{\α}}\right[z_{\mathrm{\α}}-E\left(z_{\mathrm{\α}}\right)\left]\right\}\left\{\underset{\mathrm{\β}=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ρ}}_{\mathrm{\φ}}\right[z_{\mathrm{\β}}-E\left(z_{\mathrm{\β}}\right)\left]\right\}=\end{array}$

$\underset{\mathrm{\α}=1}{\overset{n}{\mathrm{\Σ}}}\underset{\mathrm{\β}=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ρ}}_{\mathrm{\α}}{\mathrm{\ρ}}_{\mathrm{\β}}E[z_{\mathrm{\α}}-E\left(z_{\mathrm{\α}}\right)][z_{\mathrm{\β}}-E\left(z_{\mathrm{\β}}\right)]=\underset{\mathrm{\α}=1}{\overset{n}{\mathrm{\Σ}}}\underset{\mathrm{\β}=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ρ}}_{\mathrm{\α}}{\mathrm{\ρ}}_{\mathrm{\β}}\mathrm{Cov}(z_{\mathrm{\α}},z_{\mathrm{\β}})$

若将Cov(z_α，z_β)简记为C(x_α，x_β)则有： 

${\mathrm{\α}}_E^2=\underset{\mathrm{\α}=1}{\overset{n}{\mathrm{\Σ}}}\underset{\mathrm{\β}=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ρ}}_{\mathrm{\α}}{\mathrm{\ρ}}_{\mathrm{\β}}C(x_{\mathrm{\α}},x_{\mathrm{\β}})---\left(7\right)$

在用式(6)所表示的无偏性条件组的约束下，使式(7)所表示的估计方差达到极小，就要利用求条件极值的拉格朗日乘数法，即要求使 

$Q=\underset{\mathrm{\α}=1}{\overset{n}{\mathrm{\Σ}}}\underset{\mathrm{\β}=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ρ}}_{\mathrm{\α}}{\mathrm{\ρ}}_{\mathrm{\β}}C(x_{\mathrm{\α}},x_{\mathrm{\β}})-2\underset{l=0}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\μ}}_l[\underset{\mathrm{\α}=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ρ}}_{\mathrm{\α}}{f}_1\left(x_{\mathrm{\α}}\right)-f_l\left(x\right)]=\min$   (极小) 

时的ρ_α和μ_l(α＝1，2，…，n，l＝0，1，2，…，k)。 

具体方法是要求出 $\frac{\∂Q}{\∂{\mathrm{\ρ}}_{\mathrm{\α}}}{\mathrm{\μ}}_l(\mathrm{\α}=\mathrm{1,2},...,n),$ 并令其为零：同时也要求出 $\frac{\∂Q}{\∂{\mathrm{\μ}}_l}(l=\mathrm{0,1,2},...,k)$ 并令其为零。这样就可导出估计区域异常m(x)的方程组。 

③估计区域异常的方程组 

$\left\{\begin{array}{cc}\frac{\∂Q}{\∂{\mathrm{\ρ}}_{\mathrm{\α}}}=2\underset{\mathrm{\β}=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ρ}}_{\mathrm{\β}}C(x_{\mathrm{\α}},x_{\mathrm{\β}})-2\underset{l=0}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\μ}}_l{f}_l\left(x_{\mathrm{\α}}\right)=0& (\mathrm{\α}=\mathrm{1,2},...,n)\\ \frac{\∂Q}{\∂{\mathrm{\ρ}}_{\mathrm{\α}}}=-2[\underset{\mathrm{\α}=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ρ}}_{\mathrm{\α}}{f}_l\left(x_{\mathrm{\α}}\right)-f_l\left(x\right)]& (l=\mathrm{0,1},...,k)\end{array}\right.$

或即 

$\left\{\begin{array}{cc}\underset{\mathrm{\β}=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ρ}}_{\mathrm{\β}}C(x_{\mathrm{\α}},x_{\mathrm{\β}})-\underset{l=0}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\μ}}_l{f}_l\left(x_{\mathrm{\α}}\right)=0& (\mathrm{\α}=\mathrm{1,2},...,n)\\ \underset{\mathrm{\α}=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ρ}}_{\mathrm{\α}}{f}_l\left(x_{\mathrm{\α}}\right)=f_l\left(x\right)& (l=\mathrm{0,1},...,k)\end{array}\right.---\left(8\right)$

这就是估计区域异常m(x)的方程组，共有n+k+1个方程。其中，协方差矩阵为 当没有两点x_α和x_β完全重合时，永远是严格正定的。 

上述方程组可写成分块矩阵形式： $\left[\begin{array}{cc}C& f^{\′}\\ f& 0\end{array}\right]=\left[\begin{array}{c}\mathrm{\ρ}\\ -\mathrm{\μ}\end{array}\right]=\left[\begin{array}{c}0\\ f_x\end{array}\right]$

其中， $C=\left[\begin{array}{ccc}C(x_1,x_1)& ...& C(x_1,x_n)\\ .& .& .\\ .& .& .\\ .& .& .\\ C(x,x_1)& ...& C(x_n,x_n)\end{array}\right],$ $f=\left[\begin{array}{ccc}{f}_0\left(x_1\right)& ...& f_0\left(x_n\right)\\ f_1\left(x_1\right)& ...& f_1\left(x_n\right)\\ .& .& .\\ .& .& .\\ .& .& .\\ f_k\left(x_1\right)& ...& f_k\left(x_n\right)\end{array}\right]$

$\mathrm{\ρ}=\left[\begin{array}{c}{\mathrm{\ρ}}_1\\ {\mathrm{\ρ}}_2\\ .\\ .\\ .\\ {\mathrm{\ρ}}_n\end{array}\right],$ $\mathrm{\μ}\left[\begin{array}{c}{\mathrm{\μ}}_0\\ {\mathrm{\μ}}_1\\ .\\ .\\ .\\ {\mathrm{\μ}}_k\end{array}\right],$ $f_x=\left[\begin{array}{c}{f}_0\left(x\right)\\ f_1\left(x\right)\\ .\\ .\\ .\\ f_k\left(x\right)\end{array}\right]$

可以证明n+k+1个未知数(n个权系数ρ_α和k+1个拉格朗日参数μ_l)，n+k+1个方程的方程组(8)，当且仅当诸基函数f_l(x)(例如在一维时，可为1，x，x²等)在信息样品点集s上为代数独立时才有解，并有唯一解。 

诸f_l(x)在s上代数独立(即线性无关)是指：若存在诸s_l使得 

$\underset{l=0}{\overset{k}{\mathrm{\Σ}}}{c}_l{f}_l\left(x_{\mathrm{\α}}\right)=0,(\mathrm{\α}=\mathrm{1,2},...,n),$ 则必有c_l＝0(l＝0，1，2，…，k) 

由于在协方差函数存在的情况下，结构函数也存在，故方程组(8)还可用结构函数写为如下形式： 

$\left\{\begin{array}{cc}\underset{\mathrm{\β}=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ρ}}_{\mathrm{\β}}\mathrm{\γ}(x_{\mathrm{\α}},x_{\mathrm{\β}})+\underset{l=0}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\μ}}_l{f}_l\left(x_{\mathrm{\α}}\right)=C\left(0\right)& (\mathrm{\α}=\mathrm{1,2},...,n)\\ \underset{\mathrm{\α}=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ρ}}_{\mathrm{\α}}{f}_l\left(x_{\mathrm{\α}}\right)=f_l\left(x\right)& (l=\mathrm{0,1},...,k)\end{array}\right.---\left(9\right)$

2结构函数法的优点 

趋势(面)分析方法是结构函数法的一个特例。 

当局部异常R(x)相当于独立随机变量时，也可用最小二乘法估计区域异常。当局部异常R(x)相当于独立随机变量时，局部异常的结构函数γ_R(h))为： 

可以证明物探重力观测数据z(x)的结构函数就等于其局部异常的结构函数γ_R(h)(限于篇 幅，证明从略)。 

此时，式(9)变为： 

$\left\{\begin{array}{cc}{c}_0\underset{\mathrm{\β}\≠\mathrm{\α}}{\overset{n}{\underset{\mathrm{\β}=1}{\mathrm{\Σ}}}}{\mathrm{\ρ}}_{\mathrm{\β}}+\underset{l=0}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\μ}}_l{f}_l\left(x_{\mathrm{\α}}\right)=0& (\mathrm{\α}=\mathrm{1,2},...,n)\\ \underset{\mathrm{\α}=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ρ}}_{\mathrm{\α}}{f}_l\left(x_{\mathrm{\α}}\right)=f_l\left(x\right)& (l=\mathrm{0,1},...,k)\end{array}\right.---\left(10\right)$

这是一个与结构函数无关的(n+k+1)个方程，(n+k+1)个未知数ρ_α、μ_l(α＝1，2，…，n，l＝0，1，2，…，k)的线性方程组。解这个方程组就相当于用最小二乘法解正规方程组。故在这种特殊的情况下，用上述结构函数法计算出的区域异常值完全等同于在趋势(面)分析方法中用最小二乘法计算出的区域异常值换句话说，在趋势(面)分析方法中，用最小二乘法估计区域异常m(x)，只是当局部异常R(x)为独立随机变量时用结构函数法估计区域异常m(x)的一个特例而已。但是，在大多数实际情况中，局部异常R(x)是完全独立的随机变量的情况是极少见的。因此，应用最小二乘法来估计区域异常就有很大的局限性了。从这里我们也可看出结构函数法比趋势(面)分析方法适应性更大，因而也更优越。 

3实际资料处理结果 

为了检验结构函数法提取断层区域异常的效果，本文对某工区的一幅实际物探重力异常图做了实验。实验结果在附图中，实验结果表明，用结构函数法提取出的反映断层的重力区域异常图在边部和里边都是比较合理的，在边部没有突变。用结构函数法提取出的反映断层(粗线)的区域异常走向和工区主体断层异常走向吻合得较好。分离提取的断层没有多余的(假断层)，也没有漏掉真断层。 

四.附图说明

图1物探重力观测数据位置图。 

在图1中，已知平面上S₀，S_PS₁，S₂，S₃，S₄六处的物探重力观测数据分别为z₀，z_P，z₁，z₂，z₃，z₄，要估计出S₀和S_P点处的重力区域异常与局部异常。对具有某种结构的物探重力观测数据来说，S₀点待估的重力区域异常值为局部异常为而对S_P点，则有一组不同的系数λ₁，λ₂，λ₃，λ₄与其对应，使得区域异常局部异常为应用结构函数法，诸系数λ₁，λ₂，λ₃，λ₄可以利用求条件极值的拉格朗日乘数法求得。 

为了检验结构函数法提取断层区域异常的效果，本文对某工区的一幅实际物探重力异常图做了实验。图2是某工区的实际物探重力异常图，是重力异常等值线图，它是区域异常和局部异常的叠加。图3是用本文给出的结构函数法，从图2中提取出的反映断层的重力区域异常图(细线为等值线图)和断层(粗线为断层，共7条)。 

从图3可见，用结构函数法提取出的反映断层的重力区域异常图(细线为等值线图)在边部和里边都是比较合理的，在边部没有突变。图3中的反映断层(粗线)的区域异常走向和工区主体断层异常走向吻合得较好。分离提取的断层没有多余的(假断层)，也没有漏掉真断层。 

图2某工区的实际物探重力异常图(等值线图，是区域异常和局部异常的叠加) 

图3结构函数法从图2中提取的反映断层的重力区域异常图(细线为等值线图)和断层(粗线为断层，共7条)。 

具体实施方式

下面结合附图，对本发明进一步说明。为了检验本发明------结构函数法提取断层区域异常的效果，本文对某工区的一幅实际物探重力异常图做了实验。图2是某工区的实际物探重力异常图，是重力异常等值线图，它是区域异常和局部异常的叠加。图3是用本文给出的结构函数法，从图2中提取出的反映断层的重力区域异常图(细线为等值线图)和断层(粗线为断层，共7条)。 

从图3可见，用结构函数法提取出的反映断层的重力区域异常图(细线为等值线图)在边部和里边都是比较合理的，在边部没有突变。图3中的反映断层(粗线)的区域异常走向和工区主体断层异常走向吻合得较好。分离提取的断层没有多余的(假断层)，也没有漏掉真断层。 

Claims (3)

1.断层表现为重力异常梯度带，这种断层重力异常具有明显的空间结构性，即是线性的，具有很长的走向，在其走向方向上物探重力数据变化缓慢且是相关的，但在其垂直方向上数据变化剧烈且相关性很小，定量描绘这种断层结构性的“结构函数”的定义。

2.结构函数法的基本思想：在考虑了求反映断层的区域异常所用物探重力观测数据的形状、大小及其与待估区域异常值的点相互之间的空间分布位置等几何特征，以及重力观测数据的空间结构信息后，为了达到线性、无偏和最小估计方差的估计，而对每个求区域异常所用重力观测数据分别赋予一定的权系数，最后进行加权平均来估计出待估区域异常点处的重力区域异常值。

3.设反映断层的重力区域异常m(x)的估计量m^*(x)是n个已知物探重力观测数据z_α的线性组合，即其中ρ_α为权系数，用结构函数法求区域异常的问题就归结为求权系数ρ_α(α＝1，2，…，n)了，本发明是这样来选择这些权系数ρ_α的，即，在满足无偏性(也称为泛性)的条件下，求满足最小方差性(也称最优性)条件的权系数ρ_α，

①无偏性(泛性)条件：即E[m^*(x)]＝m(x)，E[z(x)]为z(x)的数学期望，

②最小估计方差(或最优性)条件：在k+1个无偏性条件下，用m^*(x)估计m(x)的估计方差为：

α²＝E[m^*(x)-m(x)]²。