层析成像的理论基础

如题所述

1.Radon变换和地震层析成像

设二维域中有一未知的连续函数f(x,y),令Oxy坐标系逆时针旋转φ,形成OPS坐标系,如图13-1-1所示。将f(x,y)沿平行于s轴的L做线积分,设为

物探数字信号分析与处理技术

式中直线L1依下述方程定义

物探数字信号分析与处理技术

有许多种地震观测都可以表示成形如(13-1-1)所示的线积分,例如f(x,y)是慢度(速度的倒数)u,那么 就表示走时,即

物探数字信号分析与处理技术

式(13-1-1)可用δ函数表示:

物探数字信号分析与处理技术

式中 是在P轴正方向的单位矢量, 为r(x,y)在P轴上的投影。对于 ,一般情况 。综合(13-1-1)式和(13-1-4)式有

物探数字信号分析与处理技术

为f(x,y)在角度为φ时沿射线L1的投影值,沿许多射线L1,L2,…,Lt投影,构成投影函数,如图13-1-1所示。

物探数字信号分析与处理技术

图13-1-1 直射线投影中的几何分布与变量

改变φ角便得到一系列投影函数 。我们把f(x,y)称为模型或目标函数,式(13-1-5)称为拉冬正变换。当用 来推断“模型”,即重复(13-1-5)式时,就称为成像或图像重建。有好几个公式可以对(13-1-5)式做反变换,也就是由 求得f(x,y)。虽然从解析上看这些公式都是等价的,但对离散数据来讲,它们的数值结果并不相同。正是这些差异反映出数值反演的困难。

让投影函数与某一核函数褶积,可求出褶积后的投影函数:

物探数字信号分析与处理技术

将一系列Gφ(t)作反投影,得目标函数:

物探数字信号分析与处理技术

(13-1-7)式的运算可看作对每一个投影角φ作反投影,即滤波后的数据要沿着射线方向被投影回去。对所有的投影而言,要对这些函数积分(或求和)。我们称(13-1-7)式为Radon递变换(简记成IRT)。

反投影成像时,为消除“月晕效应”,需做适当的滤波,滤波反投影形成的Radon递变换为

物探数字信号分析与处理技术

实际应用时,需将(13-1-8)式离散化、有限化。投影装置除平行射线束外,实际用得更多的是扇形射线束,图13-1-2是扇形射线束投影示意图。

图13-1-2 扇形射线投影

图13-1-3 坐标系示意图

2.投影切片定理

设目标函数为f(x,y),投影函数为 满足Radon变换,在图13-1-3坐标系中,投影函数的一维傅氏变换为

物探数字信号分析与处理技术

在OUV坐标系中,f(r,s)的二维傅氏变换为

物探数字信号分析与处理技术

在图13-1-3中,Oxy为坐标原点,角度为φ的直线上,u=0,所以

物探数字信号分析与处理技术

式(13-1-11)称为投影切片定理,即角度为φ的投影函数的一维傅氏变换,等于目标函数的二维傅氏变换的角度φ所做的切片。

式(13-1-11)的投影切片定理给出一种图像重建的方法。对不同投影的函数作一维傅氏变换,可获得波数域的极坐标表示的f*(v,φ),即对于一个固定的投影角φ而言,投影数据的一维傅氏变换可以给出模型的傅立叶变换在一通过波数域原点的切片上的值。在对式(13-1-11)离散并进行数值反演时,必须对极坐标网格上的数据进行内插,得到所需的域中矩形网格上的数值f*(u,v),最后用二维傅立叶变换获得目标函数f(x,y)的数字化图。

在地震勘探中,地震记录就是由震源激发的波对穿过介质区域(速度场)的投影。通过对来自不同方向的投影(地震记录)进行傅氏变换,便得到研究区的坐标函数(速度场)的二维傅氏变换,再求傅氏二维逆变换,就得到研究区域的速度分布。

3.地震波的射线理论

由于实际地震波的射线并非直线,故弹性波的观测旅行时间t可表示为沿弹性波射线的线积分,即

物探数字信号分析与处理技术

式中v为弹性波的传播速度,ds为波射线弧元。如果引入慢度 ,则(13-1-12)式可表示为

物探数字信号分析与处理技术

用线性变换表示上式:

物探数字信号分析与处理技术

R为与射线路径有关的算子,若波射线为直线,称为Radon算子。Radon算子R的逆变换写作R-。如果波线为直线,在物体内任意点对来自所有方向的波线,都能得到一个相应的弹性波旅行时间t,且没有误差,则拉冬算子的逆变换R-可写作为

物探数字信号分析与处理技术

式中的B为投影算子,D是微分算子,H为希尔伯特变换算子。

式(13-1-15)就是CT技术的理论基础,若将其引用到地震层析技术中,需要寻求数据处理方法。如对(13-1-14)式离散化,则局部范围内曲线的弧元可作直线处理,其旅行时间不会引起很大误差。另外,如果在物体周围布置的发射点和接收点足够密,则对所有方向的波的旅行时间均可得到。这在CT中容易得到,而在地震层析中,对不同观测系统、不同数据量对层析成像的影响可预先进行评价。

进行地震层析数据处理时,需将研究区分割成很多网格(像元),设第j个网格的慢度为uj,对实测波至时间也离散化,设射线i在第j个网格上的路径长度为δij,则式(13-1-14)可离散化为

物探数字信号分析与处理技术

式中ti表示第i条波线从震源到接收点的旅行时间。设有I条波线,J个格子,则(13-1-14)中的R可理解为向量t与U之间的线性变换的矩阵:

物探数字信号分析与处理技术

如果波线已知,矩阵R就可确定,对(13-1-14)进行反演,就归结为求R的逆矩阵R-。当然,由于观测数据量很大,而且观测值带有误差,就使形成的矩阵大型稀疏、强超定或欠定和非一致,因此不能直接求取R-。这就需要开发适用于地震层析问题的近似解法,如迭代处理法。由于模型和数据都存在着误差,通常人们并不关注那种高精度的解,而实用的迭代方法经几何运算能给出一个可接受的解。

实际工作中,波射线一般是未知的。迭代处理时,一般涉及行处理法(RA)或代数重建法(ART)。如能给出充分接近真实慢度分布的Ui,由Ui通过波射线法求Ri+1,再由上述反演求Ui+1,将地基模型化,反复迭代,使由慢度计算的走时与观测走时t充分接近。只要初始分布U0充分接近真分布U,这种迭代方法就能够求解。

温馨提示:答案为网友推荐,仅供参考
相似回答