褶积合成地震记录

目的：
通过编写和运行相关程序，利用褶积理论实现水平层状模型自激自收合成地震记录。
内容及要求：
1、计算水平层状模型垂直入射反射系数序列r(t)
已知双层模型，界面水平，第一层厚度100m，纵波速度1500m/s，密度2000kg/m3，第二层纵波速度2200m/s，密度2100kg/m3。
根据已知两层介质、水平界面模型，计算界面垂直入射反射系数R。
根据第一层厚度及纵波速度，计算界面反射波双程旅行时t0，计算反射系数序列r(t).
2、褶积合成地震记录
设置子波为50Hz的雷克子波，利用褶积理论合成地震记录。要求：
（1）合成的自激自收地震道，子波为零相位子波（子波最大振幅和反射系数位置对应）。
（2）将反射系数序列r(t)、子波w(t)和地震道s(t)绘制到一副图中，要求纵坐标为时间。

实验原理：
计算反射系数计算公式：，和自激自收时间：，然后计算反射系数序列。子波取雷克子波，其函数为:
线性褶积模型公式为：，式中t为反射波自激自收双程旅行时，w(t)为时间域子波序列，r(t)为时间域反射系数序列。褶积时需要对子波进行相位校正，从而实现了波零相位化，使得地震道中子波振幅最大值对应界面的反射系数。校正方法为褶积运算后地震道“向上移动”半个子波序列长度。
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
代码如下
.
.
.
1.反射系数序列r(t)Python代码

/**
 仅供学习，如有纰漏，望请指正 --by wxw
*/
1.import numpy as np
2.import matplotlib.pyplot as plt
3.import scipy
4.import math
5.
6.
7.ρ1 = 2000
8.ρ2 = 2100
9.h = 100
10.v1 = 1500
11.v2 = 2200
12.R = (ρ2*v2 - ρ1*v1)/(ρ2*v2 + ρ1*v1)
13.t0 = 2*h/v1
14.
15.r = np.linspace(-0.3, 0.3, 200)
16.x1 = 0
17.x2 = R
18.y1 = t0
19.y2 = t0
20.f1 = 50
21.t = np.linspace(-0.055, 0.055, 500)
22.rick1 = (1-2*(math.pi*f1*t)**2)*math.e**(-1*(math.pi*f1*t)**2)
23.
24.
25.plt.figure(figsize=(10,6))
26.fig,ax = plt.subplots()
27.plt.ylim(0.3, 0)
28.plt.xlim(-0.3, 0.3)
29.ax.set_aspect(4)
30.plt.axvline(0, color='b')
31.ax.plot([x1, x2], [y1, y2], linewidth=1)
32.ax.set_xlabel("R")
33.ax.set_ylabel("Time")
34.ax.set_title('r(t)')
35.ax.legend()
36.plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39

2(1).合成的自激自收地震道，子波为零相位子波Python代码

/**
 仅供学习，如有纰漏，望请指正 --by wxw
*/
1.import numpy as np
2.import matplotlib.pyplot as plt
3.import scipy
4.import math
5.
6.f1 = 50
7.
8.t = np.linspace(-0.055, 0.055, 500)
9.rick1 = (1-2*(math.pi*f1*t)**2)*math.e**(-1*(math.pi*f1*t)**2)
10.
11.plt.figure(figsize=(10,6))
12.fig,ax = plt.subplots()
13.plt.ylim(0.02, -0.02)
14.plt.xlim(-0.5, 1)
15.ax.set_aspect(80)
16.ax.plot(rick1, t, '--', label="50Hz")
17.ax.set_xlabel("Time")
18.ax.set_ylabel("Amp")
19.ax.set_title('rick')
20.ax.legend()
21.plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

2(2).反射系数序列r(t)、子波w(t)和地震道s(t)Python代码

/**
 仅供学习，如有纰漏，望请指正 --by wxw
*/
1.import numpy as np
2.import matplotlib.pyplot as plt
3.import scipy
4.import math
5.
6.
7.h = 100
8.ρ1 = 2000
9.ρ2 = 2100
10.v1 = 1500
11.v2 = 2200
12.R = (ρ2*v2 - ρ1*v1)/(ρ2*v2 + ρ1*v1)
13.t0 = 2*h/v1
14.f = 50
15.dt = 0.0001
16.t_max = 0.3
17.ts_max = 0.1
18.
19.# r(t)
20.r = np.linspace(-0.3, 0.3, 200)
21.x1 = 0
22.x2 = R
23.y1 = t0
24.y2 = t0
25.
26.# w(t)
27.t = np.linspace(-0.055, 0.055, 500)
28.rick = (1-2*(math.pi*f*t)**2)*math.e**(-1*(math.pi*f*t)**2)
29.
30.# s(t)
31.r = np.zeros((round(t_max/dt), 1))
32.rd = round(t0/dt)
33.r[rd] = R
34.r = np.array(r).flatten()
35.ts = np.linspace(-0.15, 0.33, 3000)
36.rick1 = (1-2*(math.pi*f*ts)**2)*math.e**(-1*(math.pi*f*ts)**2)
37.nts = round(ts_max/dt)+1
38.trace = np.zeros((round(t_max/dt), 1))
39.trace_conv = np.convolve(rick1, r, mode='full')
40.trace = trace_conv[np.arange(round(ts_max/2/dt)+1, round(ts_max/2/dt)+round(t_max/dt)+1, 1)]
41.
42.
43.plt.figure(figsize=(12, 6))
44.fig, ax = plt.subplots(1, 3)
45.ax1 = ax[0]
46.ax2 = ax[1]
47.ax3 = ax[2]
48.
49.ax1.set_ylim(0.3, 0)
50.ax1.set_xlim(-0.2, 0.3)
51.ax1.set_aspect(4)
52.ax1.axvline(0, color='b')
53.ax1.plot([x1, x2], [y1, y2], linewidth=2, color='g')
54.ax1.set_xlabel("R")
55.ax1.set_ylabel("Time")
56.ax1.set_title('r(t)')
57.ax1.legend()
58.
59.ax2.set_ylim(-0.02, 0.02)
60.ax2.set_xlim(-0.5, 1)
61.ax2.set_aspect(90)
62.ax2.plot(rick, t)
63.ax2.set_xlabel("Amp")
64.ax2.set_ylabel("Time")
65.ax2.set_title('w(t)')
66.ax2.legend()
67.
68.ax3.set_ylim(0.3, 0)
69.ax3.set_xlim(-0.2, 0.3)
70.ax3.set_aspect(4)
71.ax3.plot(trace, ts)
72.ax3.set_xlabel("Amp")
73.ax3.set_ylabel("Time")
74.ax3.set_title('s(t)')
75.ax3.legend()
76.
77.plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80