HVSR 及反演

生活
发布日期:   2023-06-14
阅读次数:  

尽管瑞利波椭圆度数据的主要特征是2.4 Hz处的显著波谷，但在1.5至8 Hz范围内的恒定椭圆值和缺乏其他峰值和/或波谷也约束着模型。 比如用1.5~3Hz,1.5~4Hz的有限频带进行了额外的反演测试。得到的速度剖面类似于使用全带宽的剖面，但总体不确定性增加，特别是在浅部（<20米深度）

InSight HVSR

output

Fig1. InSight Sol S0423d VEL HVSR

output

Fig2. InSight Sol S0423d VEL 1.5Hz~8.0Hz HVSR

output

Fig2. InSight 2.4Hz B HVSR

​

Model

output

Fig3. Possive Data HVSR

Table1

image-20230529220617637

Inversion

image-20230522153005250

Fig4. Possive Data HVSR Inversion

截屏2023-05-22 12.01.56

Fig4. Possive Data HVSR Inversion

image-20230522153701679

层号 层厚 [m] S波速 [m/s]
1 100 200
2 650 600 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
vp = zeros(150, 300);
vs = zeros(150, 300);
vp(1:20, :) = 520; vp(21:150, :) = 640;
vs(1:20, :) = 480; vs(21:150, :) = 600;

dx = 10.;

dt = 1e-3;

dtx=dt/dx;

den=1.8*ones(nz,nx);

dz=dx;

image-20230522154317635

image-20230523173125998

image-20230526171337579

image-20230526171426682

image-20230526171446452

image-20230530105026247

image-20230530105244352

Asd

image-20230531125952026

image-20230601140446748

image-20230601113054419

image-20230601113136995

image-20230601135520329

image-20230601140252303

image-20230601142422473

code

Fig1. InSight Sol S0423d VEL HVSR

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
81
82
83
84
85
86
87
88
89
90
91
import time
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

import hvsrpy
from hvsrpy import utils

file_name = '../InSight-seismic-data-downloader/DATA/2.4_HZ/C/S0423d/S0423d_VEL.mseed'

windowlength = 60

filter_bool = False       
filter_flow = 0.1                   
filter_fhigh = 10                  
filter_order = 5

# Width of cosine taper {0. - 1.}. Geopsy default of 0.05 is equal to 0.1 -> 0.1 is recommended
width = 0.1


# Konno and Ohmachi smoothing constant. 40 is recommended.
bandwidth = 40
# Minimum frequency after resampling
resample_fmin = 0.1
# Maximum frequency after resampling
resample_fmax = 10
# Number of frequencies after resampling 重采样后的频率点数
resample_fnum = 1000
# Type of resampling {'log', 'linear'}
resample_type = 'linear'
# Upper and lower frequency limits to restrict peak selection. To use the entire range use `None`.
peak_f_lower = None
peak_f_upper = None


# Method for combining horizontal components {"squared-average", "geometric-mean", "single-azimuth"}.
# Geopsy's default is "squared-average" -> "geometric-mean" is recommended.
method = "geometric-mean"
# If method="single-azimuth", set azimuth in degree clock-wise from north. If method!="single-azimuth", value is ignored.
azimuth = 0

# Boolean to control whether frequency domain rejection proposed by Cox et al. (2020) is applied.
# Geopsy does not offer this functionality.
rejection_bool = True
# Number of standard deviations to consider during rejection. Smaller values will reject more windows -> 2 is recommended.
n = 2
# Maximum number of iterations to perform for rejection -> 50 is recommended.
max_iterations = 50

# Distribution of f0 {"lognormal", "normal"}. Geopsy default "normal" -> "lognormal" is recommended.
distribution_f0 = "lognormal"
# Distribution of mean curve {"lognormal", "normal"}. Geopsy default "lognormal" -> "lognormal" is recommended.
distribution_mc = "lognormal"


# Manually set the ylimits of the HVSR figures. Default is None so limits will be set automatically.
ymin, ymax = 0, 10


fig = plt.figure(figsize=(8, 4), dpi=150)
ax = fig.add_subplot(1, 1, 1)

start = time.time()
sensor = hvsrpy.Sensor3c.from_mseed(file_name)
bp_filter = {"flag":filter_bool, "flow":filter_flow, "fhigh":filter_fhigh, "order":filter_order}
resampling = {"minf":resample_fmin, "maxf":resample_fmax, "nf":resample_fnum, "res_type":resample_type}
hv = sensor.hv(windowlength, bp_filter, width, bandwidth, resampling, method, f_low=peak_f_lower, f_high=peak_f_upper, azimuth=azimuth)
end = time.time()
print(f"Elapsed Time: {str(end-start)[0:4]} seconds")

individual_width = 0.3
median_width = 1.3
title="Before Rejection"

# Mean Curve 
ax.plot(hv.frq, hv.mean_curve(distribution_mc), color='k', linewidth=median_width)

# Mean +/- Curve
ax.plot(hv.frq, hv.nstd_curve(-1, distribution_mc),
        color='r', linestyle='--', linewidth=median_width, zorder=3)
ax.plot(hv.frq, hv.nstd_curve(+1, distribution_mc),
        color='g', linestyle='--', linewidth=median_width, zorder=3)


#ax.set_xlim(1.5, 8)
ax.set_xscale('log')
ax.set_yscale('log')

ax.set_xlabel("Frequency (Hz)")
ax.set_ylabel("HVSR Amplitude")

Fig2. InSight Sol S0423d VEL 1.5Hz~8.0Hz HVSR

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
# Mean Curve 
ax.plot(hv.frq, hv.mean_curve(distribution_mc), color='k', linewidth=median_width)

# Mean +/- Curve
ax.plot(hv.frq, hv.nstd_curve(-1, distribution_mc),
        color='r', linestyle='--', linewidth=median_width, zorder=3)
ax.plot(hv.frq, hv.nstd_curve(+1, distribution_mc),
        color='g', linestyle='--', linewidth=median_width, zorder=3)


ax.set_xlim(1.5, 8)
ax.set_xscale('log')
ax.set_yscale('log')

ax.set_xlabel("Frequency (Hz)")
ax.set_ylabel("HVSR Amplitude")

Fig3. Possive Data HVSR

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
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
import time
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

import hvsrpy
from hvsrpy import utils


## read syn data from .txt

data_h = np.loadtxt('x.txt')
data_v = np.loadtxt('z.txt')
dt = 1.e-3
sampling_rate = 1/dt

traceX = obspy.Trace(data=np.array(data_h[:]), header={'network': 'a', 'station': 'a', 'location': 'a', 'channel': 'BHE', 'sampling_rate': sampling_rate})
traceY = obspy.Trace(data=np.array(data_h[:]), header={'network': 'a', 'station': 'a', 'location': 'a', 'channel': 'BHN', 'sampling_rate': sampling_rate})
traceZ = obspy.Trace(data=np.array(data_v[:]), header={'network': 'a', 'station': 'a', 'location': 'a', 'channel': 'BHZ', 'sampling_rate': sampling_rate})   

st = obspy.Stream([traceX, traceY, traceZ])

#将数据保存到文件
st.write('data_copy.mseed', format='MSEED')

file_name = 'data_copy.mseed'
sensor = hvsrpy.Sensor3c.from_mseed(file_name)

#file_name = '../InSight-seismic-data-downloader/DATA/2.4_HZ/C/S0423d/S0423d_VEL.mseed'
file_name = 'data_copy.mseed'

windowlength = 60

filter_bool = False       
filter_flow = 0.1                   
filter_fhigh = 10                  
filter_order = 5

# Width of cosine taper {0. - 1.}. Geopsy default of 0.05 is equal to 0.1 -> 0.1 is recommended
width = 0.1

# Konno and Ohmachi smoothing constant. 40 is recommended.
bandwidth = 40
# Minimum frequency after resampling
resample_fmin = 0.1
# Maximum frequency after resampling
resample_fmax = 100
# Number of frequencies after resampling 重采样后的频率点数
resample_fnum = 1000
# Type of resampling {'log', 'linear'}
resample_type = 'linear'
# Upper and lower frequency limits to restrict peak selection. To use the entire range use `None`.
peak_f_lower = None
peak_f_upper = None

# Method for combining horizontal components {"squared-average", "geometric-mean", "single-azimuth"}.
# Geopsy's default is "squared-average" -> "geometric-mean" is recommended.
method = "geometric-mean"
# If method="single-azimuth", set azimuth in degree clock-wise from north. If method!="single-azimuth", value is ignored.
azimuth = 0

# Boolean to control whether frequency domain rejection proposed by Cox et al. (2020) is applied.
# Geopsy does not offer this functionality.
rejection_bool = True
# Number of standard deviations to consider during rejection. Smaller values will reject more windows -> 2 is recommended.
n = 2
# Maximum number of iterations to perform for rejection -> 50 is recommended.
max_iterations = 50

# Distribution of f0 {"lognormal", "normal"}. Geopsy default "normal" -> "lognormal" is recommended.
distribution_f0 = "lognormal"
# Distribution of mean curve {"lognormal", "normal"}. Geopsy default "lognormal" -> "lognormal" is recommended.
distribution_mc = "lognormal"

fig = plt.figure(figsize=(8, 4), dpi=150)
ax = fig.add_subplot(1, 1, 1)

start = time.time()
sensor = hvsrpy.Sensor3c.from_mseed(file_name)
bp_filter = {"flag":filter_bool, "flow":filter_flow, "fhigh":filter_fhigh, "order":filter_order}
resampling = {"minf":resample_fmin, "maxf":resample_fmax, "nf":resample_fnum, "res_type":resample_type}
hv = sensor.hv(windowlength, bp_filter, width, bandwidth, resampling, method, f_low=peak_f_lower, f_high=peak_f_upper, azimuth=azimuth)
end = time.time()
print(f"Elapsed Time: {str(end-start)[0:4]} seconds")

individual_width = 0.3
median_width = 1.3
title="Before Rejection"

# Mean Curve 
ax.plot(hv.frq, hv.mean_curve(distribution_mc), color='k', linewidth=median_width)

# # Mean +/- Curve
# ax.plot(hv.frq, hv.nstd_curve(-1, distribution_mc),
#         color='r', linestyle='--', linewidth=median_width, zorder=3)
# ax.plot(hv.frq, hv.nstd_curve(+1, distribution_mc),
#         color='g', linestyle='--', linewidth=median_width, zorder=3)


#ax.set_xlim(10., 100.)
ax.set_xscale('log')
ax.set_yscale('log')

ax.set_xlabel("Frequency (Hz)")
ax.set_ylabel("HVSR Amplitude")
文章作者: ZY
文章链接: http://erbiaoger.github.io/2023/06/14/2023-06-14-hvsr-ji-fan-yan/
版权声明: 本博客所有文章除特別声明外，均采用 CC BY 4.0 许可协议。转载请注明来源 ZY !
生活
赏
评论
  目录