生成界面坐标文件

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
# 界面1
# x z, 第一列为x坐标, 第二列为z坐标
Flow('bvel1.asc',None,
     '''
     echo
     0    700
     500  900
     1000 1000
     1500 900
     2000 700
     2500 900
     3000 1000
     3500 900
     4000 700
     n1=2 n2=9 in=$TARGET
     data_format=ascii_float
     ''')

可以按照上述,生成多个界面的坐标信息文件

将界面坐标进行插值

1
2
3
4
5
6
# n1 表示x轴
Flow('lay1','bvel1.asc',
     '''
     dd form=native |
     spline n1=401 d1=10 o1=0
     ''')

插值的目的是为了得到一个规则网格的界面坐标信息

建模

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
# 建模前首先将多个界面信息文件拼接成一个总的文件
Flow('lays','lay1 lay2','cat axis=2 ${SOURCES[1:2]}')

# 采用sfunif2来建模,和su中的unif2类似
nz = 201; dz = 10; oz = 0 # spatial sampling in z
Flow('vel1','lays',
     '''
     unif2 n1=%d d1=%g o1=%g v00=2000,2500,3000 |
     put label1=Depth unit1=m label2=Distance unit2=m
     label=Velocity unit=m/s
     ''' % (nz,dz,oz) )

将上述各部分合并成一个文件,并进行适当的优化,新的脚本如下

 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
# velocity model construction
nz = 201; dz = 10; oz = 0 # spatial sampling in z
nx = 401; dx = 10; ox = 0 # spatial sampling in x

# 注意,每个界面至少有4个坐标,与三次样条插值sfspline有关
layers = (
     ((0,700),(500,900),(1000,1000),(1500,900),(2000,700),(2500,900),(3000,1000),(3500,900),(4000,700)), # 第一个界面
     ((0,1500),(1000,1500),(2000,1500),(3000,1500),(4000, 1500)) # 第二个界面
)

def arr2str(array,sep=' '):
    return string.join(map(str,array),sep)

nlays = len(layers)

# 将生成界面坐标文件与插值合并
for i in range(nlays):
    inp = 'inp%d' % i
    Flow(inp+'.asc',None,
          '''
          echo %s in=$TARGET
          data_format=ascii_float n1=2 n2=%d
          ''' % \
          (' '.join(map(lambda x: ' '.join(map(str,x)),
                           layers[i])),len(layers[i])))
    lay = 'lay%d' % i
    Flow(lay,inp+'.asc','dd form=native | spline n1=%d d1=%g o1=%g' %(nx,dx,ox))

Flow('lays','lay0 lay1 ',
     'cat axis=2 ${SOURCES[1:2]}')

Flow('vel','lays',
    '''
    unif2 n1=%d d1=%g o1=%g v00=1500,2000 |
    put label1=Depth unit1=m label2=Distance unit2=m
    label=Velocity unit=m/s
    ''' % (nz,dz,oz) )

运行上述脚本,得到的结果图如下:

vel1.png

参考

修订历史

  • 2018-05-21:初稿;
  • 2018-12-18:修复文档中的错误信息。