CN112270064B - 一种岩层产状的计算方法及系统 - Google Patents

一种岩层产状的计算方法及系统 Download PDF

Info

Publication number
CN112270064B
CN112270064B CN202010848655.XA CN202010848655A CN112270064B CN 112270064 B CN112270064 B CN 112270064B CN 202010848655 A CN202010848655 A CN 202010848655A CN 112270064 B CN112270064 B CN 112270064B
Authority
CN
China
Prior art keywords
attitude
projection
calculating
occurrence
stratum
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202010848655.XA
Other languages
English (en)
Other versions
CN112270064A (zh
Inventor
刘修国
花卫华
赵亚博
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China University of Geosciences
Original Assignee
China University of Geosciences
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by China University of Geosciences filed Critical China University of Geosciences
Priority to CN202010848655.XA priority Critical patent/CN112270064B/zh
Publication of CN112270064A publication Critical patent/CN112270064A/zh
Application granted granted Critical
Publication of CN112270064B publication Critical patent/CN112270064B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/004—Artificial life, i.e. computing arrangements simulating life
    • G06N3/006—Artificial life, i.e. computing arrangements simulating life based on simulated virtual individual or collective life forms, e.g. social simulations or particle swarm optimisation [PSO]

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Data Mining & Analysis (AREA)
  • General Engineering & Computer Science (AREA)
  • Pure & Applied Mathematics (AREA)
  • Evolutionary Computation (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Computational Mathematics (AREA)
  • Software Systems (AREA)
  • Computing Systems (AREA)
  • Artificial Intelligence (AREA)
  • General Health & Medical Sciences (AREA)
  • Molecular Biology (AREA)
  • Computational Linguistics (AREA)
  • Computer Hardware Design (AREA)
  • Biophysics (AREA)
  • Biomedical Technology (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Health & Medical Sciences (AREA)
  • Geometry (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本发明公开了一种岩层产状的计算方法及系统，该方法包括获取实测产状及其对应的权重系数；获取地质平面图及剖面、地球物理平面图及剖面解释图，并从中提取出地层延伸方向及其对应的权重系数；在提取的地层延伸方向处设置初始产状；利用插值计算算法，计算实测产状点处及地层延伸方向点处初始产状的交叉验证产状；使用参数角θ对地层延伸方向点处初始产状的交叉验证产状进行表达，利用投影计算方法，计算参数角θ对应的投影产状；计算实测产状与实测产状点处交叉验证产状之间，以及地层延伸方向点处的交叉验证产状与其对应的投影产状之间的夹角；结合所得的两个夹角，构建目标优化函数，最后针对目标优化函数进行迭代优化，进行最优投影产状的计算。

Description

一种岩层产状的计算方法及系统
技术领域
本发明涉及地质学构造解释领域，更具体地说，涉及一种岩层产状的计算方法及系统。
背景技术
岩层产状是地质学家了解地质构造的重要依据，广泛应用于二维和三维地质填图、定性和定量地质解释等方面。然而，受制于覆盖层、地形、资金等因素，实测产状往往空间分布不均且稀疏，在地表覆盖区和地下中更甚。如何自动添加大量、分布广泛且符合地质构造的产状一直是进行三维地质建模或精细构造解释所面临的重要问题。
从计算角度上，现有添加产状的方法大致可以分为几何方法和插值方法两大类。几何方法通过对位于同一岩层界面的点、线、面等进行几何运算求取产状。三点法(Assaadet al.,2013)是最为常见的几何方法，方法使用同一岩层界面中邻近三个点的高程计算三点所在岩层平面的产状。最初，三点法用于在钻孔密集处利用邻近钻孔求取地下产状(Bucher,1943)。随着DTM精度的提高，Chorowicz等(Chorowicz et al.,1991)、Banerjee等(Banerjee and Mitra,2004)和Yeh等(Yeh et al.,2014)分别采用地质平面图、遥感图像和LiDAR测量结合高精度DTM在地表非覆盖区中使用三点法计算产状。除了三点法之外，也可以从地震信号中通过几何方法求取产状，例如Brown等(Brown et al.,1971)从交叉地震剖面中使用几何方法求取地壳产状，Dalley等(Dalley et al.,2007)通过对三维地震数据进行曲率运算获取了地下的产状场。然而，由于几何方法对资料要求苛刻，导致可添加产状的数目和空间分布范围均有限，例如无法在地表覆盖区添加产状，在实际中稀少的钻孔或地震资料导致可添加的地下产状数目较少。此外，受制于采样分辨率、采样质量等因素，导致几何方法普遍存在精度有限的问题。
基于产状之间空间相关性的插值方法(Chilès and Delfiner,1999)可以在任意位置处添加产状，多用于构建构造可视化和地质随机模拟建模中需用的产状场(Hillieret al.,2013；Xu,1996)。由于无法避免角度数据的周期性引发的算术平均值不连续问题(Grancher et al.,2012)，在对产状数据进行插值时常使用三维单位矢量表达产状(Gumiaux et al.,2003；Upton and Fingleton,1989)，而后利用三维矢量插值方法(如反距离加权(Hillier et al.,2013)、克里金(Young,1987)、径向基插值(Dodu and Rabut,2002)等)进行插值求取产状矢量。尽管插值方法在理论上可以添加任意位置处的产状，但是在构造复杂或实测产状稀疏的区域，插值方法求取的产状难以符合实际地质情况。
发明内容
本发明要解决的技术问题在于，针对现有技术在构造复杂或实测产状稀疏的区域，插值方法求取的产状难以符合实际地质情况的缺陷，提供一种岩层产状的计算方法及系统。
本发明解决其技术问题所采用的技术方案是：构造一种岩层产状的计算方法，该方法包括：
S1、获取实测产状
Figure GDA0003582677010000021
及由实测产状
Figure GDA0003582677010000022
的不确定性所确定的权重系数ωunc,M；获取地质平面图及剖面、地球物理平面图及剖面解释图，并从中提取出地层延伸方向
Figure GDA0003582677010000031
及由地层延伸方向
Figure GDA0003582677010000032
的不确定性所确定的权重系数ωunc,E；
S2、在提取的地层延伸方向
Figure GDA0003582677010000033
处设置初始产状
Figure GDA0003582677010000034
S3、利用插值计算算法，计算实测产状点处的交叉验证产状
Figure GDA0003582677010000035
及地层延伸方向点处初始产状
Figure GDA0003582677010000036
的交叉验证产状
Figure GDA0003582677010000037
S4、使用第一参数角θ对地层延伸方向点处初始产状
Figure GDA0003582677010000038
的交叉验证产状
Figure GDA0003582677010000039
进行表达后，利用投影计算方法，计算第一参数角θ对应的投影产状
Figure GDA00035826770100000310
S5、计算实测产状
Figure GDA00035826770100000311
与实测产状点处的交叉验证产状
Figure GDA00035826770100000312
之间的夹角，得到第一夹角
Figure GDA00035826770100000313
计算地层延伸方向点处的交叉验证产状
Figure GDA00035826770100000314
与其对应的投影产状
Figure GDA00035826770100000315
之间的夹角，得到第二夹角
Figure GDA00035826770100000316
S6、结合第一夹角和第二夹角，构建目标优化函数ObjFunc：
Figure GDA00035826770100000317
S7、针对目标优化函数ObjFunc进行迭代优化过程，每次的迭代优化过程中包括：
进行目标值ObjFunci的计算，ObjFunci为第i次迭代优化过程中计算得到的目标值，i≥1；
设定收敛条件，当第i次迭代优化过程中计算所得目标值ObjFunci满足收敛条件时，则停止迭代，并输出当前的投影产状
Figure GDA00035826770100000318
在目标值ObjFunci不满足收敛条件时，则在更新投影产状
Figure GDA00035826770100000319
后，将地层延伸方向
Figure GDA00035826770100000320
处设置为更新后的投影产状
Figure GDA00035826770100000321
并重复步骤S3-步骤S7；直到满足设定收敛条件时，输出当前的投影产状
Figure GDA00035826770100000322
本发明另一方面公开的一种岩层产状的计算系统，该系统包括：
数据获取单元，用于获取实测产状
Figure GDA00035826770100000323
及由实测产状
Figure GDA00035826770100000324
的不确定性所确定的权重系数ωunc,E；获取地质平面图及剖面、地球物理平面图及剖面解释图，并从中提取出地层延伸方向
Figure GDA0003582677010000041
及由地层延伸方向
Figure GDA0003582677010000042
的不确定性所确定的权重系数ωunc,E；
初始产状设置单元，用于在提取的地层延伸方向
Figure GDA0003582677010000043
处设置初始产状
Figure GDA0003582677010000044
交叉验证产状计算单元，用于利用插值计算算法，计算实测产状点处的交叉验证产状
Figure GDA0003582677010000045
及地层延伸方向点处初始产状
Figure GDA0003582677010000046
的交叉验证产状
Figure GDA0003582677010000047
投影产状计算单元，用于使用第一参数角θ对地层延伸方向点处初始产状
Figure GDA0003582677010000048
的交叉验证产状
Figure GDA0003582677010000049
进行表达后，利用投影计算方法，计算第一参数角θ对应的投影产状
Figure GDA00035826770100000410
夹角计算单元，用于计算实测产状
Figure GDA00035826770100000411
与实测产状点处的交叉验证产状
Figure GDA00035826770100000412
之间的夹角，得到第一夹角
Figure GDA00035826770100000413
计算地层延伸方向点处的交叉验证产状
Figure GDA00035826770100000414
与其对应的投影产状
Figure GDA00035826770100000415
之间的夹角，得到第二夹角
Figure GDA00035826770100000416
目标函数构建单元，用于结合第一夹角和第二夹角，构建目标优化函数ObjFunc：
Figure GDA00035826770100000417
迭代优化单元，用于针对目标优化函数ObjFunc进行迭代优化过程，而，每次的迭代优化过程中包括：
进行目标值ObjFunci的计算，ObjFunci为第i次迭代优化过程中计算得到的目标值，i≥1；设定收敛条件，当第i次迭代优化过程中计算所得目标值ObjFunci满足收敛条件时，则停止迭代，并输出当前的投影产状
Figure GDA00035826770100000418
在目标值ObjFunci不满足收敛条件时，则在更新投影产状
Figure GDA00035826770100000419
后，将地层延伸方向
Figure GDA00035826770100000420
处设置为更新后的投影产状
Figure GDA00035826770100000421
并重复步骤S3-步骤S7；直到满足设定收敛条件时，输出当前的投影产状
Figure GDA00035826770100000422
在本发明所述的一种岩层产状的计算方法及系统中，利用地层延伸方向与产状方向正交这一几何性质，并基于交叉验证方法构建了最优化目标函数，从而实现了自动求取地层延伸方向处产状。且，本发明中为了快速并精确的优化目标函数，采用角度参数化和拟梯度下降法改进了传统非线性优化算法。通过上述方式，本发明能够自动计算大量广泛分布且符合地质构造的产状。
附图说明
下面将结合附图及实施例对本发明作进一步说明，附图中：
图1是本发明一种岩层产状的计算方法的第一个实施例的流程图；
图2是本发明一种岩层产状的计算方法的第二个实施例的流程图；
图3是本发明一种岩层产状的计算系统的第一个实施例的结构示意图；
图4是本发明一种岩层产状的计算系统的第二个实施例的结构示意图；
图5为第一实施例的构造盆地褶皱模拟模型；
图6为第一实施例的产状及地层延伸方向的抽样；
图7是传统粒子群算法目标函数收敛过程图；
图8是本专利方法目标函数收敛过程图；
图9是真实产状与直接矢量插值产状(倾角)对比图；
图10是真实产状与直接矢量插值产状(倾向)对比图；
图11是真实产状与本专利方法求取产状(倾角)对比图；
图12是真实产状与本专利方法求取产状(倾向)对比图。
具体实施方式
为了对本发明的技术特征、目的和效果有更加清楚的理解，现对照附图详细说明本发明的具体实施方式。
实施例1：
请参考图1，其为本发明一种岩层产状的计算方法的第一个实施例的流程图，该方法包括：
S100、获取实测产状
Figure GDA0003582677010000061
及由实测产状
Figure GDA0003582677010000062
的不确定性所确定的权重系数ωunc,M；获取地质平面图及剖面、地球物理平面图及剖面解释图，并从中提取出地层延伸方向
Figure GDA0003582677010000063
及由地层延伸方向
Figure GDA0003582677010000064
的不确定性所确定的权重系数ωunc,E；
当前步骤下，第一方面需要说明的是，实测产状的不确定性可以依据野外采集时的估测、制图标准等进行估计。不确定性对应的权重设定范围为0≤ωunc,M≤1，其中，估计到的不确定性越大，对应的权重就越小。
当前步骤下，第二方面需要说明的是，在提取地层延伸方向时，也可以根据专家知识、数学方法给出其对应的不确定性权重ωunc,E(0≤ωunc,E≤1)。以惯性矩方法为例，可以依据下式定义地层延伸方向
Figure GDA0003582677010000065
的不确定性权重：
Figure GDA0003582677010000066
其中，Ixx和Ixy分别是地层延伸主方向与地层延伸侧方向处的惯性矩阵。
S200、在提取的地层延伸方向
Figure GDA0003582677010000067
处设置初始产状
Figure GDA0003582677010000068
S300、利用插值计算算法，计算实测产状点处的交叉验证产状
Figure GDA0003582677010000069
及地层延伸方向点处初始产状
Figure GDA00035826770100000610
的交叉验证产状
Figure GDA00035826770100000611
具体的，本步骤中，所利用插值计算算法得到实测产状点处的交叉验证产状
Figure GDA00035826770100000612
及地层延伸方向点处的交叉验证产状
Figure GDA00035826770100000613
具体为：
利用公式(1)分别进行实测产状点处以及地层延伸方向点处的交叉验证产状的计算：
Figure GDA00035826770100000614
其中：产状集合
Figure GDA00035826770100000615
为由实测产状
Figure GDA00035826770100000616
和初始产状
Figure GDA00035826770100000617
构成的集合，共M个产状；
Figure GDA00035826770100000618
为留一产状集合，
Figure GDA00035826770100000619
为留一产状集合中的第i个子集，
Figure GDA00035826770100000620
为
Figure GDA00035826770100000621
对应的表示留一产状不确定性的权重系数；interpolater(*)表示“*”进行三维矢量插值的计算。
S400、使用第一参数角θ对地层延伸方向点处初始产状
Figure GDA0003582677010000071
的交叉验证产状
Figure GDA0003582677010000072
进行表达后，利用投影计算方法，计算第一参数角θ对应的投影产状
Figure GDA0003582677010000073
需要说明的是，由于在地层延伸方向点处，地层的延伸方向与产状呈正交关系，因此，本步骤下，使用第一参数角θ对地层延伸方向点处初始产状
Figure GDA0003582677010000074
的交叉验证产状
Figure GDA0003582677010000075
进行表达的表达公式具体由公式(2)所示：
Figure GDA0003582677010000076
其中，θ＝[θ1,θ2,...,θj],j∈[1,2,...,N]，θj为地层延伸方向第j个点处待估的第一参数角，θ为地层延伸方向点处由j个待估第一参数角组成的集合，N为地层延伸方向点个数；
Figure GDA0003582677010000077
为初始产状矢量，
Figure GDA0003582677010000078
为待求产状矢量，
Figure GDA0003582677010000079
为待求产状矢量在z方向上的分量，
Figure GDA00035826770100000710
为待求产状矢量在x方向上的分量；
Figure GDA00035826770100000711
为已知的地层延伸沿三维坐标轴x,y,z三个方向形成的单位矢量。
式(2)中使用参数角θ来替代地层延伸方向点处产状，其具体形式可参考：
Figure GDA00035826770100000712
Figure GDA00035826770100000713
其中，
Figure GDA00035826770100000714
为第j个投影产状，
Figure GDA00035826770100000715
为所求的j个投影产状构成的集合；
Figure GDA00035826770100000716
为初始产状矢量
Figure GDA00035826770100000717
的转置矩阵，
Figure GDA00035826770100000718
为旋转矩阵；
Figure GDA00035826770100000719
由公式(3)定义为：
Figure GDA0003582677010000081
S500、计算实测产状
Figure GDA0003582677010000082
与实测产状点处的交叉验证产状
Figure GDA0003582677010000083
之间的夹角，得到第一夹角
Figure GDA0003582677010000084
计算地层延伸方向点处的交叉验证产状
Figure GDA0003582677010000085
与其对应的投影产状
Figure GDA0003582677010000086
之间的夹角，得到第二夹角
Figure GDA0003582677010000087
具体的，第一夹角
Figure GDA0003582677010000088
具体是通过公式(4)计算得到：
Figure GDA0003582677010000089
其中，arccos(*)为对“*”进行反余弦的计算，“·”为对两个向量进行点乘计算；
具体的，第二夹角
Figure GDA00035826770100000810
具体是通过公式(5)计算得到：
Figure GDA00035826770100000811
S600、结合第一夹角和第二夹角，构建目标优化函数ObjFunc：
Figure GDA00035826770100000812
其中，ωubc,M为由实测产状
Figure GDA00035826770100000813
的不确定性所确定的权重系数。
S700、针对目标优化函数ObjFunc进行迭代优化过程，而，每次的迭代优化过程中包括：
进行目标值ObjFunci的计算，ObjFunci为第i次迭代优化过程中计算得到的目标值，i≥1；
设定收敛条件，当第i次迭代优化过程中计算所得目标值ObjFunci满足收敛条件时，则停止迭代，并输出当前的投影产状
Figure GDA00035826770100000814
在目标值ObjFunci不满足收敛条件时，则在更新投影产状
Figure GDA00035826770100000815
后，将地层延伸方向
Figure GDA00035826770100000816
处设置为更新后的投影产状
Figure GDA0003582677010000091
并重复步骤S3-步骤S7；直到满足设定收敛条件时，输出当前的投影产状
Figure GDA0003582677010000092
需要说明的是，迭代收敛条件可以包括三类：
1、达到迭代次数上限时，停止收敛；
2、目标函数收敛到某一数值以下时，停止收敛；
3、连续一定次数目标函数减小幅值小于某一数值时，停止收敛。
重复步骤S3-步骤S7的过程为拟梯度下降法，最终输出结果即为步骤S400中求取到的投影产状
Figure GDA0003582677010000093
实施例2：
为了保证输出数据的准确性，请参考图2，其为本发明一种岩层产状的计算方法的第二个实施例的流程图，结合实施例1，在步骤S700之后还包括：
S800、基于输出投影产状的精度值，判断是否采用步骤S700的投影产状输出结果。
具体的，步骤S800下包括以下子步骤：
S810、当所求岩层投影产状的精度满足预置要求时，则采用步骤S700的输出结果；不满足的情况下则执行步骤S820；
S820、采用传统非线性算法，以θ为待优化变量，以步骤S700中求得的投影产状
Figure GDA0003582677010000094
对应的参数角为初始解，优化步骤S600中构建的目标优化函数ObjFunc；非线性优化算法完成后，输出当前的投影产状
Figure GDA0003582677010000095
需要说明的是，步骤S820中采用的传统非线性算法包括粒子群算法和遗传算法。
结合实施例1-2，在本发明所述的一种岩层产状的计算方法中，利用地层延伸方向与产状方向正交这一几何性质，并基于交叉验证方法构建了最优化目标函数，从而实现了自动求取地层延伸方向处产状。且，本发明中为了快速并精确的优化目标函数，采用角度参数化和拟梯度下降法改进了传统非线性优化算法。通过上述方式，本发明能够自动计算大量广泛分布且符合地质构造的产状。
实施例3：
请参考图3，其为本发明一种岩层产状的计算系统的第一个实施例的结构示意图；该系统包括：
数据获取单元10，用于获取实测产状
Figure GDA0003582677010000101
及由实测产状
Figure GDA0003582677010000102
的不确定性所确定的权重系数ωunc,M；获取地质平面图及剖面、地球物理平面图及剖面解释图，并从中提取出地层延伸方向
Figure GDA0003582677010000103
及由地层延伸方向
Figure GDA0003582677010000104
的不确定性所确定的权重系数ωunc,E；
初始产状设置单元20，用于在提取的地层延伸方向
Figure GDA0003582677010000105
处设置初始产状
Figure GDA0003582677010000106
交叉验证产状计算单元30，用于利用插值计算算法，计算实测产状点处的交叉验证产状
Figure GDA0003582677010000107
及地层延伸方向点处初始产状
Figure GDA0003582677010000108
的交叉验证产状
Figure GDA0003582677010000109
投影产状计算单元40，用于使用第一参数角θ对地层延伸方向点处初始产状
Figure GDA00035826770100001010
的交叉验证产状
Figure GDA00035826770100001011
进行表达后，利用投影计算方法，计算第一参数角θ对应的投影产状
Figure GDA00035826770100001012
夹角计算单元50，用于计算实测产状
Figure GDA00035826770100001013
与实测产状点处的交叉验证产状
Figure GDA00035826770100001014
之间的夹角，得到第一夹角
Figure GDA00035826770100001015
计算地层延伸方向点处的交叉验证产状
Figure GDA00035826770100001016
与其对应的投影产状
Figure GDA00035826770100001017
之间的夹角，得到第二夹角
Figure GDA00035826770100001018
目标函数构建单元60，用于结合第一夹角和第二夹角，构建目标优化函数ObjFunc：
Figure GDA00035826770100001019
迭代优化单元70，用于针对目标优化函数ObjFunc进行迭代优化过程，而，每次的迭代优化过程中包括：
进行目标值ObjFunci的计算，ObjFunci为第i次迭代优化过程中计算得到的目标值，i≥1；
设定收敛条件，当第i次迭代优化过程中计算所得目标值ObjFunci满足收敛条件时，则停止迭代，并输出当前的投影产状
Figure GDA0003582677010000111
在目标值ObjFunci不满足收敛条件时，则在更新投影产状
Figure GDA0003582677010000112
后，将地层延伸方向
Figure GDA0003582677010000113
处设置为更新后的投影产状
Figure GDA0003582677010000114
并重复步骤S3-步骤S7；直到满足设定收敛条件时，输出当前的投影产状
Figure GDA0003582677010000115
实施例4：
请参考图4，其为本发明一种岩层产状的计算系统的第二个实施例的结构示意图，结合实施例3，为了提高输出数据的准确性，该系统还包括：
投影产状精度判断单元80，用于基于输出投影产状的精度值，判断是否采用迭代优化单元的投影产状输出结果。其中，所述投影产状精度判断单元80下包括：
投影产状精度比较单元81，用于当所求岩层投影产状的精度满足预置要求时，则采用迭代优化单元的输出结果；不满足的情况下则执行传统非线性算法运用单元；
传统非线性算法运用单元82，用于采用传统非线性算法，以θ为待优化变量，以迭代优化单元中求得的投影产状
Figure GDA0003582677010000116
对应的参数角为初始解，优化目标函数构建单元中构建的目标优化函数ObjFunc；非线性优化算法完成后，输出输出当前的投影产状
Figure GDA0003582677010000117
本发明较原有技术有较大改进，使其完全符合地质勘探科研及生产的需要并产生了好用及实用的效果。使用该方法及系统增加了添加产状的数目和分布范围，添加的产状符合地质构造，以添加产状为基础可以提高地质解释和地质建模的精度，节省了勘查、勘探成本，降低了勘探风险。
实施例采用的模拟模型为一个构造盆地褶皱的半球形模型，模拟的“沉积”时间函数为
Figure GDA0003582677010000121
模型范围设置及模拟的3个地层界面(等时界面)见图5，在模型空间中有一条垂直剖面AA′。为了模拟真实地质场景，在z＝0m的“地表”上进行随机采样获取10个产状点
Figure GDA0003582677010000122
产状点处的产状矢量
Figure GDA0003582677010000123
由式
Figure GDA0003582677010000124
求取，同时在模型空间中随机采样获得50个地层延伸方向点
Figure GDA0003582677010000125
(图6)。
在角度参数化基础上，分别利用传统粒子群优化(PSO)方法和本专利方法改进的PSO方法求取地层延伸方向处产状，目标函数的收敛过程分别见图7和图8。
由图7和图8可知，其横坐标指的是迭代次数，纵坐标指的是在相应的迭代次数下目标函数的输出值。由图7和图8可知，在非线性算法迭代1000次后，传统粒子群优化(PSO)方法和本专利方法改进的PSO方法求取的最终结果相近，但以本申请公开的方法最优，原因是：本申请方法提前300次迭代就达到了PSO算法的最终优化结果，显然本申请公开的方法不仅提高了PSO算法的收敛速度，也抑制了PSO算法的早熟收敛问题。且，本申请中，在采用了拟梯度下降法进行迭代优化后，在迭代步数达到6步后，目标函数就由最初的12.08°快速的减小到8.01°，在迭代步数达到16步后，目标函数就由最初的12.08°收敛到7.91°，该方法可以快速收敛到近于最优的结果。然而，采用传统粒子群优化算法需要228次迭代才可以收敛到7.91°(图5)，其计算消耗为拟梯度下降法的228*100/16＝1425倍。
请参考图9-图12，横纵坐标分别为真实产状与直接矢量插值产状(倾角)对比图、真实产状与直接矢量插值产状(倾向)对比图、真实产状与本专利方法求取产状(倾角)对比图、真实产状与本专利方法求取产状(倾向)对比图。由图9和图10可知，在地层延伸方向处利用直接矢量插值得到的产状与真实产状之间的平均角度偏差为24.3°，两种产状之间的相关性很小，尤其是倾角的Pearson相关性弱至0.3327。而，由图11和图12可知，利用本申请公开的方法求取的产状与真实产状之间的平均角度偏差为9.3°，计算得到的两种产状之间为强线性相关，倾角和倾向的Pearson相关性分别为0.99和0.98。显示出本实施例中真实产状与本发明计算的产状相近，明显提高了直接插值产状的精度。
上面结合附图对本发明的实施例进行了描述，但是本发明并不局限于上述的具体实施方式，上述的具体实施方式仅仅是示意性的，而不是限制性的，本领域的普通技术人员在本发明的启示下，在不脱离本发明宗旨和权利要求所保护的范围情况下，还可做出很多形式，这些均属于本发明的保护之内。

Claims (11)

1.一种岩层产状的计算方法，其特征在于，该方法包括：
S1、将岩层产状用地层面法向单位矢量表达；获取实测产状
Figure FDA0003590471560000011
及由实测产状
Figure FDA0003590471560000012
的不确定性所确定的权重系数ωunc,M；获取地质平面图及剖面、地球物理平面图及剖面解释图，并从中提取出地层延伸方向
Figure FDA0003590471560000013
及由地层延伸方向
Figure FDA0003590471560000014
的不确定性所确定的权重系数ωunc,E；
S2、在提取的地层延伸方向
Figure FDA0003590471560000015
处设置初始产状
Figure FDA0003590471560000016
S3、利用插值计算算法，计算实测产状点处的交叉验证产状
Figure FDA0003590471560000017
及地层延伸方向点处初始产状
Figure FDA0003590471560000018
的交叉验证产状
Figure FDA0003590471560000019
S4、使用第一参数角θ对地层延伸方向点处初始产状
Figure FDA00035904715600000110
的交叉验证产状
Figure FDA00035904715600000111
进行表达后，利用投影计算方法，计算第一参数角θ对应的投影产状
Figure FDA00035904715600000112
S5、计算实测产状
Figure FDA00035904715600000113
与实测产状点处的交叉验证产状
Figure FDA00035904715600000114
之间的夹角，得到第一夹角
Figure FDA00035904715600000115
计算地层延伸方向点处的交叉验证产状
Figure FDA00035904715600000116
与其对应的投影产状
Figure FDA00035904715600000117
之间的夹角，得到第二夹角
Figure FDA00035904715600000118
S6、结合第一夹角和第二夹角，构建目标优化函数ObjFunc：
Figure FDA00035904715600000119
S7、针对目标优化函数ObjFunc进行拟梯度下降迭代优化过程，每次的迭代优化过程中包括：
进行目标值ObjFunci的计算，ObjFunci为第i次迭代优化过程中计算得到的目标值，i≥1；
设定收敛条件，当第i次迭代优化过程中计算所得目标值ObjFunci满足收敛条件时，则停止迭代，并输出当前的投影产状
Figure FDA00035904715600000120
在目标值ObjFunci不满足收敛条件时，则在更新投影产状
Figure FDA0003590471560000021
后，将地层延伸方向
Figure FDA0003590471560000022
处设置为更新后的投影产状
Figure FDA0003590471560000023
并重复步骤S3-步骤S7；直到满足设定收敛条件时，输出当前的投影产状
Figure FDA0003590471560000024
2.根据权利要求1所述的一种岩层产状的计算方法，其特征在于，步骤S3中，所述利用插值计算算法，得到实测产状点处的交叉验证产状
Figure FDA0003590471560000025
及地层延伸方向点处的交叉验证产状
Figure FDA0003590471560000026
具体为：
利用公式(1)分别进行实测产状点处以及地层延伸方向点处的交叉验证产状的计算：
Figure FDA0003590471560000027
其中：产状集合
Figure FDA0003590471560000028
为由实测产状
Figure FDA0003590471560000029
和初始产状
Figure FDA00035904715600000210
构成的集合，共M个产状；
Figure FDA00035904715600000211
为留一产状集合，
Figure FDA00035904715600000212
Figure FDA00035904715600000213
为留一产状集合中的第i个子集，
Figure FDA00035904715600000214
为
Figure FDA00035904715600000215
对应的表示留一产状不确定性的权重系数；interpolater(*)表示“*”进行三维矢量插值的计算。
3.根据权利要求1所述的一种岩层产状的计算方法，其特征在于，步骤S4，所述使用第一参数角θ对地层延伸方向点处初始产状
Figure FDA00035904715600000216
的交叉验证产状
Figure FDA00035904715600000217
进行表达的表达公式具体由公式(2)所示：
Figure FDA00035904715600000218
其中，θ＝[θ1,θ2,...,θj],j∈[1,2,…,N]，θj为地层延伸方向第j个点处待估的第一参数角，θ为地层延伸方向点处由j个待估第一参数角组成的集合，N为地层延伸方向点总个数；
Figure FDA0003590471560000031
为初始产状矢量，
Figure FDA0003590471560000032
为待求产状矢量，
Figure FDA0003590471560000033
为待求产状矢量在z方向上的分量，
Figure FDA0003590471560000034
为待求产状矢量在x方向上的分量；
Figure FDA0003590471560000035
为已知的地层延伸沿三维坐标轴x,y,z三个方向形成的单位矢量。
4.根据权利要求3所述的一种岩层产状的计算方法，其特征在于，步骤S4中，所述利用投影计算方法，计算第一参数角θ对应的投影产状
Figure FDA0003590471560000036
具体为：
S41、基于θj，利用投影计算方法，计算不同角度参数θj下对应的投影产状
Figure FDA0003590471560000037
计算公式为：
Figure FDA0003590471560000038
其中，
Figure FDA0003590471560000039
Figure FDA00035904715600000310
为第j个投影产状，
Figure FDA00035904715600000311
为所求的j个投影产状构成的集合；
Figure FDA00035904715600000312
为初始产状矢量
Figure FDA00035904715600000313
的转置矩阵，
Figure FDA00035904715600000314
为旋转矩阵；
Figure FDA00035904715600000315
由公式(3)定义为：
Figure FDA00035904715600000316
5.根据权利要求1所述的一种岩层产状的计算方法，其特征在于，步骤S5中，第一夹角
Figure FDA00035904715600000317
具体是通过公式(4)计算得到：
Figure FDA00035904715600000318
其中，arccos(*)为对“*”进行反余弦的计算，“.”为对两个向量进行点乘计算；
所述第二夹角
Figure FDA00035904715600000319
具体是通过公式(5)计算得到：
Figure FDA00035904715600000320
6.根据权利要求1所述的一种岩层产状的计算方法，其特征在于，还包括以下步骤：
S8、基于输出投影产状的精度值，判断是否采用步骤S7的投影产状输出结果。
7.根据权利要求6所述的一种岩层产状的计算方法，其特征在于，步骤S8中，所述基于输出投影产状的精度值，判断是否采用步骤S7的投影产状输出结果具体为：
S81、当所求岩层投影产状的精度满足预置要求时，则采用步骤S7的输出结果；不满足的情况下则执行步骤S82；
S82、采用传统非线性算法，以θ为待优化变量，以步骤S7中求得的投影产状
Figure FDA0003590471560000041
对应的参数角为初始解，优化步骤S6中构建的目标优化函数ObjFunc；非线性优化算法完成后，输出当前的投影产状
Figure FDA0003590471560000042
8.根据权利要求7所述的一种岩层产状的计算方法，其特征在于，所述传统非线性算法包括粒子群算法和遗传算法。
9.一种岩层产状的计算系统，其特征在于，该系统包括：
数据获取单元，用于获取实测产状
Figure FDA0003590471560000043
及由实测产状
Figure FDA0003590471560000044
的不确定性所确定的权重系数ωunc,E；获取地质平面图及剖面、地球物理平面图及剖面解释图，并从中提取出地层延伸方向
Figure FDA0003590471560000045
及由地层延伸方向
Figure FDA0003590471560000046
的不确定性所确定的权重系数ωunc,E；
初始产状设置单元，用于在提取的地层延伸方向
Figure FDA0003590471560000047
处设置初始产状
Figure FDA0003590471560000048
交叉验证产状计算单元，用于利用插值计算算法，计算实测产状点处的交叉验证产状
Figure FDA0003590471560000049
及地层延伸方向点处初始产状
Figure FDA00035904715600000410
的交叉验证产状
Figure FDA00035904715600000411
投影产状计算单元，用于使用第一参数角θ对地层延伸方向点处初始产状
Figure FDA0003590471560000051
的交叉验证产状
Figure FDA0003590471560000052
进行表达后，利用投影计算方法，计算第一参数角θ对应的投影产状
Figure FDA0003590471560000053
夹角计算单元，用于计算实测产状
Figure FDA0003590471560000054
与实测产状点处的交叉验证产状
Figure FDA0003590471560000055
之间的夹角，得到第一夹角
Figure FDA0003590471560000056
计算地层延伸方向点处的交叉验证产状
Figure FDA0003590471560000057
与其对应的投影产状
Figure FDA0003590471560000058
之间的夹角，得到第二夹角
Figure FDA0003590471560000059
目标函数构建单元，用于结合第一夹角和第二夹角，构建目标优化函数ObjFunc：
Figure FDA00035904715600000510
迭代优化单元，用于针对目标优化函数ObjFunc进行迭代优化过程，每次的迭代优化过程中包括：
进行目标值ObjFunci的计算，ObjFunci为第i次迭代优化过程中计算得到的目标值，i≥1；
设定收敛条件，当第i次迭代优化过程中计算所得目标值ObjFunci满足收敛条件时，则停止迭代，并输出当前的投影产状
Figure FDA00035904715600000511
在目标值ObjFunci不满足收敛条件时，则在更新投影产状
Figure FDA00035904715600000512
后，将地层延伸方向
Figure FDA00035904715600000513
处设置为更新后的投影产状
Figure FDA00035904715600000514
并重复步骤S3-步骤S7；直到满足设定收敛条件时，输出当前的投影产状
Figure FDA00035904715600000515
10.根据权利要求9述的一种岩层产状的计算系统，其特征在于，该系统还包括：
投影产状精度判断单元，用于基于输出投影产状的精度值，判断是否采用迭代优化单元的投影产状输出结果。
11.根据权利要求10所述的一种岩层产状的计算系统，其特征在于，所述投影产状精度判断单元下包括：
投影产状精度比较单元，用于当所求岩层投影产状的精度满足预置要求时，则采用迭代优化单元的输出结果；不满足的情况下则执行传统非线性算法运用单元；
传统非线性算法运用单元，用于采用传统非线性算法，以θ为待优化变量，以迭代优化单元中求得的投影产状
Figure FDA0003590471560000061
对应的参数角为初始解，优化目标函数构建单元中构建的目标优化函数ObjFunc；非线性优化算法完成后，输出当前的投影产状
Figure FDA0003590471560000062
CN202010848655.XA 2020-08-21 2020-08-21 一种岩层产状的计算方法及系统 Active CN112270064B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010848655.XA CN112270064B (zh) 2020-08-21 2020-08-21 一种岩层产状的计算方法及系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010848655.XA CN112270064B (zh) 2020-08-21 2020-08-21 一种岩层产状的计算方法及系统

Publications (2)

Publication Number Publication Date
CN112270064A CN112270064A (zh) 2021-01-26
CN112270064B true CN112270064B (zh) 2022-07-26

Family

ID=74348814

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010848655.XA Active CN112270064B (zh) 2020-08-21 2020-08-21 一种岩层产状的计算方法及系统

Country Status (1)

Country Link
CN (1) CN112270064B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114581556B (zh) * 2022-03-10 2022-12-27 青海省地质调查院 一种区域地质调查中的数字填图方法

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107045154A (zh) * 2017-02-08 2017-08-15 中国海洋石油总公司 一种水平井环境中的识别地层产状的方法和装置

Family Cites Families (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8275593B2 (en) * 2009-07-16 2012-09-25 University Of Regina Reservoir modeling method
CN103148833B (zh) * 2013-02-01 2015-03-04 淮南矿业(集团)有限责任公司 岩层产状参数获取和计算方法
CN104200039B (zh) * 2014-09-17 2017-12-29 中国石油大学(华东) 一种构造裂缝产状定量预测方法
CN105426612B (zh) * 2015-11-18 2019-01-18 中国石油天然气股份有限公司 一种地层组分最优化确定方法及装置
CN109961513A (zh) * 2019-03-22 2019-07-02 长江岩土工程总公司（武汉） 利用三维模型快速求解岩层产状的方法
CN110826215B (zh) * 2019-10-31 2024-05-10 中国地质大学（武汉） 实现高精度产状分布估计的最小夹角和最小样本容量算法

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107045154A (zh) * 2017-02-08 2017-08-15 中国海洋石油总公司 一种水平井环境中的识别地层产状的方法和装置

Also Published As

Publication number Publication date
CN112270064A (zh) 2021-01-26

Similar Documents

Publication Publication Date Title
CN108072892B (zh) 一种自动化的地质构造约束层析反演方法
CN105277978B (zh) 一种确定近地表速度模型的方法及装置
AU2008239658B2 (en) Inverse-vector method for smoothing dips and azimuths
CN106646645B (zh) 一种重力正演加速方法
CN106483559B (zh) 一种地下速度模型的构建方法
EA018815B1 (ru) Формирование геологической модели
CN105353405B (zh) 一种全波形反演方法和系统
CN111305834B (zh) 基于多探测模式电阻率测井的三维反演初始模型构建方法
CN109188519B (zh) 一种极坐标下的弹性波纵横波速度反演系统及方法
CN110187382B (zh) 一种回折波和反射波波动方程旅行时反演方法
CN110569624A (zh) 适用于PolInSAR反演的森林三层散射模型的确定及分析方法
CN110133644B (zh) 基于插值尺度函数法的探地雷达三维正演方法
CN112270064B (zh) 一种岩层产状的计算方法及系统
CN114460654A (zh) 基于l1l2混合范数的半航空瞬变电磁数据反演方法及装置
CN108303736B (zh) 各向异性ti介质最短路径射线追踪正演方法
CA2226090C (en) Trendform gridding method using distance transformations
CN109444973B (zh) 一种球坐标系下重力正演加速方法
CN113671570B (zh) 一种地震面波走时和重力异常联合反演方法与系统
CN113466933B (zh) 基于深度加权的地震斜率层析成像方法
CN111158059B (zh) 基于三次b样条函数的重力反演方法
CN111980662B (zh) 一种斜井各向异性地层阵列侧向测井资料快速处理方法
CN115730424B (zh) 基于多源大地测量数据的有限断层反演方法、装置及终端
CN116449305A (zh) 基于可控变分自编码器的稠密时变阵列构建方法及系统
CN116187168A (zh) 基于神经网络-重力信息小波分解提高水深反演精度方法
CN110794469B (zh) 基于最小地质特征单元约束的重力反演方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant