CN114139819A - 基于地统计加权随机森林的地球化学变量空间预测方法 - Google Patents

基于地统计加权随机森林的地球化学变量空间预测方法 Download PDF

Info

Publication number
CN114139819A
CN114139819A CN202111483212.6A CN202111483212A CN114139819A CN 114139819 A CN114139819 A CN 114139819A CN 202111483212 A CN202111483212 A CN 202111483212A CN 114139819 A CN114139819 A CN 114139819A
Authority
CN
China
Prior art keywords
geochemical
data
random forest
training
variable
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202111483212.6A
Other languages
English (en)
Other versions
CN114139819B (zh
Inventor
王健
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Chengdu Univeristy of Technology
Original Assignee
Chengdu Univeristy of Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Chengdu Univeristy of Technology filed Critical Chengdu Univeristy of Technology
Priority to CN202111483212.6A priority Critical patent/CN114139819B/zh
Publication of CN114139819A publication Critical patent/CN114139819A/zh
Application granted granted Critical
Publication of CN114139819B publication Critical patent/CN114139819B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q10/00—Administration; Management
    • G06Q10/04—Forecasting or optimisation specially adapted for administrative or management purposes, e.g. linear programming or "cutting stock problem"
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G06F17/12—Simultaneous equations, e.g. systems of linear equations
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/21—Design or setup of recognition systems or techniques; Extraction of features in feature space; Blind source separation
    • G06F18/214—Generating training patterns; Bootstrap methods, e.g. bagging or boosting
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/24—Classification techniques
    • G06F18/243—Classification techniques relating to the number of classes
    • G06F18/24323—Tree-organised classifiers
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N5/00—Computing arrangements using knowledge-based models
    • G06N5/01—Dynamic search techniques; Heuristics; Dynamic trees; Branch-and-bound
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Mathematical Physics (AREA)
  • Mathematical Analysis (AREA)
  • Computational Mathematics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • General Engineering & Computer Science (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Business, Economics & Management (AREA)
  • Evolutionary Computation (AREA)
  • Artificial Intelligence (AREA)
  • Software Systems (AREA)
  • Evolutionary Biology (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Operations Research (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Economics (AREA)
  • Human Resources & Organizations (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Strategic Management (AREA)
  • Probability & Statistics with Applications (AREA)
  • Development Economics (AREA)
  • Game Theory and Decision Science (AREA)
  • Computing Systems (AREA)
  • Computational Linguistics (AREA)
  • Entrepreneurship & Innovation (AREA)
  • Marketing (AREA)
  • Quality & Reliability (AREA)
  • Tourism & Hospitality (AREA)
  • General Business, Economics & Management (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了一种基于地统计加权随机森林的地球化学变量空间预测方法，包括以下步骤：S1、获取训练随机森林模型所需的训练数据，训练数据包括输入变量和模型输出变量，输入变量包括地质要素和遥感要素，模型输出变量为已知地球化学元素含量观测数据；S2、根据输入变量与已知地球化学元素含量观测数据之间的相关系数及输入变量之间的相关性筛选输入变量作为预测因子，将已知地球化学元素含量观测数据与预测因子作为训练样本；S3、基于训练样本训练随机森林模型，并采用训练好的随机森林模型对待预测位置的地球化学元素含量进行预测。本发明应用时同时考虑地学现象的空间相关性和异质性，能提升预测精度。

Description

基于地统计加权随机森林的地球化学变量空间预测方法
技术领域
本发明涉及地质勘查领域，具体是基于地统计加权随机森林的地球化学变量空间预测方法。
背景技术
勘查地球化学手段在过去几十年的找矿工作中发挥了举足轻重的作用。随着露头矿和浅表矿逐渐被发现殆尽，找矿难度日益增大，如何在复杂地质条件下挖掘有用的地球化学信息，已经成为缓解目前矿产资源紧缺形势的重要途径。地球化学分布模式是原生和次生地质过程、人类活动等长期共同作用的结果。受区域地质背景、矿体埋深、矿化类型及规模、风化强度以及地球化学景观等的影响，地球化学元素在表生介质中呈现十分复杂的变异特征，导致地球化学异常识别结果存在很大的不确定性，从而增加了勘查风险，限制了勘查地球化学手段的找矿成效。因此，准确模拟地球化学元素含量的空间分布，并刻画其不确定性显得尤为重要。
随机森林由于具有模型简单、准确度高、可解释性强等优点，在地球化学元素含量空间分布模拟中具有较大应用潜力。然而，经典随机森林算法潜在假设样本是独立同分布的，并未考虑元素含量分布的空间异质性和空间相关性特征，从而增加了模拟结果的不确定性。针对上述问题，不少国内外研究尝试改进随机森林算法。如：国外一些学者提出采用考虑空间依赖性的方法，其将观测点的欧氏距离图层作为预测因子之一，该方法操作简单、间接考虑了观测值之间的空间关联，但没有直接考虑空间异质性和观测值分布的空间结构特征。中国专利公布号CN112163731A于2021年01月01日公开了发明创造名称为“一种基于加权随机森林的专变用户电费回收风险识别方法”的发明专利申请，其针对非空间科学问题，“加权”是对随机森林中决策树的贡献进行加权，并不涉及地理位置的概念。中国专利公布号CN104391970A于2015年03月04日公开了发明创造名称为“一种属性子空间加权的随机森林数据处理方法”的发明专利申请，其通过加权属性子集构建决策树，再将其整合形成随机森林模型，该方法针对超高维数据有效处理问题，采用信息增益法为属性加权，并不涉及观测变量在地理空间上的相关性。中国专利公布号CN113033599A于2021年06月25日公开了发明创造名称为“面向露头地质体岩层分层的空间随机森林算法”的发明专利申请，其同时考虑岩层体元的空间特征和属性特征构建空间决策树模型，该处的“空间特征”主要包括地质体体元中心坐标及其产状信息，并未涉及地学变量的空间相关性特征。
综上所述，现有研究没有同时考虑地学现象的空间相关性和异质性提出行之有效的改进方案，采用现有方案模拟地球化学元素含量空间分布的精度较低。
发明内容
本发明的目的在于解决现有模拟地球化学元素含量空间分布精度较低的问题，提供了一种基于地统计加权随机森林的地球化学变量空间预测方法，其应用时同时考虑地学现象的空间相关性和异质性，能提升预测精度。
本发明的目的主要通过以下技术方案实现：基于地统计加权随机森林的地球化学变量空间预测方法，包括以下步骤：
S1、获取训练随机森林模型所需的训练数据，所述训练数据包括输入变量和模型输出变量，输入变量包括地质要素和遥感要素，地质要素包括从区域地质数据提取的岩性类别数据和从区域地质数据提取的线环构造矢量数据转换的断裂密度，遥感要素包括地形坡度和植被覆盖度，模型输出变量为已知地球化学元素含量观测数据；其中，输入变量和输出变量采用栅格格式存储；
S2、根据输入变量与已知地球化学元素含量观测数据之间的相关系数及输入变量之间的相关性筛选输入变量作为预测因子，将已知地球化学元素含量观测数据与预测因子作为训练样本；
S3、基于训练样本训练随机森林模型，并采用训练好的随机森林模型对待预测位置的地球化学元素含量进行预测。
进一步的，所述步骤S1中将线环构造矢量数据转换为断裂密度的转换公式为：
Figure BDA0003396269870000021
其中，i表示栅格的行，j表示栅格的列，ρij为第i行第j列栅格的断裂密度，Lij为对应栅格内线性断裂和环形断裂的总长度，S为栅格的面积；
所述步骤S1中岩性类别数据为将区域地质数据中的岩性分布数据转换的类别变量，即
cij∈{Ck,k＝1,2,…n}
其中，n为岩性类别总数，Ck表示n种岩性类别中的第k种岩性，cij代表栅格位置(i,j)处的岩性类别。
进一步的，所述步骤S1还包括对遥感数据进行辐射校正处理，并利用以下公式计算归一化植被指数：
Figure BDA0003396269870000022
其中，i表示栅格的行，j表示栅格的列，NDVIij为栅格位置(i,j)处的植被指数，
Figure BDA0003396269870000023
为近红外波段处的反射率，
Figure BDA0003396269870000031
为红光波段处的反射率。
进一步的，所述步骤S1还包括以下已知地球化学元素含量观测数据处理步骤：
迭代剔除位于区间[μm-3σm,μm+3σm]之外的已知观测值，其中，μm为第m个元素的均值，σm为迭代过程中第m个元素的标准差；
对低于检测限的观测值用相应元素检测限的一半代替；
对缺失值用对应元素的均值补全；
采用中心对数比变换法对已知多元地球化学元素含量观测数据进行处理，处理公式如下：
Figure BDA0003396269870000032
其中，xm和zm分别代表第m个元素变换前后的值，D代表元素个数。本发明采用对数比变换法，能消除闭合效应。
进一步的，所述步骤S1在获取训练随机森林模型所需的训练数据后还包括以下步骤：
采用待预测区域掩膜文件对栅格数据进行裁剪，并统一数据的坐标系。
进一步的，所述步骤S2中根据输入变量与已知地球化学元素含量观测数据之间的相关系数及输入变量之间的相关性筛选输入变量包括以下步骤：
计算各个输入变量与已知地球化学元素含量观测数据之间的相关系数并剔除弱相关性输入变量，再计算各个输入变量之间的相关性，任意两个输入变量之间相关性大于设定阈值时保留其一。
进一步的，所述计算各个输入变量与已知地球化学元素含量观测数据之间的相关系数、计算各个输入变量之间的相关性均包括以下步骤：对于连续型数据，计算Pearson相关系数，对于离散型数据，计算Spearman秩相关系数；
计算各个输入变量与已知地球化学元素含量观测数据之间的相关系数时若相关系数绝对值|r|<0.3，则剔除该属性特征。
进一步的，所述步骤S3包括以下步骤：
S31、定义模拟栅格图层和滑动窗口半径序列，以及确定随机森林模型中决策树的数目和决策树节点划分时属性的数目；
S32、按照从左到右，从上到下的顺序遍历栅格单元，对于当前需预测地球化学元素含量值的栅格位置，利用紧凑度指数确定最佳滑动窗口半径；
S33、计算空间邻域观测点的空间权重：提取最佳窗口内的训练样本，通过全局变差函数模型和克立格方程组求解训练样本的权重；
S34、训练地统计加权随机森林模型：采用地统计加权的均方根误差函数作为改进的目标函数，训练随机森林回归模型；
S35、预测地球化学元素含量未知点处的取值：将待预测位置s处的预测因子值输入步骤S34训练好的随机森林模型，预测地球化学元素含量值；
S36、重复步骤S32～步骤S35，遍历结束即得完整的地球化学元素含量空间分布图。本发明采用滑动窗口技术表征地球化学元素含量分布的空间异质性，定义滑动窗口半径序列，用于之后借助紧凑度指数动态确定最佳滑动窗口尺寸。如此，本发明将当前栅格位置处的预测因子作为随机森林模型的输入，即可预测地球化学元素含量值。重复步骤S32～步骤S35，遍历结束，即得完整的地球化学元素含量空间分布图。
进一步的，所述步骤S31中定义滑动窗口半径序列为ε1≤ε2≤…εT，滑动窗口采用方形窗口，当窗口半径为εt(1≤t≤T)时，滑动窗口的大小为(2εt+1)×(2εt+1)；
所述步骤S32中利用紧凑度指数确定最佳滑动窗口半径包括以下步骤：
对于当前需预测地球化学元素含量值的栅格位置s，计算滑动窗口半径为εt时的紧凑度指数，紧凑度指数计算公式如下：
Figure BDA0003396269870000041
其中，Ct代表落入半径为εt窗口内的地球化学观测样本集合，Nt是样本数目,Yj′和Yk′代表地球化学样本；若相邻两个窗口内样本几乎来自同一个统计分布总体，则Dt和Dt+1差别不大；若样本明显来自两个不同统计总体，则Di和Di+1在整个变化趋势中会出现突变，计算相邻紧凑度指数Dt(t＝1,2,…T)之差Dt+1-Dt，设差值最大时滑动窗口对应的半径为εt，则优选出的最佳半径为ε*(s)＝εt；
所述步骤S33计算空间邻域观测点的空间权重包括以下步骤：
获取位置s处最佳窗口内的训练样本
{(X1(sl′),X2(sl′),…,XK(sl′)；Y(sl′)),l′＝1,2,…L′{(L′<L)，计算各个样本所处位置sl′(l′＝1,2,…l′)与位置s的欧氏距离，计算公式为：
Figure BDA0003396269870000042
将计算得到的欧氏距离代入利用全局样本拟合的变差函数模型γ(h)；
求解克立格方程组，公式为：
Figure BDA0003396269870000051
上式共包含L′+1个方程，式中未知数为λl′(l′＝1,2,…L′)和μ共L′+1个，其中μ为拉格朗日算子，求解该方程组得到的λl′即为邻域观测点的权重；
所述步骤S34训练地统计加权随机森林模型包括以下步骤：
对于每个回归决策树，采用改进的目标损失函数进行训练，公式为：
Figure BDA0003396269870000052
其中，
Figure BDA0003396269870000053
Figure BDA0003396269870000054
Figure BDA0003396269870000055
X代表构建决策树的切分属性或因子，v代表构建决策树相应属性的切分点，Yleft和Yright分别代表左子节点和右子节点所包含的地球化学观测样本点取值集合，Ya和Yb分别代表左右子节点集合中的样本值，Nleft和Nright分别代表左子节点和右子节点中观测样本的数目。
进一步的，所述步骤S2与步骤S3之间还包括以下步骤：将训练样本分为训练数据集和验证数据集，所述步骤S3基于训练数据集训练随机森林模型；所述步骤S3之后还包括以下步骤：提取验证数据集中各验证样本所处位置处的地球化学变量预测值，建立验证样本观测值与预测值数据对(Yg,eval,Yg,obs)(g＝1,2,…Neval)的线性回归模型，得到决定系数R2，据此定量评价随机森林模型的预测效果；其中，Yg,eval代表第g个点处地球化学变量的模型预测值，Yg,obs代表第g个点处地球化学变量的参考值，Neval代表验证点的个数。如此，本发明基于构建好的训练数据集训练随机森林模型，并借助训练好的模型对未采样栅格的地球化学元素含量进行预测。借助测试数据集对得到的元素含量预测结果进行验证分析，并将模型输出的特征重要性进行空间制图，可作为后续分析的重要参考。
综上所述，本发明与现有技术相比具有以下有益效果：(1)现有经典随机森林回归算法在构建决策树时，预测值的计算方式是求取叶子节点内所有样本观测值的平均值。由于地学变量观测样本存在一定程度的空间相关性，本发明通过地统计方法为叶子节点内样本进行加权，预测得到的目标位置变量取值更加准确，元素含量分布结果更能反映地质地球化学过程的空间变异特征。
(2)现有经典随机森林算法只在研究区构建一个全局随机森林模型，难以顾及地球化学元素含量空间分布的变异特征。本发明将滑动窗口技术与随机森林模型相结合，同时借助紧凑度指数确定任意局部位置的最佳滑动窗口尺寸，能够更好地捕捉局部地球化学分布模式的尺度。
(3)本发明同时考虑空间异质性和空间相关性特征，对地学变量的空间分布进行预测，并提供度量其不确定性的指标，有助于提高元素含量的预测精度。
附图说明
此处所说明的附图用来提供对本发明实施例的进一步理解，构成本申请的一部分，并不构成对本发明实施例的限定。在附图中：
图1为本发明一个具体实施例的流程图；
图2为本发明一个具体实施例的预测方法示意图。
具体实施方式
为使本发明的目的、技术方案和优点更加清楚明白，下面结合实施例和附图，对本发明作进一步的详细说明，本发明的示意性实施方式及其说明仅用于解释本发明，并不作为对本发明的限定。
实施例：
如图1所示，基于地统计加权随机森林的地球化学变量空间预测方法，包括以下步骤：S1、获取训练随机森林模型所需的训练数据，其中，训练数据包括输入变量和模型输出变量，输入变量和输出变量采用栅格格式存储，输入变量包括地质要素和遥感要素，地质要素包括从区域地质数据提取的岩性类别数据和从区域地质数据提取的线环构造矢量数据转换的断裂密度，遥感要素包括地形坡度和植被覆盖度，模型输出变量为已知地球化学元素含量观测数据；S2、根据输入变量与已知地球化学元素含量观测数据之间的相关系数及输入变量之间的相关性筛选输入变量作为预测因子，将已知地球化学元素含量观测数据与预测因子作为训练样本；S3、基于训练样本训练随机森林模型，并采用训练好的随机森林模型对待预测位置的地球化学元素含量进行预测。
本实施例中输入变量即为模型训练所需的输入要素，模型输出变量即为模型的响应变量，本实施例针对地球化学元素含量预测问题，区域地质数据可从全国地质资料馆获取，收集1:5万或1:20万水系沉积物区域勘查地球化学数据，该数据多以原始数据表格(如Excel)或矢量格式存储(如MapGIS)，基本包含常量元素化合物及微量元素含量数据。本实施例统一以栅格数据格式(如TIFF)存储，从区域地质数据提取的线环构造矢量数据转换的断裂密度为连续型变量，从区域地质数据中提取岩性类别数据，用无序类别变量表示。遥感数据可从地理空间数据云网站获取，提取地形坡度、植被覆盖度，用连续型变量表示。本实施例的输出变量是具有完整空间分布的已知地球化学元素含量观测数据，模型训练时输入，其只有部分位置有值，地球化学元素含量未知点(即待预测位置)处的取值采用训练好的随机森林模型预测。
本实施例在具体实施时，为了进一步提升精度，输入变量不局限于地质要素和遥感要素，还可提取其它与地球化学元素含量有关的地质、地球物理等专题数据，并对这些数据进行预处理，如可对航磁数据进行航磁化极处理。本实施例的地质要素不局限于岩性类别和断裂密度，遥感要素不局限于地形坡度和植被覆盖度。
输入变量和输出变量属于类别变量或连续变量，因此，需要对收集的数据进行适当处理和变换。为便于存储与进一步建模分析，数据统一采用栅格格式存放。为此，需构建能够覆盖研究区的矩形栅格，根据所收集数据的精度和研究目的，设置合适的栅格单元大小。
本实施例的步骤S1在获取训练随机森林模型所需的训练数据后还包括以下步骤：离群值剔除、缺失值插补和成分数据变换等预处理，再采用待预测区域掩膜文件对栅格数据进行裁剪，并统一数据的坐标系。其中，掩膜文件为以矢量或栅格格式存储的研究区边界文件。本实施例应用时，若存在数据的坐标系不一致问题，采用现有的ArcGIS软件将其进行变换统一，通过统一数据的坐标系，便于不同空间数据的对比和运算。
本实施例步骤S1中将线环构造矢量数据转换为断裂密度的转换公式为：
Figure BDA0003396269870000071
其中，i表示栅格的行，j表示栅格的列，ρij为第i行第j列栅格的断裂密度，Lij为对应栅格内线性断裂和环形断裂的总长度，S为栅格的面积；线环构造矢量数据为区域地质数据中的断裂矢量数据，本实施例经转换后可得到地质岩性类别图。
本实施例步骤S1中岩性类别数据为将区域地质数据中的岩性分布数据转换的类别变量，即
cij∈{Ck,k＝1,2,…n}
其中，n为岩性类别总数，Ck表示n种岩性类别中的第k种岩性，cij代表栅格位置(i,j)处的岩性类别。如此，本实施例经转换后可得到断裂密度分布图。
所述步骤S1还包括对遥感数据进行辐射校正处理，并利用以下公式计算归一化植被指数：
Figure BDA0003396269870000081
其中，i表示栅格的行，j表示栅格的列，NDVIij为栅格位置(i,j)处的植被指数，
Figure BDA0003396269870000082
为近红外波段处的反射率，
Figure BDA0003396269870000083
为红光波段处的反射率。
本实施例的步骤S1还包括以下已知地球化学元素含量观测数据处理步骤：针对已知地球化学元素含量观测数据，首先采用阈值法对目标元素含量的离群值进行剔除，本实施例的剔除方式为迭代剔除位于区间[μm-3σm,μm+3σm]之外的已知观测值，其中，μm为第m个元素的均值，σm为迭代过程中第m个元素的标准差；
对低于检测限的观测值用相应元素检测限的一半代替；
对缺失值用对应元素的均值补全；
由于地球化学元素含量数据属于闭合数据，即所有元素或化合物的含量值的和为一常数，不能够直接用于欧氏距离等的计算，因此，需采用对数比变换方法消除闭合效应，因此，本实施例采用中心对数比变换法对已知多元地球化学元素含量观测数据进行处理，处理公式如下：
Figure BDA0003396269870000084
其中，xm和zm分别代表第m个元素变换前后的值，D代表元素个数。
本实施例步骤S2中根据输入变量与已知地球化学元素含量观测数据之间的相关系数及输入变量之间的相关性筛选输入变量包括以下步骤：计算各个输入变量与已知地球化学元素含量观测数据之间的相关系数并剔除弱相关性输入变量，在基本不影响模型精度的同时使模型尽可能简单；再计算各个输入变量之间的相关性，任意两个输入变量之间相关性大于设定阈值时保留其一，以避免由于存在强的多重共线性导致模型结果不可靠。保留下来的输入要素和已知地球化学元素含量观测数据作为模型的最终输入数据集，即作为训练样本。本实施例计算各个输入变量与已知地球化学元素含量观测数据之间的相关系数、计算各个输入变量之间的相关性均包括以下步骤：对于连续型数据，计算Pearson相关系数，对于离散型数据，计算Spearman秩相关系数。计算各个输入变量与已知地球化学元素含量观测数据之间的相关系数时若相关系数绝对值|r|<0.3，则剔除该属性特征。本实施例应用时，若存在某两个或多个输入变量之间的相关系数大于设定阈值(本实施例设为0.85)，则剔除部分输入变量，保证各个输入变量之间的相关性不高于设定阈值，以避免多重共线性问题导致模型结果不可靠。
本实施例的骤S2与步骤S3之间还包括以下步骤：将训练样本分为训练数据集和验证数据集，步骤S3基于训练数据集训练随机森林模型。本实施例的步骤S3之后还包括以下步骤：提取验证数据集中各验证样本所处位置处的地球化学变量预测值，建立验证样本观测值与预测值数据对(Yg,eval,Yg,obs)(g＝1,2,…Neval)的线性回归模型，得到决定系数R2，据此定量评价随机森林模型的预测效果；其中，Yg,eval代表第g个点处地球化学变量的模型预测值，Yg,obs代表第g个点处地球化学变量的参考值，Neval代表验证点的个数。
本实施例采用随机抽样方式划分为训练集和验证数据集，训练数据集用于模型训练，验证数据集用于模型性能验证，具体方式为：设经过上述处理得到的预测因子集合为
Figure BDA0003396269870000091
K为预测因子数目，s＝(i″,j″)代表空间位置，Xk(s)代表第k(k＝1,2,…K)个预测因子在位置s上的取值，G代表研究区域，
Figure BDA0003396269870000092
代表二维实数空间；地球化学元素含量值为{Y(s1),Y(s2),…Y(sL)},其中sl(l＝1,2,…L)代表地球化学样本的空间位置,提取位置sl处的预测因子，构成数据集D＝{(X1(sl),X2(sl),…,XK(sl)；Y(sl)),l＝1,2,…L。采用随机抽样方式按照7:3的比例将D划分为训练数据集和验证数据集D＝Dtrain∪Deval。
本实施例的步骤S3包括以下步骤：S31、定义模拟栅格图层和滑动窗口半径序列，以及确定随机森林模型中决策树的数目和决策树节点划分时属性的数目；S32、按照从左到右，从上到下的顺序遍历栅格单元，对于当前需预测地球化学元素含量值的栅格位置，利用紧凑度指数确定最佳滑动窗口半径；S33、计算空间邻域观测点的空间权重：提取最佳窗口内的训练样本，通过全局变差函数模型和克立格方程组求解训练样本的权重；S34、训练地统计加权随机森林模型：采用地统计加权的均方根误差函数作为改进的目标函数，训练随机森林回归模型；S35、预测地球化学元素含量未知点处的取值：将待预测位置s处的预测因子值输入步骤S34训练好的随机森林模型，预测地球化学元素含量值；S36、重复步骤S32～步骤S35，遍历结束即得完整的地球化学元素含量空间分布图。本实施例定义模拟栅格图层时，栅格单元大小与前述构建能够覆盖研究区的矩形栅格保持一致，并将已知地球化学观测点数据赋值给最邻近栅格单元，作为后续模拟的条件数据。这里采用滑动窗口技术表征地球化学元素含量分布的空间异质性。因此，需事先定义滑动窗口半径序列，合适的窗口尺寸应该能够反映地球化学元素含量变异的空间范围，同时窗口内数据量应足够用于统计推断。由于随机森林模型的性能依赖于两种随机机制：一是随机森林中每棵决策树所使用的训练数据是从原始训练数据集随机抽样获得的一个子集，二是构建决策树时用于节点划分的属性来源于从属性集合随机抽取的子集，因此，需兼顾随机森林性能及运算时间复杂度，定义合适的决策树数目和决策树节点划分时属性的数目。
本实施例采用滑动窗口技术表征地球化学变量分布的空间异质性，定义滑动窗口半径序列为ε1≤ε2≤…εT，T为滑动窗口的个数，之后将借助紧凑度指数动态确定最佳滑动窗口尺寸，滑动窗口采用方形窗口，滑动窗口半径以栅格单元为单位。如图2所示，当窗口半径为εt(1≤t≤T)时，滑动窗口的大小为(2εt+1)×(2εt+1)。本实施例步骤S32中利用紧凑度指数确定最佳滑动窗口半径包括以下步骤：按照从左到右，从上到下的顺序遍历栅格单元，对于当前需预测地球化学元素含量值的栅格位置s，计算滑动窗口半径为εt时的紧凑度指数，紧凑度指数计算公式如下：
Figure BDA0003396269870000101
其中，Ct代表落入半径为εt窗口内的地球化学观测样本集合，Nt是样本数目,Yj′和Yk′代表地球化学样本；随着窗口半径不断变化，Dt取值也在不断发生变化，若相邻两个窗口内样本几乎来自同一个统计分布总体，则Dt和Dt+1差别不大；若样本明显来自两个不同统计总体，则Di和Di+1在整个变化趋势中会出现突变，计算相邻紧凑度指数Dt(t＝1,2,…T)之差Dt+1-Dt，设差值最大时滑动窗口对应的半径为εt，则优选出的最佳半径为ε*(s)＝εt。本实施例步骤S33计算空间邻域观测点的空间权重包括以下步骤：获取位置s处最佳窗口内的训练样本{(X1(sl′),X2(sl′),…,XK(sl′)；Y(sl′)),l′＝1,2,…L′}(L′<L)，计算各个样本所处位置sl′(l′＝1,2,…l′)与位置s的欧氏距离，计算公式为：
Figure BDA0003396269870000102
将计算得到的欧氏距离代入利用全局样本拟合的变差函数模型γ(h)；
求解克立格方程组，公式为：
Figure BDA0003396269870000103
上式共包含L′+1个方程，式中未知数为λl′(l′＝1,2,…L′)和μ共L′+1个，其中μ为拉格朗日算子，求解该方程组得到的λl′即为邻域观测点的权重。
本实施例步骤S34训练地统计加权随机森林模型包括以下步骤：对于每个回归决策树，采用改进的目标损失函数进行训练，公式为：
Figure BDA0003396269870000104
其中，
Figure BDA0003396269870000111
Figure BDA0003396269870000112
Figure BDA0003396269870000113
X代表构建决策树的切分属性或因子，v代表构建决策树相应属性的切分点，Yleft和Yright分别代表左子节点和右子节点所包含的地球化学观测样本点取值集合，由划分属性Xk及该属性切分点vkw所得，Ya和Yb分别代表左右子节点集合中的样本值，Nleft和Nright分别代表左子节点和右子节点中观测样本的数目。该步骤和经典随机森林算法的不同之处在于：
Figure BDA0003396269870000114
和
Figure BDA0003396269870000115
的计算不再是简单地求取算数平均数，而是地统计加权的平均数，权重反映了地学变量空间分布的相关性特征。
本实施例采用地统计方法为随机森林中各决策树叶子节点样本进行加权，改进随机森林的损失函数，更好地反映地学变量的空间相关性特征。将滑动窗口技术与随机森林相结合，同时借助紧凑度指数动态确定最佳滑动窗口尺寸，更好地反映地学变量的空间变异性特征，能提升预测精度。
以上所述的具体实施方式，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施方式而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.基于地统计加权随机森林的地球化学变量空间预测方法，其特征在于，包括以下步骤：
S1、获取训练随机森林模型所需的训练数据，所述训练数据包括输入变量和模型输出变量，输入变量包括地质要素和遥感要素，地质要素包括从区域地质数据提取的岩性类别数据和从区域地质数据提取的线环构造矢量数据转换的断裂密度，遥感要素包括地形坡度和植被覆盖度，模型输出变量为已知地球化学元素含量观测数据；其中，输入变量和输出变量采用栅格格式存储；
S2、根据输入变量与已知地球化学元素含量观测数据之间的相关系数及输入变量之间的相关性筛选输入变量作为预测因子，将已知地球化学元素含量观测数据与预测因子作为训练样本；
S3、基于训练样本训练随机森林模型，并采用训练好的随机森林模型对待预测位置的地球化学元素含量进行预测。
2.根据权利要求1所述的基于地统计加权随机森林的地球化学变量空间预测方法，其特征在于，所述步骤S1中将线环构造矢量数据转换为断裂密度的转换公式为：
Figure FDA0003396269860000011
其中，i表示栅格的行，j表示栅格的列，ρij为第i行第j列栅格的断裂密度，Lij为对应栅格内线性断裂和环形断裂的总长度，S为栅格的面积；
所述步骤S1中岩性类别数据为将区域地质数据中的岩性分布数据转换的类别变量，即
cij∈{Ck，k＝1，2，…n}
其中，n为岩性类别总数，Ck表示n种岩性类别中的第k种岩性，cij代表栅格位置(i，j)处的岩性类别。
3.根据权利要求1所述的基于地统计加权随机森林的地球化学变量空间预测方法，其特征在于，所述步骤S1还包括对遥感数据进行辐射校正处理，并利用以下公式计算归一化植被指数：
Figure FDA0003396269860000012
其中，i表示栅格的行，j表示栅格的列，NDVIij为栅格位置(i，j)处的植被指数，
Figure FDA0003396269860000013
为近红外波段处的反射率，
Figure FDA0003396269860000014
为红光波段处的反射率。
4.根据权利要求1所述的基于地统计加权随机森林的地球化学变量空间预测方法，其特征在于，所述步骤S1还包括以下已知地球化学元素含量观测数据处理步骤：
迭代剔除位于区间[μm-3σm，μm+3σm]之外的已知观测值，其中，μm为第m个元素的均值，σm为迭代过程中第m个元素的标准差；
对低于检测限的观测值用相应元素检测限的一半代替；
对缺失值用对应元素的均值补全；
采用中心对数比变换法对已知多元地球化学元素含量观测数据进行处理，处理公式如下：
Figure FDA0003396269860000021
其中，xm和zm分别代表第m个元素变换前后的值，D代表元素个数。
5.根据权利要求1所述的基于地统计加权随机森林的地球化学变量空间预测方法，其特征在于，所述步骤S1在获取训练随机森林模型所需的训练数据后还包括以下步骤：
采用待预测区域掩膜文件对栅格数据进行裁剪，并统一数据的坐标系。
6.根据权利要求1所述的基于地统计加权随机森林的地球化学变量空间预测方法，其特征在于，所述步骤S2中根据输入变量与已知地球化学元素含量观测数据之间的相关系数及输入变量之间的相关性筛选输入变量包括以下步骤：
计算各个输入变量与已知地球化学元素含量观测数据之间的相关系数并剔除弱相关性输入变量，再计算各个输入变量之间的相关性，任意两个输入变量之间相关性大于设定阈值时保留其一。
7.根据权利要求6所述的基于地统计加权随机森林的地球化学变量空间预测方法，其特征在于，所述计算各个输入变量与已知地球化学元素含量观测数据之间的相关系数、计算各个输入变量之间的相关性均包括以下步骤：对于连续型数据，计算Pearson相关系数，对于离散型数据，计算Spearman秩相关系数；
计算各个输入变量与已知地球化学元素含量观测数据之间的相关系数时若相关系数绝对值|r|＜0.3，则剔除该属性特征。
8.根据权利要求1所述的基于地统计加权随机森林的地球化学变量空间预测方法，其特征在于，所述步骤S3包括以下步骤：
S31、定义模拟栅格图层和滑动窗口半径序列，以及确定随机森林模型中决策树的数目和决策树节点划分时属性的数目；
S32、按照从左到右，从上到下的顺序遍历栅格单元，对于当前需预测地球化学元素含量值的栅格位置，利用紧凑度指数确定最佳滑动窗口半径；
S33、计算空间邻域观测点的空间权重：提取最佳窗口内的训练样本，通过全局变差函数模型和克立格方程组求解训练样本的权重；
S34、训练地统计加权随机森林模型：采用地统计加权的均方根误差函数作为改进的目标函数，训练随机森林回归模型；
S35、预测地球化学元素含量未知点处的取值：将待预测位置s处的预测因子值输入步骤S34训练好的随机森林模型，预测地球化学元素含量值；
S36、重复步骤S32～步骤S35，遍历结束即得完整的地球化学元素含量空间分布图。
9.根据权利要求8所述的基于地统计加权随机森林的地球化学变量空间预测方法，其特征在于，所述步骤S31中定义滑动窗口半径序列为ε1≤ε2≤…εT，滑动窗口采用方形窗口，当窗口半径为εt(1≤t≤T)时，滑动窗口的大小为(2εt+1)×(2εt+1)；
所述步骤S32中利用紧凑度指数确定最佳滑动窗口半径包括以下步骤：
对于当前需预测地球化学元素含量值的栅格位置s，计算滑动窗口半径为εt时的紧凑度指数，紧凑度指数计算公式如下：
Figure FDA0003396269860000031
其中，Ct代表落入半径为εt窗口内的地球化学观测样本集合，Nt是样本数目，Yj′和Yk′代表地球化学样本；若相邻两个窗口内样本几乎来自同一个统计分布总体，则Dt和Dt+1差别不大；若样本明显来自两个不同统计总体，则Di和Di+1在整个变化趋势中会出现突变，计算相邻紧凑度指数Dt(t＝1，2，…T)之差Dt+1-Dt，设差值最大时滑动窗口对应的半径为εt，则优选出的最佳半径为ε*(s)＝εt；
所述步骤S33计算空间邻域观测点的空间权重包括以下步骤：
获取位置s处最佳窗口内的训练样本
{(X1(sl′)，X2(sl′)，…，xK(sl′)；Y(sl′))，l′＝1，2，…L′}(L′＜L)，计算各个样本所处位置sl′(l′＝1，2，…l′)与位置s的欧氏距离，计算公式为：
Figure FDA0003396269860000032
将计算得到的欧氏距离代入利用全局样本拟合的变差函数模型γ(h)；
求解克立格方程组，公式为：
Figure FDA0003396269860000033
上式共包含L′+1个方程，式中未知数为λl′(l′＝1，2，…L′)和μ共L′+1个，其中μ为拉格朗日算子，求解该方程组得到的λl′即为邻域观测点的权重；
所述步骤S34训练地统计加权随机森林模型包括以下步骤：
对于每个回归决策树，采用改进的目标损失函数进行训练，公式为：
Figure FDA0003396269860000041
其中，
Figure FDA0003396269860000042
Figure FDA0003396269860000043
Figure FDA0003396269860000044
X代表构建决策树的切分属性或因子，v代表构建决策树相应属性的切分点，Yleft和Yright分别代表左子节点和右子节点所包含的地球化学观测样本点取值集合，Ya和Yb分别代表左右子节点集合中的样本值，Nleft和Nright分别代表左子节点和右子节点中观测样本的数目。
10.根据权利要求1～9中任意一项所述的基于地统计加权随机森林的地球化学变量空间预测方法，其特征在于，所述步骤S2与步骤S3之间还包括以下步骤：将训练样本分为训练数据集和验证数据集，所述步骤S3基于训练数据集训练随机森林模型；所述步骤S3之后还包括以下步骤：提取验证数据集中各验证样本所处位置处的地球化学变量预测值，建立验证样本观测值与预测值数据对(Yg，eval，Yg，obs)(g＝1，2，…Neval)的线性回归模型，得到决定系数R2，据此定量评价随机森林模型的预测效果；其中，Yg，eval代表第g个点处地球化学变量的模型预测值，Yg，obs代表第g个点处地球化学变量的参考值，Neval代表验证点的个数。
CN202111483212.6A 2021-12-07 2021-12-07 基于地统计加权随机森林的地球化学变量空间预测方法 Active CN114139819B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111483212.6A CN114139819B (zh) 2021-12-07 2021-12-07 基于地统计加权随机森林的地球化学变量空间预测方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111483212.6A CN114139819B (zh) 2021-12-07 2021-12-07 基于地统计加权随机森林的地球化学变量空间预测方法

Publications (2)

Publication Number Publication Date
CN114139819A true CN114139819A (zh) 2022-03-04
CN114139819B CN114139819B (zh) 2023-05-23

Family

ID=80384286

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111483212.6A Active CN114139819B (zh) 2021-12-07 2021-12-07 基于地统计加权随机森林的地球化学变量空间预测方法

Country Status (1)

Country Link
CN (1) CN114139819B (zh)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114596498A (zh) * 2022-05-10 2022-06-07 湖北省地质调查院 一种地球化学采样盲区的赋值方法、系统及存储介质
CN116821816A (zh) * 2023-05-17 2023-09-29 河南农业大学 一种基于加权随机森林的干热风预测方法
CN117932456A (zh) * 2024-03-22 2024-04-26 中国科学院地理科学与资源研究所 一种顾及空间异质性的集成空间预测方法
CN117932456B (zh) * 2024-03-22 2024-06-07 中国科学院地理科学与资源研究所 一种顾及空间异质性的集成空间预测方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109142679A (zh) * 2018-08-13 2019-01-04 中国科学院东北地理与农业生态研究所 基于人工神经网络克里金插值的森林土壤养分的空间预测方法
CN109711597A (zh) * 2018-11-14 2019-05-03 东莞理工学院 一种基于分层随机森林模型的铜镍硫化物矿床成矿预测方法
CN111401644A (zh) * 2020-03-19 2020-07-10 南京国准数据有限责任公司 一种基于神经网络的降水降尺度空间预测方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109142679A (zh) * 2018-08-13 2019-01-04 中国科学院东北地理与农业生态研究所 基于人工神经网络克里金插值的森林土壤养分的空间预测方法
CN109711597A (zh) * 2018-11-14 2019-05-03 东莞理工学院 一种基于分层随机森林模型的铜镍硫化物矿床成矿预测方法
CN111401644A (zh) * 2020-03-19 2020-07-10 南京国准数据有限责任公司 一种基于神经网络的降水降尺度空间预测方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
齐雁冰;王茵茵;陈洋;刘姣姣;张亮亮;: "基于遥感与随机森林算法的陕西省土壤有机质空间预测" *

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114596498A (zh) * 2022-05-10 2022-06-07 湖北省地质调查院 一种地球化学采样盲区的赋值方法、系统及存储介质
CN114596498B (zh) * 2022-05-10 2022-07-29 湖北省地质调查院 一种地球化学采样盲区的赋值方法、系统及存储介质
CN116821816A (zh) * 2023-05-17 2023-09-29 河南农业大学 一种基于加权随机森林的干热风预测方法
CN116821816B (zh) * 2023-05-17 2024-05-28 河南农业大学 一种基于加权随机森林的干热风预测方法
CN117932456A (zh) * 2024-03-22 2024-04-26 中国科学院地理科学与资源研究所 一种顾及空间异质性的集成空间预测方法
CN117932456B (zh) * 2024-03-22 2024-06-07 中国科学院地理科学与资源研究所 一种顾及空间异质性的集成空间预测方法

Also Published As

Publication number Publication date
CN114139819B (zh) 2023-05-23

Similar Documents

Publication Publication Date Title
Zhu et al. Intelligent logging lithological interpretation with convolution neural networks
DeFries et al. Multiple criteria for evaluating machine learning algorithms for land cover classification from satellite data
CN109709603B (zh) 地震层位识别与追踪方法、系统
Ercanoglu Landslide susceptibility assessment of SE Bartin (West Black Sea region, Turkey) by artificial neural networks
Lemoine et al. Testing the applicability of correlations between topographic slope and VS 30 for Europe
Alnahwi et al. Mineralogical composition and total organic carbon quantification using x-ray fluorescence data from the Upper Cretaceous Eagle Ford Group in southern Texas
CN114139819B (zh) 基于地统计加权随机森林的地球化学变量空间预测方法
Zhang et al. Extracting dispersion curves from ambient noise correlations using deep learning
CN116665067B (zh) 基于图神经网络的找矿靶区优选系统及其方法
Eberle et al. Integrated data analysis for mineral exploration: A case study of clustering satellite imagery, airborne gamma-ray, and regional geochemical data suites
CN110927793B (zh) 一种基于序贯随机模糊模拟的储层预测方法及系统
CN110222832A (zh) 长江口盐沼湿地大型底栖动物栖息地模拟预测方法
CN111582387A (zh) 一种岩石光谱特征融合分类方法及系统
Talebi et al. Towards geostatistical learning for the geosciences: A case study in improving the spatial awareness of spectral clustering
CN114997501A (zh) 基于样本失衡的深度学习矿产资源分类预测方法及系统
Süzen Data driven landslide hazard assessment using geographical information systems and remote sensing
Phelps et al. Exploring viable geologic interpretations of gravity models using distance-based global sensitivity analysis and kernel methods
Maiti et al. A deep CNN-LSTM model for predicting interface depth from gravity data over thrust and fold belts of North East India
CN115511214A (zh) 基于多尺度样本不均的矿产预测方法及系统
CN114998719A (zh) 一种基于深度学习和多源遥感数据的林火预测方法
Paasche Post‐inversion integration of Disparate Tomographic Models by Model Structure Analyses
Yan et al. Depth-to-bedrock map of China at a spatial resolution of 100 meters
Agbor et al. Comparative analysis of non-linear artificial neural networks and maximum likelihood algorithms in forest cover studies
Guo et al. Exploring the potential of HySpex hyperspectral imagery for extraction of copper content
Burengengwa Comparison of approaches for spatial interpolation of weather data on a specific date

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant