发布时间 : 星期四 文章地球物理反演与层析成像-结业论文与程序设计 - 图文更新完毕开始阅读1981507bdd36a32d7375819c
********班 **** 学号:*********
基于惠更斯原理和费马原理求取地震波走时及其反射波射线路径的新方法。该方法具有原理简单、易于实现、能适应较为复杂地质模型以及易于将其推广到各向异性介质等优点，克服基本算法速度较慢的缺陷，是一种地震波走时和反射波射线路径计算的改进方法。在保证精度的条件下，该改进算法的计算速度显著提高。
地震波走时层析成像方法--交错网格法，即利用高密度网格进行射线追踪，以适应实测射线数的剖分网格进行成像计算，并采用任务并行化逐炮点进行射线追踪。交错网格法的成像网格单元内射线为曲线，具有较高的成像精度。方法原理简单，易于实现。
（2）波形拟和法
基于波动方程的层析成像一般有理论地震图法和接收函数法。由于波动方程数值模拟实质是求解地震波波动方程，因此模拟的地震波场含了地震波的所有信息，但由于基于波动方程的层析成像方法需要超大规模的三维数值计算，所以计算速度相对于几何射线法要慢，且易引进干扰波，目前还有许多困难问题没有解决。但波动方程包含了地震波场的全部信息，比仅利用走时资料仅用于模拟波的运动学特征的射线追踪层析成像更能客观地反映地下结构的信息，因此对于研究复杂条件下的各种波场最为有效，具有广阔的发展前景。目前常用的方法有伪谱法、有限元法、有限差分方法。伪谱法处理边界灵活，是有限差分法近似阶数趋于无限时的极限，它用快速傅氏变换来计算空间导数，计算精度要高于有限差分法。但是和有限差分方法一样计算量大，效率较低：有限元法由于剖分的任意性及它所依据的变分原理，对含有多种介质和自然边界条件的处理非常方便有效，已成为解决地震波传播数值模拟的一种重要方法。它是目前为止最精确的一种正演模拟方法，但计算量大。有限元法的主要优点是适宜于模拟任意地质体形态，可以任意三角形逼近地层界面，保证复杂地层形态模拟的逼真性。有限差分法和有限元法的主要缺点在于对高频分辨的限制，对地震勘探中典型的速度和频率，计算中需要大量的网格点，而伪谱法则相对更有效。
1．3反演及图像重建技术
层析成像中的反演方法可分为线性方法和非线性方法两种。目前非线性反演方法主要有速传算法、模拟退火法和神经网络法等。在体波层析成像中，使用线性反演方法的较多，如奇异值分解法、共扼梯度法和最小二乘法等。最小二乘法中加上阻尼得到了阻尼最小二乘法，其实质就是用三角矩阵分解法加上阻尼最小二乘法的超定线性方程组求解。地球物理问题大多是高度非线性问题，如层析成像问题就是一个典型的非线性问题。非线性反演方法全局搜索，且不依赖于初始模型，适用于对被研究区域的初始信息了解较少的研究，可以把其反演的结果作为初始模型进行局部最优图像重建。虽然计算速度较慢，但是对于复杂的非线性反演问题效果显著。而线性反演方法大多是人为的使非线性问题线性化，存在较大的不稳定性，且依赖于初始模型，但是其计算速度较快。因此，在层析成像中无论是线性反演方法还是非线性反演方法，都有较多的应用。
1.4分辨率分析（反演结果的评价）
层析成像解的评价是层析成像研究的重要组成部分。通过对解的分析，我们可以了解结果的可靠性、分辨率及误差等重要信息，常用的方法有：
（1）射线密度法。通过衡量每个节点附近的射线数量作为解的可靠性的一
********班 **** 学号:*********
种评价。但这种评价方法仅给出解的可靠性的初步度量，对解的分辨率的评价还需要进一步分析。
（2）线性反演理论方法。该方法是用广义线性反演理论，利用模型分辨率矩阵、数据分辨率矩阵和协方差矩阵来描述解的评价方法。
（3）尖峰试验法。该方法是通过使用合成数据去获得分辨率矩阵的列矢量，以测试方程组的病态对解的歪曲效应。但是，它只能估计解对单个参数点的分辨能力，无法评价整个解的可靠性。该方法主要目的是研究已知资料的某种形状的异常是否可以分辨。这种试验可以提供有关短波异常图像的成像能力，从而有助于区分垂向分辨率和横向分辨率的优劣，还可以对不同大小、不同形状异常体进行尖峰试验，以便来检验算法和数据对这种异常体的成像能力。
（4）棋盘分辨率试验法。该方法的基本原理是，首先用一个人工合成数据集代替已有的观测数据集。合成数据集由在一个特定的三维速度模型下，应用真实的射线分布计算得到的理论走时值构成。这一特定三维网格模型的速度分布（即棋盘格）是在初始一维速度模型基础上，加上规则分布的扰动值所构成（即各节点的扰动值大小相同，但正负相间地顺序排列，这样做的目的是便于分析）。然后，对合成数据集作反演计算，并把反演结果的三维速度结构与检测板的相似程度作为解的可靠性和分辨率的估计。由于棋盘格实验方法直观、实用，便于分析对比，现今大部分层析成像的结果都应用此方法进行评价。
（5）恢复分辨率实验法。其基本原理是，应用反演得到的结果作为人工合成模型，在此模型中进行射线追踪，计算走时，同时在这个数据中加入与真实数据误差相同数量级的随机误差，得到人工合成数据集。反演此数据集，得到的结果再与真实结果进行对比，以便分析图像恢复的情况。棋盘格实验设置的规则扰动不能模拟复杂模型情况下的图像恢复情况，相对而言，恢复分辨率实验要比棋盘格实验更接近真实情况，但不易分析和对比。
2． 地震层析成像技术的未来发展趋势
在以往对中国及邻近区域的地壳上地慢三维结构的研究中，人们比较注重地壳上地慢三维速度结构的研究。一方面是因为用速度结构可以解释地球内部三维结构某些方面的特征，另一方面是因为计算速度参数相对而言涉及的因素比较少。但是，从总的方面来看，速度参数只是利用了地震波的运动学信息，而忽略了地震波动力学的特征信息。地震层析成像研究大多根据地震记录局部上的单一观测值反演单一物理量，方法各自独立，表现出全局性和系统性的不足，在相当大的程度上阻碍了地震层析成像方法的使用。地震波动是岩石物性的综合响应，岩石物性参数间存在着必然联系，单一参数反演方法对多参数多分量综合响应的分析显得非常无力，因此，多参数多分量同步反演方法的研究非常必要。在以往的地震层析成像方法中，走时信息和振幅信息分别反演岩石的波速与衰减。不论走时和振幅都只能代表地震波动的局部特征，而地震波形的整体则是地下地质状况和岩石物性的综合响应。多分量波动层析成像突破以往只提取局部信息的作法，利用地震波动的整体信息，通过对理论地震波动数据与实测地震波动数据之间波形残差的量化分析，修改弹性波方程中的物理参数，进而修改地质模型的岩石物性参数。通过这样的研究，地震层析成像研究方法将发生实质性变化，它将使以地震波局部特征为主的研究转变为以地震波场整体动力学特征为主的研究。
********班 **** 学号:*********
地震走时层析成像
内容： 程序解释
（1）正演程序
#include\#include\struct data { double begin; double end; double slope; double length;
double transform[12][9]; double transform2[108]; int n;
double time; struct point {
double x; double y; }point[100]; }ray[144];
结构体ray[i]为第i条射线的信息: begin 射线的起点纵坐标 end 射线的终点纵坐标 slope 射线的斜率 length 射线的总长度
transform【i】为射线在每个单元格内大小 共有两种形式 二维数组为矩阵形式，一维数组是为了方便计算而把二维数组转化成的一位数组
n 计算中的中间数据，物理意义是射线所穿透的横轴的个数，为正是向下穿，为负的是向上穿
time 整个射线的走时
point 射线与纵横轴的交点（X，Y）
struct ww {
double x; double y;
}d[100]; 结构体 d
为了给交点排序的中间数据 main()
********班 **** 学号:*********
{ FILE *fp,*fp0,*fp1,*fp2,*fp3,*fp4,*fp5,*fp6; int i,j,y,m,b,e,n,q,c,k,a;
double mm,temp1,temp2; temp1 temp2 为为了给点排序的中间数据 double x[108],xc[12][9];
for(i=0;i<108;i++)x[i]=1.0/3000.0; x[20]=x[29]=1.0/5000.0; x[77]=x[78]=1.0/2000.0;
/*----------------------speedmatrix-------------------------------*/ 得出慢度矩阵
for(i=0;i<12;i++)
for(j=0;j<9;j++)xc[i][j]=x[9*i+j]; fp0=fopen(\ for(i=0;i<12;i++) {
for(j=0;j<9;j++)
fprintf(fp0,\ }
/*-------------------the144rays-----------------------------------*/ 144条射线的 起点 终点 长度 斜率 的计算 for(i=0;i<12;i++) {
for(j=0;j<12;j++) {
ray[12*i+j].begin=1.5+3.0*i; ray[12*i+j].end=1.5+3.0*j;
ray[12*i+j].slope=(ray[12*i+j].begin-ray[12*i+j].end)/45;
ray[12*i+j].length=sqrt(45*45+pow(fabs(ray[12*i+j].begin-ray[12*i+j].end),2)); } }
/*-------------------dataout(ray)---------------------------------*/输出以上数据 便于验证
fp=fopen(\射线.txt\ for(i=0;i<144;i++) {
fprintf(fp,\e,ray[i].length); }
/*---------------------dataout(speed)-----------------------------*/ 输出慢度 与所给验证
fp1=fopen(\慢度.txt\