地质对象由多个不同的地质单元构成,且具有空间连续性,因此,地质体模型是一个连续实体网络. 如图1,假设将地质结构曲面S定义为由欧氏三维空间上的曲线多边形面单元或连续小面片的集合,则该曲面无需闭合或相连.

设定离散拓扑模型A(Ω,N),Ω=Ω(S)是描述曲面S的所有点的集合,N=│Ω│是Ω集合中点的数量. 在此基础上,设定离散模型M(A(Ω,N),φ,λ), 其中,φ(α)={φ^x(α),φ^y(α),φ^z(α),…,φ^v(α),…}, α∈Ω, λ={λ₁,λ₂,……}, {λ_i}∑_α∈Ω∑_vA^v_i(α)φ^v(α)=b_i, 系数A和b可根据各种不同的约束条件确定.

图1 离散拓扑地质曲面示意图
Fig.1 Geological surface model by DSI

本研究采用离散光滑插值技术来进行高精度三维曲面的构建. 设φ(k)是定义在所有节点k∈Ω上的m个属性函数集合,φ(k)=(φ₁(k),φ₂(k),…,φ_m(k)),(φ₁(k),φ₂(k),φ₃(k))表示k点的几何坐标, φ₄(k)表示k点上的曲面法向矢量, φ_i(k)(4<i≤m)则是k点定义在曲面S上的标量或张量属性.

任意曲面S均可以由一个M×N的矩阵来表达其空间位置和属性分布,用一个m维行矩阵φ(k)=[φ₁(k),φ₂(k),…,φ_m(k)] 表达Ω上某节点k∈Ω的所有属性,用一个N维列矩阵φ_i=[φ_i(1),φ_i(2),…,φ_i(N)]^T表达定义在Ω集合中所有点的某一个属性i(1≤i≤m). 对于某一个定义在Ω集合中所有点的函数φ_i, 部分节点的函数值已知,其余节点函数值是未知的,需基于离散光滑插值方法来求解确定.

定义 Y为节点y∈Ω的集合, φ_i(y)已知; U为节点u∈Ω的集合, φ_i(u)未知; U=Ω-Y, φ_i =[φ_i(1),φ_i(2),…,φ_i(N)]^T=[φ_i(Y),φ_i(U)]^T, 因为网格节点的编号方式与排序均不影响拟合插值,故可对φ矩阵中元素的排列顺序进行交换,从而将φ矩阵表示成任意2个已知子矩阵与未知子矩阵的和,形成式(1)的简洁数学表达.

φ=[φ_I

φ_L]{φ_I=[φ(i₁)



φ(i_n)], i_α∈I

φ_L=[φ(l₁)



φ(l_m)], l_β∈L(1)

在进行三维建模时,不同可靠度约束条件的处理方法不同. 必须精确表达测量、钻孔和平硐记录等原始的空间位置及产状数据,对此类型数据处理时,需增加曲面描述的节点数量(即增加Ω空间容量). 基于地质对象的复杂性,在对地质人员的地质解释、物探数据解释等模糊性数据进行处理时,可对曲面S建立一系列不同的约束条件,并根据约束条件、地质可靠性的差异设置不同的置信因子,从而达到控制约束条件对曲面S影响的目的.

因此,可以用一个N维列矩阵来表达定义在Ω集合中所有点的某一个属性i(1≤i≤m), 若属性i有n个地质解释约束条件t(t=1,2,…,n), 则可以用一个N维行矩阵和一个常数来定义该约束条件Atφ_i=Bt, 其中, A_t=[A_t(1), A_t(2),…,A_t(N)], B_t为第t项约束条件对应常数.

约束条件对Ω空间容量无影响,故属性函数不完全满足约束关系,即|A_t(k)φ_i(k)-B_t|²>0. 需要根据每个约束条件的可靠性设置一个正数置信因子ω²_t, 并考虑所有的约束条件

ρ(φ_i(k))=∑_tω²_t|A_t(k)φ_i(k)-B_t|²(2)

在进行拟合插值时,ρ(φ_i)越小,越能准确地吻合约束条件; 可靠度大的约束条件,ω²_t的值也大,达到相同ρ(φ_i)时,要求|A_tφ_i-B_t|²越小越好,即严格遵守了约束条件.

接下来则需要确定出能使全局粗糙度与线性约束之和R^*(φ)最小的φ(k)函数,从而用φ(k)对曲面S进行估计. 设R(φ_i|k)是定义在节点k∈Ω上的属性函数φ_i的局部粗糙度, v_α(k)是定义在节点k邻域N(k)上所有点的权函数, u(k)是定义在节点k∈Ω上所有点的权函数,则有

R(φ_i|k)=|∑_α∈N(k)v_α(k)φ_i(α)|²=

|v_k(k)φ_i(k)+∑_{(α∈N(k))/(α≠k)}v_α(k)φ_i(α)|²(3)

v_k(k)=-∑_{(α∈N(k))/(α≠k)}v_α(k)(4)

R(φ_i|k)=(v_k(k))²×

|φ_i(k)-1/(∑_{(α∈N(k))/(α≠k)}v_α(k))[∑_{(α∈N(k))/(α≠k)}v_α(k)φ_i(α)/( / )]|²(5)

φ^-_i(k)=1/(∑_{(α∈N(k))/(α≠k)}v_α(k))[∑_{(α∈N(k))/(α≠k)}v_α(k)φ_i(α)/( / )](6)

式(6)中, φ^-_i(k)是定义在点k的邻域N(k)上,不包括k点的属性函数,按照该点的权函数v_α(k)(α∈N(k), α≠k)得到的平均函数值为

R(φ_i|k)=(v_k(k))²|φ_i(k)-φ^-_i(k)|²(7)

φ_i(k)与φ^-_i(k)差值的绝对值越小,节点k在曲面S上越光滑.

累积所有Ω空间上属性函数φ_i的局部粗糙度,即可求得全局粗糙度

R(φ_i(k)=∑^Nu(k)R(φ_i|k)(8)

DSI的目标是求解φ,使全局粗糙度R^～(φ)最小. 可将R(φ)近似看作R^～(φ), 即 R(φ)=R^～(φ), 定义ρ(φ|λ_i)=|∑_α∈ΩA_i(α)φ(α)-b_i|², 则加上所有约束条件的全局粗糙度R^*(φ)=∑_α∈Ωμ(k)R(φ|k)+∑_t(-overω)_iρ(φ|λ_i),要使R^*(φ)达到最优值(最小值),则要求

(R^*(φ))/(φ(1))=…=(R^*(φ))/(φ(α))=…=

(R^*(φ))/(φ(M))=0(9)

进而推导出一个半正定的稀疏线性方程组,采用LDL^T法或预条件共轭梯度法可解出方程组的未知量φ(k)={φ^x(k),φ^y(k),φ^z(k),…}. 求得R^*(φ_i(k))=Cφ²_i(α)+Dφ_i(α)+E, 系数C、 D和E不依赖于φ_i(α),当φ_i(α)=-D/(2C)时,可得最小度R^*(φ_i), 说明离散光滑插值方法是收敛的.

DSI技术的要点是确定出能够使R^*(φ)达到最小值的φ(k)函数,在理解该方法数学原理的基础上,可以进一步对其拟合精度进行实验验证.

为了检验所提出的模型,本研究用某水电站的测量、钻孔和物探等实际工程数据进行实验研究. 该水电站是红水河上一座梯级电站,1期装机容量为1 210 MW,2期扩建工程装机总量为2×300 MW,在1期厂坝工程的右岸山体内布置2期厂房,包括主厂房洞、主变室洞、引水洞、尾水主洞、尾水支洞、尾水闸门廊道、交通洞、厂房运输洞、尾水闸门运输洞及施工洞、排风洞和防渗帷幕灌浆洞等. 本实验选取了该电站1期和2期扩建工程的全部钻孔数据作为试验对象,含钻孔数量403个,按1:500地形测量点数据,以及全部钻孔柱状图信息.

基于TIN网来构建三维地表模型是多数GIS软件选择的方法,但因为TIN构模通常是基于原始数据采样点,构建的曲面不进行插值拟合,所以构建的曲面主要取决于原始采样点的质量. 在对TIN进行改进基础上形成的基于约束的DSI算法,不仅能够在TIN网模型构建时实现拟合插值,还能通过模糊信息约束、节点约束的设置等方法来达到曲面形状变形控制的目的. 离散光滑插值方法能构建任意复杂的曲面结构,有利于大倾角、非层状等结构形式的面模型构建,并具有很强的适应能力.

在三维地质曲面建模过程中,对于不同的地质对象,如捕掳体、逆掩断层和溶洞等,需要提供不同的约束方式,以达到对应的拟合目的,从而为构造出这些不同特征的地质体对象提供可编辑、可编程的数学逻辑. 而要达到能够处理不同可靠度约束条件的目的,需要提供不同的权函数,且在进行权函数的选择时要能够完全自由. 常用的权函数为

1)v_α(k)调和权

v_α(k)=1, α∈N(k)且α≠k(10)

v_k(k)=-∑_{(α∈N(k))/(α≠k)}v_α(k)(11)

2)v_α(k)高斯权

当属性函数φ_i表示定义在节点k的3×3张量 D(如应变)时,定义λ_αk为点α指向点k的距离矢量 λ_αk=[(x_k-x_α),(y_k-y_α),(z_k-z_α)], d_αk=λ_αkDλ^T_αk, v_α(k)=exp(-d_αk), α∈N(k)且α≠k, v_k(k)=-∑_{(α∈N(k))/(α≠k)}v_α(k), u(k)用于调整插值的局部光滑程度,一般取u(k)=1.

本实验对象的厂区硐室主要是辉绿岩,硐室上游侧主要为Ⅱ类围岩,硐室下游侧主要为Ⅲ类围岩,以陡立断层F₄₈断层为分界,主要有2组很发育的陡倾角节理形成不利地质因素. 另有产状为80°/NW∠20°～60°的逆掩断层F₄₄,其影响带范围为10.5～11.5 m. 主厂房硐室和引水隧道等均位于大坝的右岸山体内,且处于F₄₈断层的下盘,岩体的完整性较好. 断层F₄₄和F₄₈之间的尾水洞则裂隙发育、岩体破碎,从而导致成洞的地质条件差,且在施工期间易受地下水的影响.

在运用DSI技术进行高精度三维曲面的构建过程中,可以让多个约束条件一起生效,且各个约束条件之间,能够通过置信因子进行加权,可针对每个指定约束条件设置不同的权重,从而构建出满足多个不同约束条件的高精度的三维地质曲面模型. 基于DSI技术对地质对象几何、物理等特性进行模拟表达过程中,可把已知节点位置信息及相关地质特征信息等转化为约束条件,并纳入地质曲面建模的全过程. DSI技术还能进一步实现模糊化的向量约束、模糊化的控制节点等. 为了验证该方法在进行地质曲面三维建模时的效果,本研究基于实验数据分别模拟了地质对象中所包含的2个捕掳体的实体几何结构,分别如图2和图3.

图2 某工程中捕掳体1的几何结构
Fig.2 The geometry of the captive exile in igneous rock body 1

图3 某工程中捕掳体2的几何结构
Fig.3 The geometry of the captive exile in igneous rock body 2

在人工参与的地质解释过程中,对捕掳体施加了不同置信度的约束条件,在进行离散光滑插值时,置信度较高的约束条件将对应产生更严格的变形要求. 高精度地质曲面建模方法可以对网格模型进行自由选择和自动调整,还可以进行实时交互操作、对不确定数据进行过滤处理等,非常便于三维可视化和模型的修改更新.

本实验验证过程中,相关参数的取值为η₁=2/3, η₂=η₃=1/6, R₀=0.03. 为对预测结果进行误差评定,引入最为常用的3个精度评价指标:平均绝对误差(mean absolute error,MAE)、平均绝对相对误差(mean absolute percentage error, MAPE)和均方根误差(root-mean-square error, RMSE).为对置信区间进行评价,引用了预测区间平均覆盖概率(mean prediction interval coverage probability, MPICP). MPICP是指在一定置信度下,实际观测值包含在相对应的置信区间内的平均概率. 在进行预测精度分析的过程中,直接用钻孔的实际观测值作为真值.各误差指标计算如式(12)至式(16),其中,平均绝对误差为

MAE=1/n∑ⁿ_t=1|Y_t-X^⌒_t|(12)

平均绝对相对误差为

MAPE=1/n∑ⁿ_t=1|(Y_t-X^⌒_t)/(Y_t)|(13)

均方根误差为

RMSE=((∑ⁿ_t=1(Y_t-X^⌒_t)²)/n)^1/2(14)

区间平均覆盖率为

MPICP=1/n∑ⁿ_j=11/m∑^m_i=1c_i(15)

残差为

Z_res=Z_dat-Z_grd(16)

这里, Y_t为钻孔的实际测量值; X_t为钻孔之间插值拟合的预测值.

设预测区间为[L_i, U_i], 如果实际钻孔值x_i∈[L_i, U_i], 则c_i=1, 否则c_i=0.

一般情况下,置信度增加,置信区间的长度增加. 置信度越高,提供的信息越可靠,而置信区间越长,可能范围就越大,因此需要在两者间权衡. 通常认为当置信度达到80%以上时,比较可靠. 因此,本研究以置信度为80%时的预测区间来表现分布的预测结果.

本研究选择了10个钻孔数据来检验地质拟合插值的可靠性. 利用原始钻孔数据的某一类地层进行插值,运用式(16),计算原始钻孔数据与插值拟合结果之间的残差,并进行比较,见表1.

表1 几种插值拟合方法的残差评定
Table 1 Residuals assessment of several interpolation method

对10组被选择的钻孔资料进行不同的插值拟合效果做评价,发现DSI离散光滑插值残差比MS、Kringing和MinC等拟合方法的残差都要小. 基于DSI技术的高精度地质曲面建模方法既考虑了连续性地质曲面的插值,还能对三角网格法等方法所不能够插值的不连续地质曲面进行构建,且插值拟合后的光滑度效果更好. 同时,还能够将地质专家的经验解释纳入建模过程,在基本模型构建完成后,通过修改相关的权系数,即可得到对应的插值效果,从而极大地方便了模型的修改、更新和完善. 更关键的是,由于对原始采样点进行了硬约束,基于DSI的曲面建模方法可以保证钻孔等原始数据的精确,后期插值过程中产生的点对已有的采样数据无影响,从而有利于对整体模型的精度控制.