2. 中国地质科学院地质研究所, 北京 100037
2. Institute of Geology, Chinese Academy of Geological Sciences, Beijing 100037, China
在地面观测的重力异常,是来自不同深度和不同尺度源体的位场的叠加.因此,如果不对位场进行适当的分解,则难以定量地反演计算场源的物性及几何参数.几十年来,人们已经认识到重力异常包含着尺度很大的成分,称为区域场.因此,采用圆滑、低通滤波或向上延拓等方法估算区域场,同时把重力场与区域场之差称为剩余重力异常.这种二分重力异常的程式一直沿用至今.但是,这种人为地二分重力异常的方法并不全面,“区域场”和“剩余异常”并没有严格的定义,只是相对的概念,从数学上也不存在重力场二分的根据.另一方面,由于重力场包含不同深度场源引力作用的叠加,用地表高精度重力观测资料研究区域密度结构非常困难.区域密度结构指不同深度水平面上的密度横向不均匀性,而不同深度场源引起的重力异常与地质体尺度大小有关.
自从1995年引入重力异常小波多尺度分解以来(侯遵泽和杨文采,1995;Yang and Hou, 1996,1997),重力场的分解结果有了明显的改观,大量实践结果显示了其在重力资料处理中应用的优越性.对中国大陆布格重力异常多尺度分析的结果表明(侯遵泽和杨文采,1997; 侯遵泽等,1998),中国重力场的四阶小波逼近主要反映Moho面的起伏和上地幔岩石密度的变化,四阶小波变换细节在青藏高原反映地壳下层岩石质量的不均匀性,而在东部反映Moho面的上隆,而一阶小波变换细节主要反映最上地壳密度的不均匀性,根据1~4阶重力异常小波变换细节,可以了解中国大陆各地区地壳密度自上而下的变化趋势和密度差异的空间分布.在国家高技术研究发展计划(863)海洋领域820主题的“东海地壳结构重磁地震综合反演”研究中,采用二维小波多尺度分解技术,对东海及其邻域自由空间重力异常进行了分解(高德章等,2000),结果揭示了浅层岩石密度的不均匀性和沉积层的厚度,反映了东海陆架沉积盆地和邻区沉积基底面以及Moho面的起伏;在国家自然科学基金“中国东部大陆及其边缘深部壳幔过程”项目研究中,应用小波重力场多尺度分解方法,对布格异常采用小波变换场分离技术得到深部重力异常,利用其他地球物理方法获得较可靠的先验信息作为约束,并采用球谐函数的级数解形式对中国东部及其邻域的深部重力异常进行联合反演,得到了中国东部及其邻域岩石圈底界面的反演结果及其展布图像(方盛明等,2001);在对澳洲大陆边缘某地区的布格重力异常的小波多尺度分析中(杨文采等,2001),发现小尺度与中尺度的异常成分具有不同的走向,根据位场频谱理论,可以计算出小波变换细节及逼近所对应的场源深度,得到反演的结果;在国家自然科学基金项目“华北岩石层密度时变的综合反演”中,利用小波多尺度分解程序对首都圈的重力时变资料进行了空间域的多尺度分解,从分解结果看,可有效地划分地表干扰因素(地下水活动)、地震孕育过程和大区域性影响的重力分量;在中国大陆科学钻探场地的重力异常多尺度分解中,也同样取得了重力场来源的更多的信息,其中,在大 别苏鲁地区岩石圈分区等项研究中,效果显著.
在研究重力场的二维离散小波变换多尺度分析方法及应用的基础上,笔者又结合功率谱分析和广义线性密度反演,将这三项方法集成,形成了多尺度密度反演方法,用以描述地壳的密度结构(侯遵泽和杨文采,2011).实践证明这种集成的思路是正确的,所形成的新的方法步骤对实际资料处理和解释是有效的.
本文根据华北地区的地质特征,利用多尺度密度反演方法反演出了对应三个埋深地层中心平面上的密度扰动,建立了区域的三维密度结构,经对比分析,浅层密度分布与华北地区油气藏的分布有着很好的对应关系. 2 重力异常多尺度密度反演方法 2.1 多尺度分解
设重力异常 Δg(x,y)=f(x,y)∈V20∈L2(R2) 根据小波多尺度分析原理(Mallat,1989; Daubechies,1990; 李世雄和刘家琦,1994),可以有以下分解:
这里, 其中若假定φ(·)和Φ(·)分别为V0的尺度函数和小波函数,则V20的尺度函数为 小波函数为 而系数C1,m1,m2,d11,m1,m2,d21,m1,m2和d31,m1,m2由下面的式子计算: 其中,(2)式可继续分解为
最终(1)式可写成 这里,p是正整数, 2.2 多层密度反演在完成重力异常小波变换多尺度分解之后,首先将重力资料小波细节进行适当的组合,再用频谱分析方法,求取它们的平均功率谱,以确定各层的平均深度,最后,运用广义线性反演方法求取目的层的密度分布(杨文采,1997).其中深度的求法是根据功率谱直线段斜率或它的渐近线在横轴上的切点计算出来的(杨文采等, 1978,1979a,1979b; 杨文采, 1985,1986; Cianciara and Marcar, 1975; Sybery,1972; Bhimasankakam et al., 1977).
由重力异常小波细节频谱分析得到的平均埋深只是对应“等效”地层的中心深度,为满足地质解释的需要,应该求出每一层的横向密度变化.
位场正演问题是求下述Poisson方程边值问题的解:
其中v(x,y,z=0)为地表位场值,ρ(x,y,z)为场源项,u(x,y,z)为下半空间任意点位场函数.由数学物理方法知(17)式的解(杨文采,1997) 在地表有反演问题就是给出v(x,y),由(19)式求等效层体密度 ρ(x,y,z).(19)式是第一类 Fredholm 积分方程,求解算法很不稳定.而在本文中z′为等效层中心埋深,已由功率谱斜率计算求出,所以,根据(19)式,用广义线性反演方法就比较容易求出目的层的密度分布ρ(x,y,z′).此后,假定各层的平均密度,求出该等效层的密度扰动. 3 华北重力场多尺度反演与密度结构
在国土资源部有关地球物理勘察院、中石化和中石油等单位的大力协助下,我们汇编了各种比例尺的高精度重力测量数据.华北区域重力及密度扰 动图的地理范围为102°E—123°E,32°30′N—42°30′N. 采用等角割圆锥投影,中央经线为110°,第一标准纬线为34°,第二标准纬线为44°,零点纬线为赤道(0°纬线).网格数据坐标范围:-743~1062 km,3944~5094 km,坐标网格间距为10 km×10 km.由于测量是在不同时间、不同仪器和不同精度要求下进行的,根据《重力测量规范》“五统一”的要求,对全区1 ∶ 100万比例尺编图取得的重力数据进行了整理.布格重力异常数据改算执行《区域重力调查规 范》DZ/T0082-2006.布格重力异常总精度是±2.0×10-5m·s-2. 布格重力异常示于图 1.图中的坐标是等角割圆锥投影坐标,单位为mm,1 mm代表实际地面上的2.5 km;图 3—6雷同.
我们认为重力场本身的特征决定重力异常应该 分解成几个层次,并根据其他地质资料了解每个层 次对应着什么地质含义.图 2为华北区域布格重力 异常的典型功率谱曲线,在此曲线上可切出反映不 同深度层的四个斜率直线段,表明重力异常最少应该分解成4个层次,即反映沉积盆地的异常、反映沉积盆地基底的异常、反映中地壳的异常和反映莫霍面起伏的异常.然后,探讨将小波变换多尺度分析和密度反演方法应用于华北区域布格重力异常的分解和盆地三维结构密度研究,并根据钻探了解每个层次对应着的地质含义.
首先进行重力异常的小波多尺度分解,分为一至五阶小波细节异常和五阶逼近. Δg=D1f+D2f+D3f+D4f+D5f+A5f. 之后,用重力异常小波细节(D1+D2)、(D3+D4)、D5和五阶小波逼近A5(见图 3—6),作二维Fourie变换,求它们的功率谱,进而求得各自对应的场源平均深度,结果示于表 1—4.表中的参数K、L为某个方向功率谱曲线的两个不同的波数,GfK、GfL为对应K、L的对数功率谱.根据功率谱直线段斜率计算出场源的深度H之后,统计分析求得均值与方差,便取得等效层的平均深度参数.
对比钻探资料和区域深反射地震剖面可知,图 3中一阶和二阶小波细节揭示浅层的场源,对应沉积盆地中的地层以及造山带或隆起区近地表侵入体.图 4中三阶和四阶小波细节揭示中层的场源,主要对应上地壳的结晶基底.图 5中五阶小波细节揭示深层的场源,主要对应上地壳底部与中地壳岩石密度的变化.最后,五阶小波逼近反映的主要是莫霍面起伏以及下地壳岩石密度变化(图 6).小波多尺 度分析提供了不同深度平面上重力分量具有不同走向和极性的平面分布图像.
为了进一步了解华北重力多尺度分析的物理含义,我们对小波细节及5阶小波逼近进行了等效层密度反演,结果参见图 7—10.密度反演采用的参数分述如下.
(1)浅层密度反演
利用图 3数据和频谱分析方法求得的平均深度3.13 km;该层平均密度取为2.60 g·cm-3,反演出浅层横向相对密度变化(为了作图方便,结果中的数据乘了1000,数值单位为:10-3 g·cm-3,下同).
(2)中层密度反演
利用图 4数据和频谱分析方法求得的平均深度16.80 km;该层平均密度取为2.67 g·cm-3,反演出中层横向相对密度变化.
(3)深层密度反演
利用图 5数据和频谱分析方法求得的平均深度23.05 km;该层平均密度取为2.80 g·cm-3,反演出深层横向相对密度变化.
(4)最深层密度反演
利用图 6数据和频谱分析方法求得的平均深度39.83 km;该层平均密度取为3.10 g·cm-3,反演出最深层横向相对密度变化.
通过等效层密度反演可见,华北区域浅层揭示 的是沉积盆地中的地层以及造山带近地表侵入体的密度扰动信息,中层小波细节和密度扰动揭示的是造山带或隆起区结晶基底岩石的密度扰动信息,深层小波细节揭示的是中地壳岩石的密度扰动信息,最深层小波逼近反映莫霍面起伏以及下地壳岩石到上地幔密度扰动的信息.由表 1可见,由频谱求 平均深度的方差较大,这是由于莫霍面起伏大造成的. 4 华北区密度结构地质含义的初步分析
不同深度等效层上的密度水平差异反映了区域地壳岩性和结构的横向不均匀性.为便于解释,我们把上述反演密度结果图与地面地质图(马丽芳等2006)相拼接,得到一组图件见图 7—10.
由图 7可见,华北区浅层(0~6 km)的层平均密度扰动在横向有明显的变化,变化最大幅度为-0.23~+0.27 g·cm-3.层平均密度扰动剧烈变化出现在秦岭和燕山—阴山造山带,造山带中的强低密度异常主要反映沿造山带展布的花岗岩体.二叠纪苏鲁碰撞造山带密度扰动并未出现明显的大幅度降低,说明这个地段已被克拉通化.在华北克拉通内部,浅层整体呈现为中高密度.最高密度带沿北京—保定—邯郸展布,与唐山—邢台地震带相吻合,可能与华北平原下方沿裂谷带分布的基性岩体有关.兰州—天水地块的浅层整体呈现为高密度,说明它是一个比较坚硬的构造单元.
图 8为华北区中层(6~20 km)的层平均密度扰动平面图,中层小波细节和密度扰动揭示的结晶基底岩石的密度扰动信息.由图可见,中层平均密度 强低密度异常也出现在秦岭和燕山—阴山造山带,反映沿造山带展布的花岗岩体,变化最大幅度为-0.08~+0.1 g·cm-3.除了秦岭和燕山—阴山带 之外,强低密度异常还在太行山脉及鄂尔多斯盆地腹心出现.鄂尔多斯盆地腹心的结晶基底出现低密度,也许与这里盆地沉积层很深和天然气聚集有关.在华北克拉通其他地区,中层结晶基底整体呈现为中高密度,最高密度带沿北京—保定一带展布,推测与早第三纪大陆裂谷作用基性岩侵入有关.二叠纪苏鲁碰撞造山带密度扰动接近于零.兰州—天水地块的中层整体呈现为高密度,也说明它是一个比较坚硬的构造单元.
图 9为华北区深层(20~30 km)的层平均密度 扰动平面图,深层小波细节揭示的是中地壳岩石 的密度扰动信息.由图可见,深层平均密度扰动在横向有一定变化,变化最大幅度为-0.03~+0.04 g·cm-3,说明中地壳岩石密度横向变化远小于上地壳岩石密度扰动.深层平均密度大幅度降低出现在秦岭和燕山—阴山带,以及太行山脉及鄂尔多斯盆地腹心.兰州—天水地块的高密度区已经缩小到几乎尖灭,说明它不是一个穿透全地壳的构造单元.
图 10为华北区最深层(>30 km)的层平均密度扰动平面图,最深层小波逼近反映莫霍面起伏以及下地壳到上地幔岩石密度扰动的信息.莫霍面从西向东抬升造成地壳减薄,最薄处在渤海湾,密度扰动幅度为-0.1 g·cm-3左右.
根据以往对地质构造与油气探测的理论与实践(杨文采,2009;杨文采等,2008)的认识,初步的对比分析表明,利用小波多尺度反演提取的密度信息对资源勘探有一定意义.例如,在图 11 中对比了华北 东部区浅层密度扰动与油气坳陷.华北东部的八个 主要的油气坳陷,包括渤海、冀东、塘沽、黄骅、任丘、东营、济阳和濮阳,都表现为低幅度的低密度异常区.当然,位于造山带中的低密度异常与油气坳陷无关,主要反映沿造山带展布的花岗岩体,上面己经说 明.从图 11不难看出,位于华北盆地中的低密度异常与油气坳陷对应很好,除上述八个主要的油气坳陷表现为低幅度的低密度异常以外,华北盆地中的临清坳陷也反映为低密度异常.临清坳陷虽然是重点勘探区,但尚未发现油气藏.对比分析表明,利用小波多尺度反演提取的密度信息对油气勘探的部署有一定指导意义.要说明的是,并非所有盆地中浅层的低密度异常都与油气坳陷对应,每一个油气盆地的密度结构都有其特点,在华北平原的地壳密度结构与鄂尔多斯的地壳密度结构就很不一样.关于鄂尔多斯的地壳密度结构和油气藏的关系以后将专门 拟文讨论.
对华北克拉通的地壳三维密度反演的实践,表明了利用小波多尺度分解和功率谱方法估算地层场源深度进而进行密度反演方法的有效性;对华北地区密度结构的地质分析,对比华北东部浅层密度扰动与油气藏分布,可见二者有一定的对应关系,这可为油气勘探提供有价值的参考.
[1] | Bhimasankakam V L S, Nagendra R, Seshagiri Rao S V. 1997. Interpretation of gravity anomalies due to finite inclined dikes using Fourier transformation. Geophysics, 42(1): 51-59. |
[2] | Cianciara B, Marcak H. 1975. Interpretation of gravity anomalies by means of local power spectra. Geophysical Prospecting, 24(2): 273-286. |
[3] | Daubechies I. 1990. The wavelet transform, time frequency localization and signal analysis. IEEE Transactions on Information Theory, 36(5): 961-1006. |
[4] | Fang S M, Yu Q F, Zhang X K. 2001. Characteristics of the lithospheric bottom interface and seismicity in eastern China and its vicinity. Chinese Journal of Geophysics (in Chinese), 44(1): 48-53. |
[5] | Gao D Z, Hou Z Z, Tang J. 2000. Multi-scale analysis of gravity anomalies on East China Sea and adjacent regions. Chinese Journal of Geophysics (in Chinese), 43(6): 842-849. |
[6] | Hou Z Z, Yang W C. 1995. An operational research on the wavelet analysis. Computing Techniques for Geophysical and Geochemical Exploration (in Chinese), 17(3): 1-9. |
[7] | Hou Z Z, Yang W C. 1997. Two-dimensional wavelet transform and multiscale analysis of the gravity field of China. Chinese Journal of Geophysics (in Chinese), 40(1): 85-95. |
[8] | Hou Z Z, Yang W C. 2011. Multi-scale inversion of density structure from gravity anomalies in Tarim Basin. Science China Earth Sciences, 54(3): 399-409. |
[9] | Hou Z Z, Yang W C, Liu J Q. 1998. Multi-scale inversion of density distribution of the Chinese crust. Chinese Journal of Geophysics (in Chinese), 41(5): 642-651. |
[10] | Li S X, Liu J Q. 1994. Wavelet Transform and Mathematics of Inversion (in Chinese). Beijing: Geological Publishing House. |
[11] | Ma L F, Ding X Z, Ju Y J. 2006. Chinese Geological Map (in Chinese). Beijing: Geological Publishing House. |
[12] | Mallat S G. 1989. Multifrequency Channel decomposition and wavelet models. IEEE Transactions on Acoustics, 37(12): 2091-211. |
[13] | Sybery F J R. 1972. A Fourier method for the regional-residual problem of potential fields. Geophysical Prospecting, 20(1): 47-75. |
[14] | Yang W C. 1985. The power spectrum analysis of the field data (A). Technol. Geophys. Geochem. Prospect (in Chinese), 7(3): 188-196. |
[15] | Yang W C. 1986. The power spectrum analysis of the field data (B). Technol. Geophys. Geochem. Prospect (in Chinese), 8(1): 14-22. |
[16] | Yang W C. 1987. Generalized inversion applying in field data calculation. Chinese Journal of Geophysics (in Chinese), 30:283-291. |
[17] | Yang W C. 1997. Theory and Methods in Geophysical Inversion (in Chinese). Beijing: Geological Publishing House. |
[18] | Yang W C. 2009. Tectonophysics in Paleo-Tethys domain. Chapter 7 (in Chinese). Beijing: Petroleum Industry Press. |
[19] | Yang W C, Guo A Y, Xie Y Q, et al. 1979a. Interpretation of gravity anomalies in frequency domain (A). Technol. Geophys. Geochem. Prospect (in Chinese), 79(1): 1-16. |
[20] | Yang W C, Guo A Y, Xie Y Q, et al. 1979b. Interpretation of gravity anomalies in frequency domain (B). Technol. Geophys. Geochem. Prospect (in Chinese), 79(2): 1-14. |
[21] | Yang W C, Guo A Y, Zhao J X, et al. 1978. Theory and methods for interpretation of gravity and magnetic anomalies in the frequency domain. Bull. Inst. Geophys. Geochem. Prospect (in Chinese), 78(2): 134-178. |
[22] | Yang W C, Hou Z Z. 1996. Partition of crustal density structure of Chinese continent. Continental Dynamics, 1(2):7-11. |
[23] | Yang W C, Hou Z Z. 1997. Geophysical analysis of main minerogenetic belts in China.//Proc. 30th Int'l. Geol. Congr., 20: 112-120. |
[24] | Yang W C, Hou Z Z. 1997. Multi-scale analysis of the gravity field and partition of crustal density in China. Academic Periodical Abstracts of China (in Chinese), 3(3):56-57. |
[25] | Yang W C, Shi Z Q, Hou Z Z, et al. 2001. Discrete wavelet transform for multiple decomposition of gravity anomalies. Chinese Journal of Geophysics (in Chinese), 44(4): 534-541. |
[26] | 高德章, 侯遵泽, 唐建. 2000. 东海及邻区重力异常多尺度分解. 地球物理学报, 43(6): 842-849. |
[27] | 方盛明, 余钦范, 张先康. 2001. 中国东部及其邻域岩石圈底界面特征及地震活动性. 地球物理学报, 44(1): 48-53. |
[28] | 侯遵泽, 杨文采. 1995. 小波分析应用研究. 物探化探计算技术, 7(3): 1-9. |
[29] | 侯遵泽, 杨文采. 1997. 中国重力异常的小波变换与多尺度分析. 地球物理学报, 40(1): 85-95. |
[30] | 侯遵泽, 杨文采, 刘家琦. 1998. 中国大陆地壳密度差异多尺度反演. 地球物理学报, 41(5): 642-651. |
[31] | 侯遵泽, 杨文采. 2011. 塔里木盆地多尺度重力场反演与密度结构. 中国科学: 地球科学, 41(1): 29-39. |
[32] | 李世雄, 刘家琦. 1994. 小波变换与反演数学基础. 北京: 地质出版社. |
[33] | 马丽芳, 丁孝忠, 剧远景. 2006. 中国地质图集. 北京: 地质出版社. |
[34] | 杨文采. 1985. 位场数据的频谱分析技术(上). 物化探计算技术, 7(3): 188-196. |
[35] | 杨文采. 1986. 位场数据的频谱分析技术(下). 物化探计算技术, 8(1): 14-22. |
[36] | 杨文采. 1997. 地球物理反演的理论与方法. 北京: 地质出版社. |
[37] | 杨文采. 2009. 古特提斯域大地构造物理学(第七章). 北京: 石油工业出版社. |
[38] | 杨文采, 郭爱缨,赵敬洗等. 1978. 重磁异常频率域解释的理论与方法. 物化探研究报导, 78(2): 134-178. |
[39] | 杨文采, 郭爱缨, 谢玉清等. 1979a. 重磁异常在频率域的解释方法(上). 物探化探计算技术, 79(1): 1-16. |
[40] | 杨文采, 郭爱缨, 谢玉清等. 1979b. 重磁异常在频率域的解释方法(下). 物化探电子计算技术, 79(2): 1-14. |
[41] | 杨文采,侯遵泽. 1997. 中国重力场的多尺度分析和地壳密度分区. 中国学术期刊文摘, 3(3):56-57. |
[42] | 杨文采, 施志群, 侯遵泽等. 2001. 离散小波变换与重力异常多重分解. 地球物理学报, 44(4): 534-541. |