曹亚琴1,王永志1,2,*,卢鹏羽1
吉林大学地球探测科学与技术学院
吉林大学综合信息矿产预测研究所
关键词 机器学习;多重分形;能谱密度面积(S-A);线性回归拟合;阈值
为了弥补统计方法识别异常背景值受原始数据等因素制约而难以刻画地球系统非线性特征的不足,科学家将分形理论引入地球科学领域,用其科学描述成矿作用过程形成的自相似或广义自相似特征的多重分形场。通过元素浓度-面积(C-A,Concentration-Area)和能谱密度面积(S-A,Spatial-Area)等分形模型的广泛应用,分形/多重分形方法已成为矿产资源预测领域的背景异常场分离的重要方法,在各种地学数据异常背景分离等过程中成效显著。直接影响成矿背景异常场分离效果的异常下限常用GeoDAS软件计算并从双对数图中获得。
GeoDAS是20世纪基于MapObjects、采用VisualBasic开发的32位专业应用系统,后续虽有完善但尚未升级至64位,难以直接利用其他在线的计算资源(如机器学习库),亦无法满足地学大数据的计算需求。随着人工智能技术在地学领域的落地应用,支持机器学习库、Python语言、原生64位的ArcGISPro软件平台越来越受地学计算重视,但其缺少用于多重分形计算工具。本文探索采用机器学习方法按需自动计算分形/多重分形的一个或多个阈值,直接在ArcGISPro环境下对地学大数据进行背景场和异常场的快速分解,可拓宽非线性地学大数据成矿预测计算的适用范围。
自然界的地球化学元素(如Au、Cu等)局部高度富集可形成矿体,这种成矿过程不断叠加即可形成多重分形场,元素浓度-面积(C-A)分形模型擅于处理此种非线性场。元素含量(浓度)大于某一含量阈值的区域面积表示为:
(1)式进行变形并取对数得到:
式(3)中元素浓度大于某一含量时的面积值与浓度本身在双对数坐标中明显呈线性关系,斜率即为分形维数。在多成矿元素地球化学异常综合中,最小二乘拟合会出现多个阈值,斜率即为分形维数。以不同阈值为分界点可确定异常下限,从而有效圈定对应地球化学异常的范围。
能谱密度-面积(S-A)是C-A模型在频率域的扩展,能在傅立叶能谱空间中度量地球物理和地球化学异常所对应的各向异性的广义自相似性,根据识别不同的广义自相似性构造不同的滤波器,进而利用反傅里叶变换对异常和背景进行分解。此法圈定的物化探异常具有异常强度、大小范围和不同背景等多个显著特征,其在频率域中表现出与背景场的不同自相似性。能谱S(ωx,ωy)可用傅里叶变换结果的实虚部绝对值的平方计算。对于一维时间序列,能谱与频率的关系为:
对于二维情况,引入了同样的函数来描述平均能谱与频率之间的关系。多重分形理论证明了具有自相似性地质体分形场的指数关系为:
常用线性回归、多项式回归、神经网络等方法计算双对数图的阈值。文章采用机器学习的线性回归方法,通过分段拟合自动计算阈值(一个、两个、三个或多个),而线性回归拟合算法的本质是最小二乘法。设原始数据集为{(xi,yi)},则其线性拟合直线的斜率和截距为:
线性回归拟合从一组无序的数据中找到一条或几条直线的趋势,从而消除离散数据所带来的局部扰动。线性回归模型可拟合一条或多条直线,拟合的多条直线交点即为阈值。图1为基于机器学习的背景异常分解技术中关键参数计算的流程,数据拟合按需采用线性回归、多项式回归、神经网络等方法,可采用单一方法亦可使用多种混合方法,本文以线性回归为例。
图1 基于机器学习的背景异常分解关键参数计算框架
对于地球物理数据或地球化学数据等进行反距离权重(IDW,Inverse Distance Weight)插值和多重分形处理,对其进行C-A分析或S-A分析得到双对数图数据,利用机器学习方法(如线性回归、神经网络等)分段拟合双对数图数据,自动计算最小拟合差之和,得到优化后的最佳阈值。具体计算步骤如下:
(1)记分形方法C-A、S-A异常分解得到的双对数图数据总数为n,每个数据点的坐标为xi,yi(0≤i≤n-1);
(2)根据设置阈值数量不同(如1、2、3),采用机器学习(如Scikit-learn)的拟合算法(如线性回归)进行分段拟合,分别计算阈值数对应的数据区域;
(3)利用各拟合直线lj的斜率a和截距b,代入xi,求出新值yi;
(4)计算每种情况各拟合直线l1,l2,…lj(2≤j≤4)的最小均方差之和Min(Error(l1)+Error(l2)+…Error(lj))=Min(∑jt=2∑(̂yti-yti)2)),从而确定各拟合直线l1,l2,…lj的最佳函数关系;
(5)自动计算最佳函数关系的各个拟合直线l1,l2,…lj的交点的横坐标即为阈值。
基于Scikit-learn机器学习框架,采用Python语言,在64位的ArcGISPro软件中实现本文提出的自动求取阈值算法,并将其封装成一个带有图形界面的工具嵌入到ArcGISPro中。
实验数据采用加拿大新斯科舍省南部Nova地区的地质数据,主要包括湖底沉积物的地球化学样品数据(CSV,共671个)、地质图数据(SHP)等(如图2)。湖底沉积物的地球化学数据主要提取经度、纬度、砷(As)元素浓度三列值用于对比实验。
图2 加拿大新斯科舍省南部 Nova 地区 As 数据
(a)岩性图及As测点位置;(b)As数据IDW插值结果。
3.2 实验效果分析
为了检验阈值(异常下限)求取方法的正确性,以同一实验数据作为输入,分别采用GeoDAS软件、机器学习方法计算阈值,对计算结果进行对比。表1列出了两种方法在阈值数为1、2和3时三种不同能谱值计算结果。GeoDAS(LogS)和线性回归拟合(LogS)均为计算能谱阈值的对数,GeoDAS(S)和线性回归拟合(S)为求取出的能谱阈值,差为两种方法求取阈值之差,平均差为每种情况差的平均值。
当阈值数为1时,双对数图数据被分成两段(图3)。利用GeoDAS软件估计的能谱阈值为8.28,利用机器学习估计的能谱阈值为8.94。以GeoDAS软件阈值估计结果为基准,阈值相差7.97%。图3a和图3b的拟合效果显示,前段GeoDAS软件处理吻合效果较好,后段机器学习方法处理吻合效果较好。但从图3a中可明显看出,GeoDAS两段求解的拟合直线并未求交点,可能是产生二者两段不同的主要原因。利用机器学习中线性回归算法拟合求阈值可达到或优于GeoDAS处理效果。
当阈值数为2时,双对数图数据分成3段(图4)。利用GeoDAS软件、机器学习方法估计的阈值分别为7.49和8.94、7.88和9.15。同样以GeoDAS软件阈值估计结果为基准进行对比,第一和第二个阈值分别相差5.21%和2.35%,平均差为3.78%。从图4a和图4b中可以明显看出,两种方法拟合双对数图均取得较好的效果且估计的两个阈值结果也相差较小,验证了机器学习方法计算的正确性。但从图4a同样可清晰地看出,GeoDAS软件的三条拟合直接并未在交点处相交,而机器方法得到的三条拟合直线在两个交点处吻合,可能是造成误差的主要原因。
(b)利用机器学习估计的阈值。
当阈值数为3时,双对数图数据分成三段(图5)。利用GeoDAS软件估计的阈值为6.18、7.60和8.94,而机器学习方法估计的3个阈值是7.03、8.17和9.15。对两种处理方法对比可知,第一、二和三个阈值分别相差13.75%、7.5%和2.35%,平均差为7.86%。从双对数图拟合效果及导出数据查看(如图5a、b),GeoDAS软件在拟合过程中忽略了第一个点,只对后面数据进行较好地拟合。而机器学习方法对第一段所有数据进行完整拟合(含第一个点)。通过对阈值结果分析,两种方法虽有差距但相对较小,特别是后两个阈值估计误差较小。造成这种误差的原因(第一个阈值相差较大)可能是对第一个数据点的处理方式不同。
(a)利用GeoDAS软件估计的阈值;
(b)利用机器学习估计的阈值。
S-A绘制的是一个二维散点图,横轴代表能谱密度取自然对数值,纵轴代表能谱密度大于阈值S是面积取自然对数值。随着能谱值增加,所圈闭的面积越来越小,故双对数图总呈下降趋势。不同的直线遵循不同地幂律关系,从而将异常场从背景中分离。从两种方法的双对数图拟合及阈值估算结果对比可知,线性回归拟合的S-A散点图计算的阈值和GeoDAS求得的相差较小。在估计三个阈值情况时,仅左起第一个点(阈值1)相差较大,说明在斜率变化不大时,线性回归拟合方法仍能很好地确定阈值的位置。以阈值数为2时的情况为例,把能谱划分成3段圈定异常背景场。
图6列出了两种方法相同分段得到不同阈值求得的背景场,(a)和(b)总体趋势差异较小,(c)和(d)虽在右上聚集区域数量不同但总体空间分布趋势一致,(e)和(f)几乎完全相同。由于在S-A分形模型中,阈值对结果影响相对较小,故虽然机器学习方法分离出的背景场区域稍大,但和GeoDAS处理结果基本一致。图7对两种方法按不同分段求得的异常场进行对比可见,虽局部稍有差异,但机器学习方法分离得到的异常场和GeoDAS软件处理得到的异常场总体趋势基本相同,空间分布规律一致。说明了采用机器学习(如线性回归)拟合求阈值方法可行。
(f)(GeoDAS)全域能谱情况下分离出的背景场。
相比机器学习,人工的经验还是更胜一筹
4 结论与探讨
4.1 结论
(1)为了弥补分形软件利用新型在线计算资源的不足,本文以成矿背景异常分解技术为切入点,提出一种基于机器学习(如线性回归)的自动化计算能谱阈值方法。
(2)在64位软件ArcGISPro中,采用Python语言编程实现了该自动化算法,融入工具箱的算法可同时调用在线、本地的计算资源和数据资源。
(3)以GeoDAS自带示例数据作为输入,对GeoDAS和机器学习方法分别进行能谱阈值计算实验,经各分段阈值下限、分离异常背景场的结果对比,二者总体趋势和空间分布基本保持一致,说明了利用机器学习自动求取阈值方法的有效性。
4.2 探讨
(1)GeoDAS分段拟合过程中,发现部分拟合分段首点未参与计算、两条拟合直线未相交等现象,可能是产生其与机器学习方法阈值误差的原因,有必要进一步验证。