APP下载

基于matlab的地震资料反演模拟

2020-07-08孙福玉

中国科技纵横 2020年6期
关键词:最小二乘法

孙福玉

摘 要:本文对地球物理学反演的意义做了简要说明。以地球物理反问题中的地震波褶积模型为例,通过理论分析和数学推导给出超定问题的最小二乘解法、先验信息的应用、线性反问题的广义反演法等多种常用方法,并通过matlab编程实现做最终成图对比分析。

关键词:超定问题;matlab编程;最小二乘法;先验信息;广义反演

中图分类号:P631.4 文献标识码:A 文章编号:1671-2064(2020)06-0186-02

0 引言

在自然界中,存在大量的客观事物,这些客观事物在传统物理学中,常常可以看成是一个系统。这些事物包罗万象,小到原子大到宇宙。一般来说,所有的客观事物都能用一定的物理量来描述,但不是所有的物理系统所包含的信息都能够被直接测量,有一部分物理系统的物理量只能通过间接测量。在地球物理学领域的研究中,研究对象绝大多数都位于地表以下。例如,研究地下地层的构造特征、寻找含油储集层、寻找矿产资源、寻找含水层。将地下介质全部挖开去直接研究这些地质信息显然是不可能的。为了全面的了解一个物理系统,往往希望获得尽可能多的描述系统的物理量,往往能够直接获取的物理量十分有限。作为一种演绎性的研究手段,地球物理反演就是通过不断地推演系统中那些可以被直接观测的物理量与不可直接测量的物理量之间的联系。达到通过已知信息去探索未知信息的一种方法。

反演问题一直以来都是地球物理学中的核心理论问题。由于人类自身的局限性,获取深地信息的来源大都来源于地球表面,有极小一部分资料来自于对地球的直接探测,例如钻井信息。所以人类对于地球的认识和研究终究还是要依靠有限的地球物理信息来实现,这些信息如何应用到实际研究中,就是反演问题的科研本质。本文通过将理论方法用matlab编程实现,以反演地震子波为例,模拟地球物理反演过程。

1 超定问题的常规解法[1]

通常,用物理模型去求解物理參数的过程叫做正演,而给定一系列物理模型的参数信息去求解物理模型的过程叫做反演。反问题的研究前提是正演问题得到合理的解释。正演问题大体可分为两方面的研究,一是通过实验的方法寻找出经验公式,利用此公式建立从模型到数据的联系。二是通过严密的物理推导,通过给出的模型参数预测观测数据。由此可见,只要正演问题不能被正确的解决,反演问题的研究就不能展开。此外,地下介质环境异常复杂,及时正演问题即物理理论已经被完善的解决,反演问题还是不能得出正确答案,需要考虑多种假设条件,通常会添加一些限制条件得出某种意义下的最佳解。

对于线性离散问题,总可以描述为d=GM的形式,可以理解为一系列的线性方程组的通式,其中为已知的观测数据为待求的模型参数。 为数据核矩阵,与模型有关。

设观测数据的数量为M,待定的模型参数数量是为N,G为M×N阶矩阵,其秩为r,当M>N=r时,观测数据提供了多于模型参数数目的信息,此问题称为超定问题,求解不同类型的线性反演问题,所采用的方法是有不同的。在本文中着重探讨一下,超定问题的不同解法。

超定问题是观测数据提供了多于模型参数数目的信息,因此求得的每一个G×M都会和真实观测数据存在误差。要想寻求超定问题真解是不可能的,只能寻求其在不同条件下的最佳解。在这里通常是去构造最小误差解,也就是说寻找到这样一组模型参数,它和G核矩阵相乘所得的计算值和真正的观测数据之间的差异最小。这个差异通常用误差的二范数平方来表征。二范数平方的最小值解就是常用的最小二乘解。

所以对于超定问题的求解方案可以这样描述:期望求得一组与观测数据误差平方和为最小的预测数据所对应的模型参数[3]。

为了不失普适性,我们先假设观测数据d为M维向量,模型参数m为N维向量,且有M>N,则上式的求解可转化为一个线性方程组的求解[3],即:

通过目标函数对模型参数求一阶导函数可求得它的极值点,不难验证这个极值点即为极小值点,也就是最小值点。解得:

下面看一个地球物理学上的实际问题,地震波褶积模型。由于地震子波是延续时间非常短的信号,地震记录是数据相对丰富的信号。那么,利用反射系数(测井资料)和地震记录(井旁道地震资料)求地震子波这一反演问题,就成为经典的超定问题。应用matlab编程[2]实现这一实际问题的步骤如下:

(1)定义模型。设计一个地质模型,定义速度、密度、最小偏移距、采样间隔、道数等物理量。通过数学公式求出雷克子波,一定要注意雷克子波是一个延续时间非常短的信号,也就是十分短的向量。

(2)模型正演。因为反演问题的求解都是建立在正演问题充分解决的基础上的。所以要先求出正演数据。通过此模型生成反射系数序列,利用反射系数序列求出反射系数矩阵。正演实质上是反射系数与雷克子波向量的褶积。所以在此问题中的G核矩阵实际上就是反射系数矩阵。

(3)模型反演。利用最小二乘法解出地震子波并将其与雷克子波向量做对比。如图1为引入加权因子的反演结果对比。

2 先验信息的应用[1]

上文中的解决方案,是在观测数据基础上去反推模型参数。实际工作过程中要想全面确定模型参数大多数情况下都是信息不足的,为了使反演问题的解更切合实际情况,我们需要添加一些在观测数据中未包含的信息,这种附加给反演问题的信息叫先验信息。先验信息实际上是对观测数据的一种补充。为了使得上文中超定问题的最小二乘解更加真实,则可以通过添加先验信息来进行优化。先验信息的类型有许多种,但对于不同的形式,都要使得目标函数以定量的形式出现,并且是不依赖于实际数据的。

在本文中介绍超定问题常用的先验信息,引入加权因子。通过对观测数据的预测误差进行加权来度量。在实际工作中,经常会碰到某些观测数据的精度比另外一些观测数据的精度高的情况发生。在这种情况下,通常的解决方式是在总误差的定量表示中,在比较精确的观测数据所相应的预测误差上加的权大于不精确的观测数值所加的权。

为了完成这样的加权,通常要定义一个广义预测误差,为观测数据的加权因子,确定了每个单独的误差对总预测误差的相对贡献。加权因子是一个对角矩阵,对角线上的元素,就是相应观测数据误差所加权重的大小度量。因此可以通过广义预测误差最小来估计此时的模型参数。通过和上文类似的导函数推导即可求得解的形式。

加权因子的引入是提高某些数据的可信度,在地震波褶积模型问题中因为所有的数据都足够准确,在此基础上反演得到的模型不会因为加权因子的改变而改变,不管引入什么样的加权因子都可以反演得到雷克子波。因此在本文中为了方便读者对加权因子的引入有更加深入的理解,在matlab编程[2]时更改地震记录(观测数据)后引入加权因子。引入如图2所示的加权因子,目的是削弱错误数据的权重。

We=eye(600);We=(200,200)=0.2;We=(290,290)=0.1;

如圖2所示,引入加权因子后削弱了错误数据的权重,可以明显地看出反演出来的模型更加接近雷克子波。

3 线性反问题的广义反演[1]

首先本文先介绍什么是广义逆,设矩阵存在矩阵满足一部分的矩阵称为G的广义逆。可见,G的广义逆不唯一共有15种。同时满足所有条件的广义逆,称为Moore-Penrose逆,记为G+。对于任意矩阵的Moore-Penrose广义逆G+必然存在,并且是唯一的。这样的性质才使得我们讨论问题有意义。

本文着重讨论广义逆的求法之一,奇异值分解方法。任意一个实数矩阵A都可以写成A×A转置的特征向量构成的列向量乘上A×A转置特征值的平方根构成的矩阵乘上A转置×A的特征向量构成的列向量三个向量乘积的形式。

图中深色部分是奇异值已经它所对应的向量,由奇异值分解可以求得它的广义M-P逆,它的形式如图3。

通过公式推导不难得出超定、欠定、混定三种问题解的统一的表达式,也从另一个方面证明了解的统一性,实际上这都是基于2-范数。通过matlab编程[2]实现模型构建、G核矩阵的求取,以上两个步骤在上文已经给出,接下来运用matlab预定义的SVD函数求取A矩阵的奇异值分解,进而求出A矩阵的广义M-P逆矩阵。运用广义M-P逆的核矩阵带入通用公式即可求解。

如图4所示,运用广义M-P逆反演出的子波与原雷克子波模型完全一致。

参考文献

[1] 姚姚.地球物理反演基本理论与应用方法[M].武汉:中国地质大学出版社,2002.

[2] 童孝忠,柳建新.MATLAB程序设计及在地球物理中的应用[M].长沙:中南大学出版社,2013.

[3] 谢万学.井地联合反演方法研究[D].青岛:中国石油大学,2007.

猜你喜欢

最小二乘法
基于压电陶瓷的纳米定位与扫描平台模型辨识算法研究
基于惯导角度量测的轨道平面最佳线形参数估计算法
全国各省份经济发展影响因素的实证分析
手动求解线性回归方程的方法和技巧