高精度三维地震勘探技术在长沟煤矿的应用_《环球市场信息导报》杂志投稿_天涯博客
投稿信箱：咨询***：编辑QQ：
今日访问：[$DayVisitCount$]
总访问量：3619
开博时间：
博客排名：暂无排名
　　高精度三维地震勘探技术在长沟煤矿的应用　　石君华（中国煤炭地质总局水文物测队 河北邯郸 056001）　　摘
要：该文以阳泉煤业(集团)长沟煤矿三维地震勘探工程为例，讨论了表、浅层地震地质条件复杂、地层倾角较大、部分地段煤层埋深较浅且巷道采空区错综复杂等特点，设计合理的观测系统，充分发挥多道仪器的优势，以达到经济技术的合理性，资料处理解释中应用叠前时间偏移和属性解释技术，对采空区边界的圈定和巷道采空区的精细解释取得了较好的地震勘探成果。　　关键词：三维地震勘探；观测系统；巷道采空区；叠前时间偏移；属性解释　　作者简介：石君华(1972-)男，四川简阳人，汉族，物探高级工程师，1996年毕业于中国矿业大学勘查地球物理专业，一直从事地球物理勘探工作。　　阳泉煤业(集团) 长沟煤矿是整合矿井，勘探区为一典型的山区地貌，沟谷纵横，梁垣陡立，高差变化较大树林密布，最大相对高差约383.18m，不利于机械化施工。本区浅层地层结构就岩性而言，主要分为二种类型：黄土覆盖区，面积约占全区的30%，主要分布于山脊两侧的平台、斜坡及沟谷，厚度一般小于20m，岩性以含砂粘土和粘土为主，不含水，波速极低；基岩出露区，分布于测区的大部，占全区面积约70%，岩性主要为砂岩、泥岩等，个别岩石比较坚硬。本次勘探的主要目的层15#煤层埋藏深度变化较大，最大埋深500m，最小埋深80m，一般埋深在200m-300m；15#煤层平均厚度5.2m ，沉积稳定，结构简单，倾角变化较大（5°～25°），一般为10°左右。测区巷道采空区较多、东北部煤层埋深较浅。表、浅、深层地震地质条件复杂。　　根据测区的地震地质条件和采空区特点，选用合适的观测系统、野外数据采集、资料处理、资料解释技术获得了质量较高的地震地质成果。　　1. 本区技术难点　　①复杂的地形条件严重影响激发、接收点的到位率；②多变的浅层岩性，特别是黄土覆盖地段给激发带来一定的困难，影响高分辨率、高信噪比资料的获得；③井下作业的强噪音可能对提高资料信噪比带来一定困难；④东、北部煤层埋藏浅，对观测系统设计及施工不利，难以发挥多道仪器的作用。　　2. 技术对策 　　对偏离设计点位的激发、接收点位置进行二次测量；增加叠加次数以提高资料信噪比；合理选取激发井深、药量，提高激发主频、拓宽有效频带，提高资料信噪比。尽量统一激发岩性。在东北部煤层埋深较浅地段采用较小的炮点网格和较少的检波点数。即在东北部使用炮点网格40m×20m、单线20道的观测系统，其它地段则使用炮点网格40m×30m、单线30道的观测系统；采用16线（320--480道）接收，最大程度发挥多道优势，减少炮数，达到经济技术合理。　　3. 观测系统的设计　　针对本次三维勘探地质任务的要求和本区地震地质条件，在观测系统类型的设计上主要考虑了以下几点：　　采用较小面元，提高纵横向分辨率，确保小构造落实的精度；本区为一单斜构造，倾角较大，因此采用单边放炮；每次滚动八条检波线，在提高静校正耦合精度的同时，改善炮检距、方位角分配的均匀程度，减小采集脚印；本区埋深较浅且埋深变化较大，应保证覆盖次数均匀以及最大炮检距符合要求；充分发挥多道仪器的优势以及经济技术的合理性，采用16线接收，线距20m。通过以上分析，本次勘探采用单边放炮16线4炮制规则束状观测系统，见图1。　　
　　图1一个排列片观测系统图（左图：炮排距30m、单线30道；右图：炮排距20m、单线20道）　　观测系统主要参数：16线4炮制线束状观测系统，端点激发；接收道数：16×30=480（深部），16×20=320（浅部）；CDP网格：10m（横）×5m（纵）；覆盖次数：20次（横4×纵5）；最小炮检距：22.36m；最大炮检距：374.4m，297.0m。　　4. 地震数据采集　　由于区内浅层地震地质条件复杂，既有基岩裸露，又有黄土覆盖、还有坡积物地段。出露地层横向岩性变化大、岩性复杂，对地震波的激发接收不利。本次勘探用风镐和人工钻结合成孔，基岩出露地段井深为3m，黄土地段成孔穿过黄土到基岩内1m，坡积物地段成孔要打穿坡积物并在基岩内激发，下药后把炮孔填实；检波器布设时，必须除去地表植被，将检波器插实，实现检波器与岩土的最佳耦合。本次采用TNT成型炸药，单孔药量1Kg～2Kg；采用60HZ动圈式检波器3只串联点式接收；仪器采样率：1ms，记录长度：1s，前放增益：12dB, 全频带接收，见图2典型地震记录。　　
基岩出露区井深3m药量1Kg
黄土覆盖区井深6.3m药量1Kg　　图2
典型地震记录　　5. 资料的处理及解释　　由于本区地层倾角较大，构造复杂，采空区多且赋存情况特别复杂，随开采程度、开采方式、开采厚度、采空区年限等不同而有较大差异。测区南部15105工作面是三维地震勘探前的新近采空区，机械化程度较高、开采率较高，反射波在时间剖面上主要表现为:被采煤层反射波突然消失；而测区北部主要采用巷采，都是过去靠人工挖采造成，错综复杂，开采率较低、采空年限较长，反射波在时间剖面上主要表现为：错断下凹、振幅变弱、极性反转、波形变化无规律、凌乱等，部分地段就一条3m左右宽的巷道采过去，在地震剖面上反射波没有明显变化，部分地段采空巷道较多，反射波稍向下弯曲。这些巷道采空区给资料解释带来难度。根据矿方要求对采空区边界圈定精度要高，本次解释通过对叠前时间偏移和叠后时间偏移的对比分析以及属性体解释技术的研究，针对测区北部巷道采空区，解释主要采用叠前时间偏移数据体和主频属性解释对巷道采空区经行了圈定，通过资料验证，获得了较好的效果,见图3～4。　　 　　图3
测区北部采空区CK3在叠后时间偏移数据体（上）和叠前时间偏移数据体（下）时间剖面上的显示　　 　　图4
主频属性解释　　从图3中可以看出利用叠后时间偏移所解释的采空区CK3范围小，而利用叠前时间偏移所解释的采空区CK3范围较大，分析解释叠前、叠后两个采空区边界之差范围处的时间剖面，叠前时间剖面比叠后时间剖面反映的采空区特征更明显；图4中所圈定的范围是矿上目前已采宽度4m间距36m的3条并行巷道处，时间剖面上没有明显反映，而本区通过对属性解释技术的研究，发现在主频属性解释上有明显反映，因此根据主频属性对测区北部的巷道采空区进行了解释圈定。　　6. 地质成果及验证情况　　本次勘探查明了测区内8#、15#煤层的赋存情况。在测区内共解释了断层3条，落差均大于5m；共解释陷落柱7个；挠曲1个；基本查明了8#煤层中倾角大于15°的区段有五处，15#煤层中倾角大于15°的区段有五处；经过本次勘探对原有采空区CK3边界范围进行了修正，在测区中部和北部新发现采空区2个。本次三维地震勘探解释的3条断层（DF1、DF2、DF3），3个陷落柱（DX1、DX3、DX5）和1处采空区（CK2）现已被巷道和井工工程部分揭露，解释成果与实际验证情况基本吻合。　　通过此项目的研究发现，在成孔条件复杂、成孔成本高的地区可以充分发挥多道仪器的优势，以达到经济技术的合理性。对过去机械化程度低、开采率较低、没有规划、巷道采空区较多且错综复杂、地层倾角较大、构造复杂的三维地震勘探区，采用叠前时间偏移数据体和属性解释技术圈定的采空区更合理，解释的采空区边界更准确。　　　　参考文献　　[1]何樵登等. 地震勘探原理和方法[M].北京: 地质出版社,1986. 　　[2]李庆忠.走向精确勘探的道路--高分辩率地震勘探系统工程剖析[M].北京:石油工业出版社,1995.　　[3]关民全等. 阳泉煤业(集团) 长沟煤矿地面综合物探报告[R].中国煤炭地质总局水文物测队,2011.　　　　　　分类： |障碍物密集区三维地震勘探施工方法
&&&&2014, Vol. 29 Issue (2): 638-641
鲍五堂, 石磊军, 邱兆泰, 张灯亮, 徐忠华, 张建军. 2014. 障碍物密集区三维地震勘探施工方法[J]. 地球物理学进展, 29(2): 638-641, doi: 10.6038/pg&& BAO Wu-tang, SHI Lei-jun, QIU Zhao-tai, ZHANG Deng-liang, XU Zhong-hua, ZHANG Jian-jun. 2014. Construction method of three-dimensional seismic exploration in obstacles concentration areas. Progress in Geophysics, 29(2): 638-641, doi: 10.6038/pg &&
障碍物密集区三维地震勘探施工方法
鲍五堂, 石磊军, 邱兆泰, 张灯亮, 徐忠华, 张建军&&&&
河北省煤田地质局物测地质队, 邢台 054000
作者简介：鲍五堂，男，1974年生，高级工程师，长期从事煤炭勘探管理和技术工作.（E-mail：）
摘要：在障碍物密集区进行三维地震勘探，由于城镇、村庄、厂矿企业、铁路遍布测区，障碍物占整个施工面积的50％以上，激发点和接收点布设困难，难以获得地下资料.如何在这样的测区取得高质量地震资料是一个技术难题.需要从设计到施工的全面配合.本文以江苏某区为例说明了在障碍物密集区取得高质量地震资料的设计手段和施工措施.深入讨论了障碍物密集区三维地震施工的方法对策，并概括说明施工方法对后续处理解释的影响和结果.利用这些手段和措施在该区获得了比较满意的地震资料.从而对破解这一难题做出了有益的探索.
障碍物密集区&&&&
三维地震&&&&
处理解释&&&&
Construction method of three-dimensional seismic exploration in obstacles concentration areas
BAO Wu-tang, SHI Lei-jun, QIU Zhao-tai, ZHANG Deng-liang, XU Zhong-hua, ZHANG Jian-jun&&&&
Geophysical Prospecting & Surveying Team of Hebei Provincial Coalfield Geology Bureau, Xingtai 054000, China
Abstract: At present, the obstacles concentration areas for 3 d seismic exploration, because of the towns, villages, factories and mines enterprises, railway around the measurement area, the obstacles of the total construction area of more than 50%, shot point and acceptance point setting difficult, difficult to get information underground. How the area to get high quality seismic data in this test is a technical problem. Need full cooperation from design to construction. Taking jiangsu district as an example illustrates the high quality seismic data obtained in obstacles concentration areas design methods and construction measures. Deeply discussed the obstacles concentration areas the construction method of 3d seismic countermeasures, and summarize the construction methods for subsequent processing and interpretation and the influence of the results. Use of these means and measures obtained satisfactory seismic data in this area. Which made a beneficial exploration to solve this difficult problem. In this area has important theoretical and practical significance.
Key words:
obstacles concentration areas&&&&
construction&&&&
3D seismic&&&&
processing and interpretation&&&&
目前在障碍物密集区的三维地震勘探中如何取得高质量地震资料是一个重要技术难题，需要从设计到施工的全面配合(
在障碍物密集区如城镇、村庄、厂矿企业、铁路遍布区，特别是障碍物占整个施工面积的50％以上的测区，存在许多技术难点，主要为激发点和接收点布设困难，难以获得覆盖区的地下资料(
)其主要对策是采用特殊三维观测系统，利用计算机辅助设计软件根据障碍物的分布情况灵活地设计炮点和检波点位置，根据目的层的产状、深度和分布范围及时计算出地下目的层的覆盖次数及面元内炮检距分布情况，对不合理的分布情况进行适当调整，并制定相应的施工技术措施，从观测系统和施工质量几方面入手，进而获得在障碍物下面的高质量地震资料.
本文通过实例说明在障碍物密集区取得高质量地震资料的设计手段和施工措施，期望对破解这一重要技术难题做出探索.
1 方法原理
三维地震勘探是利用炮点网格激发，检波点网格接收，从而获得地下一定范围内均匀分布的达到一定叠加次数的数据体，来达到控制地下构造形态的勘探方法(
).特殊三维观测系统根据障碍物的分布情况灵活地设计炮点和检波点位置.所以，三维观测系统选择的正确与否，直接影响勘探精度、施工效率和勘探成本.
三维地震勘探要求地下CDP点均匀分布，且具有较高的信噪比和分辨率，因此，纵、横向应满足一定的覆盖次数.
1.1 CDP网格的确定
CDP网格的选择主要由地质任务、地震地质条件和经济效益等因素决定.为防止出现空间假频，必须满足
其中，DX≤Δx／2，DY≤Δy/2，DX和DY为纵横向样点间隔；Δx、Δy为纵横接收点间距.vR为均方根速度；fmax为有效波最高频率；θx、θy为纵、横方向上地震波入射到地面的角度，可用地层倾角Φx、Φy代替.
1.2 叠加次数选择因素
覆盖次数的选择要有助于提高资料的整体信噪比，加强构造主体部位和圈闭部位的有效覆盖次数.特别针对勘探区内干扰波的发育情况设计覆盖次数，覆盖次数的选择应能充分压制干扰(随机噪音)、增加目的层的反射能量，从而提高资料的信噪比，确保成像效果.另外，由于三维地震勘探获得的地震信息来自于地下的各个方向，为保证丰富的反射信息并 满足资料处理的要求，在覆盖次数一定的情况下，覆盖次数的选择还需考虑Crossline方向自动剩余静校正及速度分析的精度，即必须保证纵横向覆盖次数较为均匀.同时，还必须保证纵向上有足够的炮检距分布，从而提高速度分析精 度.
1.3 炮检距
三维地震勘探中，沿接收线方向的炮检距称纵向炮检距X，沿垂直接收线方向的炮检距称横向炮检距Y，最大非纵炮检距为
最大非纵炮检距Xmax的选择原则为：
(1)考虑求取速度的精度，Xmax越大越好;
(2)压制多次波的效果，Xmax越大越好;
(3)考虑动校正拉伸畸变对高频信号的影响及反射系数的变化，Xmax越小越好;
(4)Xmax一般与目的层埋深大致相当.
1.4 特观设计
三维地震勘探中，用计算机辅助设计，采用灵活多变的观测系统，在均匀覆盖的技术上，通过检波器、炮点的合理布置，加快了施工进度，保证了资料采集质量.提高精度是地震勘探的永恒追求，地震勘探的精度主要体现在成像精度和高分辨率.理想的作法不仅要增加空间采样密度，而且要采样均匀，并且要从叠前处理和偏移成像的角度来设计系统的技术方案：
(1)叠前去噪对炮点均匀性的要求.在高精度勘探中，理想的情况是各个方向都不出空间假频，这就要求空间采样不仅要密而且要均匀.有时为了弥补垂直接收线方向上单炮记录空间采样的不足，用垂直接收线的多个炮记录，按炮点到接收线的距离进行互换，模拟空间采样均匀的炮集，即所谓的正交子集.然后进行去噪，去噪后栽返回各炮记录.这也是目前国外高密度采样资料叠前去噪的新技术.但构建正交子集时用到了炮点和接收线的互换，认为各炮的干扰是相同的.这在地表不均匀，各炮干扰不相同时就有较大的局限性.通过重排，模拟出来的空间采样相对均匀的子集毕竟不是原始炮集.要达到好的去噪效果，最好使原始炮集各个方向上的空间采样密且均匀.
(2)叠后偏移对炮点均匀性的要求.偏移并非是对反射界面上二次点源子波进行收敛，实际上在地震记录或叠加剖面上并不存在孤立的二次点源子波.现行偏移的做法可简单地理解为把和绕射曲线相切处的反射波归位到绕射极小点.这样对于反射波来说，偏移实际上是不同相叠加，叠加效果与道密度和均匀性密切相关.
(3)静校正精度对炮点均匀性的要求.由于空间采样密度高并且均匀，使记录初至清晰且连续，能更细致地反映表层变化，有利于提高用初至波反演地表模型的精度和静校正的精度.
2 实际效果应用分析
2.1 工区概况
江苏某区属建筑物极密集区，地表条件极为复杂，区内河流、高压线路、矿区自营铁路、村庄密布，村庄占地总面积几乎占全区的50%以上(图 1)，其中村庄、公路连成片，给地震测线的布设、炮点的到位、检波器的安置带来较大困难.
图 1图 1工程示意图Fig. 1 Schematic diagram of engineering
2.2 解决方法
根据本区地质任务的要求，考虑到区内地表条件复杂，通过选用先进的无线遥测地震仪，应用科学合理、灵活的三维地震勘探观测系统，选用合适的观测系统参数，优化施工设计；采用常规观测系统与块状特殊观测系统相结合，同时进行现场处理，有针对性地布置有效激发点，来完善工区资料，以获取可靠的原始记录．采用特殊处理手段，经过精细处理和综合地质解释，可以获得良好的地质效果().
本次三维地震勘探激发方式采用井炮激发，借鉴邻区勘探成功经验，优选三维观测系统.我们采用高叠加次数，高频检波器接收，8线16炮制中点激发的束状观测系统施工.采用美国产BOX无线遥控数字地震仪，0.5 ms采样，1.5 s录制，CDJ-60型检波器接收.
2.3 特观设计
该区观测系统，参考邻区地震成果，以及本区的地质资料，有关地震勘探参数为：fmax=70 Hz，VR=3600 m/s，Φx=Φy≤14°.
从原理计算得出Δχ≤53.1 m、ΔY≤53.1 m，DX≤26.6、DY≤26.6.根据地质任务要求，CDP网格定为10 m×10 m，完全可以满足空间采样定理的要求.检波线间距为Δx=80 m，Δy=20 m.叠加次数24次.
由于该区公路、矿区铁路、河流及村庄已连成一片，给地震测线布设带来了极大的障碍.许多地方必须进行特观处理(图 2).根据障碍物的分布情况灵活地设计炮点和检波点位置(图 3)，并及时计算出地下目的层的覆盖次数及面元内炮检距分布情况，对不合理的分布情况进行适当调整.
特殊观测系统和叠加图Fig. 2 Special observation system and the superposed graph
村中放炮位置图Fig. 3 The shooting location map
2.4 数据处理
数据处理阶段针对该区地形起伏变化比较明显，土丘、沟坎发育，浅层速度变化较大，资料存在着较严重的静校正问题，利用美国I/O 公司的初至折射静校正软件进行静校正处理，运用分频手段分阶段进行全区去噪处理.
在处理中，进行了多参数细致、反复认真的测试，合理选择不同的流程进行资料处理.偏移速度场的建立采用“地质条件约束下建立偏移速度场”的方法. 针对本区检波点在村庄、河道的个别不到位现象及检波点的偏移记录，用检波点校正的重定位技术软件做好静校工作，确保炮、检位置关系准确.复杂地形区进行了二次精细的静校，对村庄障碍物等地段造成的检波点偏移情况均进行了精细归位，运用多种参数测试比较，确定好最佳的处理流程，最终获得了较好的资料处理成果(图 4).
障碍物下时间剖面Fig. 4 Obstacles time profile
2.5 数据解释
资料解释用GEOFRAM解释软件进行人机交互解释，纵向、横向和联井时间剖面对比，水平切片、沿层切片相结合，从时间剖面到平面反复认识，确保了解释成果的正确精细和可靠性，成果则采用CPS-3地质绘图软件输出，并配合AUTOCAD绘图.
障碍物密集区三维地震施工，需要从设计到施工全面配合.在施工中，利用特观设计软件优化三维地震观测系统，应充分利用村庄中及村庄周围空地布设炮点和检波点，并及时计算出地下目的层的覆盖次数及面元内炮检距分布情况，确保了障碍物下资料的完整性，处理中着重障碍物下资料进行精细处理和解释，从而获得可靠的地震资料.
总之，在障碍物密集区进行三维地震勘探，采用科学的野外资料采集方法，常规观测系统与特殊观测系统灵活结合，完全可以获得完整的地下地质资料.本文的施工方法对类似区域的地震勘探有指导作用.
感谢审稿专家和编辑部的大力支持。
Cao G B. 2010, Analysis of effect of offset Selection on the oblique crossing acquisition system[J]. Progress in Geophysics, 25(6): .
Feng S M, Feng X, Li Y L. 2009..
Guo Y M, et al. 2004.Building concentration areas high obliquity effect of 3 d seismic exploration[J]. Journal of China coal field geology, (1) : 51-52.
Geng L J. 2005.3D seismic wave survey for coal layer in Hami desert area of Xinjiang[J].
Geng C M, Lin N, Lin Q Z. 2012.The village populated areas 3D seismic exploration for coal mine design optimization method[J]. Shandong Coal Science and Technology..
Liu L N, Zhou L H, Zhu Z Y, et al. 2002.3D Wave Equation Prestack Depth Migration with Complex Layouts[J]. Progress in Geophysics, 17(4):658-662.
Liu B J, Hu P,Meng Y Q, et al. 2009. In the Beijing area crustal fine-structure of deep seismic reflection profile detection research[J]. Journal of geophysics,52 (9) : .
LI X J. 2011,Applications and development of city geophysical prospecting in China[J]. Progress in Geophysics, 26(6): .
Liu Y,Shen Y X,Niu X J. 2012,3D seismic prospecting in urban active fault accurate positioning of the meaning[J].Journal of earthquake, 34 (1):97-104.
Peter C Ashton,Brad Bacon,Christian Deplante,et al. 1994. 3D seismic survey design[J].
Qin G S,Cai Q X,Wang G H, et al. 2010. 3D Seismic Geometry Design Based on Pre-stack Imaging[J]. Progress in Geophysics, 25(1): 238-248.
Wu Q B,Wang X D,Ning F Z. 2005.Application and analysis of 3D seismic exploration in coal mine[J].
Wang Y J, Li G. 2009.Obstacles concentration areas of 3 d seismic observation system design and application[J].
Zhang Z, Liu C, Shao Z G. 2003. Theory and experimentation of the charge sizes in seismic-sources for seismic exploration[J].
Zhang Y H, Meng H X. 2006.High resolution 3 d seismic observation system design research[J]. China coalfield geology, S1:53-55.
曹国滨. 2010.炮检距选取对斜交观测系统的影响分析[J].
冯世民,冯鑫,李雅玲..
郭彦民. 2004.建筑物密集区高倾角三维地震勘探效果[J].
耿丽娟. 2005.新疆哈密戈壁滩区煤层三维地震勘探实践[J].
耿春明,林宁,林全章. 2012.村庄密集地区煤矿三维地震勘探设计优化方法[J].山东煤炭科技..
刘礼农, 周立宏, 朱振宇,等. 2002.复杂观测系统下的三维波动方程叠前深度偏移[J].
刘保金,胡平,孟勇奇,等.2009.北京地区地壳精细结构的深地震反射剖面探测研究[J].
李学军. 2011.我国城市物探的应用与发展[J].
骆宗强, 魏伟, 孙伟家,等. 2012.三维地震观测系统采集脚印定量分析[J].
刘影,沈月霞,牛小军.2012,三维地震勘探在城市活断层精确定位中的意义[J].
秦广胜,蔡其新,汪功怀,等. 2010.基于叠前成像的三维地震观测系统设计[J].
吴钦宝,王秀东,宁方助. 2005.煤田三维地震勘探的应用及效果分析[J].
武磊彬,徐?. 2009.叠前时间偏移技术在淮北某煤矿采区三维地震勘探中的应用与效果[J].
王玉娇,李刚. 2009.障碍物密集区三维地震观测系统的设计与应用[J].
张智,刘财,邵志刚. 2003.地震勘探中的炸药震源药量理论与实验分析[J].
张云花,孟红星.2006.高分辨三维地震勘探观测系统设计研究[J].带控制点的三维密度界面反演方法
&&&&2014, Vol. 29 Issue (6):
胡立天, 郝天珧. 2014. 带控制点的三维密度界面反演方法[J]. 地球物理学进展, 29(6): , doi: 10.6038/pg&&
HU Li-tian, HAO Tian-yao. 2014. The inversion of three-dimensional density interface with control points. Progress in Geophysics, 29(6): , doi: 10.6038/pg &&
带控制点的三维密度界面反演方法
胡立天1,2, 郝天珧1 &&&&
1. 中国科学院地质与地球物理研究所, 北京 100029;2. 中国科学院大学, 北京 100049
基金项目：国家油气重大专项(-006-30)、国土资源部海洋地质保障工程项目(GZH)、国家高技术研究发展计划(-02)和国家自然科学基金(814011)共同资助.
作者简介：胡立天, 男1989年生, 山东东营人, 2011年毕业于中国海洋大学勘查技术与工程专业, 现为中国科学院地质与地球物理研究所在读博士, 主要从事重力处理与解释.(E-mail:hulitian@mail.)
通讯作者：郝天珧, 女, 1957年生, 中国科学院地质与地球物理所研究员, 主要从事海、陆油气盆地综合地球物理研究.(E-mail:tyhao@mail.)
摘要：利用重力资料反演密度界面对于研究盆地、区域构造和深部构造有着重要意义.但是密度基准面深度和界面密度差的不准确性降低了反演结果的可靠性.本文使用已知深度的控制点,在逐步迭代中计算出最合适的密度基准面深度和界面密度差,使反演结果和控制点拟合最好.理论模型和南海地区的莫霍反演结果表明本文方法可以计算出最合适的密度基准面深度和界面密度差,取得理想的反演结果.
三维密度界面反演&&&&
已知控制点&&&&
The inversion of three-dimensional density interface with control points
HU Li-tian1,2, HAO Tian-yao1 &&&&
1. Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing 100029, C2. University of Chinese Academy of Sciences, Beijing 100049, China
Abstract: The inversion of density interface using gravity data has great values in basin research、regional tectonics and deep structure. However, the uncertainty of density datum plane and density contrast makes the inversion result unreliable. An inversion method of three-dimensional density interface with several known control points is presented in this paper. We calculate the datum plane and density contrast in each iteration to make the inversion result fit the control points best. The tests of synthetic model and south China sea gravity data prove that the method obtains suitable datum plane and density contrast and gets a good inversion result.
Key words:
three-dimensional inversion of density interface&&&&
known control points&&&&
重力异常是地下各岩性不均匀变化与各物性界面起伏的综合反映，重力密度界面反演对于研究区域构造基底起伏、莫霍面起伏以及岩石圈厚度具有非常重要的意义，一直以来都是地球物理反演中的重要问题(；；).反演过程中不同的基准面深度与界面密度差会反演出不同的结果，因此确定界面基准面以及界面密度差具有十分重要的意义.
近些年来，众多地球物理学家对如何选取界面密度差和基准面深度展开了研究，并且取得了许多研究成果.使用统计模型，对洋中脊、海洋和大陆的地壳及上地幔密度分别赋值进行反演；和在计算沉积层与地壳平均密度时，除了参照各地区在做重力反演时曾选用过的数据外，还根据地震纵波的实测值使用层速度-密度转换公式得到界面密度差；对青藏高原地区的莫霍面结构采用了变密度反演，浅层的密度由岩芯及野外标本实测得到，中深层的密度可由地震波速-密度的经验公式转换得到.然而不同的区域地质情况相差很大，对于不同的区域使用统计模型中统一的密度差是不合适的；此外在野外工作较少的地区，浅层岩芯及标本和中深层的地震资料较少，速度-密度转换公式也只是经验公式，在实际应用中会带来误差.求取重力基准面的传统方法是根据重力异常的功率谱估计界面的平均深度值，然而在功率谱的推导过程中有很多的假设和近似，在最后根据斜率求取也是根据个人经验，因此功率谱只能够大致估计界面的埋深，导致反演结果与已知深度点无法拟合的很好(；；；；；)；或者是在反演的区域选取不同参考深度进行运算，对每一种反演结果与已知深度点进行比较，选取误差较小的作为最终结果，这种算法计算量大以及所选取的不同参考深度数量有限并且在很大程度上依赖于个人经验，这无疑造成了重力反演结果的不准确性.
上述方法没有或者很少使用控制点深度约束反演过程，使最终反演结果可能与控制点相差较大.因此，本文采用一种密度界面上已知若干控制点深度的迭代反演法，通过控制点来计算最合适的界面密度差和基准面深度，逐次迭代使最终得到的界面反演深度与已知控制点相差最小，实现了三维单一密度界面起伏的模型计算，并对南海地区的实测剖面进行了反演，取得了较为理想的结果.
三维单一密度界面的正演计算
正演计算是反演过程的重要组成部分，概括起来主要分为两大类，即频率域和空间域.频率域主要是基于快速傅里叶变换的Parker反演法，计算速度很快.在空间域中对于形状不规则的地下介质，可将其分割为不同的几何单元分别求取异常，再相加得到总异常.根据分割单元方式，正演方法可以分为线元法、面元法、体元法、表面积法及三角棱柱组合法等，较为常用的为长方体组合模型.
本文采用长方体组合模型，如图 1所示，在坐标原点处产生的重力异常为
式中：l=(ε2+η2+ζ2)1/2，积分限(x1，x2)、(y1，y2)、(z1，z2)分别对应ε，η，ζ，F为万有引力常数，p为这个长方体的密度.
图 1Fig. 1
图 1 小立方体模型
Fig. 1 Single
2 带控制点的三维单一密度界面反演方法
现在的界面反演方法大都是在频率域中利用Parker-Oldenburg公式迭代反演，其优点是速度快，但存在对重力的“向下延拓”过程而产生高频噪声，为了消除高频噪声，需要采取低通滤波，在滤波的同时也把数据中有用的高频数据滤掉了(Oldenburg et al., 1974；柯小平等，2005).因此本方法结合压缩质面法与等效源法(；)，在空间域中进行反演运算，使运算过程中的误差尽量减小.
如图 2所示，h为地表，m为重力基准面，起伏的界面可以划分成许多小立方体.在反演过程中将小立方体压缩到压缩面上，根据重力异常值反演其高度，使用控制点信息调整界面密度差和基准面深度，计算出界面起伏深度，再对其正演值与实际重力异常值的偏差再次反演，直到到达一定次数或者偏差值达到指定范围为止，具体流程如下.
图 2Fig. 2
图 2 地下小立方体的划分
Fig. 2 Subsurface prism
如图 3所示，假定已知反演区域内m行n列的布格重力异常值G(i，j)，(i=1，2，…m;j=1，2，…，n)，以及反演深度内已知深度的e个控制点深度h(i)，(i=1，2，…，e)，其对应的重力异常值为Gn(i)，(i=1，2，…，e)，将地下网格也划分成m行n列，坐标与重力异常值坐标一一对应.
图 3Fig. 3
图 3 重力异常值
Fig. 3 The gravity anomaly
(1)初始基准面及界面密度差的建立
设初始基准面的深度为m0，初始密度差为p0，根据无限水平板公式：
2πFp0(h(i)－m0)=Gn(i)(i=1，2，3，…，e)
F为万有引力常数，根据这e个方程，使用最小二乘法求取m0，p0.
(2)求取初始面密度矩阵 σ 1
第一次反演将m0作为压缩面，将代表起伏界面的小立方体压缩到m0上，使其成为面密度不均匀的物质面.每一块的面密度σ1(i，j)可看为(depth1(i，j)－m0)?p，depth1(i，j)为点(i，j)处起伏界面深度，则depth1(i，j)－m0为每个小立方体的高，也就是相对于压缩面的起伏，p为界面密度差.
被压缩到m0面上的小长方形(i，j)在重力观测点(k，l)所产生的重力异常值Δg为
其中为σ1(i，j)处压缩面深度.
对于观测面每一个重力异常值G(k，l)，可得
对于地面上所有重力异常值来说，可写出下面的线性方程组：
可写成矩阵方程 G=A?σ1，其中：
其中 G为已知重力异常值，A 通过计算可以获得，因此可以解出面密度矩阵 σ 1.
(3)求取初始模型 depth 1
由于初始的基准面和密度差是由控制点使用无限水平板公式求得的，当界面起伏较大时会有较大误差，因此，需要对基准面和界面密度差进行改进.
设改进后的基准面为m1，密度差为p1，使用e个控制点所对应的的面密度Ω1(1)，Ω1(2)，…，Ω1(e)和控制点深度建立方程组：
由最小二乘法求出m1和p1，则起伏界面深度depth1(i，j)可由各个小立方体的高度和基准面深度求得：depth1(i，j)=m1+σ1(i，j)/p1，求得初始模型 depth1.
迭代的目的就是改变基准面深度和界面密度差，使结果更加逼近真实情况.使用小立方体正演depth 1得到重力异常值 g 1，将实际观测异常 G减去g1作为重力异常值G2，由于此时的反演是在 depth1上的改进，因此以depth 1为压缩界面，m1为基准面，p1为密度差，重复上述计算得到面密度 σ 2.
e个控制点对应的的面密度Ω2(1)，Ω2(2)，…，Ω2(e)，此时计算新的基准面深度时所用的面密度 Ω是第一次的面密度Ω1加上改进的面密度Ω 2，即Ω(i)=Ω1(i)+Ω2(i)，(i=1，2，…，e)，设改进后的基准面深度为m2，界面密度差为p2，则有：
通过最小二乘法，得到m2与p2，以m2为新的基准面深度，p2为新的界面密度差，σ=σ1+σ2为面密度，求得第二次的界面深度模型 depth 2.
(5)继续迭代，直到到达一定的迭代次数或者反演结果符合以下两个标准
反演结果的正演重力异常值与原始重力异常值的偏差达到指定范围.
控制点与反演结果的差值达到指定范围.
3 模型试验
下面以三维密度界面上存在若干控制点为例进行反演实验.图 4为合成的起伏界面深度模型，基准面为30 km，密度差为0.5 g/cm3，黑点为控制点位置，图 5为按照上述方法迭代五次之后反演得到的界面起伏，反演结果与模型的平均偏差为2×10-6 km，标准差为0.005 km.迭代过程中每次得到的基准面深度和密度差如图 6、图 7所示，可看出第一次使用无限板公式获取基准面深度和密度差时，所得基准面和密度差与真实值相差较大，但是在此后的迭代中会逐步逼近真实值，反演结果与模型吻合很好.
图 4Fig. 4
图 4 深度模型(黑点为控制点位置，单位：km)
Fig. 4 Moho model(black dots are control points，unit:km)
图 5Fig. 5
图 5 反演结果(单位：km)
Fig. 5 Inversion result(unit:km)
图 6Fig. 6
图 6 每次迭代得到的基准面深度
Fig. 6 Reference plane in each iteration
图 7Fig. 7
图 7 每次迭代得到的密度差
Fig. 7 Density contrast in each iteration
南海范围广阔，构造演化复杂，其莫霍面可从20多公里变化到中央海盆的10 km左右，起伏较大，前人已做了大量工作(闫贫等，2002；；；；；卫小东等，2012；丘学林等，2012).本次选以南海(图 8)为例进行带控制点的三维密度界面反演，将地震结果作为控制点，控制点点距20 km，共594个，位置如图 9所示.反演结果如图 10所示，得出的莫霍基准面深度为23.625 km，密度差为0.5402 g/cm3，控制点反演结果与地震结果的平均偏差为-2.1×10-14 km，标准差为1.56 km，对反演结果使用小立方体正演得到的重力值与原始重力异常值的平均偏差为-0.0123 mGal，标准差为0.0284 mGal.同时使用Parker法对不同平均深度和密度差进行多次实验，选出最好的反演结果(图 11)，其平均偏差为2.23 km，标准偏差为2.56 km.此外，如图 12所示，抽取不同位置的两条地震剖面比较两种反演结果，可看出带控制点的三维界面反演方法反演结果优于传统Parker法.由偏差分析和剖面比较可知，带控制点的三维界面反演方法相比于传统的Parker法与控制点拟合的更好，结果更接近于真实值，而且相比于Parker法的多次反演选最好结果，带控制点的三维界面反演方法也更加简便.
图 8Fig. 8
图 8 研究区范围
Fig. 8 The research aera
图 9Fig. 9
图 9 南海重力异常
Fig. 9 The gravity anomaly of south China sea
图 10Fig. 10
图 10 带控制点的三维界面反演结果
Fig. 10 The result of three-dimensional interface inversion with control points
图 11Fig. 11
图 11 Parker反演结果
Fig. 11 The Parker inversion result
图 12Fig. 12
图 12 两种反演结果的比较
Fig. 12 The coparision of two inversion results
本文提出的带控制点的三维密度界面反演方法使用小立方体模型，避免了传统的Parker法向下延拓时产生的误差，在每次迭代过程中使用控制点约束基准面深度和密度差，使反演结果与控制点拟合较好.从本文的模型计算和实例反演来看，已知控制点数据的控制约束作用比较明显，从而避免了用不同参数分别计算，再将不同反演结果与控制点比较来选取最优结果，减少了多解性，提高了反演效率.
感谢中国科学院地质与地球物理研究所的张丽莉副研究员、徐亚副研究员、黄松副研究员、姬莉莉博士后、吕川川博士后、胡卫剑博士后、刘丽华博士和邢健博士，他们为本文写作提出了中肯的建议，在此一并致谢.同时感谢审稿专家对本文给出的意见和建议！
Kaufman A A. 1992. Geophysical Field Theory and Method, Part A: Gravitational, Electric, and Magnetic Fields[M]: Academic Press.
Parker R L. 1973. The rapid calculation of potential anomalies[J]..
Qiu X L, Ye S Y, Wu S M, et al. 2001. Crustal structure across the Xisha Trough, northwestern South China Sea[J]. .
Telford W M, Sheriff R E. 1990. Applied geophysics[M]: Cambridge university press.
Wang T K, Chen M K, Lee C S, et al. 2006. Seismic imaging of the transitional crust across the northeastern margin of the South China Sea[J]. .
Ao W, Zhao M H, Qiu X L, et al. 2012. Crustal structure of the northwest sub-basin of the South China Sea and its tectonic implication[J]. .
Dai W M. 2010. Find depth using the power spectrum of factors that directly[J]. Journal of Jilin University(Earth Science Edition), 40(S1): 17-20.
Fang J. 1999. Global crustal and lithospheric thinkness inversion by using satellite gravity data[J]. Crustal deformation and earthquake, 19(1): 26-31.
Fang J, Xu H Z. 1999. Three dimensional distribution of lithosphere density beneath the China and its adjacent region[J]. Progress in Geophys. (in Chinese), 14(2): 88-93.
Feng J, Meng X H, Chen Z X, et al. 2014. Review of the gravity density interface inversion[J]. .
Feng R. 1985. Crustal thickness and densities in the upper mantal beneath China --the results of three dimensional gravity inversion[J]. Acta seismologica sinica, 7(2): 143-157.
Hou Z C, Li B G. 1988. Estimating depths and thickness of rock stratum by using the power spectra of potential fields and their derivatives of different order[J]. Chinese J. Geophys. (in Chinese), 31(01): 90-98.
Huang J P, Fu R S, Xu P, et al. 2006. Inversion of gravity and topography data for the crust thickness of China and its adjacency [J]. Acta Seismologica Sinica, 28(3): 250-258.
Ke X P, Wang Y, Xu H Z. 2004. Inversion of Variable Density Model of Crust from Genetic Algorithms[J].Geomatics and Information Science of Wuhan University,29(11):981-984.
Ke X P,Wang Y,Xu H Z.2006.Interface depthsinversion of potential field data with variable density model[J].Progress in Geophys.(in Chinese),21(3):825-829.
Ke X P,Wang Y,Xu H Z.2006.Moho Depths Inversion of Qinghai-Tibet Plateau with Variable Density Model[J].Geomatics and Information Science of Wuhan University,31(4):289-292.
Li C L,Xie C L.1998.Calculating top and bottom depth effect of geological body by applying potential field power spectrum[J].,17(5):45-48.
Li J.2011.Dynamics of the continental margins of South China Sea:scientific experiments and research progresses[J]..
Liu Y L,Wang Q S.1977.Inversion of gravity data by use of amethod of "compressed mass plane"to estimate crustal structure[J].Chinese J.Geophys.(in Chinese),20(01):59-69.
Liu Y L,Zheng J C,Wu C Z.1987.Inversion of the three dimensional density discontinuity by use of a method of"compressed mass plane coefficient"based on the gravitydata[J].Chinese J.Geophys.(in Chinese),30(02):186-196.
Mu S M,Ye S S.1982.Proble mon depth estimation by power spectrum[J].(2):95-102.
Qiu X L,Zhao M H,Ao W,et al.2011.OBS survey and crustal structure of the Southwest Sub-basin and Nansha Block,South China Sea[J]..
Ruan A G,Niu X W,Qiu X L,et al.2011. A wide angle Ocean Bottom Seismometer profile across Liyue Bank,the southern margin of South China Sea[J]..
Wang B B,Hao T Y.2008.The inversion of two-dimensional mono-density interface with several known controlpoints[J].Progress in Geophys.(in Chinese),23(3):834-838.
Wei X D,Ruan A G,Zhao M H,et al.2011.A wide angle OBS profile across Dongsha Uplift and Chaoshan Depression in the midnorthern South China Sea[J]..
Xu K J,Li Y S.2003.The Power Spectrum Estimation Method Basedon Continuous Wavelet Transformation[J].Journal of Applied Sciences,21(2):157-160.
Yan J P,Cai J G,Li Z Z.2008.A power spectrum method based on wavelet transformation and its application in quantitatively dividing interfaces of sedimentary unit[J].China Offshore Oil and Gas,20(2):96-98.
YAN P,LIU H.2002.Analysis on Deep Crust Sounding Results in Northern Margin of South China Sea[J].Journal of Tropical Oceanography,21(2):1-12.
Zhang S,Meng X H.2013.Constra intinterface inversion with variable density model[J]..
敖威,赵明辉,丘学林,等.2012.南海西北次海盆及其邻区地壳结构和构造意义[J].地球科学:中国地质大学学报,37(4):779-790,doi:10.3799/dqkx.2012..
曾华霖.2005.重力场与重力勘探[M]:地质出版社.
戴伟铭.2010.利用功率谱直接求埋深的影响因素[J]..
方剑.1999.利用卫星重力资料反演地壳及岩石圈厚度[J]..
方剑,许厚泽.1999.中国及邻区岩石层密度三维结构[J]..
冯娟,孟小红,陈召曦,等.2014.重力密度界面反演方法研究进展[J].地球物理学进展,29(1):223-228,doi:10..
冯锐.1985.中国地壳厚度及上地幔密度分布(三维重力反演结果)[J]..
侯重初,李保国.1988.利用位场及其高阶导数的功率谱计算岩层的深度与厚度[J]..
黄建平,傅容珊,许萍,等.2006.利用重力和地形观测反演中国及邻区地壳厚度[J]..
柯小平,王勇,许厚泽.2004.用遗传算法反演地壳的变密度模型[J]..
柯小平,王勇,许厚泽.2006.基于变密度模型的位场界面反演[J]..
柯小平,王勇,许厚泽.2006.用变密度模型反演青藏高原的莫霍面深度[J]..
李成立,谢春临.1998.利用位场功率谱计算地质体顶底深度效果[J]..
李家彪.2011.南海大陆边缘动力学:科学实验与研究进展[J]..
刘元龙,王谦身.1977.用压缩质面法反演重力资料以估算地壳构造[J]..
刘元龙,郑建昌,武传珍.1987.利用重力资料反演三维密度界面的质面系数法[J]..
穆石敏,叶水盛.1982.功率谱估算理深中的若干间题[J].长春地质学院学报(2):95-102.
丘学林,赵明辉,敖威,等.2011.南海西南次海盆与南沙地块的OBS探测和地壳结构[J]..
阮爱国,牛雄伟,丘学林,等.2011.穿越南沙礼乐滩的海底地震仪广角地震试验[J]..
王贝贝,郝天珧.2008.具有已知深度点的二维单一密度界面的反演[J]..
卫小东,阮爱国,赵明辉,等.2011.穿越东沙隆起和潮汕坳陷的???广角地震剖面[J]..
徐科军,李永三.2003.基于连续小波变换的功率谱估计方法[J]..
闫建平,蔡进功,李尊芝.2008.基于小波变换的功率谱方法及其在沉积单元界面定量划分中的应用[J]..
阎贫,刘海龄.2002.南海北部陆缘地壳结构探测结果分析[J]..
张盛,孟小红.2013.约束变密度界面反演方法[J]..