基于模糊数学的露天铁矿采空区地球物理探测评价研究

贾三石 付建飞 王恩德 郭 凯 门业凯

（1.东北大学秦皇岛分校资源与材料学院，河北秦皇岛066004；2.东北大学资源与土木工程学院，辽宁沈阳110819）

采空区作为一种矿产资源开采后的遗留产物，大量存在于金属和非金属矿山［1］，特别是在铁矿开采历史悠久的鞍本地区表现得最为典型和突出。该地区既有矿业整合开发前偷采滥挖、采富弃贫和无序开采留下的采空区，又有日伪时期掠夺式开采留下的采空区，且它们数量庞大、规模不等和赋存空间不明。但从20 世纪末以来，随着资源开发利用的科学技术水平不断提高，主要是采掘运输能力和选矿工艺水平的大幅提升，矿产资源整合露天开采成为一种趋势，具有开采量大、矿石回收率高、生产安全、效益高、成本低等优点［2］。由于历史原因和资源开采方式的转变，露天矿山多由散乱小规模井下开采变为整合大规模露天开采，导致采场多处于不明采空区的发育区。而采空区顶板厚度会随着持续的开采扰动不断减小，直至突然塌陷形成深坑，轻则造成采掘设备坠毁，重则造成生产运输车辆和人员损失。且随着露天采矿活动的持续进行，采坑深度不断加大，形成了高陡边坡作用下的矿山开采应力场不稳多变，导致矿区采空区安全问题日益突出，不断诱发坍塌、变形移位或含水空区突水等矿山次生地质灾害，成为时刻威胁露天金属矿山生产安全的主要难题［3-6］。

为解决上述遗留采空区诱发的生产安全技术难题，近年来，科研院所和生产单位采用了多种技术方法和手段对采空区进行探测，如人工调查法、地球物理探测、工程钻探、放射性氡气测量、三维激光扫描和水下声呐测量等［7］，尤以地球物理探测方法最受青睐。该类方法具有种类多、成本低、方便快捷高效、无损非接触的地质CT 探测成像等特点，属于超前的“绿色探测”范畴，适用于对面积性采场区域内的采空区异常进行快速扫面筛查和定位，确保后期高效的钻探验证和三维激光或水下声呐精准测量工作有效开展。但地球物理探测成果影像显示异常具有间接性、多解性和圈定指标的不确定性，这使得探测成果的解译和应用需要有丰富经验的地质与地球物理专家级的人才联合参与。现阶段，在露天矿山快速持续不断的平台作业生产中，每天采集的大量多方法物探影像数据难以得到具有丰富经验的专家级人才及时处理，致使采矿作业平台采空区异常得不到及时、准确圈定，造成钻探验证和采空区安全处理工作严重滞后，在制约矿山正常生产的同时威胁了矿山安全。针对该类应用地球科学问题仅依靠专家经验解决的局限性，随着计算机科学技术和数学地质建模技术的发展，大量学者提出了多种基于专家经验和统计分析驱动的数据分析处理判别模型，实现了大部分地球科学数据成果的定性和半定量化分析研究与应用普及［8-9］。特别是模糊数学理论的出现，使得地球物理异常数据的处理和定量解译成为可能［10］。采空区地球物理异常作为矿致地球物理异常之上的二次人造异常，更具多样性、变异性和复杂性，造成的不确定性具有模糊性，尤其是采空区地球物理异常的最终解译成果只能反映异常圈定区内有无采空区存在，而这类“有、“无”之间就是一种模糊关系。

为解决采空区地球物理异常圈定的不确定性和专家经验圈定的滞后性，便于及时开展钻探验证和安全处理工作，基于多年的BIF 型铁矿山采空区探测成果，特别是采空区类型和地球物理特征研究成果的总结，应用模糊数学理论和模型评价方法，基于GIS 平台构建采空区空间地球物理异常影像快速处理与定量评价模型，以期实现金属矿山不明采空区地球物理异常圈定的模型化和智能化。

1 露天铁矿采空区类型和地球物理识别特征

1.1 露天铁矿采空区类型

采空区作为一种地表以下一定深度矿产资源开采后形成的空间，目前主要存在如下类型划分思路［1］。

（1）基于资源开采时间的采空区分类。根据采空区的形成时间，采空区可以分为3类：第1类为历史时期形成的老采空区，第2类为目前正在开采形成的现采空区，第3类为未来计划开采区域形成的采空区。

（2）基于开采矿产资源类型的采空区分类。根据开采矿产资源类型的不同，采空区可以分为金属矿山采空区和非金属矿山采空区两大类，其中金属矿山采空区依据采矿方法又可分为空场法采空区、崩落法采空区和充填法采空区。

（3）基于采空区规模大小分类。可以简单分为独立矿山采空区和整体矿山采空区两类。其中，独立矿山采空区按体积可以分为0.1～1 万m3、1～3 万m3和3～10万m33类，整体矿山采空区按体积规模可以分为5类，分别为＜50万m3、50～100万m3、100～500万m3、500～1 000 万m3和＞1 000 万m3。此外，还可以划分为小规模的巷道式采空区和大规模的采场式采空区。

（4）基于采空区埋藏深度分类。根据采空区埋藏深度，可以将采空区分为两类：第1 类为浅部采空区，埋藏深度在60 m 以内，是矿山开采和工程施工必须及时探明和安全处理的采空区；第2类为深埋采空区，埋藏深度大于60 m，需要进行安全评价、监测与处理的采空区。

（5）基于采空区赋存状态分类。根据采空区的赋存状态，既可以分为已知采空区和未知采空区，也可以分为已塌陷采空区和未塌陷采空区。未知采空区，即为不明采空区或隐伏采空区，是地球物理探测和矿山安全防治的重点对象。

（6）基于采空区充填介质类型的物性分类。不论何种采空区，一旦形成，在未被人为充填和完全塌陷之前，总保持一定的形态。对于其空间内的充填介质，除了常见的围岩和矿石碎块外，对其物性影响最大的介质即为四处流动的空气和水。依据对采空区物性差异的影响程度，一般可以将不同充填介质类型的采空区划分为两大类：第1类为空气充填型采空区，其内充填介质以空气为主；第2类为水充填型采空区，其内充填介质以水为主。该分类方式也是不明采空区可以进行地球物理探测的前提条件和理论基础。

1.2 采空区地球物理识别特征

不同类型的采空区具有不同的地球物理特征，而不同地球物理特征的采空区来自于不同物性差异对采空区类型及赋存围岩类型的可区分程度，这种区分程度也是选择不同地球物理探测方法的理论基础。基于物性差异的认识，本研究将采空区分为两大类，即空气充填型采空区和水充填型采空区。

（1）空气充填型采空区。空气充填型采空区以一种极低密度绝缘体空气为充填介质。该类采空区与赋存空间的岩矿石相比，具有极低密度、高电阻率和低电导率的地球物理识别特征，可以采用高密度电阻率法、瞬变电磁法和地震映像法进行组合探测。相关研究表明，空气充填型采空区异常在高密度电阻率法和瞬变电磁法数据影像上均表现为高阻值异常圈闭空间，且阻值从异常边界到异常核心逐渐增大，呈正向递增梯度变化趋势［11-12］；在对地震映像法探测数据进行解译时，可通过抽取最佳偏移距道数据来快速定位采空区异常位置，具体表现为地震波形图同相轴的突然增减和错断，且主要以同相轴的错断增多为采空区异常存在特征［13-14］。

（2）水充填型采空区。水充填型采空区以一种来源于矿山水循环系统、富含多种矿物质元素的水为充填介质，对外显示为低密度的良导体。该类采空区与赋存空间的岩矿石相比，具有低密度、低电阻率和高电导率的地球物理识别特征，也可以选用分辨率高的高密度电阻率法、瞬变电磁法和地震映像法进行组合探测。相关研究表明，水充填型采空区异常在高密度电阻率法和瞬变电磁法数据影像上均表现为低阻值异常圈闭空间，且阻值从异常边界到异常核心逐渐减小，呈负向递减梯度变化趋势［7，11］；对于地震映像法探测数据，水充填型采空区与空气充填型采空区的地球物理异常识别特征一致［13-14］。

2 模糊数学理论评价方法

模糊数学评价方法作为一种综合评价方法，主要对受到多个不确定性因素影响的目标进行整体性评价［15-16］，这就需要首先确定评价对象的影响因素。该类因素对评价结果的影响程度不同，而需要确定该类因素对上层评价目标的权重，并经过模糊变换得到该类因素对评价总目标的隶属度，最终完成对研究目标的综合评价［17］。假设2 个有限集W=，且W 与K 的模糊关系可表示为K=WY。其中，Y 是W 到K 的一个模糊变换；W 为影响因素权重集；K 为模糊综合评价成果集，表示最终综合评价结果；Y 为模糊判别矩阵，描述因素权重集W和评价成果集K之间的模糊关系。

本研究考虑到物探数据的异常不确定性和间接性，为最大限度减少人为干预对评价结果的影响，故采用全程模糊综合评价模型。首先，依据前期研究成果选出采空区异常评价指标因素集，通过数据重构提取每个评价指标因素，形成评价指标数据影像集；然后，对于获取的每个评价指标数据进行自然间断分级和模糊赋权分类，形成评价指标概率分布影像集；最后，通过模糊综合评价模型判别采空区综合地球物理异常，形成采空区异常分布概率影像图。

3 露天铁矿采空区异常评价指标

根据大量的BIF 型铁矿采空区地球物理探测研究和钻探验证成果经验总结［7，11-12，18］，优化组合出集高密度电阻率法、瞬变电磁法和地震映像法于一体的采空区异常探测方法，并利用方法中采空区与周围岩矿石物性差异最大的一组物性参数组合来评价采空区的存在与否，如电法中的电阻率和电磁法中的电导率等。此外，依据探测数据中的隐含数据信息和不明采空区异常的空间赋存规律，如空气充填型采空区的高阻异常圈闭性和地表的高阻异常不可能为采空区异常等，建立了如图1所示的采空区地球物理异常评价指标体系。对于具体评价指标的数据特征和数据重构及其综合评价，本研究以空气充填型采空区为例进行分析。

3.1 高密度电阻率法数据特征和数据重构

利用高密度电阻率法探测采空区，主要是利用采空区与周围岩矿石的电阻率差异，具体表现在电阻率数据二维影像上显示为相对高阻值块区和相对低阻值块区（图2（a）和图2（b））。对于空气充填型采空区，高密度电阻率法数据特征主要表现为高阻值异常圈闭区，从异常区的边界到核心，阻值逐渐增高，表现为正向梯度变化。考虑到空气充填型采空区具有正向梯度变化特征，可以采用离散傅里叶变换垂向求导获得正向梯度值圈闭区域［19］（图2（c）），即评价指标体系中的采空区异常梯度矢量特征G2T。对于电阻率高低的评价指标G1R，本研究按照传统地球科学异常评价方法，在无遗漏采空区异常标准指导下，以每幅影像的电阻率平均值为背景值，高于背景值的区块为正异常→高阻值圈闭区→空气充填型采空区异常，低于背景值的区块为负异常→低阻值圈闭区→水充填型采空区异常。此外，探测区域的地表接地条件、不明采空区的赋存位置和探测方法数据的空间分布等均可对采空区异常圈定的可靠性产生极大影响，如地表的正向梯度异常和高阻值圈闭难以成为采空区异常。为此，本研究提出了采空区异常埋深条件G3H作为评价指标，评价电阻率法探测数据点的空间分布可信度。

上述3 个评价指标G1R、G2T和G3H，由于数据格式和数据取值范围存在很大差异，需要进行标准化的数据重构处理，主要是对各个评标指标涉及的参数值区间进行分类赋值。为使各个类所涵盖的数值区间差异最大化，并体现出固有的自然分组圈闭特征，本研究采用了自然间断点10 级分类法，形成了高密度电阻率法的3 个采空区异常评价指标影像分级图（图2（d）、图2（e）和图2（f））。

3.2 瞬变电磁法数据特征和数据重构

瞬变电磁法通过测量地下介质由于一次发射脉冲电磁场产生的二次感应电磁场随时间变化的衰减特征，来获取地下不均匀体的导电性能和位置（图3（a）），其反演结果经过处理后如图3（b）所示，但由于瞬变电磁法存在关断时间效应，其浅部为探测盲区。空气充填型采空区在瞬变电磁法成果图中也表现为高阻值异常圈闭特征，从圈闭异常的边缘到核心，电阻率逐渐增高，同样可以采用离散傅里叶变换垂向求导获得正向梯度值圈闭区域（图3（c）），即评价指标体系中的采空区异常梯度矢量特征S2T。对于电阻率评价指标S1R，可采用探测成果电阻率的平均值为背景值，高于背景值的区块即可判定为空气充填型采空区异常。由于瞬变电磁法存在探测盲区，且随着探测信号的衰减造成深部成果的可信度较低，故本研究提出采空区异常埋深条件S3H来综合评价探测数据点的空间异常可信度。

为使瞬变电磁法的探测成果可与高密度电阻率法的探测成果进行模糊综合评价，需对其3个评价指标（S1R、S2T和S3H）进行标准化的数据重构处理。数据重构处理同样采用自然间断点10 级分类法，形成了瞬变电磁法的3 个评价指标影像分级图（图3（d）、图3（e）和图3（f））。

3.3 地震映像法数据特征和数据重构

基于最佳偏移距的地震映像法主要用于快速确定采空区异常存在与否和空间位置，其埋深和规模可通过与高密度电阻率法和瞬变电磁法综合评价获得。采空区在地震映像法中主要表现为同相轴的错断和增多（图4（a）），个别显示为同相轴的错断和减少。为此，本研究将空气充填型采空区在地震映像法中的数据异常特征分为3 个类别进行位置区间赋值重构，无采空区异常位置区段赋值为1，采空区异常显示一般区段赋值为3，采空区异常显示强区段赋值为5，得到地震映像法的采空区异常评价指标D1W（图4（b））。

4 基于模糊数学模型的露天铁矿采空区探测地球物理异常评价

依据高密度电阻率法、瞬变电磁法和地震映像法的探测成果，提取了7个采空区地球物理异常评价指标，分别为G1R、G2T、G3H、S1R、S2T、S3H和D1W，且当它们每个指标在一个指定区段内都达到权重概率预定值以上，方可综合判定为采空区异常。同时，需在保证无遗漏采空区的条件下缩小探测数据异常圈定范围，以降低钻探验证成本。因此，有必要对7 个不确定性评价指标合理提取关键参数并进行赋权，确保评价模型的可靠性和拓展性，从而可适用于复杂多变的矿山地球物理探测评价，达到采空区地球物理异常评价的定量化和智能化。本研究采用已验证区域的空气充填型采空区的震—电—磁探测数据来开展研究（图5），同时，为了确保一致性和系统性，文中用到的震—电—磁地球物理探测成果数据为同一组数据。

4.1 评价指标模糊分类

7 个采空区异常评价指标可以分为4 类：第1 类为探测数据值指标，为G1R和S1R；第2类为探测数据提取指标，为G2T和S2T；第3 类为有效探测区间指标，为G3H和S3H；第4类为采空区探测异常位置指标，为D1W。对于上述每类评价指标的模糊分类，可使已进行数据重构的评价指标归一化和概率化，即将采空区异常评价指标分类数据转换到0～1 的概率区间内。该区间内，“0”分配给确定不属于指定采空区异常集合的分类值，“1”分配给确定属于指定采空区异常集合的分类值，0～1 区间内的概率值分配给可能为采空区异常的分类值。

4.1.1 第一类采空区异常评价指标模糊分类

G1R代表高密度电阻率法的直接电阻率测量值，对于空气充填型采空区，在无遗漏采空区异常标准指导下，探测的电阻率值若大于背景值或探测影像的电阻率平均值，即可认定为采空区异常，可以采用模糊较大值计算模型进行分析：

式中，x为输入计算值；f1为定义过渡区形状和特征的展开参数，本研究取3；f2为定义的起算中点，用于确定交叉点（起算中点分配概率值为0.5，大于起算中点的值为采空区异常的概率较高，小于起算中点的值为采空区异常的概率较低）。

该模型以评价指标G1R中探测影像电阻率平均值所代表的分类值（图2（e））为起算中点f2，其起算概率为0.5，将栅格数据G1R（图2（e））图像值作为输入值进行模糊计算，得到采空区异常评价指标G1R指示采空区异常存在的概率分布图（图6（a）），其概率值越大，属于采空区的可能性便越大。同理，可以采用模糊较大值计算模型对瞬变电磁法采空区异常评价指标S1R（图3（e））进行采空区异常概率分布计算，其起算中点f2为探测影像电阻率均值所代表的分类值，可以得到采空区异常评价指标S1R指示采空区异常存在的概率分布图（图6（b））。

4.1.2 第2类采空区异常评价指标模糊分类

第2 类采空区异常评价指标是高密度电阻率法的G2T和瞬变电磁法的S2T。G2T为电阻率的空间变化梯度值，正梯度值为高阻值异常圈闭区显示，代表空气充填型采空区异常。该评价指标可以采用模糊高斯模型进行计算：

式中，x为输入计算值；f1为定义过渡区形状和特征的展开参数，本研究取3；f2为定义的起算中点，用于确定交叉点（起算中点分配概率值为1，大于起算中点的值或小于起算中点的值对应为采空区异常的概率均较低）。

模糊高斯模型以指标G2T平均值所代表的分类值为起算中点f2，起算概率为1，小于或大于起算中点的概率均较低，该模型可以有效消除地表接触不良引起的正梯度值超大“采空区异常”，也可以规避将正梯度值过低的采空区异常边界过渡区确定为采空区异常。本研究将栅格数据G2T（图2（d））图像值作为输入计算值进行模糊计算，得到了采空区异常评价指标G2T指示采空区异常存在的概率分布图（图6（c））。

对于瞬变电磁法的采空区异常评价指标S2T，由于存在关断时间导致的浅部探测盲区，可有效消除因地表接触不良造成的“高阻采空区异常圈闭”，其正梯度值（图3（d））变化越大，显示为采空区异常的概率便越大，也可采用模糊较大值计算模型进行计算，其起算中点f2为正梯度值平均值所代表的分类值，计算得到的采空区异常评价指标S2T指示采空区异常存在的概率分布如图6（d）所示。

4.1.3 第3类采空区异常评价指标模糊分类

第3 类采空区异常评价指标为有效探测区间指标，为G3H和S3H（图2（f）和图3（f））。该类评价指标可以有效筛除浅部接触不良导致的异常数据、可信度低的深边部探测数据和关断时间导致的探测盲区数据，以降低采空区异常误判概率。该类指标可采用模糊高斯模型（式（2））进行计算，并依据前期探测和验证成果，将有效探测数据深度区间的中心深度值所代表的分类值作为起算中点f2，赋权概率为1，而浅部或深部异常数据为采空区的概率均较低，得到G3H和S3H的采空区异常概率分布图（图6（e）和6（f））。

4.1.4 第4类采空区异常评价指标模糊分类

地震映像法的采空区异常位置指标D1W，可依据其异常分类特征采用模糊较大值模型（式（1））计算，f1取值为2，其起算中点f2为3 级异常分类中的1 级和2 级采空区异常的中间值2（图4（b）），起算概率为0.5，图像计算得到的采空区异常评价指标D1W的采空区异常概率分布图如图6（g）所示。

4.2 评价模型构建

本研究采用数据重构方法和模糊数学计算模型提取了7 个采空区异常评价指标G1R、G2T、G3H、S1R、S2T、S3H和D1W，综合评价采空区异常存在与否。鉴于上述每个采空区异常评价指标均来自不确定性判别数据集或依据专家经验总结得出，且评价指标模糊分类中的概率赋值也选用了最低可能性，以免遗漏采空区异常，具体表现在有效探测数据区间内，空区存在的区域必有采空区异常，而有采空区异常的区域未必全是空区。因此，考虑到指标体系的不确定性，本研究空区异常评价采用模糊综合评价模型（图6（h））［20］，即模糊叠加模型，具体可为模糊And 函数评价模型：

式中，arg1,arg2,arg3,…,argn 为不同的评价指标取值，本研究n ∈[ ]1,7，Fuzzy And 叠加类型将返回单元位置所属集合的最小值。

模糊And叠加模型会创建一个输出栅格，为像元位置所属的每个集合的每个像元值给定了最小模糊值。本研究中，采空区异常评价模型是7个指标条件中的每个条件都已相对于其属于合适集合的隶属度被模糊化，模糊And模型会识别多个条件内像元属于一个合适集合的最小可能性。对于采空区地球物理异常评价，主要依据无遗漏采空区异常标准和模糊分类的概率赋值原则，只需要选择符合所有指标中概率值大于或等于0.5 的具有采空区异常的空间位置集合，叠加生成采空区异常指标综合评价概率分布图，图中概率值大于0.5 的区块圈闭即为采空区异常，需要进一步进行钻探验证。

4.3 评价结果

基于上述综合评价结果，采用ArcGIS 软件中的空间分析模块进行弓长岭铁矿区多采空区异常评价指标图像矩阵运算，计算模型为Fuzzy And，结果见图7。图7中测线起点以高密度电阻率法G8测线起点为0 点，与其平行的瞬变电磁法和地震映像法测线则以相对于电法起点叠加距离作为起点计算空间坐标。在图7 中，根据0.5 及其以上概率值作为异常圈定标准，圈定出了一个采空区异常区（即空气充填型采空区）。在对异常区进行验证时，可将采空区异常区上部边界与验证得出的采空区顶板埋深进行对比分析。图7 中采空区异常区圈定的顶板埋深为35.7 m，经过钻探验证埋深为33.8 m，两者相差1.9 m，相比于传统经验法的圈定误差（±5.0 m 以内），精度提高了3.0 m以上。

通过对空气充填型采空区的地球物理探测数据进行指标化提取，并采用本研究提出的综合评价模型，在ArcGIS 平台上定量圈定了采空区异常边界，通过钻探成果相比，可以在无遗漏采空区的前提下，大大提高采空区圈定效率和准确度。在上述评价方法指导下，对弓长岭铁矿区水充填型采空区进行了综合评价，圈定的采空区异常顶板埋深18.3 m，钻探验证结果为20.1 m，两者相差在2 m 以内。可见，本研究选择的采空区异常圈定指标具有合理性，即对于不确定性的地球物理异常可以采用模糊数学方法进行定量研究，实现地球物理探测数据评价的指标化、定量化和智能化。此外，本研究评价方法可以通过ArcGIS 平台对采集到的大量探测数据进行流程式的智能化快速综合评价，并可依据评价结果进行快速钻探验证，确保矿山生产安全。

5 结 论

（1）划分了多种采空区类型，依据采空区充填介质和物性差异的不同，将采空区分为空气充填型采空区和水充填型采空区两大类，并分别建立了相应的地球物理识别特征和震—电—磁探测方法组合。

（2）根据高密度电阻率法、瞬变电磁法和地震映像法的采空区异常数据识别特征，构建并提取了7个采空区异常评价指标G1R、G2T、G3H、S1R、S2T、S3H和D1W。

（3）通过对7个采空区异常指标进行模糊分类和模糊数学模型综合评价能够定量圈定出采空区异常边界，可使探测精度保持在±2.0 m以内，为采空区异常快速定量化和智能化评价提供了技术手段。