HTML
-
巴哈马滩是国内外研究碳酸盐岩地层沉积的热点地区,在其西缘斜坡有几口科学钻探井和实测的地震剖面,前人在该地区开展了大量研究[22⁃24]。台地整体形态、初始地形、构造沉降速度、沉积演化过程等都有比较清楚的认识,也有沉积过程定量模拟的相关研究,但都侧重于沉积过程的控制因素分析,缺少精细的沉积过程建模研究。结合地质分析和前人研究成果[25⁃26],用地层沉积正演的方法检测地区地震解释剖面,初步获取该地区碳酸盐岩海平面变化、构造沉降、碳酸盐岩产率等参数,搭建了该地区的地层沉积过程模拟基本框架。沉积模拟范围为32 km×32 km,时间跨度为5.3 Ma,平面网格数为129×129,总时间步数为100。
通过正演模拟的参数尝试,获取了主要参数的分布区间,为生成大量样本和代理模型构建奠定了基础。地层沉积正演模拟要输入的参数有初始地形、构造沉降、海平面变化曲线、碳酸盐岩产率相关参数、沉积物搬运参数等34个。根据对该地区前期的资料分析认为初始地形和构造沉降的不确定性较低,而海平面曲线、碳酸盐岩产率和沉积物搬运的相关参数不确定性较大,建立反演模拟系统时确定性较强的参数给了较小的变化区间,不确定性较大的参数给了较大的变化区间(表1)。这34个参数中,第1~4个参数与初始地形有关(沿着斜坡方向均匀分布的个4点,用来插值生成初始地形),第5~11个参数与海平面曲线有关(分别是表示三级和四级海平面旋回的正弦函数的振幅、周期和相位,以及海平面基准值),第12~16个参数与构造沉降有关(沿着斜坡方向均匀分布的4个点,用来插值生成构造沉降面),第17~20个参数与沉积物搬运有关(分别是势能和动能在X和Y方向的搬运系数),第21~26个参数与产率有关(分别是势能产率最大幅度、势能产率递减系数、透光带厚度,动能产率系数,动能产率基准值,势能产率权重),第27~32个参数与水体动能有关(动能幅度、波浪能下降系数、地形消浪能系数、浪基面、风能系数),第33个参数为地层挠曲系数,第34个参数为生物群落系数,各个参数更详细的介绍可以参考文献[18]。
序号 参数名字 最小值 最大值 序号 参数名字 最小值 最大值 1 InitTopoElev1 -267.00 -67.08 18 PDY1 30 000.00 99 972.00 2 InitTopoElev2 -341.00 -141.08 19 KDX1 200.00 599.84 3 InitTopoElev3 -467.50 -267.08 20 KDY1 200.00 599.84 4 InitTopoElev4 -632.00 -450.07 21 Aprod1 1.30 5.00 5 AmpSeaL0 40.00 79.98 22 Kprod1 0.04 0.12 6 PerSeaL0 230.00 699.81 23 Wprod1 15.00 44.99 7 PhaseSeaL0 -0.50 1.00 24 kEng1 0.50 1.50 8 AmpSeaL1 10.00 49.98 25 bEng1 0.20 0.60 9 PerSeaL1 10.00 99.96 26 pPr1 0.10 0.90 10 PhaseSeaL1 -1.00 1.00 27 kV1 10.00 39.99 11 SeaLevelDatum -50.00 39.96 28 kW1 0.02 0.20 12 SubsRate1 -0.30 0 29 kR1 0.01 0.06 13 SubsRate2 -0.30 -0.05 30 kH1 0.05 0.17 14 SubsRate3 -0.55 -0.01 31 WaveBase1 20.00 59.98 15 SubsRate4 -0.65 -0.05 32 kL1 0 0 16 SubsRate5 -0.87 -0.05 33 DFlex 1E+24 4.998E+26 17 PDX1 30 000.00 79 980.00 34 rkM1 0 0.10 -
纳入全部34个可调整的物理参数,并通过随机抽样生成了大量的样本,构建了深度学习所需的数据库。样本为沉积相剖面相模型及相对应的34个输入参数,沉积相取值0~7,其中0表示背景值,数字1~7对应图2中图例中的a1、b1、b2、c1、c2、c3、c4,表示不同的相类型,该相类型的划分方案同时考虑了沉积物的水体深度和水体能量。模型以二维矩阵形式表示,矩阵尺寸为100×129,建立了地层沉积过程正演模型样本库。通过人工和自动筛选,去除模拟结果明显不真实不合理的样本后,该数据库共计纳入训练样本416 079个,约占总样本的83%;测试样本83 167个,约占总样本的17%。训练样本和测试样本相互独立,保证深度学习模型不会陷入过拟合。从训练样本中随机抽取几个样本(图2),可以看出它们具备不同的沉积相结构样式,体现了样本包含的沉积模式的多样性,这9个样本采用的正演模拟参数具体见附表1。
-
通过测试,选择了改进后的卷积对抗神经网络GLS-GAN[27]。原始GAN要最小化生成分布与真实分布的KL散度,同时要最大化两者的JS散度,在数值上会导致梯度不稳定(通常存在梯度消失和梯度爆炸的可能),并且KL散度具有不对称性,降低了生成器的多样性,会造成模式坍缩问题。GLS-GAN具备损失敏感的自适应损失函数,在训练过程中自动调整损失函数,如果某些区域的真实样本和生成样本已经很接近,则生成器的优化重点就转移到真实样本和生成样本依然差异很大的区域,给GAN提供“按需分配”的建模能力,解决了原版GAN及其衍生算法的梯度消失和模式坍缩问题。
为了保持物理参数和正演结果的一一对应关系,我们在训练生成器时还需要考虑逐个像素的均方误差(Mean Square Error, MSE)。整个网络均由卷积Conv2D单元和转置卷积ConvTranspose2d单元组成,不包含全连接层,每个单元配套一个批量标准化层BatchNorm2d和非线性化层ReLU,通过大量测试推荐采用5~8个卷积层,生成网络和判别网络的具体结构如图3。
利用上述416 079个训练样本,进行56轮次的迭代训练,取得了比较理想的效果。将不同训练阶段的网络应用在测试集的同一套参数上,可以看出不断优化的模拟效果,相序关系变得更加合理(图4)。把最终训练后的网络应用在测试集,随机抽取几个模型,都具备合理的沉积相分布模式(图5),这6个样本的输入参数见附表2。可见相同输入参数下LS-GAN模型的输出与地层沉积过程模型的输出一致性很强,虽然没有达到完全相同,但是基本实现了合理的相带组合特征。通过34个参数生成符合条件的相模型本身是个难度很大问题,GLS-GAN很大程度上改善了模式坍缩问题,但训练难度大、收敛过程不稳定是GAN类方法自身的属性,无法完全避免。另外,网络结构、样本数量,样本质量等也会影响最终的效果。目前认为形成的替代模型可以用于下一步的反演模拟。
Figure 5. Comparison between generator simulation (false) and conventional forward simulation (true) on randomly selected 6 samples from test set
参数名称 样本A 样本B 样本C 样本D 样本E 样本F InitTopoElev1 -235.872 -211.007 -240.939 -226.946 -216.507 -134.497 InitTopoElev2 -294.878 -317.912 -305.656 -283.186 -314.672 -294.664 InitTopoElev3 -440.537 -404.218 -405.011 -363.598 -402.461 -429.075 InitTopoElev4 -611.515 -604.388 -572.029 -611.646 -603.645 -602.346 AmpSeaL0 45.907 7 44.462 8 54.435 1 44.457 3 45.539 2 45.516 2 PerSeaL0 302.895 297.319 283.025 327.798 305.036 285.534 PhaseSeaL0 -0.313 27 -0.268 81 -0.295 00 -0.265 74 -0.125 55 -0.318 12 AmpSeaL1 31.055 7 16.873 3 14.874 9 16.619 7 24.722 8 16.481 8 PerSeaL1 22.732 1 24.420 0 21.763 3 42.429 3 51.685 6 38.696 4 PhaseSeaL1 -0.751 93 -0.711 99 -0.647 79 -0.664 00 -0.199 81 -0.670 46 SeaLevelDatum -23.082 0 -29.752 5 -33.811 0 -35.430 4 -37.047 0 0.349 351 SubsRate1 -0.236 21 -0.122 30 -0.224 62 -0.231 28 -0.232 31 -0.252 31 SubsRate2 -0.248 49 -0.230 61 -0.260 81 -0.253 01 -0.245 30 -0.262 11 SubsRate3 -0.450 76 -0.489 44 -0.455 16 -0.483 23 -0.454 18 -0.485 23 SubsRate4 -0.572 49 -0.528 56 -0.519 46 -0.463 54 -0.536 64 -0.538 34 SubsRate5 -0.802 85 -0.600 26 -0.794 73 -0.775 51 -0.636 48 -0.770 71 PDX1 39 367.2 51 846.1 50 318.1 48 834.0 40 938.0 37 031.7 PDY1 49 483.2 38 700.8 41 724.9 46 108.6 40 312.3 40 762.6 KDX1 303.702 319.104 317.759 346.898 267.878 261.639 KDY1 262.673 258.522 300.591 355.141 252.191 270.105 Aprod1 2.488 89 1.720 79 2.116 18 1.935 73 1.895 37 1.905 39 Kprod1 0.078 654 0.057 464 0.058 098 0.071 617 0.050 401 0.054 544 Wprod1 18.519 7 19.082 2 21.298 2 27.134 6 19.764 0 18.834 8 kEng1 0.682 237 0.835 305 0.670 850 0.660 123 0.678 857 0.826 163 bEng1 0.348 771 0.268 214 0.268 690 0.302 373 0.271 550 0.259 995 pPr1 0.345 856 0.399 625 0.372 746 0.233 566 0.392 512 0.267 308 kV1 16.119 3 17.201 6 19.910 2 14.708 3 16.862 4 15.730 3 kW1 0.063 514 0.042 187 0.049 663 0.048 818 0.040 613 0.063 422 kR1 0.016 335 0.018 669 0.019 181 0.035 421 0.023 036 0.015 796 kH1 0.074 975 0.070 531 0.065 382 0.063 510 0.086 667 0.063 957 WaveBase1 26.208 4 36.326 0 34.947 2 24.898 5 26.646 7 33.491 0 kL1 0.000 105 0.000 138 8.20E-05 4.97E-05 0.000 124 4.43E-05 DFlex 1.34E+26 1.17E+26 9.29E+25 6.42E+25 8.00E+25 1.58E+26 rkM1 0.039 695 0.012 137 0.030 700 0.020 077 0.014 142 0.022 487 -
在开始真实观测数据反演模拟之前,先从测试集中随机选择样本进行参数反演。该样本的参数及其沉积过程模型模拟结果都是已知的,可以用于验证参数反演是否有效。在水平方向等间距设置了6个虚拟测井,作为反演系统的观测数据。用观测数据与模拟结果的均方根误差作为目标函数,优化器采用上述提到的SCE-UA算法[19]。
用SCU-UA进行了500 000次迭代优化,得到了预期的效果(图6)。图6a展示了参数的真实值(灰色虚线)和SCE-UA算法得到的最优参数值(蓝色实线),图6b展示了SCE-UA优化过程中目标函数的变化,图6c展示了替代模型LS-GAN给出的地层沉积相,其中虚拟测井位置上的沉积相用真实值(即沉积过程模型模拟值)代替。模拟结果表明,随着优化过程的进行,目标函数稳定下降并不再降低,表明SCE-UA算法已经收敛。然而反演所得到的参数值与参数的真实值有较大差异,且互不相同。这一结果表明地层沉积过程的参数反演是一个高度复杂的问题,存在多个局部极小值。
采用类似的方法进行真实观测数据的反演。将代理模型加入地层沉积反演模拟系统,替换传统的地层正演模拟器,实现参数反演。选择4口科考钻井的人工解释的沉积数据为观测数据,将钻井穿过网格的像素均方根误差为目标函数,SCE-UA为反演优化算法,进行参数反演。对于真实数据,我们并不知道测井以外区域的沉积相,也不知道真实的物理参数的值,因此只能靠专家经验判断结果是否合理。从反演结果可以看出,使用替代模型反演得到的地层沉积相的层序结构明显,分布较为合理(图7)。与基于正演模拟器的反演方法相比,结果一致性较强,效率有大幅提升(图8)。
2.1. 巴哈马滩反演模型基本框架搭建
2.2. 基于沉积正演模拟的大规模样本库建立
2.3. 基于生成对抗网络的代理模型构建
2.4. 基于代理模型的反演模拟
-
(1) 地层沉积过程反演模拟系统是在优化器的驱动下,不断调整正演模拟输入参数,提高模拟结果与观测数据的吻合程度,大大增加了沉积过程模拟建模方法的实用性,但是存在反演系统非线性强、耗时长、效率较低的问题。
(2) 生成对抗网络具有较强的模式生成能力,在一定条件下,利用该方法建立的碳酸盐岩地层沉积过程模拟系统的代理模型,具有较强的泛化能力。
(3) 将代理模型带入沉积模拟反演系统,替换正演模拟器,可以大大提高反演的效率。通过巴哈马滩西缘斜坡现代沉积的反演模型实例,验证了该方法的可行性。
参数名称 样本A 样本B 样本C 样本D 样本E 样本F 样本G 样本H 样本I 样本J InitTopoElev1 -241.740 -239.795 -240.203 -189.967 -210.110 -237.965 -212.366 -176.158 -213.844 -228.109 InitTopoElev2 -304.687 -308.444 -288.532 -302.394 -315.569 -316.927 -311.806 -313.322 -299.977 -290.627 InitTopoElev3 -443.803 -420.896 -440.692 -372.690 -426.142 -429.549 -419.037 -445.090 -444.184 -438.012 InitTopoElev4 -608.549 -594.739 -606.268 -608.710 -611.325 -600.513 -599.063 -603.759 -597.293 -608.271 AmpSeaL0 47.993 0 49.966 8 49.034 4 52.207 7 46.511 8 47.745 1 51.338 5 63.144 5 46.858 5 50.374 2 PerSeaL0 415.966 333.573 350.251 319.018 288.357 427.918 302.105 338.303 314.502 308.967 PhaseSeaL0 -0.054 27 -0.013 71 -0.303 06 -0.270 38 -0.078 13 -0.331 04 -0.128 83 0.177 19 -0.295 01 -0.309 59 AmpSeaL1 19.470 5 15.821 2 15.208 0 14.819 3 18.243 8 19.648 7 29.210 7 15.980 5 15.315 4 17.789 8 PerSeaL1 36.966 0 50.087 9 29.274 7 25.920 3 28.986 7 29.643 4 35.675 9 30.717 8 37.475 2 27.578 2 PhaseSeaL1 -0.440 86 -0.645 22 -0.632 37 -0.568 44 -0.768 05 -0.692 32 -0.678 35 -0.655 46 -0.690 46 -0.746 04 SeaLevelDatum -39.560 5 -38.872 9 -37.119 8 -34.170 6 -34.217 7 -38.279 3 -39.223 3 -37.920 3 -39.168 8 -27.961 8 SubsRate1 -0.243 89 -0.247 15 -0.253 77 -0.258 66 -0.258 47 -0.228 18 -0.243 00 -0.236 90 -0.242 21 -0.195 93 SubsRate2 -0.260 98 -0.219 83 -0.251 05 -0.240 68 -0.228 05 -0.246 30 -0.163 37 -0.237 58 -0.260 05 -0.249 80 SubsRate3 -0.370 92 -0.417 32 -0.476 49 -0.453 03 -0.408 79 -0.467 76 -0.427 09 -0.454 19 -0.458 59 -0.487 05 SubsRate4 -0.563 53 -0.496 77 -0.527 97 -0.487 82 -0.444 86 -0.581 18 -0.539 59 -0.519 66 -0.311 60 -0.544 13 SubsRate5 -0.752 94 -0.781 76 -0.671 17 -0.763 93 -0.687 47 -0.667 36 -0.593 09 -0.723 55 -0.767 86 -0.653 68 PDX1 37 175.7 43 489.0 39 140.9 35 975.7 37 643.6 36 044.0 36 179.2 40 122.3 40 944.5 52 077.4 PDY1 57 366.0 54 775.4 44 088.4 48 277.3 44 436.4 40 488.0 39 890.3 40 007.2 38 728.0 38 879.8 KDX1 275.378 266.399 245.489 256.212 251.925 276.658 245.710 255.417 396.307 275.302 KDY1 262.749 249.093 278.053 312.542 309.070 284.470 270.730 259.157 435.309 438.423 Aprod1 2.617 08 2.159 54 2.729 96 2.279 54 2.192 18 1.835 66 1.766 42 3.566 48 2.022 24 1.734 97 Kprod1 0.072 577 0.049 217 0.049 603 0.051 510 0.052 695 0.053 192 0.049 730 0.059 071 0.049 713 0.049 448 Wprod1 18.598 8 19.906 0 22.985 8 22.905 6 26.301 7 42.079 8 23.508 3 24.068 4 23.353 3 18.509 9 kEng1 0.671 643 0.738 103 0.727 784 0.666 155 0.819 835 0.799 826 0.620 340 0.714 996 0.923 437 0.725 410 bEng1 0.271 113 0.282 144 0.299 423 0.336 546 0.261 021 0.317 633 0.246 501 0.282 42 0.257 869 0.391 195 pPr1 0.206 796 0.196 949 0.268 695 0.204 356 0.305 793 0.231 975 0.261 185 0.364 455 0.210 016 0.258 421 kV1 28.296 3 15.356 4 14.039 5 13.641 2 17.244 1 15.175 5 15.189 6 15.865 20.368 9 16.541 7 kW1 0.058 621 0.051 851 0.061 631 0.044 944 0.053 190 0.125 966 0.043 592 0.060 993 0.050 890 0.041 966 kR1 0.022 455 0.024 601 0.025 026 0.034 426 0.015 812 0.016 974 0.016 930 0.017 06 0.017 712 0.016 261 kH1 0.063 897 0.065 008 0.120 705 0.069 272 0.063 445 0.077 070 0.137 212 0.083 297 0.077 088 0.076 075 WaveBase1 28.713 3 27.629 2 30.370 3 25.878 6 31.746 9 28.013 1 27.847 9 24.5619 30.966 8 31.201 4 kL1 0.000 107 4.65E-05 4.53E-05 6.79E-05 5.09E-05 4.26E-05 6.15E-05 4.73E-05 5.87E-05 4.86E-05 DFlex 1.19E+26 7.95E+25 7.63E+25 6.58E+25 1.13E+26 5.80E+25 6.33E+25 5.81E+25 7.15E+25 1.94E+26 rkM1 0.016 568 0.027 596 0.031 403 0.012 224 0.031 896 0.017 919 0.015 168 0.018 852 0.014 158 0.034 389