层析反演实习内容

1.SVD 方法求解方程组 （来自：geophysical inversion_book.pdf ）
⎪⎪⎭⎫ ⎝⎛=⎪⎪⎪⎭
⎫ ⎝⎛⎥⎦⎤⎢⎣⎡11321100011m m m SVD 求解结果：⎥⎥⎥⎥⎥⎥⎥⎥⎦
⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=10021021 A
2.最速下降法
1.滤波反投影重建方法
2.Fourier 重建方法
3.代数重建方法
4.共轭梯度方法（来自：geophysical inversion_book.pdf---p164、数值分析与实验(薛毅-2006，p102)）
5.投影方法-跨孔地震层析成像研究_硕士论文
6. SVD 方法
7.射线追踪
1. Matlab 初识与概述
2.Radon 变换
计算下列函数的radon 变换结果：
⎩⎨⎧≤≤=其他 01,0 1),(y x y x f
3.Radon 逆变换
对1进行Radon 逆变换
4.射线追踪（跨孔地震层析成像研究_硕士论文）
模型1为一三层模型，范围100m ×100m 。

该模型含有三个水平层，上下两个层速度为3000m/s ，中间为一高速层，速度为4500m ／s ，高速层厚度范围是50～60m 。

激发点有19个，坐标依次为是E l (0,5)、E 2(0,10)、E 3(0,15)，…，E 10(0,95)接收点有21个，坐标分别为R l (100，0)、R 2 (100,5)、R 3(100,10)，…，R 21(100,100)。

将此模型网格化为10×10个单元格．根据LTI 法先求出10×21条射线的旅行时。

5.滤波反投影重建方法
第一步：利用matlab 里的phantom 变换生成头模型；
phantom 功能:
产生一个头部幻影图像.
语法:
P = phantom(def,n)
P = phantom(E,n)
[P,E] = phantom(...)
举例
P = phantom('Modified Shepp-Logan',200);
imshow(P)
相关命令:
radon, iradon
第二步：利用randon 函数生成0o 到180o 19个方向上的投影数据；
第三步：编写radon 逆变换程序。

6.Fourier 重建方法
7.代数重建方法
实验题：CT 图像的代数重建问题（线性方程组的应用）
X 射线透视可以得到3维对象在2维平面上的投影，CT 则通过不同角度的X 射线得到3维对象的多个2维投影，并以此重建对象内部的3维图像。

代数重建方法就是从这些2维投影出发，通过求解超定线性方程组，获得对象内部3维图像的方法。

这里我们考虑一个更简单的模型，从2维图像的1维投影重建原先的2维图像。

一个长方形图像可以用一个横竖均匀划分的离散网格来覆盖，每个网格对应一个像素，它是该网格上各点像素的均值。

这样一个图像就可以用一个矩阵表示，其元素就是图像在一点的灰度值（黑白图像）。

下面我们以33 图像为例来说明。

每个网格中的数字代表其灰度值，范围在]1,0[内，记0表示白，1表示黑，0.5为中间的灰色，沿某个方向的投影就将该方向上的灰度值相加（见上图所示）。

如果我们不知道网格中的数值，只知道沿竖直方向和水平方向的投影。

设网格按第1列、第2列、第3列的顺序排列，为了确定网格中的灰度值，可以建立线性方程组：
1 1
1.5 1 1 0
⎪⎪⎪⎪⎪⎪⎪⎪⎭
⎫ ⎝⎛=⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛5.1115.15.05.1100100100010010010001001001111000000000111000
000000111X 显然该方程组的解是不唯一的，为了重建图像，必须增加投影数量。

如我们增加从右上到左下的投影，则方程组将增加5个方程，成为超定方程组。

考虑到测量误差，可以将超定方程组的近似解作为重建的图像数据。

问题：给定一个33⨯图像的2个投影，沿左上到右下，投影数据依次为0.8，1.2，1.7，0.2，0.3；从右上到左下的投影数据为0.6，0.2，1.6，1.2，0.6。

求：建立可以确定网格数据的线性方程组，并用MATLAB 求解；
8.共轭梯度方法
用共轭梯度发求解下列方程组：
⎪⎩
⎪⎨⎧-=-+=-+=++53.048.104.368.203.130.193.048.105.053.448.151.2321321321x x x x x x x x x
9.投影方法 （LSQR 方法）（跨孔地震层析成像研究_硕士论文）
模型1为一三层模型，范围100m ×100m 。

该模型含有三个水平层，上下两个层速度为3000m/s ，中间为一高速层，速度为4500m ／s ，高速层厚度范围是50～60m 。

激发点有19个，坐标依次为是E l (0,5)、E 2(0,10)、E 3(0,15)，…，E 10(0,95)接收点有21个，坐标分别为R l (100，0)、R 2 (100,5)、R 3(100,10)，…，R 21(100,100)。

将此模型网格化为10×10个单元格．根据LTI 法先求出10×21条射线的旅行时．构成求解方程min||Ax-b||2时b 时间向量，初始速度向量x (0)=3000m/s 。

通过反复调用射线追踪程序与反演算法给定选代误差5m/s ，迭代13次程序停止。

反演速度平面图如图5.4和反演速度三维图5.5。

10. SVD 方法
SVD 方法求解方程组 （来自：geophysical inversion_book.pdf ）
⎪⎪⎭⎫ ⎝⎛=⎪⎪⎪⎭⎫ ⎝⎛⎥⎦⎤⎢⎣⎡113
21100011m m m SVD 求解结果：⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢
⎢⎢⎣⎡=10021
02
1
A。