离散光滑插值方法是由Mallet教授[15-16]提出的,其思想是基于对目标体的离散化,类似有限元法,用一系列具有物体的几何和物理特性的相互连结的节点来模拟地质体. 已知节点和地质学的典型信息被转化为线性约束,被引入到模型生成的全过程.利用此技术来构建地质曲面,对采样数据少的区域,可以添加一系列的约束条件,控制曲面的插值过程,并自动调整插值点,以满足约束条件,使建立的曲面能较好地拟合采样点. 同时,还可考虑到地质变量本身的特性,因此对地质变量的统计分析能达到很好的效果[17-18].
本研究采用DSI技术来构造三维地质体网格模型.设离散拓扑模型为A(Ω, N),Ω=Ω(S)是描述曲面S的所有点的集合,N=|Ω|是Ω集合中点的数量.在此基础上,设离散模型为M(A(Ω, N),φ(α),λ), 其中,φ(α)={φx(α),φy(α), φz(α),…,φv(α),…},α∈Ω,λ={λ1, λ2,…},{λi} ∑α∈Ω∑vAvi(α)φv(α)=bi, 可根据各种不同的约束条件确定系数A和b的值.
设φ(k)是定义在节点k∈Ω上的m个属性函数的集合,φ(k)=(φ1(k), φ2(k), …, φm(k)),(φ1(k), φ2(k), φ3(k))代表k点的几何坐标,φ4(k)用于表达在k点上曲面的法向矢量,φi(k)(4﹤i≤m)表示在k点用标量或张量定义在曲面S上的属性.
对于一个给定的曲面S, 其空间位置和曲面上的属性分布,可以通过一个m×n的矩阵来表达. Ω集合中节点k∈Ω的所有属性,可以用一个m维行矩阵表达, 即 φ(k)=[φ1(k), φ2(k), …, φm(k)]. 对于定义在Ω集合中所有点的某一个属性i(1≤i≤m), 可以用一个n维列矩阵表达, 即 φi =[φi(1), φi(2), …, φi(n)]T, 对于某一个定义在Ω集合中所有点的函数,有部分节点的函数值已知,其余节点的函数值未知,需通过离散光滑插值的方法来确定.
定义Y为节点y∈Ω的集合, φi(y)已知; U为节点u∈Ω的集合, φi(u)未知; U=Ω-Y, φi =[φi(1),φi(2),…,φi(n)]T=[φi(Y),φi(U)]T, 由于对网格节点的编号顺序和方式并不影响插值问题的解决,可以假设交换φ矩阵中元素的排列顺序,即把φ矩阵分成任意两个子矩阵 φI和 φL,其中, φI为已知矩阵, φL为未知矩阵.从而使符号更简洁.
φ=[φI
φL]{φI=[φ(i1)
φ(i2)
φ(in)], iα∈I
φL=[φ(l1)
φ(l2)
φ(lm)], lβ∈L(1)
在地质建模的过程中,需要处理不同可靠度的约束条件.对于诸如钻孔平硐记录这样的空间位置数据和产状数据,必须严格遵守.对于该类型数据,一般都需要通过增加对曲面描述的节点数量,即增加Ω空间的容量来实现.对于另外一些约束条件,如地质人员的地质解释和物探数据的解释等,这些解释会形成对曲面S的一系列不同的约束条件.由于地质体的高度复杂性,这些通过钻探或物探得到的解释条件会存在相互矛盾的状况.根据对获得这些约束条件的地质可靠性的判断,对每个约束条件设置一个置信因子来控制约束条件对曲面S的影响.
如上所述,对于定义在Ω集合中所有点的某一个属性i(1≤i≤m), 可以用一个n维列矩阵 φi=[φi(1), φi(2), …, φi(n)]T来表达.对于该属性,如果存在一系列的地质解释约束条件t(t=1, 2, 3, …), 那么对于这些约束条件对曲面S的影响,可以用一个n维的行矩阵和一个常数来定义该约束条件Atφi=Bt, 其中At=[At(1), At(2), …, At(n)], Bt为第t项约束条件对应常数.
由于约束条件不增加Ω空间的容量,属性函数φi不能完全满足约束关系,即|At(k)φi(k)-Bt|2>0. 对于每个约束条件,根据其可靠性设置一个正数置信因子ω2t并考虑所有的约束条件ρ(φi(k))=∑tω2t|At(k)φi(k)-Bt|2.
在插值过程中, ρ(φi)越小,则对约束条件的吻合就越准确.同时可以看出,对于可靠度比较大的约束条件ω2t的值也较大,达到相同ρ(φi)时要求|Atφi-Bt|2越小,即对约束条件的遵守就越严格.
接下来需要解决的问题,是要确定出能使全局粗糙度与线性约束的和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(α)|2=
|vk(k)φi(k)+∑(α∈N(k))/(α≠k)vα(k)φi(α)|2(2)
vk(k)=-∑(α∈N(k))/(α≠k)vα(k)(3)
则有
R(φi|k)=(vk(k))2|φi(k)-1/(∑(α∈N(k))/(α≠k)vα(k))×
( )[∑(α∈N(k))/(α≠k)vα(k)φi(α)]|2(4)
φ^-i(k)=1/(∑(α∈N(k))/(α≠k)vα(k))[∑(α∈N(k))/(α≠k)vα(k)φi(α)](5)
可以看出,φ^-i(k)就是定义在点k的邻域N(k)上不包括k点的属性函数φi按照该点的权函数vα(k)(α∈N(k), α≠k)得到的平均函数值,则式(4)为
R(φi|k)=(vk(k))2|φi(k)-φ^-i(k)|2(6)
当定义在节点k的属性函数φi(k)和定义在k点邻域的平均函数φ^-i(k)的差值越小,则曲面S在节点k上就越光滑.
将定义在Ω空间上所有节点的属性函数φi的局部粗糙度累积,可得到全局粗糙度
R(φi(k))=∑Nu(k)R(φi|k)(7)
离散光滑插值方法的目标就是求解φ, 使得全局粗糙度R～(φ)最小. 可以将R(φ)近似地看作为全局粗糙度R～(φ), 即R(φ)=R～(φ), 定义局部约束度为w^-i, 可得局部粗糙度函数为ρ(φi|λi)=|∑α∈ΩAi(α)φi(α)-bi|2, 则加上所有约束条件的全局粗糙度为R*(φ)=∑α∈Ωu(k)×R(φi|k)+∑i(-overω)iρ(φi|λi),要使R*(φ)达到最优值(最小值),则要求
(R*(φ))/(φ(1))=…=(R*(φ))/(φ(α))=…=(R*(φ))/(φ(m))(8)
由此可推导得出一个半正定的稀疏线性方程组,采用预条件共轭梯度法或LDLT法,求解方程组的未知量φ(k)={φx(k), φy(k),φz(k)}. 得R*(φi(k))=Cφ2i(α)+Dφi(α)+E, 其中,C、 D和E为系数,不依赖于φi(α),当φi(α)=-D/(2C)时,可得最小粗糙度R*(φi), 从而说明DSI法收敛.
