利用航空磁测资料研究结晶基底表面起伏

古老地台结晶基底顶面起伏的信息，对于解决区域性和勘探性的多种地质问题十分重要，其中包括在基底本身及沉积盖层的继承性构造中寻找油气藏。

基底顶面起伏图通常依据钻探和地震资料编制，有时加入电测深信息。然而，拥有较可靠基底深度信息的点数通常十分有限，且分布极不均匀。在这种情况下，纯粹插值的结果极不可靠，因此往往借助研究区重力场资料填补空缺。但重力场仅部分由基底构造决定，相当一部分来自沉积盖层，这自然会给最终图件带来误差。因此，在由多位研究者编制此类图件的地区，图上基底高程通常相差100—300 m，且往往大得多。

引入航空磁测提供的补充信息，可显著提高结晶基底表面填图的效果，这一点几乎在航空磁测出现后立即就被认识到。异常磁场能量的很大部分——通常超过99%——正是来自结晶基底上部，因此在多数地台地区，航空磁测实际上是其直接填图方法。问题仅在于开发有效提取所需信息的方法，长期以来对此投入了大量研究。

目前实际应用中确定基底表面深度的主要方向，是分析描述异常场的函数奇异点。1956年，V. K. Ivanov首次提出根据振幅谱确定最上部奇异点深度[2]。他从理论上证明，高频段振幅谱渐近于指数函数，其衰减系数等于最上部奇异点的深度。在此基础上，开发并广泛应用了多个计算机程序[4]。国外A. Spector和F. Grant于1970年针对多块体模型提出了几乎相同的思路[12]，此后广泛用于解决多种实际问题。上世纪90年代，一系列研究指出，实际介质中异质块体的共同作用可产生其他频谱类型的异常场。一方面，J. Garcia-Abdeslem和J. Ness于1994年证明[7]，块体呈均匀分布和正态分布的介质，可能使Spector—Grant方法确定的深度系统性偏小。另一方面，M. Fedi、T. Quarta和A. De Santis于1997年指出，如果实际介质具有分形特征，则此类方法确定的深度反而可能系统性偏大[6]。需要强调，这类计算通常使用化极场，未考虑测量高程变化。

另一种广泛应用的方法称为欧拉反褶积，其基础是齐次函数的欧拉方程。美国的D. Thompson于1982年首次提出剖面形式的方法[13]；1990年，A. Reid领导的英国地球物理学家将其推广至三维情形[11]，随后纳入所有主要解释系统。该方法在滑动窗口内求解欧拉方程的线性反问题，确定一个指定类型、代表窗口内大部分场能量的奇异点坐标。其结果为点云，在平面上通常对应异常场的梯度带。现代解释系统手册将欧拉反褶积确定奇异点深度的精度估计为15%。尽管该方法会在地图上留下完全不含奇异点的大面积空白，欧拉反褶积仍常被推荐为快速估算异常源参数的主要方法[3]。

与此方法在一定程度上相近的是所谓Werner反褶积。S. Werner于1953年首次提出该方法以确定薄板参数[14]，此后经多位研究者修改。它实际上也是在滑动窗口内用有理分式近似场，归结为求解线性反问题。最近，M. Nabighian和R. Hansen指出，可将这两种方法统一用于三维模型[10]。
本文讨论的ROMGAS技术综合了该领域以往成果，是SIGMA-3D解释软件包的一部分。它依据滑动窗口内计算的振幅谱，确定描述磁异常的函数奇异点，具有以下主要特点：

1. 在正方形滑动窗口内确定奇异点深度，窗口尺寸应比预期基底深度大约一个数量级。计算直接使用异常场，无须预先化极。解释人员得到规则网格上的奇异点高程估计值，其集合实际上描述了奇异点集的包络面，通常称为主磁性界面（MMAS）。MMAS大多与结晶基底表面重合，但自然也可能偏离。一方面，若结晶基底上部由几乎无磁性的岩石组成，MMAS位于更深处，此时反映的是基底内磁性岩石的埋深。另一方面，若沉积盖层内存在高磁化强度的侵入体或厚层火山喷出—沉积岩系，MMAS则浅于基底表面。该技术的另一限制是无法确定狭窄坳拉谷内的深度。显然，当滑动窗口的平面尺寸远大于坳拉谷本身时，定位的最上部奇异点位于两侧，而不是底部。尽管有上述限制，ROMGAS技术在多种地质条件下仍显示出很高的效果。
2.若原始场给定于平面正方形网格节点，则按常规采用快速傅里叶变换计算二维频谱。
3. 若原始场给定于不规则网格，并考虑观测点相对大地水准面的高度，则用相应二维傅里叶级数的截断形式近似计算二维频谱。
根据不规则三维网格观测近似计算频谱的思路，似乎最早由R. Henderson和L. Cordell于1971年发表于广泛可得的文献[8]；他们指出，该方法是对W. James于1966年更早工作的扩展[9]。此后，多位地球物理学家在不同程序中采用了这一思路，包括本文一位作者开发的Singular集成系统[1]。目前，V. N. Strakhov和I. A. Kerimov正在积极开发基于积分表示法的原创逼近频谱分析变体[5]。
ROMGAS近似计算的谐波数量自然远少于观测点数。通常，30×30 km滑动窗口内包含1:50,000航空磁测的数万个观测点，而稳定确定表征MMAS的奇异点深度，仅需约一千个初始谐波。在确定这些谐波时，采用正则化方法求解线性反问题。

图1. 在钻孔上方30×30 km窗口内确定主磁性界面深度：а）异常磁场DT；б）地表起伏；в）归算至平均观测高度的场之振幅谱模对数与空间频率的交会图，以及斜率表征深度的回归线。

4. 可通过两种方式直接确定表征MMAS的奇异点深度。第一种方式先对计算的二维频谱作径向平均，再取对数，然后在程序自动选择的区段上用直线近似，其斜率即为所求深度的估计。该方法应用广泛，例如著名的Geosoft集成系统也采用它。
第二种方式直接拟合振幅谱模对数值及对应空间频率，并利用所有低频建立回归线。以下实际例子可说明这种方法。图1а和1б为莫斯科向斜盆地南部一口钻孔上方30×30 km窗口内的异常磁场DT及地表起伏。钻探资料表明，此处结晶基底顶面高程为-586 m。图1в显示归算至平均观测高度的场之振幅谱模对数与空间频率的交会图，所用谐波总数为900（30×30）。回归线斜率对应的奇异点高程为-599 m。需要注意，频谱最低频部分的点相对较少，因此不会像第一种预先径向平均的方法那样，对所求奇异点深度的确定造成显著干扰。
无论采用哪种方式，所得信息通常被视为深度的相关量，而非真实深度，因为研究区异质块体的共同影响及基底磁化强度的分形维数一般事先未知。然而，如果完全没有所研究表面的先验信息，按上述方法计算的深度本身可作为基底顶面形态的估计。
4. 按该方法完成全区计算后（自然不包括紧邻边界、宽度为半个滑动窗口的区域），为校正块体共同的集合影响，采用线性回归将结果与研究区现有的表面先验信息联系起来，这些信息主要来自钻探。由此可针对研究区特征作最终校准，包括校准基底模型的分形维数。

图2. 莫斯科向斜盆地内ROMGAS与钻探资料所得基底顶面高程的比较：а）394口钻孔资料与ROMGAS结果的交会图；б）钻探与ROMGAS确定的深度差值直方图。

ROMGAS程序已应用于多种地质环境，并显示出很高的效果。实际上，它已用于研究俄罗斯联邦约20%领土的结晶基底顶面起伏。MMAS确定精度通常不劣于10%，在有利条件下明显更高。
图2а为莫斯科向斜盆地394口钻孔资料与ROMGAS结果的交会图，图2б为相应深度差值的直方图。分析表明，除坳拉谷发育区外，该技术均获得了相当可靠的结果。在这些地区，50×50 km滑动窗口内最上部奇异点位于两侧，因而不能反映坳拉谷底部的真实深度。尽管如此，最终确定莫斯科向斜盆地结晶基底表面深度的均方根误差为267.5 m，钻探与ROMGAS深度的相关系数为0.940。需要指出，采用该技术的盆地大部分区域覆盖了1:200,000航空磁测。对于现代详细、高精度测量资料，MMAS深度确定精度还会更高。

图3. 根据鞑靼斯坦伏尔加河以西地区1:25,000航空磁测资料确定的结晶基底顶面起伏：а）标有深钻位置的DT等磁力线图；б）根据ROMGAS资料绘制的主磁性界面图；в）Geosoft系统采用结构指数0进行欧拉反褶积所得奇异点图；г）钻探与ROMGAS所得基底顶面高程交会图。

例如，利用1:25,000航空磁测资料研究了鞑靼斯坦伏尔加河以西地区的MMAS。图3а为标有深钻位置的原始DT等磁力线图，图3б为ROMGAS所得最终主磁性界面图。在该区，结晶基底顶面平均深度约2 km，深度确定的均方根误差为26.4 m。图3г比较了ROMGAS结果与钻探资料，两者相关系数为0.805。作为对比，图3в为同区采用结构指数0、即针对岩石接触面的欧拉反褶积结果，计算使用Geosoft系统。图中可见，欧拉反褶积所得奇异点分布在300 m至2.5 km的深度范围内，据此很难判断基底顶面的真实形态。实际上，在本例中，欧拉反褶积的结果远逊于ROMGAS。
因此，ROMGAS技术已在多种地质条件下证明其有效性，可推荐作为编制基底表面起伏图时对钻探和地震资料插值的基本方法。即使没有先验信息，ROMGAS也能对所求深度给出足够准确的估计。

参考文献
1. Blokh Yu.I., Kaplun D.V., Konyaev O.N. 在Singular集成系统中利用奇异点方法解释位场的能力 // Izvestiya vuzov. Geologiya i razvedka. 1993. 第6期. 第123-127页.
2. Ivanov V.K. 关于位势奇异点的分布 // Uspekhi matematicheskikh nauk. 1956. 第11卷. 第5(71)期. 第67-70页.
3. 1:200 000比例尺地质填图地球物理保障工作方法建议. 圣彼得堡: VIRG-Rudgeofizika. 2000. 240页.
4. Serkerov S.A. 重力勘探与磁法勘探中的谱分析. 莫斯科: Nedra. 1991. 279页.
5. Strakhov V.N., Kerimov I.A. 重力测量和磁力测量中谱分析的逼近实现 // 重力与磁异常解释理论的主要问题. 莫斯科: OIFZ RAN. 1999. 第183-206页.
6. Fedi M., Quarta T., De Santis A. Inherent power-law behavior of magnetic field power spectra from a Spector and Grant ensemble // Geophysics. 1997. v. 62. N 4. pp. 1143-1150.
7. Garcia-Abdeslem J.G., Ness G.E. Inversion of the power spectrum from magnetic anomalies // Geophysics. 1994. v. 59. N 3. pp.391-401.
8. Henderson R.G., Cordell L. Reduction of unevenly spaced potential field data to a horizontal plane by means of finite harmonic series // Geophysics.1971. v. 36. N 5. pp. 856-866.
9. James W.R. Fortran IV program using double Fourier series for surface fitting of irregularly spaced data // Kanzas Geol. Survey. Computer Cdontr. 1966. 5. 19 p.
10. Nabighian M.N., Hansen R.O. Unification of Euler and Werner deconvolution in three dimensions via the generalized Hilbert transform // Geophysics. 2001. v. 66. N6. pp. 1805-1810.
11. Reid A.B., Allsop J.M., Granser H., Millett A.J., Somerton I.W. Magnetic interpretation in three dimensions using Euler deconvolution // Geophysics. 1990. v. 55. N 1. pp. 80-91.
12. Spector A., Grant F.S. Statistical models for interpreting aeromagnetic data // Geophysics. 1970. v. 35. N 2. pp. 293-302.
13. Thompson D.T. EULDPH: A new technique for making computer-assisted depth estimates from magnetic data // Geophysics. 1982. v. 47. N1. pp. 31-37.
14. Werner S. Interpretation of magnetic anomalies at sheet-like bodies // Sveriges Geologiska Undersok. 1953. Ser. C. C. Arsbok 43. No. 6.