各向异性VTI介质快速扫描算法FSM走时计算及其改进算法

已于 2024-06-23 22:04:03 修改
阅读量190 收藏
点赞数 2
分类专栏： 地球物理 计算物理 计算数学 文章标签： 算法 c++ c语言
于 2023-10-18 20:38:16 首次发布
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/weixin_42262783/article/details/133914237
版权

(PDF) A shortest‐path‐aided fast‐sweeping method to improve the accuracy of traveltime calculation in vertically transverse isotropic media VTI介质最短路径辅助的快速扫描法提高走时计算精度 (researchgate.net)

代码：

#include<stdio.h>
#include<iostream>
#include<fstream>
#include<math.h>
#include<stdlib.h>
#include<malloc.h>
#include<string.h>
using namespace std;
float pi=3.1415926;
int nx =100;
int nz =100;
float dx= 10;
float dz =10;
float T_control(float a)
{
	if(a>0) return  a;
	else return 0;

}
float min(float a,float b)
{
	if (a>b) return b;
	else return a;
}
float max(float a,float b)
{
	if (a>b) return a;
	else return b;
}



void main()
{
	int i,j;
	int x0 = nx/2; int z0 = nz/2;
	float **s=new float*[nx],**T_old=new float*[nx],**T_old0=new float*[nx],**T_new=new float*[nx],**T=new float*[nx];
	float **Eps=new float*[nx],**Del=new float*[nx],**Eta=new float*[nx],**Vn=new float*[nx],**V0=new float*[nx];
	for(int i=0;i<nx;i++)
	{
		s[i]=new float[nz];T_old[i]=new float[nz];T_old0[i]=new float[nz];T_new[i]=new float[nz];T[i]=new float[nz];
		Eps[i]=new float[nz];Del[i]=new float[nz];Eta[i]=new float[nz];Vn[i]=new float[nz];V0[i]=new float[nz];
		for(int j=0;j<nz;j++)
		{	s[i][j]=0;T_old[i][j]=0;T_old0[i][j]=0;T_new[i][j]=0;T[i][j]=0;Eps[i][j]=0;Del[i][j]=0;Eta[i][j]=0;Vn[i][j]=0;}
	}
	/*ifstream Vnelin;
	  Vnelin.open("mar_140_501_h15");
	  if(!Vnelin)
	  {
	  cerr<< "open file failed" << endl;
	  exit(1);
	  }
	  else 
	  {
	  for(i=0;i<nx;i++)
	  for(j=0;j<nz;j++)
	  Vnelin>>Vn[i][j];
	  }
	  Vnelin.close();
	  */

	FILE *fp;FILE *filename;
	/*	FILE *fp;FILE *filename;
		filename=fopen("mar_140_501_h15","rb");
		for(i=0;i<nx;i++)
		for(j=0;j<nz;j++)
		fread(&Vn[i][j], sizeof(float), 1, filename);
		fclose(filename);
		*/	
	for(i=0;i<nx;i++)   
		for(j=0;j<nz;j++)   
		{   
			Eps[i][j]=0.0;
			Del[i][j]=0.0;
			Eta[i][j]=(Eps[i][j]-Del[i][j])/(1+2*Del[i][j]);
			V0[i][j]=1000;
			Vn[i][j]=1000*sqrt(1.0+2*Del[i][j]);
			s[i][j]=(1.0/Vn[i][j]);
			T_old[i][j]=10.0;
			T[i][j]=22.0;
		}  
	T_old[x0][z0]=0.0;

	float T_xmin=33.0; float T_zmin=33.0;float A=0.0;float B=0.0;float C=0.0,D=0.0;
	float T_xmin0=33.0; float T_zmin0=33.0;
	int loop=0,Max_loop=50;
	int sx,sz;
	int xa,xb,za,zb;

	for(loop=0;loop<Max_loop;loop++)
	{
		for(i=0;i<nx;i++)   
			for(j=0;j<nz;j++)   
			{   
				T_old[i][j]=1e5;
			}  
		T_old[x0][z0]=0.0;


		for(i=0;i<nx;i++)
			for(j=0;j<nz;j++)
			{
				xa=i+1;xb=i-1;za=j+1;zb=j-1;
				if(xa==nx) xa=i;
				if(xb==-1) xb=i;
				if(za==nz) za=j;
				if(zb==-1) zb=j;
				T_xmin=min(T_old[xa][j],T_old[xb][j]);
				T_zmin=min(T_old[i][za],T_old[i][zb]);
				T_xmin0=min(T_old0[xa][j],T_old0[xb][j]);
				T_zmin0=min(T_old0[i][za],T_old0[i][zb]);

				A=Vn[i][j]*Vn[i][j]*(1.0+2*Eta[i][j]);
				B=Vn[i][j]*Vn[i][j]/(1.0+2*Del[i][j]);
				C=2*Eta[i][j]*Vn[i][j]*Vn[i][j];
				if (loop==0) D=1.0;
				else
					D=1+B*C*(T_old0[i][j]-T_xmin0)*(T_old0[i][j]-T_xmin0)*(T_old0[i][j]-T_zmin0)*(T_old0[i][j]-T_zmin0)/dx/dx/dz/dz;
				//cout<<"D="<<D<<endl;
				//D=1+B*C*(T_old0[i][j]-T_old[xa][j])*(T_old0[i][j]-T_old[xa][j])*(T_old0[i][j]-T_old0[i][za])*(T_old0[i][j]-T_old0[i][za])/dx/dx/dz/dz;
				//if(T_xmin-T_zmin>=sqrt((D*dx*dx)/B))   T_old[i][j]=min(T_old[i][j],T_zmin+sqrt((D*dx*dx)/B));
				//else if(T_zmin-T_xmin>=sqrt((D*dx*dx)/A))   T_old[i][j]=min(T_old[i][j],T_xmin+sqrt((D*dx*dx)/A));
				T_old[i][j]=min(T_old[i][j],T_zmin+sqrt((D*dx*dx)/B));
				T_old[i][j]=min(T_old[i][j],T_xmin+sqrt((D*dx*dx)/A));
				if(T_old[i][j]>T_xmin&&T_old[i][j]>T_zmin)
				{
					T_old[i][j]=min(T_old[i][j],(A*T_xmin+B*T_zmin+sqrt((A*T_xmin+B*T_zmin)*(A*T_xmin+B*T_zmin)-(A+B)*(A*T_xmin*T_xmin+B*T_zmin*T_zmin-D*dx*dx)))/(A+B));
				}
				float	fi=pi/4;
				float	Q=1.0+2*(Eps[xb][zb]-Del[xb][zb])/(1.0+2*Del[xb][zb]);
				float	Ag=1.0/(V0[xb][zb]*V0[xb][zb])/(1.0+2*Eps[xb][zb]);
				float	Cg=1.0/(V0[xb][zb]*V0[xb][zb]);
				float	Eg=Ag*pow(sin(fi),2)+Cg*pow(cos(fi),2);
				float	Dg=sqrt(Eg*Eg+4*(Q*Q-1.0)*Ag*Cg*pow(sin(fi),2)*pow(cos(fi),2));
				float	Sg=sqrt((1.0+2*Q)/2.0/(1.0+Q)*Eg+1.0/2.0/(1.0+Q)*Dg);
				T_old[i][j]=min(T_old[xb][zb]+Sg*sqrt(dx*dx+dz*dz),T_old[i][j]);
			}


		for(i=nx-1;i>=0;i--)
			for(j=0;j<nz;j++)
			{
				xa=i+1;xb=i-1;za=j+1;zb=j-1;
				if(xa==nx) xa=i;
				if(xb==-1) xb=i;
				if(za==nz) za=j;
				if(zb==-1) zb=j;
				T_xmin=min(T_old[xa][j],T_old[xb][j]);
				T_zmin=min(T_old[i][za],T_old[i][zb]);
				T_xmin0=min(T_old0[xa][j],T_old0[xb][j]);
				T_zmin0=min(T_old0[i][za],T_old0[i][zb]);

				A=Vn[i][j]*Vn[i][j]*(1.0+2*Eta[i][j]);
				B=Vn[i][j]*Vn[i][j]/(1.0+2*Del[i][j]);
				C=2*Eta[i][j]*Vn[i][j]*Vn[i][j];
				if (loop==0) D=1.0;
				else
					D=1+B*C*(T_old0[i][j]-T_xmin0)*(T_old0[i][j]-T_xmin0)*(T_old0[i][j]-T_zmin0)*(T_old0[i][j]-T_zmin0)/dx/dx/dz/dz;
				//D=1+B*C*(T_old0[i][j]-T_old[xa][j])*(T_old0[i][j]-T_old[xa][j])*(T_old0[i][j]-T_old0[i][za])*(T_old0[i][j]-T_old0[i][za])/dx/dx/dz/dz;
				//if(T_xmin-T_zmin>=sqrt((D*dx*dx)/B))   T_old[i][j]=min(T_old[i][j],T_zmin+sqrt((D*dx*dx)/B));
				//else if(T_zmin-T_xmin>=sqrt((D*dx*dx)/A))   T_old[i][j]=min(T_old[i][j],T_xmin+sqrt((D*dx*dx)/A));
				T_old[i][j]=min(T_old[i][j],T_zmin+sqrt((D*dx*dx)/B));
				T_old[i][j]=min(T_old[i][j],T_xmin+sqrt((D*dx*dx)/A));
				if(T_old[i][j]>T_xmin&&T_old[i][j]>T_zmin)
				{
					T_old[i][j]=min(T_old[i][j],(A*T_xmin+B*T_zmin+sqrt((A*T_xmin+B*T_zmin)*(A*T_xmin+B*T_zmin)-(A+B)*(A*T_xmin*T_xmin+B*T_zmin*T_zmin-D*dx*dx)))/(A+B));
				}
				float	fi=pi/4;
				float	Q=1.0+2*(Eps[xa][zb]-Del[xa][zb])/(1.0+2*Del[xa][zb]);
				float	Ag=1.0/(V0[xa][zb]*V0[xa][zb])/(1.0+2*Eps[xa][zb]);
				float	Cg=1.0/(V0[xa][zb]*V0[xa][zb]);
				float	Eg=Ag*pow(sin(fi),2)+Cg*pow(cos(fi),2);
				float	Dg=sqrt(Eg*Eg+4*(Q*Q-1.0)*Ag*Cg*pow(sin(fi),2)*pow(cos(fi),2));
				float	Sg=sqrt((1.0+2*Q)/2.0/(1.0+Q)*Eg+1.0/2.0/(1.0+Q)*Dg);
				T_old[i][j]=min(T_old[xa][zb]+Sg*sqrt(dx*dx+dz*dz),T_old[i][j]);
			}

		for(i=nx-1;i>=0;i--)
			for(j=nz-1;j>=0;j--)
			{
				xa=i+1;xb=i-1;za=j+1;zb=j-1;
				if(xa==nx) xa=i;
				if(xb==-1) xb=i;
				if(za==nz) za=j;
				if(zb==-1) zb=j;
				T_xmin=min(T_old[xa][j],T_old[xb][j]);
				T_zmin=min(T_old[i][za],T_old[i][zb]);
				T_xmin0=min(T_old0[xa][j],T_old0[xb][j]);
				T_zmin0=min(T_old0[i][za],T_old0[i][zb]);

				A=Vn[i][j]*Vn[i][j]*(1.0+2*Eta[i][j]);
				B=Vn[i][j]*Vn[i][j]/(1.0+2*Del[i][j]);
				C=2*Eta[i][j]*Vn[i][j]*Vn[i][j];
				if (loop==0) D=1.0;
				else
					D=1+B*C*(T_old0[i][j]-T_xmin0)*(T_old0[i][j]-T_xmin0)*(T_old0[i][j]-T_zmin0)*(T_old0[i][j]-T_zmin0)/dx/dx/dz/dz;
				//D=1+B*C*(T_old0[i][j]-T_old[xa][j])*(T_old0[i][j]-T_old[xa][j])*(T_old0[i][j]-T_old0[i][za])*(T_old0[i][j]-T_old0[i][za])/dx/dx/dz/dz;
				//if(T_xmin-T_zmin>=sqrt((D*dx*dx)/B))   T_old[i][j]=min(T_old[i][j],T_zmin+sqrt((D*dx*dx)/B));
				//else if(T_zmin-T_xmin>=sqrt((D*dx*dx)/A))   T_old[i][j]=min(T_old[i][j],T_xmin+sqrt((D*dx*dx)/A));
				T_old[i][j]=min(T_old[i][j],T_zmin+sqrt((D*dx*dx)/B));
				T_old[i][j]=min(T_old[i][j],T_xmin+sqrt((D*dx*dx)/A));
				if(T_old[i][j]>T_xmin&&T_old[i][j]>T_zmin)
				{
					T_old[i][j]=min(T_old[i][j],(A*T_xmin+B*T_zmin+sqrt((A*T_xmin+B*T_zmin)*(A*T_xmin+B*T_zmin)-(A+B)*(A*T_xmin*T_xmin+B*T_zmin*T_zmin-D*dx*dx)))/(A+B));
				}
				float	fi=pi/4;
				float	Q=1.0+2*(Eps[xa][za]-Del[xa][za])/(1.0+2*Del[xa][za]);
				float	Ag=1.0/(V0[xa][za]*V0[xa][za])/(1.0+2*Eps[xa][za]);
				float	Cg=1.0/(V0[xa][za]*V0[xa][za]);
				float	Eg=Ag*pow(sin(fi),2)+Cg*pow(cos(fi),2);
				float	Dg=sqrt(Eg*Eg+4*(Q*Q-1.0)*Ag*Cg*pow(sin(fi),2)*pow(cos(fi),2));
				float	Sg=sqrt((1.0+2*Q)/2.0/(1.0+Q)*Eg+1.0/2.0/(1.0+Q)*Dg);
				T_old[i][j]=min(T_old[xa][za]+Sg*sqrt(dx*dx+dz*dz),T_old[i][j]);
			}

		for(i=0;i<nx;i++)
			for(j=nz-1;j>=0;j--)
			{
				xa=i+1;xb=i-1;za=j+1;zb=j-1;
				if(xa==nx) xa=i;
				if(xb==-1) xb=i;
				if(za==nz) za=j;
				if(zb==-1) zb=j;
				T_xmin=min(T_old[xa][j],T_old[xb][j]);
				T_zmin=min(T_old[i][za],T_old[i][zb]);
				T_xmin0=min(T_old0[xa][j],T_old0[xb][j]);
				T_zmin0=min(T_old0[i][za],T_old0[i][zb]);

				A=Vn[i][j]*Vn[i][j]*(1.0+2*Eta[i][j]);
				B=Vn[i][j]*Vn[i][j]/(1.0+2*Del[i][j]);
				C=2*Eta[i][j]*Vn[i][j]*Vn[i][j];
				if (loop==0) D=1.0;
				else
					D=1+B*C*(T_old0[i][j]-T_xmin0)*(T_old0[i][j]-T_xmin0)*(T_old0[i][j]-T_zmin0)*(T_old0[i][j]-T_zmin0)/dx/dx/dz/dz;
				//D=1+B*C*(T_old0[i][j]-T_old[xa][j])*(T_old0[i][j]-T_old[xa][j])*(T_old0[i][j]-T_old0[i][za])*(T_old0[i][j]-T_old0[i][za])/dx/dx/dz/dz;
				//if(T_xmin-T_zmin>=sqrt((D*dx*dx)/B))   T_old[i][j]=min(T_old[i][j],T_zmin+sqrt((D*dx*dx)/B));
				//else if(T_zmin-T_xmin>=sqrt((D*dx*dx)/A))   T_old[i][j]=min(T_old[i][j],T_xmin+sqrt((D*dx*dx)/A));
				T_old[i][j]=min(T_old[i][j],T_zmin+sqrt((D*dx*dx)/B));
				T_old[i][j]=min(T_old[i][j],T_xmin+sqrt((D*dx*dx)/A));
				if(T_old[i][j]>T_xmin&&T_old[i][j]>T_zmin)
				{
					T_old[i][j]=min(T_old[i][j],(A*T_xmin+B*T_zmin+sqrt((A*T_xmin+B*T_zmin)*(A*T_xmin+B*T_zmin)-(A+B)*(A*T_xmin*T_xmin+B*T_zmin*T_zmin-D*dx*dx)))/(A+B));
				}
				float	fi=pi/4;
				float	Q=1.0+2*(Eps[xb][za]-Del[xb][za])/(1.0+2*Del[xb][za]);
				float	Ag=1.0/(V0[xb][za]*V0[xb][za])/(1.0+2*Eps[xb][za]);
				float	Cg=1.0/(V0[xb][za]*V0[xb][za]);
				float	Eg=Ag*pow(sin(fi),2)+Cg*pow(cos(fi),2);
				float	Dg=sqrt(Eg*Eg+4*(Q*Q-1.0)*Ag*Cg*pow(sin(fi),2)*pow(cos(fi),2));
				float	Sg=sqrt((1.0+2*Q)/2.0/(1.0+Q)*Eg+1.0/2.0/(1.0+Q)*Dg);
				T_old[i][j]=min(T_old[xb][za]+Sg*sqrt(dx*dx+dz*dz),T_old[i][j]);
			}

		for(i=0;i<nx;i++)
			for(j=0;j<nz;j++)
			{
				T_old0[i][j]=T_old[i][j];
			}

		if(loop==0)
		{
			for(i=0;i<nx;i++)
				for(j=0;j<nz;j++)
				{
					T_new[i][j]=T_old[i][j];
				}
		}
	}


	fp=fopen("traveltime.dat","wb+");
	for(i=0;i<nx;i++)
		for(j=0;j<nz;j++)
		{
			//if(fabs(T_old[i][j]-0.6+0.033)<0.01) T_old[i][j]=1;
			//else T_old[i][j]=0;
			T_new[i][j]=dx*sqrt(pow((i-x0),2)+pow((j-z0),2))*s[i][j];
			float dt=T_old[i][j]-T_new[i][j];
			//fwrite(&dt,sizeof(float),1,fp);
			fwrite(&T_old[i][j],sizeof(float),1,fp);
		}
	fclose(fp);

}

张建明0704
关注
专栏目录
二维VTI介质弹性波高阶有限差分正演模拟+pml边界条件
11-17
基于matlab编写的，VTI介质弹性波方程的高阶交错网格有限差分正演模拟+pml吸收边界条件。带有注释说明，添加了pml边界代码，很适合基础新手参考学习，波场模拟的结果以动画的方式展示。同时也可以进行波场快照的输出...
快速扫描法（Fast Sweeping Method，FSM）及最短路径辅助的快速扫描法SPAFSM走时计算
weixin_42262783的博客
10-18 650
代码包含快速扫描法（FSM）走时计算的基本方法（Zhao 2005）以及高阶FSM方法和对FSM方法的一点改进：添加最短路径辅助解提高计算精度（SPAFSM，Zhang 2023），仅供初学者参考，如有疑问，欢迎随时交流讨论。//i是次数，j是具体下标。
各向异性介质多参数化走时反演敏感核（辐射模式）分析
weixin_42262783的博客
04-02 390
在阅读文章的基础上，很多朋友对此研究很感兴趣，想进一步深入学习，在此，提供了敏感核分析的代码，下面是Thomsen参数化的代码，其他参数化的代码一并打包上传。（好吧，发现不能打包上传，谁有办法，评论区告诉我！（需要代码的私信我就好））说明：做该项研究当时用好几种办法推导过 结果都很相似 验证了我的推导是对的。
各向异性VTI介质地震波场模拟，速度-应力弹性波方程交错网格有限差分法，PML吸收边界（C++和MATLAB两种语言的代码都有）
weixin_42262783的博客
10-18 454
参考各向同性的程序，进一步修改为VTI介质，不懂的随时交流。
快速行进法（FMM）求解程函方程的基本原理，以及VTI介质和TTI介质中的走时计算
weixin_42262783的博客
07-25 141
【代码】快速行进法（FMM）求解程函方程的基本原理。
程函方程伴随状态法走时层析，初至及反射，VTI介质多参数
weixin_42262783的博客
10-18 337
实际资料例子：
各向异性TTI介质地震波场模拟，速度-应力弹性波方程交错网格高阶有限差分法，PML吸收边界（C++和MATLAB两种语言的代码都有）
weixin_42262783的博客
10-18 870
逐渐改进，从各向同性到VTI各向异性，再到TTI介质。模拟效果可以在走时计算的文章里看到，有对比。不懂的欢迎随时交流。
matlab绘制vti群速度,VTI介质中P波相速度与群速度的比较
weixin_28958485的博客
03-16 514
己l健绪,l月近年来,随着地展资料质量及其处理和解释质量的不断提高,地展各向异性的研究在勘探和开发中取得了长足的进步(例如:ThomsonL.,1986;Tsvankin1.,1995;AikhalifaTandTsvankin1.,1995;GreehkaV.andTsvankin11998;GreehkaV.andTsvankin1.,2000;LeStunffY.,GreehkaV.,and...
根据克里斯托弗方程（Christoffel equation）计算相速度和偏振(极化)方向，进而计算群速度
weixin_42262783的博客
03-12 988
弹性波方程是通过应力（stress）应变（strain）关系（广义胡克定律，或本构关系）+运动微分方程（动量守恒线性方程 Linearized equation of momentum conservation Navier 方程或牛顿第二定律）导出，描述了波的所有动力学和运动学传播特性。而Christoffel 方程是各向异性弹性波动方程的平面波（单相波，对应相速度）解形式，可以用来独立描述三种类型的波的极化和相速度特征，即运动学（速度）和偏振（振幅以及方向）。3.几大方程的物理意义是什么？
基于MPI并行的VTI介质逆时偏移成像与ADCIGs提取
RToax
02-01 1270
简单明了“炮并行”、“MPI”、“C语言”、“VTI介质”、“RTM”、“中间激发两边接收”、“全孔径接收”、“照明”、“拉普拉斯滤波”、“角度域共成像点道集”、“Poynting矢量”； 在此做个备份！ 直接上代码吧！头文件请到SeismicUnix中找。 //#############################################################a
行业资料-交通装置-VTI介质中地震波的走时计算方法.zip
08-17
行业资料-交通装置-VTI介质中地震波的走时计算方法.zip
算法刷题：300. 最长递增子序列、674. 最长连续递增序列、718. 最长重复子数组、1143. 最长公共子序列
qq_61936886的博客
09-12 512
如果text1[i - 1] 与 text2[j - 1]不相同，那就看看text1[0, i - 2]与text2[0, j - 1]的最长公共子序列 和 text1[0, i - 1]与text2[0, j - 2]的最长公共子序列，取最大的。但dp[i][0] 和dp[0][j]要初始值，因为 为了方便递归公式dp[i][j] = dp[i - 1][j - 1] + 1;2.递推公式：当A[i - 1] 和B[j - 1]相等的时候，dp[i][j] = dp[i - 1][j - 1] + 1;
算法知识点——常用输入输出数据的方式
5233的博客
09-12 301
除了换行之外，还会清空缓存区，而清空缓存区这件事情会耗时很久，因此为了提高速度，一般直接：cout << “\n”;2、n组输入输出(n不确定）1、n组输入输出(n确定）3、 复杂的输入输出。
输出一个三位数
最新发布
2301_77501988的博客
09-13 403
本题要求编写程序，输入一个三位数，输出其个位数字、十位数字和百位数字。
Java数据结构（十）——冒泡排序、快速排序
xiaokuer_的博客
09-11 1194
文章介绍了冒泡排序和快速排序，快速排序的内容包括：划分方式、优化策略以及非递归实现
井地联合提取火成岩VTI介质各向异性参数新方法
在新疆塔里木盆地的宽方位三维地震数据中，由于这种介质的各向异性和垂向速度变化，导致非双曲线时差现象，使得传统的地震资料处理方法如DMO反演法、移动震源VSP反演法以及多偏移距VSP初至旅行时反演法无法准确获取...
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值