CWPSU用自定义模型产生反射波法模拟数据

目录
  1. 利用CWPSU自建模型进行反射波法模拟数据
    1. 前情提要
    2. SUSYNLV
      1. Usage
      2. Main Parameters
      3. Other Parameters
      4. Codes
    3. 产生结果
      1. 运行
      2. 抽道
    4. 预告

利用CWPSU自建模型进行反射波法模拟数据

拖了很久的反射波法终于出炉了。限于时间，本期只有数据的合成，用法的解释。对于数据后续的处理，放到明天再展开，同学们可以自己先捣鼓。由于有了数据，处理也会如鱼得水。

前情提要

为什么反射波法拖了这么久，其实还是因为一个矛盾：自建模型和线性数据的矛盾。此话怎讲呢，本人看了很多网上包括前辈们的操作，大抵相似，都是使用susynlv来做，不过susynlv是一个模拟线性地层的工具，那不就先前自建的各种崎岖模型就失去意义了吗？对此，网上很多攻略忽略了这一点，有的甚至还滥竽充数，分析的头头是道，结果前后矛盾。如果说只是单纯的模拟一个在深度和方向上完全线性的地层，那么很简单，直接用susynlv就可以做到，不过显然我们不想这么咸鱼。

那么CWPSU提供了susynvxz和susynvxzcs等一系列可以输入模型的方法，这些方法参数众多，琢磨很久也没有完全摸透，索性放弃。不过之前在外网看到一个用trimodel做的模型，再用triseismic做的反射成像，图很漂亮：

cwpsureflect1

cwpsureflect2

奈何不是多次覆盖。难以调用。

我们回到最初的起点，回到susynlv

SUSYNLV

我们回到了起点，我们之前一直不想用它是因为所有攻略中它不能带入模型文件，只能做一个线性的地层，或者一个速度为常数的模型，而且那些作者也大多是这么干的。那么susynlv能不能做一个自定义的地层反射波模拟呢？我们先对它的Usage做一个介绍：

不在意原理的同学可以直接到代码块,点击跳到Codes

Usage

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
 SUSYNLV - SYNthetic seismograms for Linear Velocity function

 susynlv >outfile [optional parameters]

 Optional Parameters:
 nt=101                 number of time samples
 dt=0.04                time sampling interval (sec)
 ft=0.0                 first time (sec)
 kilounits=1            input length units are km or kilo-feet
                         =0 for m or ft
                        Note: Output (sx,gx,offset) are always m or ft
 nxo=1                  number of source-receiver offsets
 dxo=0.05               offset sampling interval (kilounits)
 fxo=0.0                first offset (kilounits, see notes below)
 xo=fxo,fxo+dxo,...     array of offsets (use only for non-uniform offsets)
 nxm=101                number of midpoints (see notes below)
 dxm=0.05               midpoint sampling interval (kilounits)
 fxm=0.0                first midpoint (kilounits)
 nxs=101                number of shotpoints (see notes below)
 dxs=0.05               shotpoint sampling interval (kilounits)
 fxs=0.0                first shotpoint (kilounits)
 x0=0.0                 distance x at which v00 is specified
 z0=0.0                 depth z at which v00 is specified
 v00=2.0                velocity at x0,z0 (kilounits/sec)
 dvdx=0.0               derivative of velocity with distance x (dv/dx)
 dvdz=0.0               derivative of velocity with depth z (dv/dz)
 fpeak=0.2/dt           peak frequency of symmetric Ricker wavelet (Hz)
 ref="1:1,2;4,2"        reflector(s):  "amplitude:x1,z1;x2,z2;x3,z3;...
 smooth=0               =1 for smooth (piecewise cubic spline) reflectors
 er=0                   =1 for exploding reflector amplitudes
 ls=0                   =1 for line source; default is point source
 ob=1                   =1 to include obliquity factors
 tmin=10.0*dt           minimum time of interest (sec)
 ndpfz=5                number of diffractors per Fresnel zone
 verbose=0              =1 to print some useful information

Notes:
Offsets are signed - may be positive or negative.  Receiver locations
are computed by adding the signed offset to the source location.

Specify either midpoint sampling or shotpoint sampling, but not both.
If neither is specified, the default is the midpoint sampling above.

More than one ref (reflector) may be specified. Do this by putting
additional ref= entries on the commandline. When obliquity factors
are included, then only the left side of each reflector (as the x,z
reflector coordinates are traversed) is reflecting.  For example, if x
coordinates increase, then the top side of a reflector is reflecting.
Note that reflectors are encoded as quoted strings, with an optional
reflector amplitude: preceding the x,z coordinates of each reflector.
Default amplitude is 1.0 if amplitude: part of the string is omitted.

东西很多，不过这里只对我们需要了解的东西做一个介绍。

Main Parameters

我们做反射波法模拟，多次覆盖法。需要重点定义：

  • 1 道间距
  • 2 道数
  • 3 偏移距
  • 4 覆盖次数

本博客做出如下定义：

紧跟上次模型，长500米，深80米。100道，5米一道，起始道位于5米处。25次覆盖（打25炮），根据公式，每炮移动2个道间距（10米），我们采用零偏。那么这些参数怎么在这里体现呢？

经过实验发现：

nxo代表道数、dxo代表道间距、fxo代表第一个检波点位置（默认模型起始为0米）

我们还需要定义炮点参数。在Notes里我们可以了解到，定义炮点参数只需定义nxm和nxs中的一个系列即可，我们选择nxs系列参数：

nxs代表炮数（覆盖数）dxs代表每炮移动距离 fxs代表起始炮位置

Other Parameters

我们还需定义一些其他参数，例如采样频率，单位制，速度分布，界面信息等。

本文定义：

采样频dt=0.01 采样次数301 起始采样0s开始

单位制为mkilounits=0

我们说过，我们想把自己的自定义模型用起来，那么susynlv如何做到？我们发现，它提供了一个v00和ref的参数，这使得我们的想法成为可能。

我们定义速度分布为v00=1000,1700,2300

界面信息ref为

ref1=”0.3:0.0,30;100.0,36;200.0,42;300,47;400,40;500,24”
ref2=”0.7:0,60;500,60”
ref3=”0.8:0,80;500,80”

第一个为反射系数，后面是界面信息。

其实这就是我们之前建立的Model。你可以前去查看：

点击前去

其实到这里你就发现了，我们的susynlv是可以做到加载自定义而非线性模型了，只是我们之前用unif2建立的文件只是用来显示给别人看的，在这里我们还需重新设置一次相同的参数，同时把线性模型参数全部设为0。

Codes

我们解释了很多参数，这里给出本文使用的shell脚本dg.sh

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
#constant parameters
dt=0.01
nt=301
ft=0.0
kilounits=0

nxo=100
dxo=5
fxo=5

nxs=25
dxs=10
fxs=5
ref1="0.3:0.0,30;100.0,36;200.0,42;300,47;400,40;500,24"
ref2="0.7:0,60;500,60"
ref3="0.8:0,80;500,80"
susynlv > dataraw.su nt=$nt dt=$dt ft=$ft dvdx=0 dvdz=0 nxo=$nxo dxo=$dxo fxo=$fxo nxs=$nxs dxs=$dxs fxs=$fxs ref=$ref1 ref=$ref2 ref=$ref3 v00=1000,1700,2300 kilounits=$kilounits verbose=1

因为我们不再使用完全线性地层，所以两个梯度全部是0。verbose只是可以显示更多信息，不用在意。
在这里我们特别说明一下，这里我们让nxs=25是为了理解多次覆盖的方便，我们将在真正数据处理的过程中将它改为49，这是为了满足覆盖的点多一些 数据处理的链接在此处
所以我们的模型可以概括如下：

在这样的一个地层中：

1571400551475

（注：我们为了使得界面更加明显，崎岖，中间凹的很厉害）

我们在上面布设了100个检波点，5米一个，起始位置5米处。25次覆盖，每次移动2个道间距（10米）。零偏。文件储存在dataraw.su中。

产生结果

运行

我们运行

1
sh ./dg.sh

注意此时不要直接用suxwigb打开，因为它是2500个道集，打开就全是黑的，且道数增倍：

1571401128981

抽道

我们需要抽出我们想要的道数，为了说明一般性，我们抽1和25炮的记录：

输入命令：

1
suwind <dataraw.su key=fldr max=1 max=1 |suxwigb &
1
suwind <dataraw.su key=fldr min=25 max=25 |suxwigb &

分别获得第1炮和第25炮的数据记录：

第1炮：

1571401456015

第25炮：

1571401426460

我们还是能看出区别的。

不过我们需要特别说明：

在第n炮的记录中，第一道并不是我们原先布设的第一道，可以理解，因为我们炮点移动了，比如第2炮的第一道就相对于第一炮的第3道了。也就是说我们的100道是往后顺延的。

这也是为什么我们后面处理需要抽道选排了。

后续的数据处理我们还需要加噪声，动校正等等，这个后续会更新，其实只要数据模拟出来了，后面的处理大同小异，相信大家能够完成，当然你也可以期待近期更新（预计19号） 已更新，这是链接

预告

下期我们将带来模拟数据的处理，噪声、选排、校正、偏移等。再会！