复杂地形条件下山区土石界面插值策略与预测方法研究

Interpolation Strategy and Prediction Method of Soil-Rock Interface in Mountainous Areas under Complex Topography

摘要

土石界面预测对山区机场建设具有重要意义.现有方法通常基于钻孔数据直接估计土石界面高程，但在山区地形起伏剧烈、界面高程差异显著的条件下，其结果准确性难以满足工程需求.提出一种基于地形高程与界面深度的插值策略，土石界面高程可表示为已知地形高程减去土石界面深度.该策略只需预测土石界面深度，因此降低了地形起伏对预测结果的影响.系统对比了该策略在多种预测方法中的表现，包括反距离加权、径向基函数核回归、高斯过程回归等核方法以及神经网络方法，如多层感知机和Kolmogorov-Arnold网络.山区机场工程案例表明，所提策略适用于不同预测方法，并显著提升土石界面预测准确性，为复杂地形条件下的机场选址、土石方工程量估算及施工方案优化提供技术支撑.

Abstract

The prediction of the soil-rock interface is crucial for airport construction in mountainous areas. Conventional methods typically rely on borehole data to directly estimate the elevation of the soil-rock interface. However, when topographic variability is pronounced and interface-elevation contrasts are large, the accuracy of such estimates is often insufficient to satisfy engineering requirements. In this paper it introduces an interpolation strategy that incorporates terrain elevation and interface depth. In this framework, the elevation of the soil-rock interface is expressed as the terrain elevation minus the depth of the interface. By shifting the prediction task to estimating interface depth alone, this approach reduces the influence of topographic variability on prediction accuracy. The performance of the proposed strategy is systematically evaluated using various prediction methods, including kernel methods such as inverse distance weighting, radial basis function kernel regression, and Gaussian process regression, as well as neural network approaches (e.g., multilayer perceptron and Kolmogorov-Arnold Networks). Case studies from airport projects in mountainous regions demonstrate that the strategy can be readily integrated with different prediction methods and substantially improves the accuracy of soil-rock interface predictions. The findings provide technical support for airport site selection, earthwork volume estimation, and construction planning in complex terrains.

Graphical abstract

关键词

山区 / 地形 / 土石界面 / 插值策略 / 预测方法 / 工程地质学.

Key words

mountainous area / topography / soil-rock interface / interpolation strategy / prediction method / engineering geology

0 引言

近年来，我国越来越多的机场选址于地形复杂、地质条件多变的山区地带（谢春庆等，2015）.受地形起伏影响，这类山区机场施工通常需进行大规模土石方开挖与填筑（王巧，2011；潘凯等，2017）.因此，确定土石界面对于地基承载力计算（仇志龙等，2025）、场地安全（苏正洋等，2025）和优化工程量估算、施工成本控制等具有关键意义（Zhang and Dasaka， 2010）.

开展山区机场工程勘察时，钻探仍是获取地层信息的主要手段.钻孔资料可用于确定土层厚度、基岩埋深及土石分界点位置，为土石界面预测提供依据（谢春庆等，2013）.工程上通常采用直线连接钻孔揭示的土石分界点以重建土石界面.但是考虑到钻探成本和地形条件，实际工程中钻孔一般稀疏布置（Bai et al.， 2024），进而使上述土石界面重建方法的预测准确性受限.因此，在有限钻孔条件下合理预测复杂山区土石界面是亟待解决的难题.

为此，文献中提出并应用了多种插值方法.Samui et al.（2015）比较了自适应神经模糊推理系统（ANFIS）、多元自适应样条回归（MARS）、普通克里金（OK）和反距离加权法（IDW）的岩层深度预测准确性；Li et al.（2016）基于条件随机场（CRF）方法预测岩层深度，可量化结果不确定性；Qi et al.（2020）开展了二维剖面内土石界线预测研究，对比了多项式回归、样条回归、样条插值法以及条件随机场法；Qi et al.（2021，2022）进一步开展了三维空间内土石界面的预测研究，对比了MARS、CRF与薄板样条插值法（TPSI）；Deng et al.（2023）采用高斯过程回归（GPR）方法预测了岩层深度.表1汇总了前述文献的关键信息，包括钻孔数量、场地高程差、均方根误差（RMSE）.由表1可见，现有研究主要针对高程差较小（≤64.1 m）的场地.土石界面预测的准确性主要受钻孔数量（或密度）与场地高程差共同影响：RMSE整体上随钻孔密度增加而降低、随高程差增大而升高.因此，在地形起伏剧烈且钻孔分布稀疏的山区机场场地，土石界面的准确预测问题尚未得到有效解决.

虽然在山区机场场地勘察过程中钻孔数量有限，但遥感与地面测量可获取高精度数字高程模型DEM，为研究区提供连续、致密的地表高程信息（张世殊等，2025）.已有研究表明岩土体界面高程与地表高程具有较强相关性（Cordonnier et al.， 2018）.除钻孔外，地表高程数据（如DEM）可作为先验信息用于土石界面预测.据此，本文将DEM与稀疏钻孔联合用于土石界面预测，通过引入地形先验约束（Huang et al.，2025），以缓解钻孔稀疏性并提升预测准确性.

本文提出了一种基于地形高程与界面深度插值策略：先预测土石界面深度，随后从DEM提取预测点高程，两者相减得到界面高程，从而在一定程度上降低地形起伏对预测的影响.为评估所提策略的预测效果，本文将其与传统基于界面高程的插值方法进行比较.在反距离加权（IDW）、径向基函数（RBF）核回归、高斯过程回归（GPR）、多层感知机（MLP）与Kolmogorov-Arnold网络（KAN）5种典型预测模型中分别应用两种策略，并基于山区机场工程数据对两种策略的预测表现进行对比分析.本文结构安排如下：先阐述所提土石界面插值策略，再介绍预测方法及对应评价指标，继而开展工程案例研究，最后总结并提出应用建议.

1 土石界面插值策略

1.1　基于界面高程的插值策略

在工程实践中，土石分界的判定通常依托钻孔资料，并结合相关规范所规定的室内试验、原位测试结果及土石材料分级标准加以确定.例如，《民用机场勘测规范》（MH/T5025）则按照开挖难易程度和工程性质对材料进行分级分类（表2）.如图1所示，传统的土石界面预测流程分为两步：首先基于土体分类规范确定钻孔位置的土石分界点高程；然后选择适当预测方法预测未知区域的土石分界点高程，直接生成密集土石分界点.本研究将其称为“基于界面高程的插值策略”.然而，在山区机场工程应用中，因钻孔稀疏且场地高程差显著，该策略的预测准确性难以保障.

1.2　基于地形高程与界面深度的插值策略

为了克服基于界面高程插值策略的局限性，本文提出了一种基于地形高程与界面深度的插值策略（图2）.所提策略分为3步：先从钻孔数据提取土石界面深度，然后预测未知区域的土石界面深度，最后以预测点的地形高程（来自DEM）减去预测深度，获得相应位置的土石界面高程，从而生成密集的界面点集.相较于传统高程插值策略，所提策略创新之处在于：（1）土石界面深度预测过程不依赖界面高程数据，一定程度降低地形起伏对预测准确性的影响；（2）充分利用地表高程信息，在稀疏钻孔条件下提升预测准确性.

2 土石界面预测方法与评估指标

为系统评估所提插值策略，本研究选取了文献中多种典型土石界面预测方法，并采用3种评估指标，开展所提插值策略和基于界面高程插值策略的对比研究.下文依次介绍插值方法与评价指标.

2.1　预测方法

考虑不同预测方法的建模范式差异，本文选取了两类方法：（1）以核方法为代表的非参数方法；（2）以神经网络为代表的参数化回归方法（张瀚等，2024）.核方法以核函数度量待预测位置与钻孔的相似性，据此计算权重，对钻孔处的土石界面高程（或深度）进行加权平均，以估计预测位置的土石界面高程（或深度）；常见实现包括反距离加权（IDW）、径向基函数（RBF）核回归和高斯过程回归（GPR）.神经网络方法则基于可学习的参数化结构，通过迭代训练不断优化权重或函数模块，从而逼近输入特征与土石界面高程（或深度）之间的复杂非线性关系；典型方法包括多层感知机（MLP）和Kolmogorov-Arnold网络（KAN）.

2.1.1　核方法

核方法用于预测任务时，对应预测方程表达式如下：

\widehat{y}={\displaystyle \sum _{i=\mathrm{1}}^n}\omega \left(\underline{\widehat{\mathit{x}}},{\underline{\mathit{x}}}_i\right)y_i，（1）

式中：\underline{\widehat{\mathit{x}}}为预测点特征，本研究特指预测点的平面坐标；{\underline{\mathit{x}}}_i为第i个观测点（钻孔）的特征，即其平面坐标；\widehat{y}为预测值，代表预测点的高程（或深度）；y_i为第i个观测点的土石界面高程（或深度）；\omega \left(\underline{\widehat{\mathit{x}}},{\underline{\mathit{x}}}_i\right)为权重函数；本文将单次预测所用的观测点数记为n：在IDW中，n=k（目标点的近邻观测点数）；在RBF核回归和GPR中，n=N（全部观测点数）.

反距离加权法（IDW）采用反距离幂次权重函数（李章林等，2015；Chen et al.，2015），记作{\omega }_I\left(\underline{\widehat{\mathit{x}}},{\underline{\mathit{x}}}_i\right)，其表达式如下：

{\omega }_I\left(\underline{\widehat{\mathit{x}}},{\underset{\_}{x}}_i\right)=\left(\frac{\mathrm{1}}{{‖\underline{\widehat{\mathit{x}}}-{\underline{\mathit{x}}}_i‖}^p}\right)/{\displaystyle \sum _{i=\mathrm{1}}^k}\left(\frac{\mathrm{1}}{{‖\underline{\widehat{\mathit{x}}}-{\underline{\mathit{x}}}_i‖}^p}\right),（2）

式中：‖\underline{\widehat{\mathit{x}}}-{\underline{\mathit{x}}}_i‖为预测点和其近邻观测点之间的欧式距离（i=1， 2， …， N）；p为幂指数，可以是任意正实数或0.

径向基函数（RBF）核回归采用径向基函数核（Murphy，2012，2022；扶金铭等，2024），其权重函数{\omega }_R\left(\underline{\widehat{\mathit{x}}},{\underline{\mathit{x}}}_i\right)为：

{\omega }_R\left(\underline{\widehat{\mathit{x}}},{\underset{\_}{x}}_i\right)=\left(\mathrm{e}\mathrm{x}\mathrm{p}\left[\frac{-{\left(\underline{\widehat{\mathit{x}}}-{\underline{\mathit{x}}}_i\right)}^{\mathrm{2}}}{\mathrm{2}{\lambda }^{\mathrm{2}}}\right]\right)/{\displaystyle \sum _{i=\mathrm{1}}^N}\left(\mathrm{e}\mathrm{x}\mathrm{p}\left[\frac{-{\left(\underline{\widehat{\mathit{x}}}-{\underline{\mathit{x}}}_i\right)}^{\mathrm{2}}}{\mathrm{2}{\lambda }^{\mathrm{2}}}\right]\right)，（3）

式中：\lambda为带宽，带宽越小，权重主要分配给距离预测点较近的观测点；其他符号的含义参考公式（2）.

高斯过程回归（GPR）作为一种贝叶斯非参数方法，不仅能够预测土石界面高程或深度，还能量化预测的不确定性（Wang，2021）.根据高斯过程（GP）定义，假设预测点高程（或深度）y^{\mathrm{*}}和所有观测点高程（或深度）向量\underline{\mathit{y}}服从如下联合高斯分布：

\left(\begin{array}{c}\underline{\mathit{y}}\\ y^{\mathrm{*}}\end{array}\right)~\mathrm{{\rm N}}\left(\left(\begin{array}{c}\underline{\mathit{\mu }}\\ {\mu }_{\mathrm{*}}\end{array}\right),\left(\begin{array}{cc}{\mathbf{{\rm K}}}_{\underline{\mathit{X}},\underline{\mathit{X}}}& {\mathbf{{\rm K}}}_{\mathrm{*},\underline{\mathit{X}}}\\ {\mathbf{{\rm K}}}_{\mathrm{*},\underline{\mathit{X}}}^T& k_{\mathrm{*}\mathrm{*}}\end{array}\right)\right)，（4）

式中：\underline{\mathit{y}}=[y_{\mathrm{1}},y_{\mathrm{2}},...,y_N{]}^T；\underset{\_}{\mu }与{\mu }_{\mathrm{*}}分别为观测点的均值向量和预测点均值，一般假设为0；\underline{\mathit{X}}=[{\underline{\mathit{x}}}_{\mathrm{1}},{\underline{\mathit{x}}}_{\mathrm{2}},...,{\underline{\mathit{x}}}_N]，{\mathbf{{\rm K}}}_{\underline{\mathit{X}},\underline{\mathit{X}}}=K\left(\underline{\mathit{X}},\underline{\mathit{X}}\right)，{\mathbf{{\rm K}}}_{\mathrm{*},\underline{\mathit{X}}}=K\left(\underline{\widehat{\mathit{x}}},\underline{\mathit{X}}\right)，k_{\mathrm{*}\mathrm{*}}=K\left(\underline{\widehat{\mathit{x}}},\underline{\widehat{\mathit{x}}}\right)，K代表正定核函数，用于计算两点间协方差.本文选取一种常用的K，即高斯核函数：

K\left({\underline{\mathit{x}}}_i,{\underline{\mathit{x}}}_j\right)=C\mathrm{e}\mathrm{x}\mathrm{p}\left[\frac{-{\left({\underline{\mathit{x}}}_i-{\underline{\mathit{x}}}_j\right)}^{\mathrm{2}}}{\mathrm{2}{l}^{\mathrm{2}}}\right]+{\sigma }_n^{\mathrm{2}}，（5）

式中：C为核函数方差，表示函数在整体上的变化幅度；{\sigma }_n^{\mathrm{2}}为噪声方差；l为长度尺度，控制K的宽度.按高斯分布条件化的标准规则（Murphy，2012），y^{\mathrm{*}}的后验分布仍是高斯分布：

p\left(y^{\mathrm{*}}|\underline{\widehat{\mathit{x}}},\underline{\mathit{X}},\underline{\mathit{y}}\right)=N\left(y^{\mathrm{*}}\left|{\mathbf{K}}_{\mathrm{*},\underline{\mathit{X}}}^T{\mathbf{{\rm K}}}_{\underline{\mathit{X}},\underline{\mathit{X}}}^{{}^{-\mathrm{1}}}\underline{\mathit{y}},k_{\mathrm{*}\mathrm{*}}-{\mathbf{{\rm K}}}_{\underline{\mathit{X}},\underline{\mathit{X}}}^{{}^{-\mathrm{1}}}{\mathbf{K}}_{\mathbf{*},\underline{\mathit{X}}}\right.\right)，（6）

因此，y^*的后验最可能值（\widehat{y}）等于其后验分布均值，即

\widehat{y}={\mathbf{K}}_{\mathrm{*},\underline{\mathit{X}}}^T{\mathbf{{\rm K}}}_{\underline{\mathit{X}},\underline{\mathit{X}}}^{{}^{-\mathrm{1}}}\underline{\mathit{y}}={\displaystyle \sum _{i=\mathrm{1}}^N}\left({\mathbf{K}}_{\mathrm{*},\underline{\mathit{X}}}^T{\mathbf{{\rm K}}}_{\underline{\mathit{X}},{\underline{\mathit{x}}}_i}^{{}^{-\mathrm{1}}}\right)y_i，（7）

式中：{\mathbf{{\rm K}}}_{\underline{\mathit{X}},{\underline{\mathit{x}}}_i}=K\left(\underline{\mathit{X}},{\underline{\mathit{x}}}_i\right).令{\omega }_G\left(\underline{\widehat{\mathit{x}}},{\underline{\mathit{x}}}_i\right)代表GPR预测方程的权重函数，则{\omega }_G\left(\underline{\widehat{\mathit{x}}},{\underline{\mathit{x}}}_i\right)={\mathbf{K}}_{\mathrm{*},\underline{\mathit{X}}}^T{\mathbf{{\rm K}}}_{\underline{\mathit{X}},{\underline{\mathit{x}}}_i}^{{}^{-\mathrm{1}}}.

2.1.2　神经网络方法

多层感知机（MLP）通过多层非线性映射和反向传播优化参数，用于刻画复杂输入输出关系（杜昕等，2020）.该网络包含输入层、输出层以及多个隐含层，以神经元为最小单位，将输入的多个节点信息非线性映射到单一输出上.f\left(x\right)代表激活函数，可以为Sigmoid、Tanh或者ReLU函数中的一种.第1层神经元的解析表达式：

z_j^{\mathrm{1}}=f\left[{\displaystyle \sum _{t=\mathrm{1}}^{N_x}}\left({\omega }_{t,j}^{\mathrm{1}}{x}_t+b_{t,j}^{\mathrm{1}}\right)\right]，（8）

式中：z_j^{\mathrm{1}}为代表第1个隐含层上的第j个神经元，N_x为输入层维度，x_t为特征\underline{\mathit{x}}中第t个分量，{\omega }_{t,j}^{\mathrm{1}}是输入层第t个分量与第1个隐含层上第j个神经元的连接权重，b_{t,j}^{\mathrm{1}}为偏置常数.向后继续传播时，以z_j^{\mathrm{1}}作为新的输入，对公式（8）依次递归，由公式（9）进一步得到输出层的解析表达式：

\widehat{y}=f\left[{\displaystyle \sum _{t=\mathrm{1}}^{N_{M-\mathrm{1}}}}\left({\omega }_{t,j}^M{z}_t^M+b_{t,j}^M\right)\right]，（9）

式中：\widehat{y}代表输出层的预测值，M为神经网络隐含层总数，N_{M-\mathrm{1}}代表第M-\mathrm{1}个隐含层包含的神经元总数.

与MLP不同，Kolmogorov-Arnold网络（KAN）将用被参数化的函数代替权重参数（Liu et al.， 2024；Somvanshi et al.， 2026）.KAN中的每一条连接可表示为单变量连续函数加和的形式，通过一组参数化的基函数（样条函数、多项式函数等）对这些单变量进行近似.KAN网络输出层的解析表达式可表示为：

\widehat{y}=f\left[{\displaystyle \sum _{j=\mathrm{1}}^{\mathrm{2}{N}_x+\mathrm{1}}}{\mathrm{\Phi }}_t\left({\displaystyle \sum _{t=\mathrm{1}}^{N_x}}{\varphi }_{t,j}\left(x_t\right)\right)\right],（10）

式中：{\varphi }_{t,j}是对输入做非线性变换的单变量函数，{\mathrm{\Phi }}_t表示在隐藏层做进一步变换的单变量函数.

2.2　评估指标

为全面、客观地评估不同预测方法的表现，本文采用引入均方根误差（RMSE）、归一化均方根误差（nRMSE）与决定系数（R²）3项评估指标.RMSE表示预测值与真实值的平均绝对误差，其计算公式为：

RMSE=\sqrt[]{\frac{\mathrm{1}}{N}{\displaystyle \sum _{i=\mathrm{1}}^N}{\left({\widehat{y}}_i-y_i\right)}^{\mathrm{2}}}，（11）

式中：N为预测点个数，\widehat{y}为预测值，y为已知值.由于山区地形起伏较大，RMSE可能难以全面反映预测准确性的相对水平.为此，本文提出归一化均方根误差（nRMSE），即

nRMSE=\sqrt[]{\frac{\mathrm{1}}{N}{\displaystyle \sum _{i=\mathrm{1}}^N}{\left({\widehat{y}}_i-y_i\right)}^{\mathrm{2}}}/\mathrm{\Delta }h，（12）

式中：\mathrm{\Delta }h为最大地形高差.nRMSE将\mathrm{\Delta }h作为归一化因子对RMSE进行标准化，从而更客观地比较高差较大场景下不同方法的误差水平.R²表示模型对观测数据变异性的解释能力，其计算公式为：

R^{\mathrm{2}}=\mathrm{1}-\frac{{\displaystyle \sum _{i=\mathrm{1}}^N}{\left(y_i-{\widehat{y}}_i\right)}^{\mathrm{2}}}{{\displaystyle \sum _{i=\mathrm{1}}^N}{\left(y_i-\overline{y}\right)}^{\mathrm{2}}}，（13）

式中：\overline{y}为观测数据的均值.R²值越接近1，说明模型拟合效果越好.

综上所述，本节简要介绍了典型土石界面预测方法与评价，下文将通过工程案例验证所提基于地形高程与界面深度的插值策略.

3 工程案例——某山区机场工程

3.1　案例简介

该机场位于我国西南山区，场地范围约2 000 m×3 500 m.图3a显示了研究区域的钻孔布置，共计473个钻孔.由于钻孔分布稀疏，仅基于钻孔资料难以准确重构土石界面.图3b为该场地数字高程模型（DEM），其显示地表高程分布在2 100~2 350 m之间，最大高差约250 m.综合钻孔和DEM，该场地土石界面高程范围为2 120~2 322 m，标准差52 m；对所有钻孔的地表高程与对应土石界面高程进行相关分析，相关系数为0.99，表明二者相关性极强.因此，本文在土石界面预测过程中引入地表高程信息，以约束预测结果.图3c为剖面1-1'（图3a）的工程地质剖面图.

3.2　二维土石界线预测

本文首先在纵剖面1-1'上进行土石界线（一维）插值，用于结果展示与方法比较.由于地形条件限制，1-1'剖面内钻孔的X坐标并不完全一致，但其最大偏差不超过10 m，而剖面沿Y方向的长度约为1 600 m.与剖面长度相比，该偏移量很小，因此，在工程上可认为这些钻孔位于同一剖面.图4展示了剖面1-1'内25个真实钻孔数据的岩性分布与土石分界点深度分布.强风化泥质粉砂岩与中等风化泥质粉砂岩的分界线代表土石分界线，其深度范围为1.9~25.0 m，标准差5.9 m.基于图4所示土石分界点数据，本研究分别采用两种插值策略，并结合不同预测方法（即IDW、RBF核回归、GPR、MLP与KAN）预测剖面1-1'内的土石分界线.

本文在预设参数网格上采用留一交叉验证（LOOCV）对各预测方法（GPR除外）进行参数选择与性能评估：以LOOCV的RMSE最小为准则确定最优参数配置，并在该最优参数配置下计算nRMSE与R².GPR为特例，其最优参数采用最大化对数边缘似然直接估计，随后在该最优参数配置下，用LOOCV计算RMSE、nRMSE与R².表3给出了各方法的最优参数配置.此外，GPR噪声方差极小，对结果无显著影响，因此未单独列出.MLP激活函数为ReLU函数；KAN学习函数为B-样条函数.神经网络方法的输入特征为平面X和Y坐标与地表高程，输出为界面高程或深度.为了防止神经网络方法过拟合，在所有全连接层加入了L2正则化项并采用早停机制.

图5分别展示了两种插值策略下5种预测方法的留一交叉验证结果，即钻孔处土石界面高程预测值.所提插值策略下，每种预测方法结果均接近对应真实值，其中未发现预测值超出地表的情况.而基于界面高程插值策略下，所有预测方法在1-1'剖面边界附近钻孔处土石界面高程预测值超出地表；此外，神经网络模型（MLP和KAN网络）只能捕捉土石界面的整体趋势，难以准确反映其局部起伏.图5结果说明了所提插值策略通过引入地形约束，提升了预测的合理性.

为了进一步定量评估预测方法，图6给出了两种插值策略下不同方法的留一验证RMSE（或nRMSE）对比结果.与传统基于界面高程的插值策略相比，采用所提策略后，各预测方法的RMSE均有所降低，降幅在38%~47%之间（视具体预测方法而定）.nRMSE的对比结果亦与RMSE保持一致，各预测方法均呈现相同的性能提升趋势.图6结果说明所提插值策略能够显著提升预测方法的准确性，其先预测土石界面深度，再由地形高程减去深度预测值估计土石界面高程，有助于降低地形起伏对预测的影响.表4给出了两种插值策略下不同方法的R²指标.所提插值策略下各方法预测结果的R²均大于基于界面高程的插值策略，表明所提插值策略能够提升各预测结果的拟合优度.

需要指出的是，不同评估指标在土石界面预测中侧重点各异.在本研究中R²指标普遍较高（>0.9），说明各策略与各方法均能较好地反映整体趋势，但无法揭示复杂地形下的绝对误差.相比之下，RMSE表示绝对误差，但其受地形起伏影响显著，在高差较大的山区场地往往数值较大，而在高差较小的场地即使预测准确度不高，RMSE也可能较小，因此有其局限性.而相对误差指标nRMSE在基于最大地形高差对RMSE进行了归一化，因此适合用于不同场地的预测结果对比.由此可见，单一指标难以全面评价预测效果，需结合多种指标并结合工程实际背景进行综合判断.

在剖面1-1'所在区域内，基于DEM选取1 734个预测点（该处界面深度未知），在所提插值策略下，分别采用IDW、RBF核回归、GPR、MLP、KAN对这些点的土石界面深度进行估计.作为对比，给定基于界面高程插值策略，在相同位置直接预测土石界面高程.图7显示了剖面1-1'土石界线预测结果.对非参数方法（IDW、RBF核回归、GPR），无论采用哪种插值策略，预测结果均能较好地刻画土石界线的空间起伏.相比之下，对参数化方法（MLP和KAN），在基于界面高程的策略下，结果整体呈平滑上升趋势，难以准确表征局部起伏，这可能与样本量较小（N=25）有关.所提插值策略将地形信息作为先验约束引入预测过程，在一定程度上缓解了数据稀缺导致的参数化方法对局部起伏欠拟合问题.

GPR方法不仅可以给出平均预测结果，还可量化预测结果的不确定性.图7展示了GPR预测结果的95%贝叶斯可信区间（95% Bayesian credible interval），即不同色带，其宽度=\mathrm{2}\times \mathrm{1.96}{\sigma }_{\mathrm{p}}表示预测结果的不确定性大小；{\sigma }_{\mathrm{p}}代表预测标准差.结果显示所提插值策略的色带宽度明显更窄，25个预测点的标准差平均值从5.9 m减小为3.9 m，由此说明所提插值策略有助于降低土石界面预测的不确定性，从而提高了预测结果的可靠性.

3.3　三维土石界面预测

为评估所提插值策略在三维场景中的适用性，本文将研究区域内473个钻孔的土石界面点作为观测数据，进行土石界面预测.这些土石界面点深度范围约0.6~42.0 m，标准差等于6.5 m.

本研究按照8︰2的比例划分了训练数据集与测试集，并在训练数据集内进行五折交叉验证.将训练数据集随机划分为5个大小相同的互斥子集.然后，每次用4个子集的并集作为训练集，余下的子集作为验证集.本研究的重点在于比较不同模型的预测性能，因此针对各模型均采用相同的训练/验证/测试集划分，并进行一次完整的交叉验证，确定最优模型参数；将该最优模型参数对应的预测值作为模型的预测结果，以此计算评估指标RMSE、nRMSE和R².预测方法IDW、RBF核回归、GPR、MLP、KAN对应的参数寻优策略参考3.2节.表5展示了不同预测方法的最优参数配置，MLP激活函数为ReLU函数，KAN学习函数为B-样条函数.

图8展示了在两种插值策略下，各预测方法RMSE（或nRMSE）的对比图.表6则列出了每种预测方法的预测性能指标R².在所提插值策略下，每种预测方法的RMSE均在6.5 m左右；这些结果均小于对应基于界面高程的插值策略，并且RMSE降幅在30%~57%之间（视方法而定）.nRMSE结果也呈现类似趋势.此外，所提插值策略下R²计算结果整体大于基于界面高程的插值策略.上述结果再次表明，所提插值策略能够提升预测结果的准确性.

为全面评估所提插值策略，本文将第3.2节的二维案例与本节的三维案例进行了对比.结果显示给定所提插值策略条件下，在三维案例中各预测方法的RMSE相较二维案例有所增大，但经由最大地形高差归一化后的nRMSE以及决定系数R²基本保持稳定.这表明所提策略能够在不同地形起伏与数据变异性条件下维持一致的预测性能，具有良好的普适性.因此，单独使用RMSE不足以全面评价插值策略，需结合nRMSE与R²等指标开展综合评估.作为对照，在基于界面高程的传统插值策略下，当研究区由二维拓展至三维时，各预测方法的结果更易出现nRMSE升高或R²降低的趋势，表明传统插值策略对最大高程差与数据变异性的敏感性更强，从而进一步凸显了所提策略的普适性.

由于所提插值策略下，IDW、RBF核回归、GPR、MLP、KAN方法的预测误差以及R²总体一致，本文仅以IDW为例展示场地三维土石界面预测结果.在区域内的DEM上选取了18 000预测点，并将IDW的参数设置为其最优参数配置（表6），IDW在不同插值策略下的预测结果如图9所示.图中红色实心点标示了界面预测结果超过地表高程的异常位置.采用所提插值策略时，IDW预测结果中未发现超出地表点，较好刻画了土石界面空间起伏.相反，采用基于界面高程的插值策略，在场地边界以及地形起伏较大区域出现众多超出地表的预测点.由此可见，所提插值策略在三维土石界面预测中提升了结果准确性以及合理性.

在该三维场景中，钻孔地表高程与对应土石界面高程的相关系数为0.99；地表最大高差和标准差分别为212.4 m和50.4 m，而界面深度的最大差值和标准差仅为41.4 m和6.5 m.由此可见，界面深度的变化幅度明显小于界面高程.所提策略在工程中取得较高准确性主要源于两点：（1）将预测对象由界面高程转为界面深度，可有效降低预测难度，从而提升预测精度；（2）界面高程由“地表高程-预测深度”得到，地表高程天然作为土石界面高程的上限约束，能确保预测界面始终位于地表以下，使结果更加合理.基于上述原因，该策略能够在复杂地形条件下保持较高的预测准确性，可为机场选址、土石方工程量估算及施工方案优化等提供可靠的技术支持.

本研究中地表高程与界面高程高度相关，在该情况下所提策略优势尤为显著.然而，该策略的适用性在一定程度上依赖于地形对界面起伏的控制作用.当地表高程与界面高程相关性减弱时（如岩溶发育区、断裂破碎带），该策略有效性和稳定性有待验证.未来研究可从以下方面展开：进一步评估该策略在地表高程与界面高程弱相关或不相关条件下的有效性，考虑引入更多反映地下结构的先验信息（如物探数据等）以增强稳定性，也可以考虑基于地质构造进行分区预测，以提升策略在复杂地质条件下的适用性和可靠性.

4 结论

本文针对山区复杂地形条件下的土石界面预测问题，提出了一种基于地形高程与界面深度的插值策略，并将该策略用于IDW、RBF核回归、GPR、MLP和KAN五种预测方法，通过RMSE、nRMSE和R²三种评估指标系统评估了该策略的应用效果.采用山区机场工程案例验证了所提策略，主要结论如下：

（1）与传统基于界面高程的插值策略相比，本文提出的基于地形高程与界面深度的插值策略先预测土石界面深度，再由地表高程减去预测深度估计土石界面高程，能够降低山区地形剧烈起伏对预测结果影响.工程案例结果表明，所提策略提升了预测准确性，预测结果准确地刻画了土石界面地空间起伏特征，为机场选址、土石方工程量估算及施工方案优化提供了技术支撑.

（2）本研究表明在小样本条件下，核方法（如IDW、RBF核回归和GPR）较神经网络方法（如MLP和KAN）的表现更为稳定；随着样本数量增加，两类方法的误差均显著降低并趋于一致.所提插值策略通过引入地表高程作为先验约束，可在小样本条件下一定程度缓解神经网络方法的欠拟合问题，因此具有更广的适用性.

（3）本文采用的nRMSE将RMSE相对于研究区最大地形高差进行归一化，提升了在不同地形起伏条件下误差指标的可比性，尤其适用于最大高差差异显著的场地.在对比不同场地（或研究区域）上的预测结果时，本研究建议联合使用多种评估指标（如nRMSE、RMSE和R²），以确保评价结果的客观性.

参考文献

[1] Bai, J., Wang, S., Xu, Q., et al., 2024. Intelligent Regional Subsurface Prediction Based on Limited Borehole Data and Interpretability Stacking Technique of Ensemble Learning. Bulletin of Engineering Geology and the Environment, 83(7): 272. https://doi.org/10.1007/s10064-024-03758-y

[2] Chen, C. F., Zhao, N., Yue, T. X., et al., 2015. A Generalization of Inverse Distance Weighting Method via Kernel Regression and Its Application to Surface Modeling. Arabian Journal of Geosciences, 8(9): 6623-6633. https://doi.org/10.1007/s12517-014-1717-z

[3] Cordonnier, G., Cani, M. P., Benes, B., et al., 2018. Sculpting Mountains: Interactive Terrain Modeling Based on Subsurface Geology. IEEE Transactions on Visualization and Computer Graphics, 24(5): 1756-1769. https://doi.org/10.1109/TVCG.2017.2689022

[4] Deng, Z. P., Pan, M., Niu, J. T., et al., 2023. Spatial Prediction of Rockhead Profile Using the Gaussian Process Regression Method. Canadian Geotechnical Journal, 60(12): 1849-1860. https://doi.org/10.1139/cgj-2022-0372

[5] Du, X., Fan, T. E., Dong, J. H., et al., 2020. Characterization of Thin Sand Reservoirs Based on a Multi-Layer Perceptron Deep Neural Network. Oil Geophysical Prospecting, 55(6): 1178-1187, 1159 (in Chinese with English abstract).

[6] Fu, J. M., Hu, M. S., Fang, F., et al., 2024. Complex Orebody 3D Modeling Using Radial Basis Function Surface Incorporating Stacking Integration Strategy. Earth Science, 49(3): 1165-1176 (in Chinese with English abstract).

[7] Huang, D. Z., Zuo, R. G., Wang, J., et al., 2025. Spatially Constrained Variational Autoencoder for Geochemical Data Denoising and Uncertainty Quantification. Journal of Earth Science, 36(5): 2317-2336. https://doi.org/10.1007/s12583-025-0180-y

[8] Li, X. Y., Zhang, L. M., Li, J. H., 2016. Using Conditioned Random Field to Characterize the Variability of Geologic Profiles. Journal of Geotechnical and Geoenvironmental Engineering, 142(4): 04015096. https://doi.org/10.1061/(asce)gt.1943-5606.0001428

[9] Li, Z. L., Wu, C. L., Zhang, X. L., et al., 2015. Uncertainty Assessment for IDW Ore Grade Estimates. Earth Science, 40(11): 1796-1801 (in Chinese with English abstract).

[10] Liu, Z. M., Wang, Y. X., Vaidya, S., et al., 2024. KAN: Kolmogorov-Arnold Networks. arXiv:

[11] Murphy, K. P., 2012. Machine Learning (4th Printing). The MIT Press, Cambridge, Mass.

[12] Murphy, K. P., 2022. Probabilistic Machine Learning. The MIT Press, Cambridge, Mass.

[13] Pan, K., Xie, C. Q., Cheng, R. Y., et al., 2017. Engineering Example of Comprehensive Prospecting of Underground Karst in Mountain Airport Area. Site Investigation Science and Technology, (1): 28-32 (in Chinese with English abstract).

[14] Qi, X. H., Pan, X. H., Chiam, K., et al., 2020. Comparative Spatial Predictions of the Locations of Soil-Rock Interface. Engineering Geology, 272: 105651. https://doi.org/10.1016/j.enggeo.2020.105651

[15] Qi, X. H., Wang, H., Chu, J., et al., 2022. Two-Dimensional Prediction of the Interface of Geological Formations: A Comparative Study. Tunnelling and Underground Space Technology, 121: 104329. https://doi.org/10.1016/j.tust.2021.104329

[16] Qi, X. H., Wang, H., Pan, X. H., et al., 2021. Prediction of Interfaces of Geological Formations Using the Multivariate Adaptive Regression Spline Method. Underground Space, 6(3): 252-266. https://doi.org/10.1016/j.undsp.2020.02.006

[17] Qiu, Z. L., Wu, J. D., Wan, P., et al., 2025. Experimental Study on the Influence of Soil-Rock Ratio on the Dynamic Compaction Reinforcement Effect of High Fill Gravel Soil Subgrade. Chinese Journal of Underground Space and Engineering, 21(1): 123-130 (in Chinese with English abstract).

[18] Samui, P., Kim, D., Viswanathan, R., 2015. Spatial Variability of Rock Depth Using Adaptive Neuro-Fuzzy Inference System (ANFIS) and Multivariate Adaptive Regression Spline (MARS). Environmental Earth Sciences, 73(8): 4265-4272. https://doi.org/10.1007/s12665-014-3711-x

[19] Somvanshi, S., Javed, S. A., Islam, M. M., et al., 2026. A Survey on Kolmogorov-Arnold Network. ACM Computing Surveys, 58(2): 1-35. https://doi.org/10.1145/3743128

[20] Su, Z. Y., Li, D. Q., Wang, S., et al., 2025. Large Deformation Analysis of 3D Soil-Rock Mixture Slopes Using SPH-DEM Method. Earth Science (in Chinese with English abstract) (in Press).

[21] Wang, J., 2021. An Intuitive Tutorial to Gaussian Processes Regression. arXiv:2009.10862

[22] Wang, Q., 2011. A Study on Quantitative Research Method about Soil-Rock Ratio in Red-Beds Airports’ Excavation Area in Sichuan (Dissertation). Chengdu University of Technology, Chengdu (in Chinese with English abstract).

[23] Xie, C. Q., Li, Z. Z., Pan, K., 2015. Study on Method Determining Soil-Rock Ratio in Large-Area Excavation Works in Southwest Karst Mountains. Subgrade Engineering, (2): 1-5, 12 (in Chinese with English abstract).

[24] Xie, C. Q., Rong, S. J., Wang, W., 2013. Study on Soil-Rock Ratio of the Projects in Red Bed Area. Subgrade Engineering, (6): 83-88 (in Chinese with English abstract).

[25] Zhang, H., Gui, L., Wang, T. F., et al., 2024. Prediction of Quaternary Cover Thickness and 3D Geological Modeling Based on BP Neural Network. Earth Science, 49(2): 550-559 (in Chinese with English abstract).

[26] Zhang, L. M., Dasaka, S. M., 2010. Uncertainties in Geologic Profiles versus Variability in Pile Founding Depth. Journal of Geotechnical and Geoenvironmental Engineering, 136(11): 1475-1488. https://doi.org/10.1061/(asce)gt.1943-5606.0000364

[27] Zhang, S. S., Li, Q. C., Li, H., et al., 2025. Intelligent Glacial Lake Identification in Complex Plateau Terrain Regions Using Multi-Source Remote Sensing Data and Mask R-CNN Deep Learning Model. Earth Science, 50(8): 3132-3143 (in Chinese with English abstract).

基金资助

国家自然科学基金项目(52278368)

四川省自然科学基金项目(24NSFSC2017)