地震层析成像的论文中很多采用Paraview来绘制三维速度模型。在此简单介绍一下如何将地震层析成像得到的速度模型用 Paraview绘制出来。

假设我们使用的层析成像代码不支持生成VTK文件，那么第一步就需要手动生成一个VTK文件。再次假定我们得到的层析成像 的结果的格式如下：

1

longitude latitude depth vp

我们可以使用Python来生成VTK文件。这里主要参考 https://blog.xumijian.me/post/paraview-model/, 由于现在 PyVista的版本已经更新，因此原博文的代码直接运行会报错，这里针对新的版本进行了修订。

安装PyVista 和 netCDF4

1
2

$ pip install pyvista
$ pip install netCDF4

我这里安装的版本如下：

1
2

netCDF4                  1.7.2
pyvista                  0.45.2

生成vtk文件

写一个python程序，命名为convertToVtk.py,具体内容如下：

 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

import numpy as np
from scipy.interpolate import griddata
from scipy import ndimage
import scipy.signal as signal
import pyvista as pv
from netCDF4 import Dataset

#https://blog.xumijian.me/post/paraview-model/

class XYZ2VTK():
    def __init__(self, lat1=21, lon1=97, lat2=34, lon2=108, maxdep=70):
        self.lat1 = lat1
        self.lat2 = lat2
        self.lon1 = lon1
        self.lon2 = lon2
        self.lats = np.arange(lat1, lat2, 0.5)
        self.lons = np.arange(lon1, lon2, 0.5)
        self.dep = np.arange(-maxdep, 0, 5)
        self.vp_file = 'vp.txt'
        # self.moho_file = '../ccp3d/good_moho_ET2020_repick.dat'
        # self.topo_file = '../figures/topo_resample.grd'

    def read_xyz(self):
        lon, lat, dep, vp = np.loadtxt(self.vp_file,  unpack=True)
        dep /= -1
        new_lon, new_lat, new_dep = np.meshgrid(self.lons, self.lats, self.dep, indexing='ij')
        grid_data = griddata((lon, lat, dep), vp, (new_lon, new_lat, new_dep))
        return grid_data

    def create_vtk(self):
        grid = pv.ImageData()
        grid_vp = self.read_xyz()
        grid.dimensions = grid_vp.shape
        grid.origin = (self.lon1, self.lat1, self.dep[0])
        grid.spacing = (0.5, 0.5, 5)
        grid.point_data["vp"] = grid_vp.flatten(order="F")
        grid.save('vp.vtk')

if __name__ == "__main__":
    converter = XYZ2VTK()
    converter.create_vtk()
    print("VTK file created successfully: vp.vtk")

我们从网上找一个川滇的公共速度模型简单测试一下

1
2
3
4
5
6
7

# 下载速度模型
$ git clone https://github.com/liuyingustc/SWChinaCVM-1.0.git
$ cd SWChinaCVM-1.0
# 只取前4列数据
$ awk '{print $1,$2,$3,$4}' SWChinaCVMv1.0.txt > vp.txt
# 将其转换成vtk文件
$ python convertToVtk.py

在Paraview中画图

有了vtk文件后，在Paraview中画图用视频教程更容易上手，这里大家可以参考B站视频教程。这里我主要讲解一下如何调整一些参数，让画图更好看。

改变坐标轴的相对比例

如图1所示，首先在红框区域点击右键，选择open可导入vtk文件，或者直接把vtk文件拖到红框区域也可以打开vtk文件。

图1

图1显示的模型垂向上有点过长，因此需要设置一些参数使其看上去更美观一些。首先先找到图2中1标记的位置，打开齿轮工具，然后 找到图2中2标记的scale，第三列对应的就是深度方向的scale，我这里将其调整为之前的0.1，可以看到对应右边的绿框中显示的深度 轴已经变成原来的十分之一了。但同时也带来了一个问题，就是对应的深度的坐标刻度也发生了变化。

图2

接下来，首先找到图3中1标记的Axes Grid点击Edit，然后找到图3中2标记的齿轮工具并点击它，接着找图3中3标记的Data scale的第3列，同样，将其值调整为0.1，然后点击下面的Apply, 你会发现速度模型图的深度坐标的刻度恢复正常了。这里我由于 已经点过了Apply，所以图3中显示的是灰色的。

图3

通过刚才的试验，可以发现只需将图2和图3的两个scale设置为相同的值，最后的模型图中的数值是和实际是一致的。

只绘制指定的几个刻度值

要改变坐标轴的刻度值，步骤和图3类似，主要步骤见图4。安装图中的顺序依次进行操作，其中前2步和图3一样，不再赘述。第三步，首先勾选Z Axis Use Custom Labels, 然后，点击4处的+,这里，我依次输入了-10，-30，-50，-70四个深度值。然后点击 5处的Apply,可以看到图中只显示了4个刻度值。这里我由于 已经点过了Apply，所以图4中显示的是灰色的。

图4

参考资料