地下煤层自燃区磁性体三维反演

（一）频率域功率谱可视化深度计算方法

在频率域中通过计算场的功率谱求取磁性体深度的方法广泛用于大面积的航磁及地磁测量资料的解释中。该方法曾成功运用于下扬子航磁资料、辽东湾海域磁测资料以及康滇大陆古裂谷带等许多地区的磁测资料深度计算，取得了较好的应用效果。这些测区大都是大面积的区域性测区，局部异常数目多达几千甚至上万个，数据运算量非常大。因此以前主要采用计算机等间隔分块（即分窗口），窗口自动滑动方法来计算。这样，窗口能否包含目标异常将无法控制。大量的模型实验结果表明，当反演目标异常位于窗口范围内时可以获得较高的计算精度。采用等间隔窗口自动滑动进行计算，对精度造成一定的影响。

研究表明，如果能根据异常实际位置人工选择窗口的位置，则将大大提高对航磁异常深度计算的精度。人工选择窗口位置若不能在可视化环境下 *** 作，计算过程的反复性将大大降低反演速度。为了实现可视化反演目的，我们对原方法和程序进行了必要的改进，并将其移植到基于MapObject（简称MO）开发的GIS平台上。利用GIS对图形数据的显示、编辑、修改和对属性数据的编辑、存储、导入、导出等功能，实现了计算过程的可视化。

1.方法原理

关于在频率域用功率谱法进行深度计算的基本原理及公式导出过程在文献（郭志宏，1988）中有详尽的介绍，这里对方法原理进行简单介绍。

1）求磁性体顶面深度的功率谱法

将窗口内地下磁性体视为有限延伸垂直棱柱体，将磁异常化到地磁极，在较高频区域求角度平均径向功率谱有如下表达式：

地下煤层自燃遥感与地球物理探测技术

式中：F（s，Ψ）为磁异常傅里叶谱；s= ，u、υ为沿x及y轴水平方向的圆频率，s为径向圆频率，Ψ为s的方位角；α为常数；zt为顶深。

这样就可以通过取数组zt（s）的平均值，得到顶深zt。

2）求磁性体质心深度的似功率谱法

在较低频区域求角度平均径向似功率谱有如下表达式：

地下煤层自燃遥感与地球物理探测技术

则有：

地下煤层自燃遥感与地球物理探测技术

式中：G（s，Ψ）= （s，Ψ）；β为常数；其他符号同顶深公式中符号。同样可以通过取数组z0（s）的平均值，得到磁性体质心深度z0。

3）求磁性体底面深度

假设磁性体顶部和质心间的距离与质心和底界面的距离相等，则磁性体底界面深度zb为

zb= 2z0-zt

应当指出的是，上述公式都是从单一均匀磁化的直立棱柱体导出的。在大范围区域中的应用是以Spector和Graft提出的统计分析理论为前提的，即多个随机分布磁性体其频谱的总体平均和某一个形体的频谱相等。这一结论拓宽了频率域计算方法的应用领域，使得前述方法照样适用于实际叠加异常区域。

2.功率谱深度计算可视化方法

1）软件研制

通过调用MO提供的有关文件管理、图形编辑修改、属性数理等一系列功能，采用VC+ +语言开发基于地图应用组件MapObject的GIS软件，形成可视化的计算平台。

修改原深度计算FORTRAN程序，分别形成顶深计算、质心计算和底深计算三个功能模块（DLL动态连接库）。通过VC+ +接口，将这些模块连接到基于MO开发的可视化计算平台上，形成基于GIS的可视化深度计算程序。

2）关键技术和难点

软件包括以下两个关键技术：

（1）网格数据的图形显示和窗口数据提取。由于原始网格数据为航遥中心内部数据格式agb格式，因此要在屏幕上显示其图形，必须研制转换接口，并生成shapefile文件。根据异常等值线位置，用计算窗口（矩形框）在屏幕上提取窗口内的网格数据，坐标一致是关键。只有在等值线坐标、矩形窗口坐标和agb网格数据坐标一致的前提下，采用正确的数据搜索方法，才能提取正确的窗口计算网格数据。

（2）人机交互控制功率谱拟合切线的移动。由磁异常数据计算得到的对数功率谱显示到屏幕上后，必须拟合一条最佳的直线，才能根据其斜率得到较好的深度计算结果。以往采用的黑匣子 *** 作计算使得计算过程不能及时控制，需根据需要多次反复，过程相当烦琐。要使深度计算过程真正可视化，研究人工控制拟合切线的移动方法非常重要。通过多次尝试，最后采取由鼠标左右键控制拟合直线移动的方法，即通过鼠标左键使拟合直线向低频方向移动，一次移动一个点，每移动一点进行一次新的拟合计算，拟合结果、深度计算结果及拟合误差即时显示在屏幕上，所见即所得，直至获得最佳深度结果，方可结束计算。同理，通过鼠标右键可以向较低频方向移动，进行上述系列计算。

几个关键技术实现以后，如何在统一的GIS平台上将一系列的可视化过程有机地联系起来是整个软件实现过程中最大的难点。只有在各部分功能之间建立正确的连接关系，才能实现整个深度计算过程的可视化人机交互。

3）软件功能介绍

图4-3-1 可视化功率谱法深度计算软件界面

图4⁃3⁃1为软件界面。软件主要包括如下功能：

（1）“文件”菜单包括文件新建、文件打开、文件保存、另存为、退出等功能。可以对ARCGIS的shape格式文件进行 *** 作。

（2）“查看”菜单可以选择打开和关闭工具栏和菜单栏。

（3）“视图”菜单可以对屏幕进行放大、缩小、移动、全屏显示。

（4）“图层”菜单可以添加tif格式的图像文件和shape格式的矢量文件到主窗口，以便与深度计算结果进行对比分析。

（5）“图形”菜单可以在主窗口画点、画线、画椭圆、画矩形以及任意多边形，还可以对这些图元进行修改、删除等 *** 作。

以上功能为深度计算的可视化搭建了平台。

（6）“深度计算”菜单包括网格数据文件的输入、顶面深度的计算、质心深度的计算，以及通过顶面深度和质心深度的计算底面深度。

这部分是深度计算的核心。包括以下一系列功能：

①以图形形式显示磁异常网格数据；②根据异常位置及形态，人机交互选择计算窗口位置；③鼠标点击窗口自动提取窗口内网格数据；④屏幕显示功率谱计算结果；⑤根据窗口显示的深度和误差数据，人机交互移动功率谱拟合切线，以获得最佳深度计算结果；⑥将深度计算结果以点的形式显示在屏幕上，并将计算值作为计算点的属性保存。

3.乌达煤田火区深度计算及其结果分析

乌达煤火区航磁异常大都由煤层燃烧形成的烧变岩或正在燃烧的热源体引起，异常形态多为单个孤立异常，异常范围大都500m×500m左右，数据网格间距为30m。我们选择计算窗口为20点×20点。根据煤矿边界以及磁性体边界，参考遥感热红外异常区范围，对全区的航磁异常进行了顶面深度和中心深度的计算。计算结果见图4⁃3⁃2和图4⁃3⁃3。

在图4⁃3⁃2和图4⁃3⁃3中，黑色边界线是煤矿区边界。这个边界是利用高分辨率Quickbird影像资料，西、南、北 以乌达17＃层煤露头为界，东以乌达城区边界为界圈出的。该区域内煤层露头出露广泛，是引起煤层火灾的危险区域。另外，本章第一节航磁异常检查结果分析表明，乌达煤田航磁异常主要由煤火和地表民用设施干扰引起，主要煤火异常均位于此区域内。

图中红色边界线是由水平梯度极大值法计算出的航磁圈出的火区边界。该边界区域是由于煤火燃烧产生的磁性体分布范围，磁性体可能是烧变岩，也可能是正在燃烧火区。具体判断应结合地面物探资料及地质调查资料来进行。图中不同颜色的小圆点为可视化功率谱法计算得到的磁性体顶深结果（图4⁃3⁃2）和磁性体中心深度（火点深度）计算结果（图4⁃3⁃3）。

结合物性特征分析结果，我们认为磁性体顶部深度就是煤火产生的高于500℃热源体的顶部深度。该热源体或者是现在已经冷却，即为烧变岩体，或者是正在燃烧，即为热源体，其内部存在着火点。

从图4⁃3⁃2可以看出，乌达煤田热源体顶深在0～70m之间，其中20～40m深度点居多，占所有计算点的51%。主要分布在乌达煤田北部，Ⅲ号火区周围。0～20m深度较浅的深度点主要位于测区东部偏南，Ⅺ号火区附近。40～70m深度点主要位于矿区西南部，Ⅷ号火区。

我们认为磁性体中心深度可能就是煤层着火点的深度。主要是基于以下考虑：假设着火点位于热源体内部，则根据热传导理论，若火点周围岩石热导率相等，等温的热源体就会围绕高温点呈辐射状分布，温度由里向外逐渐降低，着火点只能位于热源体中心或其附近，并且可以认为热源体中心即为磁性体中心。因此，采用能谱法计算磁性体质心的深度就可以是煤层着火点深度。

从图4⁃3⁃3可以看出，乌达煤田着火点深度在30～90m之间，其中50～70m深度点居多，占所有计算点的57%。着火点主要分布在乌达煤田北部，Ⅲ号火区和Ⅺ号火区周围。70～90m深度点主要位于矿区西南部，Ⅷ号火区。

通过对乌达煤田航磁资料实际处理与解释，我们认为人机交互频率域功率谱可视化深度计算方法，由于对窗口位置的选择完全在可视化人机交互状态下进行，深度计算结果在屏幕上即时显示，并以属性形式存储，大大提高了计算过程的可 *** 作性和计算速度。测区内异常范围大小毕竟存在差异，今后需要改进单一窗口大小的模式。应根据测区异常范围大小的不同，选择不同大小的窗口进行计算，将会进一步提高计算精度。

图4-3-2 内蒙古乌达煤田可视化功率谱法计算磁性体顶深结果

图4-3-3 内蒙古乌达煤田可视化功率谱法计算磁性体中心深度结果

（二）磁异常三维反演——人机联做方法

根据第二章的分析，地下煤层自燃到一定的阶段后，煤层及其围岩的磁性增强，可引起磁异常。由此可以推断磁异常主要由地下煤层自燃形成的磁性体（即燃烧体，主要为烧变岩）引起。因此，通过对磁异常的定量反演可以确定地下燃烧体（主要是烧变岩体）的埋深、大小及产状等要素，通过综合解释甚至推断煤层燃烧趋势。

作者选择乌达煤田Ⅷ号火区（局部）航磁和地磁实测资料，采用人机联做的磁异常三维反演方法对该区磁异常进行定量反演。首先根据磁异常特征，地质工作探明的该区煤层及煤火的分布特征初步设置磁源体模型，通过正演计算模型体场值，使其与实测场相拟合来确定地下燃烧体（主要是烧变岩体）的埋深、大小及产状等要素。

使用的软件界面如图4⁃3⁃4。图中右侧是乌达Ⅷ号火区地面磁异常图，矩形体为地下磁源体（烧变岩体）在地表的投影，异常图中间的水平线为当前拟合剖面位置。左下部模型体输入及编辑窗口，左上部为原始测量异常与理论异常对比窗口。

图4-3-4 三维重磁异常正反演解释软件使用界面

1.航磁异常反演

反演采用的拟合参数为：地磁场强度55000 nT；地磁场倾角58.7°；地磁场偏角-2.4°。图4⁃3⁃5、表4⁃3⁃1分别为乌达煤田Ⅷ号火区（局部）航磁异常反演地下磁性燃烧体模型和模型体参数。

图4⁃3⁃6为航磁实测异常（上图）与理论异常（下图）对比图。由图可见，实测异常与理论异常非常相似，说明拟合效果较好。

2.地磁异常反演

图4⁃3⁃7、表4⁃3⁃2分别为乌达煤田Ⅷ号火区（局部）地磁异常反演地下煤火磁性燃烧体模型和模型体参数，地面磁法测量范围见图4⁃3⁃5所标识的“地面详查区”，主要对应航磁异常反演的12、13、14号模型体分布区。图4⁃3⁃8为地磁实测异常（上图）与理论异常（下图）对比图。由图可见，实测异常与理论异常非常相似，同样说明拟合效果很好。

图4-3-5 乌达煤田Ⅷ号火区（局部）航磁异常反演模型体平面投影及编号

表4-3-1 乌达煤田Ⅷ号火区（局部）航磁异常反演模型体顶底面埋深及磁化强度

表4-3-2 乌达煤田Ⅷ号煤火区（局部）地磁异常反演模型体参数

图4-3-6 乌达Ⅷ号火区（局部）航磁实测异常（左图）与理论异常（右图）对比图

图4-3-7 乌达煤田Ⅷ号煤火区（局部）地磁异常反演模型体平面图

3.航磁异常与地磁异常反演结果分析

下面根据拟合的结果对各磁源体模型的地质意义及各磁异常反映的地下煤层燃烧情况推断解释。图4⁃3⁃9、图4⁃3⁃10分别为航磁和地磁异常三维解释立体图效果。从图中可以看出磁性体产状等分布特征，及其与磁异常之间的关系。

由磁异常反演结果与地质资料的综合分析得出以下结果。

（1）根据航磁异常反演各模型体的分布，可将模型体分为南北两部分共4个模型组合体组成，见表4⁃3⁃3。模型体的顶深为15～23.9m，厚度31.1～43m。根据该区Ⅱ号勘探线资料，该区主要分布的为9＃、10＃、12＃煤层，主要燃烧体为9＃、10＃煤层，其中9＃煤层埋深约19m，3个煤层的总厚度为35～45m。由此可见航磁异常反演的深度和厚度与地质情况非常吻合。其中航磁异常反演地下煤火燃烧体深度最大误差为25.8%，并且推断该区12＃煤层也在燃烧。

（2）地磁异常反演结果与航磁异常基本吻合，反映更细致，说明磁性体确系地下煤层燃烧形成的烧变岩引起。航磁异常反演模型体与地磁异常反演结果对比见表4⁃3⁃4。由表可见，两者计算深度误差小于20%。

图4-3-8 乌达煤田Ⅷ号煤火区（局部）实测地磁异常（上图）与理论异常（下图，对比图）

（3）航磁、地磁异常反演的磁性体顶深、底深和厚度与地质勘探确定的燃烧的地下煤层分布特征吻合。

（4）磁性体的形态和磁化强度变化特征可大致反映地下煤火的发展趋势，一般规律如下。

a.磁性体埋藏浅、厚度大、磁化强度大——表明地下煤层燃烧规模大、时间长，如1号燃烧体。

b.磁性体埋深大、厚度小、磁化强度小——表明地下煤层燃烧规模小、时间短，正在燃烧（间接反映燃烧方向），如4号燃烧体。

图4-3-9 乌达Ⅷ号火区（局部）航空物探异常反演地下煤火燃烧体可视化模型

图4-3-10 乌达煤田Ⅷ号火区（局部）地磁异常反演地下煤火燃烧体可视化模型

表4-3-3 乌达煤田Ⅷ号火区（局部）航磁推断的燃烧体特征与地质勘探结果比较

表4-3-4 乌达煤田Ⅷ号火区（局部）航磁异常与地磁异常反演燃烧体深度结果比较

（三）磁异常三维反演——三维成像方法

磁异常三维反演是使用中国地质大学安玉林教授开发研制的“复杂地形下局部重磁场源全方位成像系统”完成的。该系统是基于复坐标系（二维）和球坐标系（三维）建立起来的。可以直接对曲面上的测量数据进行反演，能定量地反演出异常体的位置（包括埋深和水平投影位置）和边界，并能给出三度体的三维立体形态。适合于中、高山区重、磁资料直接定量解释。

系统主界面及三度体磁矩中心反演功能窗口如图4⁃3⁃11。

图4-3-11 全方位成像系统主界面及三度体磁矩中心反演功能窗口

图4⁃3⁃12为乌达地区航磁ΔT10号异常等值线，图4⁃3⁃13为该异常对应的地形等高线。

从图中可以看出10号异常受旁侧异常及背景场的影响稍有变形，尤其是伴生负异常变形较大。图4⁃3⁃12中部椭圆轮廓线为本次反演场源的大致位置及范围边界线。测区中部地磁场偏角为2.6°、倾角为56.8°。测区地形较平坦，飞行高度基本与地形平行。取航测GPS数据和地形高程数据作为确定反演结果的重要约束条件。

图4⁃3⁃14是反演后提取出的主要异常等值线，与实测异常很接近。将不规则的二次背景剔除后，等值线更为规则，相对圆滑。

图4-3-12 异常区ΔT10号异常等值线图

图4-3-13 DTM等值线图的地形等值线

图4-3-14 异常区ΔT10号异常理论等值线图

图4⁃3⁃15是反演场源立体图。从反演结果可以看到场源基本上为等轴状，沿东西方向略长。经反演得到的场源中心坐标为（300m，320m，80m），磁性体规模大致为110m×90m×60m。磁性体顶深约20m。

图4-3-15 反演场源立体图

（四）水平梯度极大值法

水平梯度极大值法利用重力水平梯度异常在垂直物性边界的正上方取极大值的特点来确定边界位置。对于磁异常则不是直接应用，需要经过三步处理。

第一步，做伪重力变换。利用重力位与磁位的泊松关系，把磁异常换算为伪重力异常（磁源重力异常）。

第二步，计算伪重力异常水平梯度的模（强度）。浅层地质体产生的重力异常有较大的水平梯度，在地质体边界附近水平梯度取极值。因此，经过两步处理后得到的水平梯度模可用来确定物性边界的位置。

第三步，寻找和确定重力水平梯度模的极大值及其位置，并且在找到的极大值位置处绘出标志以表示边界位置。

图4-3-16 水平梯度极大值法圈定火区边界实例

理论及模型试验结果（郭志宏，2004）表明，水平梯度极大值法获得较好效果的前提是局部异常简单，受区域磁场干扰小；磁性体磁化方向尽可能与地磁场方向一致，满足重力位与磁位的泊松关系成立的前提条件。纵观乌达地区磁场特征特别是煤火区磁场特征，不难看出区内局部磁异常孤立，互不干扰，区域磁场极为平缓；多数正异常伴生负异常位于北侧，且负异常极大值小于正异常极大值，符合高磁纬度区局部异常特征，亦即多数磁性体磁化方向与地磁场方向一致。因此，水平梯度极大值法在乌达地区具良好的应用前提。由于方法本身局限，当物性体边界不垂直时，水平梯度模的极大值就会偏离边界位置；且对于很薄的地质体，水平梯度极大值也不能反映它们的位置。水平梯度极大值偏离地质体边界位置的距离，主要决定于边界顶端的埋深和边界的倾角。埋深越大，倾角越小，偏距也越大。因此，在使用本方法计算结果时应结合地质资料详细分析，方能准确勾绘出磁性体边界。

图4⁃3⁃16给出了水平梯度极大值法在乌达煤火区圈定磁性岩体边界应用效果。该航磁异常形态规整，正异常位于南侧，负异常位于北侧，并有向南包围南侧正异常的趋势，类似球状磁性体理论异常特征。水平梯度方法计算出的边界点围绕异常中心圈成一个封闭，并与地面磁异常检查剖面正异常范围吻合很好。从野外地质观察及火点分布情况分析，此处煤火是从南部引起，正逐渐向北发展。

（五）航磁异常深度的快速计算（切线法、外奎尔法等）

本次航磁资料测量剖面点距一般在1m左右，异常反映很精细。实际工作中，使用了Airprobe软件中自动计算深度的切线法、外奎尔法等方法，给出地下磁性体顶深初值，作为反演拟合的初始值。

（六）小结

采用频率域功率谱可视化深度计算方法快速确定煤火点埋深；采用水平梯度模方法快速圈定烧变岩体范围；可视化二维、三维磁反演方法能确定烧变岩的埋深、大小及产状；根据磁性体的形态和磁化强度变化特征可大致反映地下煤火的发展趋势。

Fortran源自于“公式翻译”（英语：FormulaTranslation）的缩写，是一种编程语言。

它是世界上最早出现的计算机高级程序设计语言，广泛应用于科学和工程计算领域。FORTRAN语言以其特有的功能在数值、科学和工程计算领域发挥着重要作用。

随着FORTRAN语言版本的不断更新和变化，语言不兼容性问题日益突出，语言标准化工作被提上了日程。

1962年5月：美国标准化协会(简称ANSI)着手进行FORTRAN语言标准化的研究工作。

1966年：ANSI正式公布了两个标准文本：美国国家标准FORTRAN(ANSI X3.9-1966)和美国国家标准基本FORTRAN(ANSI X3.10-1966)，前者相当于FORTRAN Ⅳ，后者相当于FORTRANⅡ。基本FORTRAN是美国国家标准FORTRAN的一个子集，从而实现了语言的向下兼容，初步解决了语言的兼容性问题。

---