地球物理反演与层析成像-结业论文与程序设计

地球物理反演与层析成像结业论文及地震走时层析成像程序
姓名：
学号：
班级：地球物理学*班
[在此处键入文档的摘要。

摘要通常是对文档内容的简短总结。

在此处键入文档的摘要。

摘要通常是对文档内容的简短总结。

]
概述
地学层析成像是用医学X射线CT的理论详细调查地下物性参数分布状况的物探技术。

分为地震层析成像、电磁波层析成像和电阻率层析成像。

地震层析成像就是用地震数据来反演地下结构的物质属性，并逐层剖析绘制其图像的技术。

其主要目的是确定地球内部的精细结构和局部不均匀性。

相对来说，地震层析成像较其他两种方法应用更加广泛，这是因为地震波的速度与岩石性质有比较稳定的相关性，地震波衰减程度比电磁波小，且电磁波速度快，不易测量。

地震层析成像按研究区域的尺度可分为全球层析成像、区域层析成像、局部层析成像：按所用资料的来源可分为天然地震层析成像（大尺度深部横向不均匀性研究）、人工地震测深（主要研究浅部界面分布）。

按所依据的理论基础一般分为基于射线方程的层析成像和基于波动方程的层析成像。

前者按射线追踪时所用的地震波资料的不同又可分为体波（反射波、折射波）和面波层析成像：按反演的物性参数区分，可分为利用地震波走时反演地震波速度的波速层析成像以及利用地震波振幅衰减反演地震波衰减系数的层析成像。

基于射线理论，地震波走时层析成像方法由于走时具有较高信噪比、无论是柱面波还是球面波走时的规律都相同等优点，相对来说发展较早，技术方法比较成熟，是目前地震层析成像的主要方法。

但是射线理论只适用于波速在一个波长范围内变化很小的场合，是波动方程的高频近似，因此它有一定的局限性。

而基于波动方程的层析成像方法由于需要超大规模的三维数值计算，目前还有许多问题没有解决。

但波动方程包含了地震波场的全部信息，比仅利用走时资料的射线追踪层析成像更能客观地反映地下结构的信息，因此是未来地震层析成像的主要发展方向。

层析成像技术能以图像的方式直观清晰地显示地下物质结构的属性，所以这种方法一产生就受到了极大关注，被广泛应用于内部地球物理和地球动力学、能源勘探开发、工程和灾害地质、金属矿勘探等领域。

地震层析成像技术起源于20世纪30年代，自该技术应用以来，已取得了很多重大的成果。

如以美国哈佛大学和加州理工学院为代表所做的全球三维层析成像工作，首次为人类提供了地球内部的三维结构影像图，其中最重要的结果是地震波速度成像结果与大地水准面的相关性，地球动力学对其给出了很好的解释，为板块运动的热对流学说提供了证据。

再如，用层析成像方法人们首次发现非洲超级地幔柱等大型地幔柱均起源于地幔边界。

在大洋洋脊、板块消减带、克拉通地区，地壳和上地幔中的火山、地壳和地幔顶部、造山带、断裂区和震源区等地方层析成像技术也都有大量的应用成果。

无论是能源和矿产等资源勘探，还是地球内部结构及地球动力学研究，地震层析成像技术都是有效的、重要的技术之一。

1地震层析成像技术
地震层析成像涉及3个方面：数据采集、数据处理（数据正反演计算和图像重建）、成像结果解释。

地震层析成像是采集数据的主要目的、数据解释的基础和数据处理的主要部分。

地震层析成像主要包括以下几部分：模型的参数化、正演计算地下介质属性的理论值（射线追踪、波形拟和）、反演及图像重建、反演结果的评价（分辨率分析）。

现分别就用于这4个方面的各种方法作一阐述。

1.1建模及模型的参数化
层析成像的结果是在初始模型的基础上迭代反演得来的，因此初始模型与真实地下结构接近程度直接关系到成像的结果能否准确反映客观物质属性。

如何合理、准确地描述初始模型至关重要。

早期研究一般都是假设模型为均匀层状水平各向同性介质模型，这只是一个粗略的模型，远远不能满足实际应用需要。

随着研究的深入，模型逐步过渡到三维非均匀各向异性任意界面介质模型。

国际上一些标准的模型有二维的Marmousi模型，三维的盐丘模型和逆掩模型等。

Gjoystdal（1985）提出的模型生成技术，可方便地生成任意复杂结构的壳幔模型。

块状建模的方法，对地质体的描述采取体→块→面→点→坐标的层次结构。

改变了传统的层状地层的建模方式，引进了块状结构的描述，并采用三角形面片构造块状模型的界面，可以适应非常复杂的三维介质。

新的建模方式从根本上改变了层状地层建模不能适应复杂地质结构的状况。

在地震层析成像技术中，由于最终反演的地下介质属性是通过将研究区域划分成不重叠的多个像元，依据各像元的灰度（反演得到的地下介质属性）来成图的，所以在地震层析成像中多采用网格的方法来进行模型参数化。

网格化方式也由最初网格内速度均匀分布模型发展到后来的给出节点速度值，采用插值的方法求得网格内各点的速度；由规则均匀网格发展到动态变尺度的不规则网格。

有人在2001年提出了一种模型参数化方法，把模型网格分为正演网格、地震网格（实际地下结构模型）、反演网格，各种网格划分密度和大小（根据射线分布情况）不同。

对正演网格采用精细划分，而对成像网格则采用相对较粗糙的划分，这种模型参数化方法在不增加计算复杂度的情况下，提高了成像的分辨率，是一种理想的模型参数化方法。

类似的研究有交错网格法，不同尺度的成像网格和射线追踪网格，彼此通过双曲线插值映射相互关联。

但他们采用尺寸相等的网格划分方法，对于射线分布不均匀的情况有可能造成某些网格没有射线通过的问题，虽然通过更小尺度的划分可以避免上述问题，但是这无疑增加了计算量，这是和采用不等尺寸网格划分方法相比不足的地方。

在正演数值模拟之前，还需要做的一项重要工作就是数据预处理。

地震层析成像结果的优劣除了跟初始模型的选取有关，很大程度上还取决于数据空间的完备程度。

如数据量的大小，数据的精度，射线分布的均匀程度及密度等。

这些对于人工地震资料来说，炮点和接收点是可以人为选择的，因此上述要求是可以得到满足的。

然而对于天然地震资料来说，只能通过数据预处理尽可能地提高成像精度，如震源深度校正、地震重新定位、时差校正、远震的高度校正和地球椭圆扁率的校正等。

1.2正演数值模拟正演计算
在层析成像中起着极其重要的作用。

正演计算的精度和计算速度，直接决定着成像的分辨率和可靠程度。

正演数字模拟技术分为求解偏微分方程的波动方程数值模拟和由积分方程以求解波场传播旅行时为主的射线追踪数值模拟。

下面分别就这两种数值模拟技术的主要方法进行评述。

（1）射线追踪数值模拟方法
射线追踪的方法种类较多。

经典的方法是基于初值问题的试射法和基于边值问题的弯曲法。

经典方法存在的不足有：难以处理介质中较强的速度变化，难以求出多值走时中的全局最小走时，计算效率较低。

而且，试射法不能对首波和阴影区内（射线理论不成立）的射线路径进行追踪；弯曲法对于两点距离较远的情况效率较低。

随着射线追踪方法的发展，出现了大量不同于传统方法的新型算法。

这些方法的主要特点在于不再局限于地震波的射线路径描述，而是直接从惠更斯原理或费马原理出发，采用等价的波前描述地震波场的特征。

1）有限差分求解程函方程法
Vidale基于扩张波前的思想提出了用有限差分法求解程函方程来进行射线追踪的方法，开辟了一条射线追踪的新途径，后又于1990年将该方法推广到三维。

但是他的方法仍然没有解决首波的射线追踪问题，同时当介质中存在较大的速度间断面时会出现不稳定（不满足因果性关系，出现负数开平方）。

将地震波的传播路径近似看成一条射线的情况下，波的旅行时只与沿着波射线的速度分布有关。

由于低频波的波前恢复效应，实际情况是当波穿过的介质远远小于菲涅尔带时，这种近似的影响可以忽略，然而当波穿过的介质大小和菲涅尔带相当时，旅行时和整个菲涅尔带内的速度分布情况相关，这时把波的传播路径近似成射线对旅行时所造成的误差是不可忽略的。

2）最短路径法
最短路径法的基础是费马原理及图论中的最短路径理论，是用网络节点之间的最小旅行时连线近似地震射线路径的。

这种方法可以同时计算出从震源到达空间所有点的初至走时及相应的射线路径，并且不受射线理论的约束，准确地追踪出阴影区内的折射波射线路径。

波速模型的复杂性与空间维数也不会影响算法的实现，而且所得初至走时保证了全局最小的特性。

这种算法灵活高效，实用性强并且克服了经典算法的缺陷，是一种较理想的算法。

但是最短路径法的精度和速度并不比传统方法强。

当节点较稀时，射线常常呈之"字形路径，计算出的旅行时将比实际旅行时系统偏大，且在波传播方向上节点越少，误差越大。

另一个问题是，当网络稀疏时，特别对于速度变化平缓的区域，在两个点之间常常会有几条等时最短路径，其中可能有一条能较好地近似真实射线路径，而其它的却不能，因此具有一定的不确定性。

3）解析计算法
解析的射线追踪方法本质上都是对射线方程解析求解，一般有以下几类。

基于费马原理的解析方法，即通过求解射线方程的极小值来求得射线的路径。

解析法适用范围较小，因为实际地质构造复杂，即速度分布比较复杂，而且解析法只能对少数特殊的速度分布实现射线追踪（如速度或慢度平方是常梯度，以及慢度平方是多项式的情形）。

4）传统方法的改进方法
拟弯曲法在模型中引入了复杂形状的速度间断面，使用了可以计算含有速度间断面的非均匀模型中的地震波走时和射线路径的射线追踪技术，适合复杂结构地区的地震成像。

该方法对射线路径的扰动主要分两种情况处理。

当扰动点在间断面上时，根据斯奈尔定律，采用两分法逐渐缩小范围，找出在间断面上的折射点：当扰动点在连续介质内部时，以射线路径上3个相邻点为例，先固定不相邻两点的位置，扰动中间点的位置，以寻求不相邻两点之间最小走时路径。

由于这些特点，该算法近年来在国内外获得了较广泛的应用。

5）其他方法
旅行时线性插值法是Vidale差分法的一种高级形式。

为了提高地震波旅行时的计算精度，提出了一种改进的射线追踪法——旅行时二次/线性联合插值法（QLTI），即在震源的近场采用二次插值计算旅行时，在远场仍用LTI算法。

QLTI 算法改善了近场旅行时的精度，降低了累加误差，从而提高了全场（尤其是远场）旅行时的计算精度，QLTI算法较传统的LTI明显提高了精度。

基于惠更斯原理和费马原理求取地震波走时及其反射波射线路径的新方法。

该方法具有原理简单、易于实现、能适应较为复杂地质模型以及易于将其推广到各向异性介质等优点，克服基本算法速度较慢的缺陷，是一种地震波走时和反射波射线路径计算的改进方法。

在保证精度的条件下，该改进算法的计算速度显著提高。

地震波走时层析成像方法--交错网格法，即利用高密度网格进行射线追踪，以适应实测射线数的剖分网格进行成像计算，并采用任务并行化逐炮点进行射线追踪。

交错网格法的成像网格单元内射线为曲线，具有较高的成像精度。

方法原理简单，易于实现。

（2）波形拟和法
基于波动方程的层析成像一般有理论地震图法和接收函数法。

由于波动方程数值模拟实质是求解地震波波动方程，因此模拟的地震波场含了地震波的所有信息，但由于基于波动方程的层析成像方法需要超大规模的三维数值计算，所以计算速度相对于几何射线法要慢，且易引进干扰波，目前还有许多困难问题没有解决。

但波动方程包含了地震波场的全部信息，比仅利用走时资料仅用于模拟波的运动学特征的射线追踪层析成像更能客观地反映地下结构的信息，因此对于研究复杂条件下的各种波场最为有效，具有广阔的发展前景。

目前常用的方法有伪谱法、有限元法、有限差分方法。

伪谱法处理边界灵活，是有限差分法近似阶数趋于无限时的极限，它用快速傅氏变换来计算空间导数，计算精度要高于有限差分法。

但是和有限差分方法一样计算量大，效率较低：有限元法由于剖分的任意性及它所依据的变分原理，对含有多种介质和自然边界条件的处理非常方便有效，已成为解决地震波传播数值模拟的一种重要方法。

它是目前为止最精确的一种正演模拟方法，但计算量大。

有限元法的主要优点是适宜于模拟任意地质体形态，可以任意三角形逼近地层界面，保证复杂地层形态模拟的逼真性。

有限差分法和有限元法的主要缺点在于对高频分辨的限制，对地震勘探中典型的速度和频率，计算中需要大量的网格点，而伪谱法则相对更有效。

1．3反演及图像重建技术
层析成像中的反演方法可分为线性方法和非线性方法两种。

目前非线性反演方法主要有速传算法、模拟退火法和神经网络法等。

在体波层析成像中，使用线性反演方法的较多，如奇异值分解法、共扼梯度法和最小二乘法等。

最小二乘法中加上阻尼得到了阻尼最小二乘法，其实质就是用三角矩阵分解法加上阻尼最小二乘法的超定线性方程组求解。

地球物理问题大多是高度非线性问题，如层析成像问题就是一个典型的非线性问题。

非线性反演方法全局搜索，且不依赖于初始模型，适用于对被研究区域的初始信息了解较少的研究，可以把其反演的结果作为初始模型进行局部最优图像重建。

虽然计算速度较慢，但是对于复杂的非线性反演问题效果显著。

而线性反演方法大多是人为的使非线性问题线性化，存在较大的不稳定性，且依赖于初始模型，但是其计算速度较快。

因此，在层析成像中无论是线性反演方法还是非线性反演方法，都有较多的应用。

1.4分辨率分析（反演结果的评价）
层析成像解的评价是层析成像研究的重要组成部分。

通过对解的分析，我们可以了解结果的可靠性、分辨率及误差等重要信息，常用的方法有：（1）射线密度法。

通过衡量每个节点附近的射线数量作为解的可靠性的一
种评价。

但这种评价方法仅给出解的可靠性的初步度量，对解的分辨率的评价还需要进一步分析。

（2）线性反演理论方法。

该方法是用广义线性反演理论，利用模型分辨率矩阵、数据分辨率矩阵和协方差矩阵来描述解的评价方法。

（3）尖峰试验法。

该方法是通过使用合成数据去获得分辨率矩阵的列矢量，以测试方程组的病态对解的歪曲效应。

但是，它只能估计解对单个参数点的分辨能力，无法评价整个解的可靠性。

该方法主要目的是研究已知资料的某种形状的异常是否可以分辨。

这种试验可以提供有关短波异常图像的成像能力，从而有助于区分垂向分辨率和横向分辨率的优劣，还可以对不同大小、不同形状异常体进行尖峰试验，以便来检验算法和数据对这种异常体的成像能力。

（4）棋盘分辨率试验法。

该方法的基本原理是，首先用一个人工合成数据集代替已有的观测数据集。

合成数据集由在一个特定的三维速度模型下，应用真实的射线分布计算得到的理论走时值构成。

这一特定三维网格模型的速度分布（即棋盘格）是在初始一维速度模型基础上，加上规则分布的扰动值所构成（即各节点的扰动值大小相同，但正负相间地顺序排列，这样做的目的是便于分析）。

然后，对合成数据集作反演计算，并把反演结果的三维速度结构与检测板的相似程度作为解的可靠性和分辨率的估计。

由于棋盘格实验方法直观、实用，便于分析对比，现今大部分层析成像的结果都应用此方法进行评价。

（5）恢复分辨率实验法。

其基本原理是，应用反演得到的结果作为人工合成模型，在此模型中进行射线追踪，计算走时，同时在这个数据中加入与真实数据误差相同数量级的随机误差，得到人工合成数据集。

反演此数据集，得到的结果再与真实结果进行对比，以便分析图像恢复的情况。

棋盘格实验设置的规则扰动不能模拟复杂模型情况下的图像恢复情况，相对而言，恢复分辨率实验要比棋盘格实验更接近真实情况，但不易分析和对比。

2．地震层析成像技术的未来发展趋势
在以往对中国及邻近区域的地壳上地慢三维结构的研究中，人们比较注重地壳上地慢三维速度结构的研究。

一方面是因为用速度结构可以解释地球内部三维结构某些方面的特征，另一方面是因为计算速度参数相对而言涉及的因素比较少。

但是，从总的方面来看，速度参数只是利用了地震波的运动学信息，而忽略了地震波动力学的特征信息。

地震层析成像研究大多根据地震记录局部上的单一观测值反演单一物理量，方法各自独立，表现出全局性和系统性的不足，在相当大的程度上阻碍了地震层析成像方法的使用。

地震波动是岩石物性的综合响应，岩石物性参数间存在着必然联系，单一参数反演方法对多参数多分量综合响应的分析显得非常无力，因此，多参数多分量同步反演方法的研究非常必要。

在以往的地震层析成像方法中，走时信息和振幅信息分别反演岩石的波速与衰减。

不论走时和振幅都只能代表地震波动的局部特征，而地震波形的整体则是地下地质状况和岩石物性的综合响应。

多分量波动层析成像突破以往只提取局部信息的作法，利用地震波动的整体信息，通过对理论地震波动数据与实测地震波动数据之间波形残差的量化分析，修改弹性波方程中的物理参数，进而修改地质模型的岩石物性参数。

通过这样的研究，地震层析成像研究方法将发生实质性变化，它将使以地震波局部特征为主的研究转变为以地震波场整体动力学特征为主的研究。

地震走时层析成像
内容：
程序解释
（1）正演程序
#include"stdio.h"
#include"math.h"
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("shijimoxing.txt","w");
for(i=0;i<12;i++)
{
for(j=0;j<9;j++)
fprintf(fp0,"%f",xc[i][j]);fprintf(fp0,"\n");
}
/*-------------------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","w");
for(i=0;i<144;i++)
{
fprintf(fp,"%d%f%f%f%f\n",i+1,ray[i].begin,ray[i].end,ray[i].slop e,ray[i].length);
}
/*---------------------dataout(speed)-----------------------------*/ 输出慢度与所给验证
fp1=fopen("慢度.txt","w");
for(i=0;i<12;i++)
{
for(j=0;j<9;j++)fprintf(fp1,"%f ",xc[i][j]);
fprintf(fp1,"\n");
}
/*----------------------------------------------------------------*/计算中间数据输出
fp2=fopen("中间数据.txt","w"); //
for(i=0;i<12;i++)
//
{
//
for(j=0;j<12;j++)
//
{
//
ray[12*i+j].n=j-i;//
fprintf(fp2,"%d ",ray[12*i+j].n); //
}
//
fprintf(fp2,"\n");
//
}
//
/*--------------------initialize----------------------------------*/初始化数组
for(i=0;i<144;i++)
for(j=0;j<12;j++)
for(y=0;y<9;y++)ray[i].transform[j][y]=0;
/*----------------------------------------------------------------*/
for(i=0;i<144;i++)
{
for(j=0;j<100;j++)
{
ray[i].point[j].x=0;
ray[i].point[j].y=0;
}
}
/*----------------------------------------------------------------*/
for(i=0;i<144;i++)
{
ray[i].time=0;
}
fp4=fopen("y.txt","w");
for(i=0;i<12;i++) {
for(j=0;j<12;j++)
{
/*----------------------------------------------------------------*/
for(q=0;q<100;q++)
{
d[q].x=100;
d[q].y=0;
}
/*----------------------------------------------------------------*/计算交点
for(m=0;m<10;m++)
{
d[m].x=m*5.0;
d[m].y=-d[m].x*ray[i*12+j].slope+ray[i*12+j].begin;
}
n=10;
e=ray[i*12+j].n;
if(e>0)
{for(y=1;y<=e;y++)
{
d[n].y=ray[i*12+j].begin+3.0*y-1.5;
d[n].x=(3.0*y-1.5)/fabs(ray[12*i+j].slope);
n++;
}}
if(e<0)
{for(y=-1;y>=e;y--)
{
d[n].y=ray[i*12+j].begin+3.0*y+1.5;
d[n].x=-(3.0*y+1.5)/fabs(ray[12*i+j].slope);
n++;
}}
/*----------------------------------------------------------------*/对交点进行排序
b=0;
while(d[b].x!=100)
{
b++;
}
/*----------------------------------------------------------------*/ for(c=0;c<b;c++)
for(k=0;k<b-1;k++)
if(d[k].x>d[k+1].x)
{
temp1=d[k].x;
d[k].x=d[k+1].x;
d[k+1].x=temp1;
temp2=d[k].y;
d[k].y=d[k+1].y;
d[k+1].y=temp2;
}
/*----------------------------------------------------------------*/输出交点
b=0;
while(d[b].x!=100)
{
ray[12*i+j].point[b].x=d[b].x;
ray[12*i+j].point[b].y=d[b].y;
fprintf(fp4,"%d\t%f,%f
\n",12*i+j+1,ray[12*i+j].point[b].x,ray[12*i+j].point[b].y);
b++;
}
/*----------------------------------------------------------------*/利用两点间距离计算每个射线在每个单元格上的的长度
for(a=1;a<b;a++)
{
mm=sqrt(pow((ray[12*i+j].point[a-1].x-ray[12*i+j].point[a].x),2)+ pow((ray[12*i+j].point[a-1].y-ray[12*i+j].point[a].y),2)) ;
if(ray[12*i+j].n>=0)
ray[12*i+j].transform[(int)(ray[12*i+j].point[a-1].y/3)][(int)(ra y[12*i+j].point[a-1].x/5)]=mm;
if(ray[12*i+j].n<0)
ray[12*i+j].transform[(int)(ray[12*i+j].point[a].y/3)][(int)(ray[ 12*i+j].point[a-1].x/5)]=mm;
}
}
}
/*----------------------------------------------------------------*/计算走时
for(i=0;i<144;i++)
{
for(j=0;j<12;j++)
{
for(y=0;y<9;y++)
{
ray[i].time+=ray[i].transform[j][y]*xc[j][y];
}
}
}
输出最后正验结果
fp3=fopen("矩阵A.txt","w");
for(i=0;i<144;i++)
{
for(j=0;j<12;j++)
{for(y=0;y<9;y++)
fprintf(fp3,"%f ",ray[i].transform[j][y]);
fprintf(fp3,"\n");
}
fprintf(fp3,"\n");
}
fp5=fopen("矩阵A2.txt","w");
for(i=0;i<144;i++)
{
for(j=0;j<12;j++)
{
for(y=0;y<9;y++)
fprintf(fp5,"%f ",ray[i].transform[j][y]);
}
fprintf(fp5,"\n");
}
fp6=fopen("走时.txt","w");
for(i=0;i<144;i++)
{
fprintf(fp6,"%f\n",ray[i].time);
}
}
（ 2 ) 反演程序
#include"stdio.h"
#include"math.h"
struct data
{
float A[108];
float b;
float r;
}ray[144];
结构体 ray【i】表示
b 每条射线的走时
A 每条射线在每条单元格内的长度
r 反演每条射线计算出的误差
main()
{int m[108];
float x[108],R,a1,a2; x[i]反演的每个单元格的慢度
int n=0,j,i,j1;
FILE *fp1,*fp2,*fp3;/*,*fp3,*fp4*/
/*----------------------------------------------------------------*/读取正演的结果
fp1=fopen("矩阵A2.txt","r");
fp2=fopen("走时.txt","r");
for(i=0;i<144;i++)
for(j=0;j<108;j++)
fscanf(fp1,"%f",&ray[i].A[j]);
for(i=0;i<144;i++)
fscanf(fp2,"%f",&ray[i].b);
for(i=0;i<108;i++)x[i]=0; 给第一次的慢度全部付零值
for(i=0;i<144;i++)
for(j=0;j<108;j++)
printf("%f ",ray[i].A[j]);
printf("\n");
for(j1=0;j1<108;j1++) 算出系数矩阵中每一列不为零值的数的个数
{m[j1]=0;
x[j1]=0;
for(i=0;i<144;i++) if(ray[i].A[j1]!=0)m[j1]+=1;
printf("%d ",m[j1]);
n=20; 迭代次数20次
while(n) {
误差的计算
for(i=0;i<144;i++)
{ R=0;
for(j=0;j<108;j++) R+=ray[i].A[j]*x[j];
ray[i].r=ray[i].b-R;
}
for(j1=0;j1<108;j1++)
{
{ a2=0;
/*----------------------------------------------------------------*/ 利用代数重建算法迭代慢度矩阵
for(i=0;i<144;i++)
{
a1=0.0;
for(j=0;j<108;j++) a1+=ray[i].A[j]*ray[i].A[j];
a2+=ray[i].r*ray[i].A[j1]/a1;
}
x[j1]+=a2/m[j1];
}
//printf("%d %d \n",j1,m[j1]); 输出中间数据
}n--; }
// for(j=0;j<108;j++)
//{ if(x[j]>2)x[j]=10;
// if(x[j]<0)x[j]=0;
// }
输出最后数据
fp3=fopen("out.txt","w");
n=0;
for(i=0;i<12;i++){for(j=0;j<9;j++)
{fprintf(fp3,"%f ",x[n]);n++;}
fprintf(fp3,"\n");
}
}。