python 地震数据可视化

参考源码

https://github.com/agile-geoscience/seisplot

准备工作：

在windows10下安装python3.7，下载参考源码到本地。

1. demo绘图测试

demo绘图指令

cmd> python seisplot.py --demo

问题1）缺少依赖包

 File "D:/Desktop/python可视化/seisplot/seisplot.py", line 16, in <module>import yaml
ModuleNotFoundError: No module named 'yaml'>pip3 install yamlCould not find a version that satisfies the requirement yaml (from versions: )
No matching distribution found for yaml

查看environment.yml:

channels:
- defaults
- conda-forge
dependencies:
- matplotlib
- numpy
- obspy
- pillow
- pyyaml

安装相关依赖包，如 obspy pillow pyyaml

cmd> pip3 install obspy

出现proxy error，添加选项 --proxy server:port

安装ok

问题2）配置文件的编码异常

D:\Desktop\python可视化\seisplot>python seisplot.py --demo
[91mWelcome to┌─┐┌─┐┬┌─┐┌─┐┬  ┌─┐┌┬┐└─┐├┤ │└─┐├─┘│  │ │ │└─┘└─┘┴└─┘┴  ┴─┘└─┘ ┴Good luck[0m
Traceback (most recent call last):File "seisplot.py", line 421, in <module>cfg = yaml.load(f)File "D:\ProgramFiles\Python36\lib\site-packages\yaml\__init__.py", line 70, in loadloader = Loader(stream)File "D:\ProgramFiles\Python36\lib\site-packages\yaml\loader.py", line 34, in __init__Reader.__init__(self, stream)File "D:\ProgramFiles\Python36\lib\site-packages\yaml\reader.py", line 85, in __init__self.determine_encoding()File "D:\ProgramFiles\Python36\lib\site-packages\yaml\reader.py", line 124, in determine_encodingself.update_raw()File "D:\ProgramFiles\Python36\lib\site-packages\yaml\reader.py", line 178, in update_rawdata = self.stream.read(size)
UnicodeDecodeError: 'gbk' codec can't decode byte 0x99 in position 2003: illegal multibyte sequence

分析源码中参数传递情况：

def load(stream, Loader=Loader):"""Parse the first YAML document in a streamand produce the corresponding Python object."""print(stream) ## debug by huanying03loader = Loader(stream)try:return loader.get_single_data()finally:loader.dispose()

打印结果如下：

<_io.TextIOWrapper name='config.yml' mode='r' encoding='cp936'>

编码不对？

65001 ：UTF-8代码页

936 ：默认的GBK

437 ：是美国英语

仔细检查发现config.yml中有异常字符，更正config.yml中错误字符：ok

问题3）python调试打印语句

print ("value=%d"%value)

demo数据绘图结果

如下：

2. 本地数据绘图测试

应用指令

cmd > python seisplot.py data/my.sgy

问题1）数据道的采样值异常

filename   data2/t10.sgy
seismic.py:104:dt=4000
seismic.py:104:ns=512
seismic.py:144:header=b'C 1 CLIENT'
n_traces   500
n_samples  512
dt         0.004
t_start    0
t_end      2.044
max_val    inf
min_val    -inf
clip_val   180462501764003194804057887685476352
Read data in 0.4 s++++++++++++++++++++++++++++++++++++++++++++++++++++++++
PLOTTING
plot width   17.00 in
plot height  11.00 in
Traceback (most recent call last):File "seisplot.py", line 442, in <module>main(t, cfg)File "seisplot.py", line 214, in maincfg)File "D:\Desktop\python可视化\seisplot\plotter.py", line 123, in plot_histogramy, x, _ = ax.hist(np.ravel(data), bins=int(100.0 / (clip_val / largest)),
OverflowError: cannot convert float infinity to integer

解决办法：将数据做归一化处理

问题2）数据文件的卷头异常

filename   data2/marmousi_vel.segy.hdrs.sgy.cdp.sgy
seismic.py:104:dt=10000
seismic.py:104:ns=350
seismic.py:144:header=b'\xc3@\xf1@\xc3\xd3\xc9\xc5\xd5\xe3'
Traceback (most recent call last):File "seisplot.py", line 442, in <module>main(t, cfg)File "seisplot.py", line 40, in mains = Seismic.from_segy(target, params={'ndim': cfg['ndim']})File "D:\Desktop\python可视化\seisplot\seismic.py", line 164, in from_segyreturn cls.from_obspy(stream, params=params)File "D:\Desktop\python可视化\seisplot\seismic.py", line 147, in from_obspyx =np.array(list(stream.textual_file_header.decode()))
UnicodeDecodeError: 'utf-8' codec can't decode byte 0xc3 in position 0: invalid continuation byte

解决办法：重写3200字节的segy卷头

marmousi相关数据绘图

速度模型绘图结果如下：

marmousi模型正演的炮集数据绘图结果如下：

测试备忘

1

cygwin执行指令，可以看到反馈信息为彩色，分析信息很方便。

相比之下，cmd终端执行指令后的反馈信息可读性较差。

2

本地数据出现异常后，通过cwp中的segyread, sushw, sugain, segywrite等程序进行修改，满足绘图程序要求。

1）3200字节的卷头必须规范

2）道头字中的cdpt取值必须规范（xline）

3）道采样值最好归一化处理。