1. 雷克子波
中心频率f=25,控制频带宽度的参数γ=3]
2. 设置一个均匀介质模型,整个区域的速度值可设为常数,即只有一种介质,将震源点放在模型中间,分别记录两个时刻的波前快照(即该时刻区域内所有网格点的波场值)。第一时刻为地震波还未传播到边界上的某时刻,第二时刻为地震波已经传播到边界上的某时刻,体会其人工边界反射
0.4s Tn=800
3. 给定水平层状介质模型,震源位于浅层中央处,绘制波前快照和地震记录,观测直达波、透射波、反射波
震源3,100
上层 2000m/s 0-80
中层 2500m/s 80-130
下层 3000m/s 130-200
波前快照
Tn=100
Tn=500
Tn=900
Tn=1100
Tn=1200
Tn=1300
地震记录
4. 自主设计模型和观测方式,至少观测如下一种波:多次波、绕射波、折射波,波前和记录都要有。
模型
分界面方程
i>(10*(j-100)^2+50)
震源3,60
低速2000
高速3000
Tn=500
Tn=700
Tn=800
Tn=1200
地震记录
5. 拓展:
①模拟倾斜界面、起伏界面、断层、盆地等特殊地质构造中地震波传播
②研究实现吸收边界处理
③地震波传播过程中的几何扩散效应及波形变化规律(后面节内容)
分界面方程
i>(50sin(jpi/50)+100)
低速2000
高速3000
震源3,100
波前记录
Tn=200
Tn=500
Tn=700
Tn=900
地震记录
6. 源码
6.1 C++
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
int Tn;
int Xn;
int Zn;
#define PI 3.141592653589793
#define X 3
#define Y 60
#define DH 8.0
#define DT 0.0005
#define Frequency 25
#define Round 3
int main()
{
printf("Hello this is a program to caculate acoustic wave equation\n");
printf("Please input the Xn, Zn and Tn\n");
printf("Xn is :");
if (scanf_s("%d", &Xn) != 1 || Xn <= 0)
{
printf("What youinput is nonesense! \n");
exit(1);
}
printf("Zn is :");
if (scanf_s("%d", &Zn) != 1 || Xn <= 0)
{
printf("what youinput is nonesense !\n");
exit(1);
}
printf("Tn is : ");
if (scanf_s("%d", &Tn) != 1 || Xn <= 0)
{
printf("what youinput is nonesense!\n");
exit(1);
}
double* w = (double*)malloc(sizeof(double) * Tn);
for (int i = 0; i < Tn; i++)
{
w[i] = exp(-(2 * PI * Frequency / Round) * (2 * PI * Frequency / Round) * DT * (i - 90) * (i - 90) * DT) * cos(2 * PI * Frequency * DT * (i - 90));
}
double** vel = (double**)malloc(sizeof(double*) * Xn);
for (int i = 0; i < Xn; i++)
{
vel[i] = (double*)malloc(sizeof(double) * Zn);
}
//绕射
/*
for (int i = 0; i < Xn; i++)
for (int j = 0; j < Zn; j++)
{
vel[i][j] = 3000;//均匀介质
}
*/
/*
for (int i = 0; i < Xn; i++)
for (int j = 0; j < Zn; j++)
if (i > (10 * pow((j - 100) , 2) + 50))
{
vel[i][j] = 3000;
}
else
{
vel[i][j] = 2000;
}
*/
//起伏界面
/*
for (int i = 0; i < Xn; i++)
for (int j = 0; j < Zn; j++)
if (i > (50* sin(j*PI/50)) + 100))
{
vel[i][j] = 3000;
}
else
{
vel[i][j] = 2000;
}
*/
//盆地
/*
for (int i = 0; i < Xn; i++)
{
for (int j = 0; j < Zn; j++)
{
if (j < (i-400)||(j>900-i)||(i>450))
{
vel[i][j] = 2000;
}
else
{
vel[i][j] = 3000;
}
}
}*/
//malloc the size for the grid
double** p_now = (double**)malloc(sizeof(double*) * Xn);
for (int i = 0; i < Xn; i++)
{
p_now[i] = (double*)malloc(sizeof(double) * Zn);
}
double** p_past = (double**)malloc(sizeof(double*) * Xn);
for (int i = 0; i < Xn; i++)
{
p_past[i] = (double*)malloc(sizeof(double) * Zn);
}
double** temp = (double**)malloc(sizeof(double*) * Xn);
for (int i = 0; i < Xn; i++)
{
temp[i] = (double*)malloc(sizeof(double) * Zn);
}
//波前
double** wavefront = (double**)malloc(sizeof(double*) * Tn);
for (int i = 0; i < Tn; i++)
{
wavefront[i] = (double*)malloc(sizeof(double) * Zn);
}
double** p_temp;
for (int i = 0; i < Xn; i++)
{
for (int j = 0; j < Zn; j++)
{
p_now[i][j] = 0;
p_past[i][j] = 0;
temp[i][j] = 0;
}
}
for (int i = 0; i < Tn; i++)
{
for (int j = 0; j < Zn; j++)
{
wavefront[i][j]=0;
}
}
for (int k = 0; k < Tn; k++)
{
for (int i = 2; i < Xn - 2; i++)
{
for (int j = 2; j < Zn - 2; j++)
{
//吸收边界
/*int A = vel[i][j] * DH/DT ;
p_now[2][j] = (2 - 2 * A - pow(A, 2)) * temp[2][j] + 2 * A * (1 + A) * temp[3][j] - pow(A, 2) * temp[4][j] + (2 * A - 1) * p_past[2][j] - 2 * A * p_past[3][j];//左边界
p_now[Xn - 3][j] = (2 - 2 * A - pow(A, 2)) * temp[Xn - 3][j] + 2 * A * (1 + A) * temp[Xn - 4][j] - pow(A, 2) * temp[Xn - 5][j] + (2 * A - 1) * p_past[Xn - 3][j] - 2 * A * p_past[Xn - 4][j];//右边界
p_now[i][2] = (2 - 2 * A - pow(A, 2)) * temp[i][2] + 2 * A * (1 + A) * temp[i][3] - pow(A, 2) * temp[i][4] + (2 * A - 1) * p_past[i][2] - 2 * A * p_past[i][3];//上边界
p_now[i][Zn - 3] = (2 - 2 * A - pow(A, 2)) * temp[i][Zn - 3] + 2 * A * (1 + A) * temp[i][Zn - 4] - pow(A, 2) * temp[i][Zn - 5] + (2 * A - 1) * p_past[i][Zn - 3] - 2 * A * p_past[i][Zn - 4];//下边界*/
p_now[i][j] = 2 * temp[i][j] - p_past[i][j] + (vel[i][j] * vel[i][j] * DT * DT / (DH * DH)) * (-1.0 / 12 * (temp[i - 2][j]
+ temp[i + 2][j]) + 4 / 3.0 * (temp[i - 1][j] + temp[i + 1][j]) - 5 / 2.0 * temp[i][j] - 1.0 / 12 * (temp[i][j - 2] + temp[i][j + 2])
+ 4 / 3.0 * (temp[i][j - 1] + temp[i][j + 1]) - 5 / 2.0 * temp[i][j]) + w[k] * (i == X && j == Y) ;
/*wavefront[k] [j]= p_now[2][j];//地震记录*/
}
}
p_temp = p_now;
p_now = p_past;
p_past = temp;
temp = p_temp;
}
printf("Xn is %d,Zn is %d,Tn is %d\n", Xn, Zn, Tn);
//将矩阵数据写入文件
FILE* fp;
errno_t error;
error = fopen_s(&fp, "wavefront.dat", "w+");
//波前记录
if (error != 0)
{
printf("打开失败");
}
else
{
for (int i = 0; i < Xn; i++)
{
for (int j = 0; j < Zn; j++)
{
fprintf(fp, "%lf ", temp[i][j]);
}
fprintf(fp, "\n");
}
}
//地震记录
/*
if (error != 0)
{
printf("打开失败");
}
else
{
for (int i = 0; i < Tn; i++)
{
for (int j = 0; j < Zn; j++)
{
fprintf(fp, "%lf ", wavefront[i][j]);
}
fprintf(fp, "\n");
}
}
*/
//释放指针
fclose(fp);
free(w);
for (int i = 0; i < Xn; i++)
{
free(p_now[i]);
}
for (int i = 0; i < Xn; i++)
{
free(p_past[i]);
}
for (int i = 0; i < Xn; i++)
{
free(temp[i]);
}
for (int i = 0; i < Xn; i++)
{
free(vel[i]);
}
free(temp);
free(p_now);
free(p_past);
free(vel);
return 0;
}