专利名称：：基于标准格架的快速三维重力、磁力物性反演的方法
技术领域：
：本发明涉及地球物理勘探技术，具体是一种基于标准格架的快速三维重力、磁力物性反演的方法。
背景技术：
：目前使用的地面重力、磁力勘探方法，以石油勘探行业重力、磁力勘探为例，执行的标准为《地面重力勘探技术规程(SY/T5819-2002)》和《地面磁法勘探技术规程SY/T5771-2004》，对重磁资料的处理解释技术要求均是基于平面处理、2D\2.5D(二度\二度半)剖面反演、单界面重磁界面反演而进行的，尚没有提出要求三维重磁反演处理解释的规定和成熟推广的技术。现有重磁反演技术一般采用2D\2.5D剖面反演或单界面重、磁界面反演，而现在的简单三维重磁反演仅对简单的标准几何形体(球体、方柱体、台阶)进行简单三维反演；这些反演方法仅能实现二维反演或二度半剖面反演、单个简单几何形体的反演，或平面的单一界面反演，这种反演结果在复杂三维情况下与实际地质体模型差异较大，甚至出现失真，无法模拟地下复杂三维构造、结构，无法满足石油勘探对复杂三维反演的需要。而重磁三维物性反演方法尚没有实用或推广应用技术可比较。
发明内容本发明的目的是提供一种快速、实用化的基于标准格架的快速三维重力、磁力物性反演的方法。本发明的实施步骤和方法如下1)野外实测采集重力、磁力数据，处理得到的重力异常和磁力异常平面网格数据，磁异常要求为磁力化极异常；2)对重磁数据进行网格，将平面数据网格插值转换为正方形节点网格；歩骤2)所述的网格插值采用常规的插值常用克里金算法。3)对下半空间的剖分并进行标准格架快速正演计算与存储；步骤3)所述的对下半空间的剖分是采用等间距的立方体剖分，或采用立方体的高度的一半作为新长方体的剖分方式。上述剖分的平面网格节点数与步骤2)所述的平面数据网格相等。步骤3)所述的标准格架是标准模型几何格架，水平范围大小及深度范围与反演深度有关。上述标准模型几何格架体是对称的包含大量独立几何格架单元的巨型长方体，其水平范围的半边长是反演深度的2倍，其水平面的四个边长相等，而深度范围与反演深度相当；同时，其内部剖分出的各个几何格架单元长方体各边的边长与下半空间剖分时的单元体的各个边长对应。步骤3)所述的正演与存储是计算标准几何格架俯视图第1象限上区(ls)的几何格架系数的几何系数。步骤3)所述的正演计算与存储是由标准几何格架俯视图第1象限上区(ls)的几何格架数据体对应地衍生标准几何格架俯视图第1象限下区(lx)的几何格架数据体，再利用对称性由1区直接得到标准几何格架俯视图第2、3、4象限区的几何格架数据体。步骤3)所述的标准格架正演存储是存储标准格架正演获得的格架体系数三维数据，存储时仅存储标准几何格架的第一象限的三维数据系数，存储量为标准格架的1/4。4)对标准几何格架中心柱体在计算点的格架系数求取平方根，获得面元反演系数数组；步骤4)所述的面元反演系数数组是取中心几何格架的平方根所获得的系数序列，由上至下，构成一个一维系数数组。5)利用己有实测数据进行快速外推的外扩计算；步骤5)所述的快速外推的外扩法计算公式y(i+b)=2*y(b)—y(b-i)，其中边界的坐标为b，边界坐标的外扩增加量为i。6)用波长滤波器进行平面重磁异常分离，获得三维重磁异常数据体；步骤6)所述的波长滤波器是数字滤波器，深度随波长的增加而增加。步骤6)所述的波长滤波器计算为波长等于l个网格距时，滤波器系数为gr=L02g。-0.5g2-0.0g3(15-1)而波长等于2个网格距时，滤波器系数为gr=L02g。-—0.75g2-0.0g3(15—2)而波长大于或等于3个网格距时，滤波器系数为gr=1.02g。-0.15g-0.35g21-0.5g31(15-3)其中g,为分离求取的剩余异常，g。为计算点异常，g。g2、&为第1、2、3各个环带上参与计算点的异常的平均值，3至10个网格距时，g、g21、g^分别取环带数为1、2、3，1、2、4，1、3、5，2、4、6，2、4、7，2、5、8，3、6、9，3、6、10，大于10个网格距时类推，获得不同波长对应的计算公式。步骤6)所述的波长滤波器用于三维重磁波长层离异常的分离；7)对获得的三维重磁异常进行重磁层离面元反演，采用下式计算o(x，y，n)"Agr6S(x，y，n)/K，得到面元的地层视密度、视磁化强度，其中Agm(x，y，n)为三维重磁异常中的第n层重磁异常，(x，y，n)表示第n层面元，o(x，y，n)代表面元密度，K为系数，K由步骤4)几何格架中引导得到；8)进行反演物性数值区间约束和场级约束；步骤8)所述的反演物性数值区间约束是根据实际确定的数值区间，约束反演计算的物性参数的数值不允许超过这个区间，当数值超过区间的最大值或最小值时，将物性参数数值赋值为该区间最大值或最小值。步骤8)所述的场级约束是依据重磁异常数值大小特征，建立重磁异常的数值高低与场源分布特征对应的物性参数约束数组，约束数组是将区内异常分为12个级别，重磁物性模型为纵向剖分为k层的立体区域，建立不同异常幅值级别对应的物性参数系数，k层的各层内系数的均值，第一层取0.005，目标层最高取O.6，层间递增10%5%，而层内从低至高的变化量为-20%至+20%。实现对反演模型的约束建模和特性约束。9)对反演结果模型进行正演，检验反演拟合程度，进行多次迭代反演，直至得到高精度的反演结果。步骤9)所述的三维重磁正演是利用反演获得的密度或磁化率三维数据体，采用已存储的标准几何格架系数数组，利用其快速正演计算模型的重力或磁力正演异常。上述三维重磁正演是第一步，分别计算第l、2、3、4区模型在计算点产生的重磁异常，计算范围纵向第一层至最低层，扩边增量i方向从l到Uj方向从l到L。，计算公式分别为gl=wx(m0+i，n0+j，k)氺gml(i，j，k)g2=wx(m。-i，n0+j，k)*gml(i，j，k)g3=wx(m0-i，n0_j，k)*gml(i，j，k)g4二wx(m0+i，n。-j，k)氺gml(i，j，k)第二步，计算模型在计算点产生的总的重磁异常值，g(m0，n0)=gl+g2+g3+g4，其中设计算点坐标为(m。，n。)，反演模型深度网格数为k。，设扩边边长增量为L。，g(m。，n。)为计算点正演重力异常值，wx(i，j，k)为反演所得的三维物性数组，gml(i，j，k)为标准几何格架第一象限的三维系数数组，gl、g2、g3、g4分别为第l、2、3、4象限区产生的重磁异常。步骤9)所述的迭代反演是进行多次迭代逐次逼近，在迭代过程中，对初步反演拟合差所得异常视为下次反演拟合目标，进行迭代循环，对其进行再次重磁异常层离法分离和场级约束反演，对上次的反演拟合差，重复采用6)、7)、8)、9)步骤，实现迭代循环，直至反演拟合差满足精度要求而结束反演，或者，迭代次数已满足终止要求而结束反演。本发明具有高效快速、可靠有效、实用化的效果，应用于实测重磁数据的反演处理，获得了三维反演数据体，反映了地下地质体的三维密度、磁化率分布，解决了复杂区实际地质问题，推动了3D重磁反演的技术进步。图l基于标准格架的快速三维重力、磁力物性反演法流程图。图2标准几何格架(俯视图)数据体衍生平面分区示意图。图3反对称式外扩方法剖面示意图。图4三维物性反演重力异常图。a图为观测异常，b图为三维反演模型的正演计算异常。异常单位mgal，异常等值线间距0.2mgal。图5三维重力物性反演密度切片图。a图为三维密度375m切片，b图为三维密度-125m切片，c图为三维密度-625m切片，d图为三维密度-1125m切片，e图为三维密度-1625m切片，f图为三维密度-2125m切片。密度等值线间距0.02g/cm3。具体实施例方式以下结合附图和实施例进一步说明本发明。依据本发明进行了实际资料的反演处理。以下实例为场级约束法3D重磁物性反演法在S区所做的三维重力反演的实施例及其效果。采用基于标准格架的快速三维重力、磁力物性反演法流程(图1)，依次进行。执行步骤1)，准备需要反演的野外采集的重磁异常数据，要求磁力异常为化极磁力异常；执行本发明步骤2)，网格化实测的某地S区重力异常数据，获得平面网格数据(节点102x102)，平面网格距500mx500m;执行步骤3)，剖分下半空间，下半空间剖分网格采用500mx500mx250m;标准几何格架剖分网格采用500mx500mx250m，深度取8000m，半边长取说明书第7/ll页24000m;快速正演标准几何格架，并存储第一象限几何格架系数数组。正演计算与存储是其实现步骤如下(图2):第一步计算第1象限上区(IS)标准格架系数数组，由常规的计算单个长方体的重磁公式，生成ls区的数组各系数gnus(j，i，k)，计算范围为纵向(k)为第一层至最低层，东西(i)方向从0到L。(即扩边增量)，南北(j)方向从O到i;第二步由Is区域数据生成第1象限下区(lx)数据，采用计算式gmlX(i,j，k)=gmlS(j，i，k)计算，计算范围纵向第一层至最低层，东西(i)方向从O到L。(即扩边增量)，南北(j)方向从i到U;第三步存储1区格架系数数组；第四步由1区域数据生成2区域数据，采用计算式gm2(i，-j，k)=gnu(i,j，k)计算，计算范围纵向第一层至最低层，东西(i)方向从1到L。(即扩边增量)，南北(j)方向从1到L。；第五步由1区域数据生成3区域数据，采用计算式gm3(-i，-j，k)=gIIU(i，j，k)计算，计算范围纵向第一层至最低层，东西(i)方向从1到L。(即扩边增量)，南北(j)方向从1到L。；第六步由1区域数据生成4区域数据，采用计算式gm4(-i，j，k)二gmi(i，j，k)计算，计算范围纵向第一层至最低层，东西(i)方向从1到L。(即扩边增量)，南北(j)方向从1到L。；以上L。为设扩边增量，i为东西向变量，j为南北向变量。执行步骤4)，计算面元反演系数数组。反演深度为8km，平面等距剖分为500mX500m，纵向剖分为250m，几何格架中心柱体的磁场响应系数产生的系数K的数组如下he(1)=0.9830493，he(2)=0.6616288，he(3)=0.4505335，he(4)二O.3347392，he(5)二O.2648031，he(6)二O.2185804，he(7)=0.1859178，he(8)=0.1616664，he(9)=0.1429713，he(10)=0.1281265，he(11)=0.1160728，he(12)=0.1060659，he(13)=9.7659901E-02，he(14)=9.0471692E-02，he(15)二8.4280781E-02，he(16)=7.8869708E-02，he(17)二7.4110262E-02，he(18)=6.9896609E-02，he(19)=6.6141181E-02，he(20)=6.2747397E-02，he(21)=5.9716422E-02，he(22)=5.6933377E-02，he(23)=5.4432638E-02，he(24)=5.2103728E-02，he(25)=4.9993400E-02，he(26)=4.8040159E-02，he(27)=4.6198893E-02，he(28)=4.4556264E-02，he(29)=4.2952269E-02，hc(30)二4.1550443E-02,he(31)=4.014966犯一02，he(32)=3.8917955E_02。执行步骤5)，实施数据扩边，扩边增量取24000m，采用S地区处理后的观测异常重力异常，对重力异常进行半边长取24000m的反对称式扩边(图3)，外扩计算式为y(i+b)=2*y(b)_y(b-i);执行步骤6)，进行位场波长层离，分离三维重力异常；计算公式如下所述波长等于l个网格距时，滤波器系数为gr=L02g。-0.5g2-0.0g3(15-1)而波长等于2个网格距时，滤波器系数为&=1.02g。--0.75g2-0.0g3(15-2)而波长等于3个网格距时，滤波器系数为gr=1.02g0-0.15&-0.35g2-0.5g3(15-3)波长等于4个网格距时，滤波器系数为&=1.02g。-0.15&-0.35&-0.5g4(15-4)波长等于5个网格距时，滤波器系数为gr=1.02g0-0.15g,-0.35g3-0.5g5(15-5)波长等于6个网格距时，滤波器系数为^=1.02g0-0.15&-0.35g4-0.5g6(15-6)同理，波长等于7个网格距时，滤波器系数为^=1.02g。-0.15g2-0.35g4-0.5g7(15-7)波长等于8个网格距时，滤波器系数为g,=1.02g。-0,15g2-0,35g5-0.5g8(15-8)波长等于9个网格距时，滤波器系数为&=1.02g。—0.15g3-0.35g6—0.5g9(15—9)波长等于10个网格距时，滤波器系数为gr=1.02g。-0.15g3-0.35g6-0.5gl。(15-10)大于10个网格距时，可以以此类推，获得不同波长对应的计算公式。g,为分离求取的剩余异常，g。为计算点异常，gl、g2…g,。为第l、2…10各个环带上参与计算点的异常的平均值。执行步骤7)，进行位场层离面元法反演，计算公式如下；o(x，y，n)"△gres(x，y，n)/K。其中Ag^(x，y，n)为局部重力异常，(x，y，n)表示第n层面元，o(x，y，n)代表面元密度，K为系数，K由歩骤4)几何格架中引导得到；执行步骤8)，进行重磁场级约束反演，采用极值约束为-O.12g/cm3、15-O.12g/cm3，场级约束数组见重磁场级约束系数表，目标区内的异常分为12个级别，重磁物性模型为纵向剖分为32层的立体区域，建立了不同异常幅值级别对应的物性参数系数。重磁场级约束系数表如下:<table>tableseeoriginaldocumentpage16</column></row><table>终止反演迭代次数取300。重复采用6)、7)、8)、9)步骤，实现迭代循环，直至反演结束。反演拟合重力异常与实测异常拟合较好，拟合精度达到O.10x10—5ms—2，相b图为三维反演模型的正演计0.2mgal。获得了三维密度数据，展示出三维重力反演求解的三维密度分布的不同海拔高度的水平切片，展示的水平切片的埋藏深度为3000m5000m(图5)。对精度大2%(图4)。图中a图为观测异常，算异常。异常单位mgal，异常等值线间距图中a为三维密度375m切片，b图为三维密度-125m切片，c图为三维密度-625m切片，d图为三维密度-1125m切片，e图为三维密度-1625m切片，f图为三维密度-2125m切片。密度等值线间距0.02g/cm3。由反演结果可见，从3000m到5000m，地下密度分布具有明显的变化，在研究区中部，深层密度分布显示出更加清楚的南北向构造的形态，而且从浅至深，局部密度高的水平位置也发生了变化，表现出不同深度的高点偏移。反演速度对于三维反演模型节点(102x102x24)的反演，使用的计算机主频为1.80Hz、内存为512MB，反演计算耗时仅7分钟，对于三维重磁反演来说，这是一个极快的计算速度。权利要求1、一种基于标准格架的快速三维重力、磁力物性反演的方法，其特征在于采用以下步骤1)野外实测采集重力、磁力数据，处理得到的重力异常和磁力异常平面网格数据，磁异常要求为磁力化极异常；2)对重磁数据进行网格，将平面数据网格插值转换为正方形节点网格；3)对下半空间的剖分并进行标准格架快速正演计算与存储；4)对标准几何格架中心柱体在计算点的格架系数求取平方根，获得面元反演系数数组；5)利用已有实测数据进行快速外推的外扩计算；6)用波长滤波器进行平面重磁异常分离，获得三维重磁异常数据体；7)对获得的三维重磁异常进行重磁层离面元反演，采用下式计算σ(x，y，n)≈Δgres(x，y，n)/K，得到面元的地层视密度、视磁化强度，其中Δgres(x，y，n)为三维重磁异常中的第n层重磁异常，(x，y，n)表示第n层面元，σ(x，y，n)代表面元密度，K为系数，K由步骤4)几何格架中引导得到；8)进行反演物性数值区间约束和场级约束；9)对反演结果模型进行正演，检验反演拟合程度，进行多次迭代反演，直至得到高精度的反演结果。2、根据权利要求l所述的基于标准格架的快速三维重力、磁力物性反演的方法，其特征在于步骤2)所述的网格插值采用常规的插值常用克里金算法，所述的对下半空间的剖分是采用等间距的立方体剖分，或采用立方体的高度的一半作为新长方体的剖分方式；剖分的平面网格节点数与步骤2)所述的平面数据网格相等。3、根据权利要求l所述的基于标准格架的快速三维重力、磁力物性反演的方法，其特征在于步骤3)所述的标准格架是标准模型几何格架，水平范围大小及深度范围与反演深度有关。4、根据权利要求1或3所述的基于标准格架的快速三维重力、磁力物性反演的方法，其特征在于上述标准模型几何格架体是对称的包含大量独立几何格架单元的巨型长方体，其水平范围的半边长是反演深度的2倍，其水平面的四个边长相等，而深度范围与反演深度相当；同时，其内部剖分出的各个几何格架单元长方体各边的边长与下半空间剖分时的单元体的各个边长对应。5、根据权利要求l所述的基于标准格架的快速三维重力、磁力物性反演的方法，其特征在于步骤3)所述的正演与存储是计算标准几何格架俯视图第l象限上区(ls)的几何格架系数的几何系数；所述的正演计算与存储是由标准几何格架俯视图第1象限上区(ls)的几何格架数据体对应地衍生标准几何格架俯视图第1象限下区(lx)的几何格架数据体，再利用对称性由1区直接得到标准几何格架俯视图第2、3、4象限区的几何格架数据体；所述的标准格架正演存储是存储标准格架正演获得的格架体系数三维数据，存储时仅存储标准几何格架的第一象限的三维数据系数，存储量为标准格架的1/4。6、根据权利要求l所述的基于标准格架的快速三维重力、磁力物性反演的方法，其特征在于步骤4)所述的面元反演系数数组是取中心几何格架的平方根所获得的系数序列，由上至下，构成一个一维系数数组。7、根据权利要求l所述的基于标准格架的快速三维重力、磁力物性反演的方法，其特征在于步骤5)所述的快速外推的外扩法计算公式y(i+b)=2*y(b)-y(b-i)，其中边界的坐标为b，边界坐标的外扩增加量为i。8、根据权利要求l所述的基于标准格架的快速三维重力、磁力物性反演的方法，其特征在于歩骤6)所述的波长滤波器是数字滤波器，深度随波长的增加而增加；波长滤波器计算为波长等于l个网格距时，滤波器系数为gr=1.02g。-0.5g2_0.0g3(15-1)而波长等于2个网格距时，滤波器系数为gr=1.02g。-0.25g,-0.75g2-0.0g3(15-2)而波长大于或等于3个网格距时，滤波器系数为gr=1.02g。-0.15gn-0.35g2I-0.5g31(15-3)其中g,为分离求取的剩余异常，g。为计算点异常，g。g2、&为第1、2、3各个环带上参与计算点的异常的平均值，3至10个网格距时，gll、g21、^分别取环带数为1、2、3，1、2、4，1、3、5，2、4、6，2、4、7，2、5、8，3、6、9，3、6、10，大于10个网格距时类推，获得不同波长对应的计算公式。9、根据权利要求l所述的基于标准格架的快速三维重力、磁力物性反演的方法，其特征在于步骤6)所述的波长滤波器用于三维重磁波长层离异常的分离；所述的反演物性数值区间约束是根据实际确定的数值区间，约束反演计算的物性参数的数值不允许超过这个区间，当数值超过区间的最大值或最小值时，将物性参数数值赋值为该区间最大值或最小值。10、根据权利要求1所述的基于标准格架的快速三维重力、磁力物性反演的方法，其特征在于步骤8)所述的场级约束是依据重磁异常数值大小特征，建立重磁异常的数值高低与场源分布特征对应的物性参数约束数组，约束数组是将区内异常分为12个级别，重磁物性模型为纵向剖分为k层的立体区域，建立不同异常幅值级别对应的物性参数系数，k层的各层内系数的均值，第一层取0.005，目标层最高取0.6，层间递增10%5%，而层内从低至高的变化量为-20%至+20%。实现对反演模型的约束建模和特性约束。11、根据权利要求1所述的基于标准格架的快速三维重力、磁力物性反演的方法，其特征在于步骤9)所述的三维重磁正演是利用反演获得的密度或磁化率三维数据体，采用已存储的标准几何格架系数数组，利用其快速正演计算模型的重力或磁力正演异常；上述三维重磁正演是第一歩，分别计算第l、2、3、4区模型在计算点产生的重磁异常，计算范围纵向第一层至最低层，扩边增量i方向从l到L。，j方向从l到L。，计算公式分别为gl二wx(m。+i，n。+丄k)*gml(i，j，k)g2=wx(mo-i，n0+j，k)%ml(i,j，k)g3=wx(m0_i,n0_j，k)氺gml(i，j，k)g4=wx(m0+i，n。-j，k)氺gml(i，j，k)第二步，计算模型在计算点产生的总的重磁异常值，g(m0，n0)=gl+g2+g3+g4，其中设计算点坐标为(m。，n。)，反演模型深度网格数为k。，设扩边边长增量为L，g(m。，n。)为计算点正演重力异常值，wx(i，j，k)为反演所得的三维物性数组，gml(i,j，k)为标准几何格架第一象限的三维系数数组，gl、g2、g3、g4分别为第l、2、3、4象限区产生的重磁异常。12、根据权利要求1所述的基于标准格架的快速三维重力、磁力物性反演的方法，其特征在于步骤9)所述的迭代反演是进行多次迭代逐次逼近，在迭代过程中，对初步反演拟合差所得异常视为下次反演拟合目标，进行迭代循环，对其进行再次重磁异常层离法分离和场级约束反演，对上次的反演拟合差，重复采用6)、7)、8)、9)步骤，实现迭代循环，直至反演拟合差满足精度要求而结束反演，或者，迭代次数已满足终止要求而结束反演。全文摘要本发明是地球物理勘探中基于标准格架的快速三维重力、磁力物性反演的方法。先对重磁数据进行网格，再转换为正方形节点网格，对下半空间的剖分并进行标准格架快速正演计算与存储，对架中心柱体计算点的格架系数求取平方根，获得面元反演系数数组，快速外推的外扩计算，平面重磁异常分离后获得三维重磁异常数据体，进行反演物性数值区间约束和场级约束，再正演得到高精度的反演模型。本发明可反映地下地质体的三维密度、磁化率分布，解决了复杂区实际地质问题，具有高效快速、可靠有效、实用化的效果。文档编号G01V11/00GK101661115SQ20081011946公开日2010年3月3日申请日期2008年8月29日优先权日2008年8月29日发明者刘云祥申请人:中国石油集团东方地球物理勘探有限责任公司