首页 > 专利 > 桂林理工大学 > 一种二度体重力异常计算方法专利详情

一种二度体重力异常计算方法 0 0

有效专利查看PDF

申请进展

基本信息

申请人信息

代理人信息

摘要

法律状态

权利要求

说明书

专利申请流程有哪些步骤？

申请

申请号：指国家知识产权局受理一件专利申请时给予该专利申请的一个标示号码。唯一性原则。

申请日：提出专利申请之日。

2017-01-10

申请公布

申请公布指发明专利申请经初步审查合格后，自申请日（或优先权日）起18个月期满时的公布或根据申请人的请求提前进行的公布。

申请公布号：专利申请过程中，在尚未取得专利授权之前，国家专利局《专利公报》公开专利时的编号。

申请公布日：申请公开的日期，即在专利公报上予以公开的日期。

2017-07-11

授权

授权指对发明专利申请经实质审查没有发现驳回理由，授予发明专利权；或对实用新型或外观设计专利申请经初步审查没有发现驳回理由，授予实用新型专利权或外观设计专利权。

2019-10-15

预估到期

发明专利权的期限为二十年，实用新型专利权期限为十年，外观设计专利权期限为十五年，均自申请日起计算。专利届满后法律终止保护。

2037-01-10

基本信息

有效性	有效专利	专利类型	发明专利
申请号	CN201710017147.5	申请日	2017-01-10
公开/公告号	CN106855904B	公开/公告日	2019-10-15
授权日	2019-10-15	预估到期日	2037-01-10
申请年	2017年	公开/公告年	2019年
缴费截止日
分类号	G06F17/50	主分类号	G06F17/50
是否联合申请	独立申请	文献类型号	B
独权数量	1	从权数量	0
权利要求数量	1	非专利引证数量	1
引用专利数量	4	被引证专利数量	0
非专利引证	1、赵东东等.基于高斯牛顿法的二维直流电阻率法的快速反演《.中国有色金属学报》.2015,第1662-1671页. 朱自强等.二度体的重力张量有限元正演模拟《.物探与化探》.2010,第668-671,685页.;
引用专利	CN104750983A、CN105549106A、CN105842745A、WO2012135020A2	被引证专利
专利权维持	6	专利申请国编码	CN
专利事件		事务标签	公开、实质审查、授权

申请人信息

申请人	桂林理工大学	第一申请人	桂林理工大学
专利权人	桂林理工大学	当前专利权人	桂林理工大学
发明人	陈龙伟、熊彬、陈欣、罗天涯	第一发明人	陈龙伟
地址	广西壮族自治区桂林市雁山街319号桂林理工大学逸夫楼503	邮编	541006
申请人数量	1	发明人数量	4
申请人所在省	广西壮族自治区	申请人所在市	广西壮族自治区桂林市

代理人信息

代理机构

专利代理机构是经省专利管理局审核，国家知识产权局批准设立，可以接受委托人的委托，在委托权限范围内以委托人的名义办理专利申请或其他专利事务的服务机构。

北京中济纬天专利代理有限公司

代理人

专利代理师是代理他人进行专利申请和办理其他专利事务，取得一定资格的人。

陆薇薇

摘要

本发明提出了一种二度体重力异常计算方法，本发明通过复杂二度体模型表示和矩形模型重力异常计算步骤，其中矩形模型重力异常计算包括加权系数计算、一维离散卷积计算、重力异常值合成步骤。本发明实现了重力异常计算在效率和精度上的统一。本发明解决了现有二度体重力异常计算方法不能同时保证计算效率和计算精度，无法满足大规模重力异常数据密度精细反演、人机交互建模和解释的需求的问题。

摘要附图
说明书附图：图1
说明书附图：图2
说明书附图：图3
说明书附图：图4
说明书附图：图5

法律状态

序号	法律状态公告日	法律状态	法律状态信息
1	2019-10-15	授权
2	2017-07-11	实质审查的生效	IPC(主分类): G06F 17/50 专利申请号: 201710017147.5 申请日: 2017.01.10
3	2017-06-16	公开

权利要求

权利要求书是申请文件最核心的部分，是申请人向国家申请保护他的发明创造及划定保护范围的文件。

1.一种二度体重力异常计算方法，其特征在于，包括以下步骤：
第一步复杂二度体模型表示：
确定观测点，观测点坐标(xm,z0)，基于观测点建立包含所有目标区域的矩形模型，确定矩形模型在x，z方向的起始位置，使得包含起伏地形的目标区域完全嵌入在该矩形模型中；
将该矩形模型均匀划分成若干个小矩形，确定小矩形的几何尺寸Δx，Δz；
最后，根据目标区域的密度分布，对每个小矩形密度进行赋值，每个小矩形密度为常值，不同小矩形密度取值不同，以此刻画任意截面形状、任意密度分布的二度体；将位于空气部分的小矩形的密度值设为零，以此刻画起伏地形；
第二步矩形模型重力异常计算；
第一步中给出的矩形模型，其重力异常计算公式为
式中，(xm,z0)表示观测点坐标，z0为常值；L表示矩形模型在z方向剖分小矩形的个数；M表示矩形模型在x方向剖分小矩形的个数；(ξp,ζr)表示编号为(p,r)的小矩形几何中心坐标；ρ(ξp,ζr)表示编号为(p,r)的小矩形的密度值；h(xm-ξp,z0-ζr)表示加权系数；所述矩形模型其重力异常计算公式的解算方法如下：
a，根据观测点坐标(xm,z0)和小矩形几何中心坐标(ξp,ζr)，计算加权系数h(xm-ξp,z0-ζr)，其计算公式为
式中，γ表示万有引力常数，arctan()表示反余切函数运算符，ln()表示自然对数运算符；其它符号含义如下
X1＝ξp-0.5Δx-xm，X2＝ξp+0.5Δx-xm，Z1＝ζr-0.5Δz-z0，Z2＝ζr+0.5Δz-z0Δx，Δz表示小矩形几何尺寸；
b，在矩形模型中在x方向上位于同一行的所有小矩形为一层，采用一维离散卷积快速计算方法来计算一层矩形模型重力异常，其计算公式为
式中，表示第r层矩形模型在测线z0产生的重力异常，其中r＝1,2,…,L；(xm,z0)表示观测点坐标；
其中的计算方法如下：
(1)将加权系数h(x1-ξp,z0-ζr)排列成向量t，记为
t＝[t0,t1,t2,…,tM-1,0,tM-1,tM-2,…,t2,t1]T         (5)
式中，矩阵元素ti与加权系数h(x1-ξp,z0-ζr)存在关系
ti＝h(x1-ξi+1,z0-ζr)                 (6)
(2)将第r层密度值ρ(ξp,ζr)排列成向量ρ，其中p＝1,2,…,M，向量元素ρi与密度值存在关系
ρi＝ρ(ξi,ζr)               (7)
将向量ρ补零扩展成向量ρext
式中，0M×1表示M×1零向量；
(3)计算
式中，fft()表示一维快速傅里叶变换；
(4)计算
式中，“.*”表示对应元素相乘运算；
(5)计算
式中，ifft()表示一维快速傅里叶反变换；
ext
(6)提取矩阵g 的前M行元素，构成向量g，即为一维离散卷积计算结果，向量g中元素就是所求
c，将各层矩形模型重力异常进行累加，得到整个矩形模型的重力异常，即

说明书

技术领域

[0001] 本发明涉及一种二度体重力异常计算方法，特别是任意截面形状、任意密度分布的二度体重力异常快速、高精度计算方法。

背景技术

[0002] 重力勘探以地球重力场为根据、密度差异为物理基础，是一种研究地球构造及寻找矿产资源的地球物理勘探方法，该方法效率高、操作简单、成本低、勘探深度大、实施过程没有过多条件限制。重力异常正演计算是指根据密度分布计算重力异常，反演是指根据观测重力值计算密度分布。正演是反演的基础，正演计算的效率直接影响反演计算的效率，而正演计算精度直接影响反演结果的质量。一直以来，学者们都非常重视重力异常的正演计算。实际情况中，常见这样一类线性地质体：其走向方向的尺度要远比垂直其走向方向的尺度大，诸如此类地质体，实际场源分布便可以用走向方向无限延伸的二度体代替。

[0003] 伴随测绘技术和仪器的发展，重力测量在测量精度、空间分辨率和测量范围上都有了显著提高，为重力反演提供了大规模高精度、高分辨率数据。同时，随着计算机软硬件水平的不断提高，人机交互建模、解释也日益得到人们的重视，在地球物理勘探中发挥着越来越重要作用。人机交互建模能够使人们通过直观的方式对地质体进行建模，更容易融合地下结构的先验信息。反演方法与人机交互建模、解释方法相辅相成，将极大提高人们对地球内部结构的认识。

[0004] 对于任意截面形状、任意密度分布的二度体重力异常计算，一般采用数值方法。文献(张岭,郝天珧.基于Delaunay剖分的二维非规则重力建模及重力计算.地球物理学报,2006,49(3):877-884.)公开了一种Delaunay剖分方式，将截面分为若干三角形，将二度体分解为若干三角棱柱的组合，利用解析公式计算变密度三角棱柱重力异常，并将其累加，这种方法保证了计算精度，但计算效率低；文献(朱自强,曾思红,鲁光银等.二度体的重力张量有限元正演模拟.物探与化探,2010,34(5):668-671.)公开了一种计算二度体重力梯度张量的有限元方法；文献(Reeder K,Louie J,Kent G.Efficient 2D finite element gravity modeling using convolution.SEG Technical Program Expanded Abstracts,
2014:1254-1258.)公开了一种有限元方法和卷积算法相结合的二度体重力异常计算方法。
采用有限元方法计算二度体重力异常，能够较精确刻画复杂二度体，计算精度高但计算效率很低。

[0005] 剖分方式和计算方法共同决定了重力异常计算的效率和精度。计算效率和计算精度是一对矛盾体，目前已有的方法存在的最大问题是不能同时保证计算效率和计算精度，无法满足大规模重力异常数据密度精细反演、人机交互建模和解释的需求。因此，寻找一种计算效率高、同时能保证计算精度的计算方法具有重要的现实意义。

发明内容

[0006] 本发明针对目前现有二度体重力异常计算方法不能同时保证计算效率和计算精度，无法满足大规模重力异常密度精细反演成像、人机交互建模和解释的需求问题，本发明其目的在于提出一种二度体重力异常计算方法，其是一种任意截面形状、任意密度分布的二度体重力异常快速、高精度计算方法。

[0007] 本发明的技术方案是：

[0008] 一种二度体重力异常计算方法，包括以下步骤：

[0009] 第一步复杂二度体模型表示：

[0010] 确定观测点，观测点坐标(xm,z0)，基于观测点建立包含所有目标区域的矩形模型，确定矩形模型在x，z方向的起始位置，使得包含起伏地形的目标区域完全嵌入在该矩形模型中；

[0011] 然后将该矩形模型均匀划分成若干个规则小矩形，确定小矩形的几何尺寸Δx，Δz；

[0012] 最后根据目标区域的密度分布，对每个小矩形密度进行赋值，每个小矩形密度为常值，不同矩形密度取值不同，以此刻画任意截面形状、任意密度分布的二度体；将位于空气部分的小矩形的密度值设为零，以此刻画起伏地形；

[0013] 第二步矩形模型重力异常计算；

[0014] 第一步中给出的矩形模型，其重力异常计算公式为

[0015]

[0016] 式中，(xm,z0)表示观测点坐标，z0为常值；L表示矩形模型在z方向剖分小矩形的个数；M表示矩形模型在x方向剖分小矩形的个数；(ξp,ζr)表示编号为(p,r)的小矩形几何中心坐标即在矩形模型中x方向为第p个且在z方向为第r个的小矩形的几何中心坐标；ρ(ξp,ζr)表示编号为(p,r)的小矩形的密度值；h(xm-ξp,z0-ζr)表示加权系数。

[0017] 进一步地，第二步中，矩形模型其重力异常计算公式的解算方法如下：

[0018] a，根据观测点坐标(xm,z0)和小矩形几何中心坐标(ξp,ζr)，计算加权系数h(xm-ξp,z0-ζr)，其计算公式为

[0019]

[0020] 式中，γ表示万有引力常数，arctan()表示反余切函数运算符，ln()表示自然对数运算符；其它符号含义如下

[0021] X1＝ξp-0.5Δx-xm，X2＝ξp+0.5Δx-xm，Z1＝ζr-0.5Δz-z0，Z2＝ζr+0.5Δz-z0[0022] Δx，Δz表示小矩形几何尺寸。

[0023] b，在矩形模型中在x方向上位于同一行的所有小矩形为一层。采用一维离散卷积快速计算方法来计算一层矩形模型重力异常，其计算公式为

[0024]

[0025] 式中，表示第r层(r＝1,2,…,L)矩形模型在测线z0产生的重力异常；(xm,z0)表示观测点坐标；

[0026] c，将各层矩形模型重力异常进行累加，得到整个矩形模型的重力异常，即[0027]

[0028] 在第二步中的步骤b中，采用一维离散卷积快速计算方法来计算一层矩形模型重力异常，其步骤为：

[0029] (1)将加权系数h(x1-ξp,z0-ζr)排列成向量t，记为

[0030] t＝[t0,t1,t2,…,tM-1,0,tM-1,tM-2,…,t2,t1]T (5)

[0031] 式中，矩阵元素ti与加权系数h(x1-ξp,z0-ζr)存在关系

[0032] ti＝h(x1-ξi+1,z0-ζr) (6)

[0033] (2)将第r层密度值ρ(ξp,ζr)(p＝1,2,…,M)排列成向量ρ，向量元素ρi与密度值存在关系

[0034] ρi＝ρ(ξi,ζr) (7)

[0035] 将向量ρ补零扩展成向量ρext

[0036]

[0037] 式中，0M×1表示M×1零向量；

[0038] (3)计算

[0039] 式中，fft()表示一维快速傅里叶变换；

[0040] (4)计算

[0041] 式中，“.*”表示对应元素相乘运算；

[0042] (5)计算

[0043] 式中，ifft()表示一维快速傅里叶反变换；

[0044] (6)提取矩阵gext的前M行元素，构成向量g，即为一维离散卷积计算结果，向量g中元素就是所求

[0045] 与现有技术相比，本发明具有以下优点：

[0046] (1)模型表示方法简单、灵活，很容易刻画任意截面形状、任意密度分布复杂二度体以及起伏地形；

[0047] (2)能够实现任意截面形状、任意密度情况下复杂二度体重力异常的快速、高精度计算，可以满足大规模重力密度精细反演、人机交互建模和解释的需求；

[0048] (3)大规模正演计算时，算法不但计算效率和计算精度高，并且所需计算机内存小。

实施方案

[0056] 为使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明实施方式作进一步地详细描述。

[0057] 参照图1，为本实施例一种二度体重力异常计算方法的流程图，包括以下步骤：

[0058] 1、复杂二度体模型表示：

[0059] 首先，确定观测点，观测点坐标(xm,z0)，基于观测点建立包含所有目标区域的规则矩形模型，确定矩形模型在x，z方向的起始位置，使得目标区域(包含起伏地形)完全嵌入在该矩形模型中；

[0060] 其次，根据实际问题需求即根据预定的针对目标区域确定的小矩形剖分个数，将矩形模型均匀划分成多个规则小矩形(如图2所示)，确定小矩形的几何尺寸Δx，Δz；

[0061] 最后，根据已知的目标区域的密度分布，对每个小矩形密度进行赋值，根据小矩形和二度体截面的几何关系，当小矩形中心位于二度体截面内时，赋值给所在位置的二度体密度，每个小矩形密度为常值，不同小矩形密度取值不同，以此刻画任意截面形状、任意密度分布的二度体。位于空气部分的小矩形，其密度值设为零，以此刻画起伏地形。

[0062] 2、矩形模型重力异常计算：

[0063] 步骤一中给出的矩形模型，其重力异常计算公式为

[0064]

[0065] 式中，(xm,z0)表示观测点坐标，z0为常值；L表示矩形模型在z方向剖分小矩形的个数；M表示矩形模型在x方向剖分小矩形的个数；(ξp,ζr)表示编号为(p,r)的小矩形几何中心坐标即在矩形模型中x方向为第p个且在z方向为第r个的小矩形的几何中心坐标；ρ(ξp,ζr)表示编号为(p,r)的小矩形的密度值；h(xm-ξp,z0-ζr)表示加权系数。

[0066] 实现式(1)的计算，分为三个环节：

[0067] 首先，根据观测点坐标(xm,z0)和小矩形几何中心坐标(ξp,ζr)，计算加权系数h(xm-ξp,z0-ζr)，其计算公式为

[0068]

[0069] 式中，γ表示万有引力常数，arctan()表示反余切函数运算符，ln()表示自然对数运算符；其它符号含义如下

[0070] X1＝ξp-0.5Δx-xm，X2＝ξp+0.5Δx-xm，Z1＝ζr-0.5Δz-z0，Z2＝ζr+0.5Δz-z0[0071] Δx，Δz表示小矩形几何尺寸。

[0072] 其次，在矩形模型中在x方向上位于同一行的所有小矩形为一层，采用一维离散卷积快速计算方法来计算一层矩形模型重力异常，其计算公式为

[0073]

[0074] 式中，表示第r层(r＝1,2,…,L)矩形模型在测线z0产生的重力异常；(xm,z0)表示观测点坐标；

[0075] 最后，将各层矩形模型重力异常进行累加，得到整个矩形模型的重力异常，即

[0076]

[0077] 步骤二中所述的一维离散卷积快速计算方法，其步骤为：

[0078] (1)将加权系数h(x1-ξp,z0-ζr)排列成向量t，记为

[0079] t＝[t0,t1,t2,…,tM-1,0,tM-1,tM-2,…,t2,t1]T (5)

[0080] 式中，矩阵元素ti与加权系数h(x1-ξp,z0-ζr)存在关系

[0081] ti＝h(x1-ξi+1,z0-ζr) (6)

[0082] (2)将第r层密度值ρ(ξp,ζr)(p＝1,2,…,M)排列成向量ρ，向量元素ρi与密度值存在关系

[0083] ρi＝ρ(ξi,ζr) (7)

[0084] 将向量ρ补零扩展成向量ρext

[0085]

[0086] 式中，0M×1表示M×1零向量；

[0087] (3)计算

[0088] 式中，fft()表示一维快速傅里叶变换；

[0089] (4)计算

[0090] 式中，“.*”表示对应元素相乘运算；

[0091] (5)计算

[0092] 式中，ifft()表示一维快速傅里叶反变换；

[0093] (6)提取矩阵gext的前M行元素，构成向量g，即为一维离散卷积计算结果，向量g中元素就是所求

[0094] 下面对本发明方法的效果进行检验。

[0095] 为了说明本发明所提出的方法用于计算任意截面形状、任意密度分布情况下复杂二度体重力异常时的效率和精度，设计了如下二度体组合模型(图3所示)：

[0096] 密度为常值的矩形区域内嵌一个密度均匀的圆形，圆心与矩形中心重合。矩形范围为：x方向从-1000m到1000m，z方向从0m到1000m(z轴向下为正)；圆形半径为400m。矩形密度为0g/cm3，圆形的密度为1g/cm3。将矩形剖分成10000×10000个大小相同的小矩形，计算高度为-50m测线(图3中虚线所示)上的重力异常，计算点个数为10000。

[0097] 重力异常算法利用Fortran语言编程实现，运行程序所用的个人台式机配置为：CPU为i7-2620，主频为2.7GHz，内存为32GB，四核八线程。运行所需时间约为6秒，由此可见算法效率很高。重力异常计算值和理论值如图4所示，从形态上看，两者是一致的。相对误差由理论值减去计算值得到差值的绝对值除以理论值得到(图5所示)，对相对误差进行统计，统计结果由表1给出，可知算法精度很高。

[0098] 表1重力异常理论值和计算值相对误差统计

[0099]

[0100] 本发明是一个有机整体，即在特定的模型表示方式条件下，建立矩形二度体重力异常叠加模型，根据一种特殊的加权系数计算公式，采用一维离散卷积快速计算方法，实现了重力异常计算在效率和精度上的统一。

[0101] 以上包含了本发明优选实施例的说明，这是为了详细说明本发明的技术特征，并不是想要将发明内容限制在实施例所描述的具体形式中，依据本发明内容主旨进行的其他修改和变型也受本专利保护。本发明内容的主旨是由权利要求书所界定，而非由实施例的具体描述所界定。

附图说明

[0049] 图1为本发明的计算流程图；

[0050] 图2为复杂二度体模型表示；

[0051] 图3为截面为圆形二度体模型图；

[0052] 图4为重力异常计算值和理论值对比图；

[0053] 图5为计算值与理论值的相对误差；

[0054] 图中符号说明如下：

[0055] ρ：表示密度。

1一种二度体重力异常计算方法 2一种流体通道间隔设置的重力热管 3一种利用人体重力自动开启的节能环保型门灯