工程科学学报，第37卷，第5期：556560,2015年5月 Chinese Journal of Engineering,Vol.37,No.5:556-560,May 2015 DOI:10.13374/j.issn2095-9389.2015.05.003:http://journals.ustb.edu.cn 四维空间分形插值算法在品位估算中的应用 胡乃联，陈建均品，李国清，杨桦 北京科技大学土木与环境工程学院，北京100083 ☒通信作者，E-mail:cj2012@sina.com 摘要克里金法是一种应用广泛的低通滤波性插值方法，但其无法重建原始信息中的高频、低频和局部信息.分形插值算 法可利用自相似性，在保留原始信息的基础上，克服克里金插值中低通滤波的局限性，从而提高插值的准确性.本文在传统 分形插值算法的基础上，结合地质空间信息，提出了适用于矿床品位估算的四维空间分形插值算法，并将其应用于钼矿的品 位估算.结果表明：在该钼矿的品位估算中，四维空间分形插值算法明显优于克里金法 关键词探矿：品位估算：插值算法：分形；均方误差 分类号TD166 Application of a four-dimensional space fractal interpolation algorithm in grade estimation HU Nai-lian,CHEN Jian-jun,LI Guo-qing,YANG Hua School of Civil and Environmental Engineering,University of Science and Technology Beijing,Beijing 100083,China Corresponding author,E-mail:cjj2012@sina.com ABSTRACT Kriging interpolation is a widely used low-pass filter interpolation method,but it cannot reconstruct the high-frequency, low-frequency and partial information of original information.Fractal interpolation using self-similarity,which can retain original infor- mation,overcomes the limitations of Kriging interpolation low-pass filters,thereby improving the interpolation accuracy.On the basis of the traditional fractal interpolation algorithm and in combination with geological spatial information,this paper introduces a four- dimensional space fractal interpolation algorithm suitable for ore grade estimation.The interpolation algorithm is applied to molybde- num ore grade estimation and then compared with the Kriging interpolation algorithm.The results show that the interpolation algorithm is superior to the Kriging interpolation algorithm. KEY WORDS mine exploration:grade estimation:interpolation algorithms:fractals:mean square error 由于勘探采样取得的钻孔品位数据有限，无法直 阵的方差中占的比重很小，这些信息很容易被压低 接判断矿床的开采价值，因此寻求一种有效的品位插 或消除 值方法对矿山项目投资十分重要.近年来，多维概率 地质数据存在复杂性和不规则性：但大量研究 分布模型四、支持向量机模型-、克里金插值法等品 表明，地质数据存在广泛的自相似性可，可以考虑将 位插值技术已应用于品位插值处理中，其中克里金插 分形插值应用到地质品位估算的研究中.针对克里 值法应用最为广泛. 金插值法存在的低通滤波性，本文提出基于分形理 克里金插值法是一种最佳线性无偏估算的方 论的四维空间分形插值算法，并将该方法应用于某 法，但它的求解过程决定了它是一个空间滑动平均 钼矿的钻孔品位数据插值过程中，同时与克里金插 或低通滤波过程，高频、局部与弱信息在相关系数矩 值法进行了对比. 收稿日期：2013-12-24 基金项目：国家自然科学基金资助项目(51104010)：中央高校基本科研业务费专项资金资助项目(FRF-SD-12001A)
工程科学学报，第 37 卷，第 5 期: 556--560，2015 年 5 月 Chinese Journal of Engineering，Vol． 37，No． 5: 556--560，May 2015 DOI: 10． 13374 /j． issn2095--9389． 2015． 05． 003; http: / /journals． ustb． edu． cn 四维空间分形插值算法在品位估算中的应用 胡乃联，陈建均，李国清，杨 桦 北京科技大学土木与环境工程学院，北京 100083 通信作者，E-mail: cjj2012@ sina． com 摘 要 克里金法是一种应用广泛的低通滤波性插值方法，但其无法重建原始信息中的高频、低频和局部信息． 分形插值算 法可利用自相似性，在保留原始信息的基础上，克服克里金插值中低通滤波的局限性，从而提高插值的准确性． 本文在传统 分形插值算法的基础上，结合地质空间信息，提出了适用于矿床品位估算的四维空间分形插值算法，并将其应用于钼矿的品 位估算． 结果表明: 在该钼矿的品位估算中，四维空间分形插值算法明显优于克里金法． 关键词 探矿; 品位估算; 插值算法; 分形; 均方误差 分类号 TD166 Application of a four-dimensional space fractal interpolation algorithm in grade estimation HU Nai-lian，CHEN Jian-jun ，LI Guo-qing，YANG Hua School of Civil and Environmental Engineering，University of Science and Technology Beijing，Beijing 100083，China Corresponding author，E-mail: cjj2012@ sina． com ABSTＲACT Kriging interpolation is a widely used low-pass filter interpolation method，but it cannot reconstruct the high-frequency， low-frequency and partial information of original information． Fractal interpolation using self-similarity，which can retain original information，overcomes the limitations of Kriging interpolation low-pass filters，thereby improving the interpolation accuracy． On the basis of the traditional fractal interpolation algorithm and in combination with geological spatial information，this paper introduces a fourdimensional space fractal interpolation algorithm suitable for ore grade estimation． The interpolation algorithm is applied to molybdenum ore grade estimation and then compared with the Kriging interpolation algorithm． The results show that the interpolation algorithm is superior to the Kriging interpolation algorithm． KEY WOＲDS mine exploration; grade estimation; interpolation algorithms; fractals; mean square error 收稿日期: 2013--12--24 基金项目: 国家自然科学基金资助项目( 51104010) ; 中央高校基本科研业务费专项资金资助项目( FＲF--SD--12--001A) 由于勘探采样取得的钻孔品位数据有限，无法直 接判断矿床的开采价值，因此寻求一种有效的品位插 值方法对矿山项目投资十分重要． 近年来，多维概率 分布模型［1］、支持向量机模型［2--3］、克里金插值法等品 位插值技术已应用于品位插值处理中，其中克里金插 值法应用最为广泛． 克里金 插 值 法 是 一 种 最 佳 线 性 无 偏 估 算 的 方 法，但它的求解过程决定了它是一个空间滑动平均 或低通滤波过程，高频、局部与弱信息在相关系数矩 阵的方差中占的比重很小，这些信息很容易被压低 或消除［4］． 地质数据存在复杂性和不规则性; 但 大 量 研 究 表明，地质数据存在广泛的自相似性［5］，可以考虑将 分形插值应用到地质品位估算的研究中． 针对克里 金插值法存在的低通滤波性，本文提出基于分形理 论的四维空间分形插值算法，并将该方法应用于某 钼矿的钻孔品位数据插值过程中，同时与克里金插 值法进行了对比．
胡乃联等：四维空间分形插值算法在品位估算中的应用 557 1克里金插值法的低通滤波性 这四个方程是由约束条件给出的.可见，每个变 换心：的系数中存在一个自由变量：取d,为自由变量， 克里金插值法是一种最优无偏估量的储量计算方 并称之为垂直尺度因子.限定自由变量0≤d,<1.因 法，它的插值公式为阿 为变换，能够把每条与y平行的直线变换成另一条与 2(- y平行的直线.所以，选取d,为自由变量，就能指定由 入：Z(x) (1) 变换心，所生成的垂直尺度.令d,为任意一个限定的实 式中，入(i=1,2,…,n)为权重系数，Z*(x。)为估计 数，则地，其他的系数可以表示为 值，Z(x:)为已知值. Xi-xi-1 从克里金插值理论可以看出，克里金插值法是利 a= XN -x0 用空间自相关性，在方差极小意义下估算待估点的属 XxXi-1-xoxi 性值，其实质是一种加权平均方法.由于高频信号仅 e= XN -X0 有较低的空间相关性，对加权系数贡献较小，造成高频 9.=-_4-) (5) 信号损失较大，因此克里金插值是一个明显的低通滤 XN-Xo XN-X0 波过程.为此，本文提出了基于分形理论的四维空 l-xoyd,(xyo-xoyx） 间分形插值算法，试图解决克里金插值的低通滤波性. Xx-x0 2分形插值原理和方法 可以证明此定义求取迭代函数系统总用唯一的吸 引子，且该吸引子必定是某个连续函数的图形，并同时 分形是对事物的形状、形态、结构与组织的分解、 通过各个插值点，而这个连续函数就称为分形插值 分割与分析，它具有自相似性m。自相似性是指局部 函数四 与整体在形态、功能、信息等方面具有统计意义上的相 2.2四维空间分形插值法 似性圆.分形插值是一种构造分形曲线的方法，是由 分形插值基于地质数据的自相似性，具体表现为 Barnsley在迭代函数系统基础上提出来的.用分形插 矿床在空间上是不均匀分布的，但矿床整体品位分布 值可以得到相邻两插值点之间的局部变化特征，从而 规律与其局部品位分布规律相似，即矿床在空间上呈 使插值结果更加符合实际 分形丛集分布画.通过钻孔品位数据得到整体品位 2.1分形插值原理 分布规律，并利用其对矿床进行品位插值，能够有效地 对于给定的数据点{(x,):x-1<x,i=1,2, 反映局部品位变化规律，从而提高插值精度. …,N},定义二维欧氏空间中的仿射变换为心：R2→ 相对于克里金插值法，分形插值利用拼贴原理得 R,可构造迭代函数系统(FS):{R2;w,i=1,2,, 到迭代函数系统.拼贴原理的实质是在矿床空间内通 n},使得此迭代函数系统的吸引子等于插值函数∫(x) 过寻找合适的压缩仿射变换，使得钻孔品位数据拼贴 的图形.设迭代函数系统中每个函数都是仿射变换， 成一个矿床整体品位数据集合。在拼贴过程中能够很 其构造为 好地保留包含高频、局部与弱信息在内的原始信息，从 而克服克里金插值法的低通滤波性 -(9份 (2) 由于品位估算时需要利用三维空间信息，加上地 质品位就成为四维空间，而传统的分形插值最多仅能 这时有集合函数L(x)=ax+e,插值函数F,(x, 计算三维空间1-，未充分考虑地质品位的空间信 y）=cx+dy+∫，其中a、cd、e,和f为变换函数的 息，在一定程度上影响了分形插值的推广能力，因此本 系数. 文提出四维空间分形插值算法 该迭代函数系统的求取是通过拼贴完成的，满足 令I=a,B],J=[y,δ]，K=[e,],设区域D=I 如下条件： ×J×K={(x,y,z)la≤x≤B,y≤y≤6，e≤z≤}， ）-） i=1,2,…,n.(3) 以△、△和上为步长，将D剖分为网格： ,a=x0<x1<·<xy=B 又由L:（xo)=x-1,L:（xx)=x,F:(xo,o）= y=%<为<…<yM=8, (6) y,F(xx,y)=y:,可得方程组： e=0<1<…<21=g a;xo +ei=xi-1' (1)X坐标的插值公式，其中N为x的可能取值. aixy +ei=xi, (4) cixo+diyo+f=yi ,()=1+)(x- -,n∈{1,2，…，N}. cixx +dyx +f=y (7)
胡乃联等: 四维空间分形插值算法在品位估算中的应用 1 克里金插值法的低通滤波性 克里金插值法是一种最优无偏估量的储量计算方 法，它的插值公式为［6］ Z* ( x0 ) = ∑ n i = 1 λiZ( xi ) ． ( 1) 式中，λi ( i = 1，2，…，n) 为权重系数，Z* ( x0 ) 为估计 值，Z ( xi ) 为已知值． 从克里金插值理论可以看出，克里金插值法是利 用空间自相关性，在方差极小意义下估算待估点的属 性值，其实质是一种加权平均方法． 由于高频信号仅 有较低的空间相关性，对加权系数贡献较小，造成高频 信号损失较大，因此克里金插值是一个明显的低通滤 波过程［4］． 为此，本文提出了基于分形理论的四维空 间分形插值算法，试图解决克里金插值的低通滤波性． 2 分形插值原理和方法 分形是对事物的形状、形态、结构与组织的分解、 分割与分析，它具有自相似性［7］． 自相似性是指局部 与整体在形态、功能、信息等方面具有统计意义上的相 似性［8］． 分形插值是一种构造分形曲线的方法，是由 Barnsley 在迭代函数系统基础上提出来的． 用分形插 值可以得到相邻两插值点之间的局部变化特征，从而 使插值结果更加符合实际． 2. 1 分形插值原理 对于给定的数据点{ ( xi，yi ) ; xi － 1 ＜ xi，i = 1，2， …，N } ，定义二维欧氏空间中的仿射变换为 w: Ｒ2 → Ｒ2 ，可构造迭代函数系统( IFS) : { Ｒ2 ; wi，i = 1，2，…， n} ，使得此迭代函数系统的吸引子等于插值函数 f( x) 的图形． 设迭代函数系统中每个函数都是仿射变换， 其构造为 wi( ) x y = ai 0 ( c d )i ( ) x y + ei ( f )i ． ( 2) 这时有集合函数 Li ( x) = ai x + ei，插值函数 Fi ( x， y) = ci x + di y + fi，其中 ai、ci、di、ei和 fi为变换函数的 系数． 该迭代函数系统的求取是通过拼贴完成的，满足 如下条件: wi x0 ( y ) 0 = xi － 1 yi ( ) － 1 ，wi xn ( y ) n = xi ( y )i ，i = 1，2，…，n． ( 3) 又由 Li ( x0 ) = xi － 1，Li ( xN ) = xi，Fi ( x0，y0 ) = yi － 1，Fi ( xN，yN ) = yi，可得方程组: aix0 + ei = xi － 1， aixN + ei = xi， cix0 + diy0 + fi = yi － 1， cixN + diyN + fi = yi ． ( 4) 这四个方程是由约束条件给出的． 可见，每个变 换 wi的系数中存在一个自由变量; 取 di 为自由变量， 并称之为垂直尺度因子． 限定自由变量 0≤di ＜ 1． 因 为变换 wi能够把每条与 y 平行的直线变换成另一条与 y 平行的直线． 所以，选取 di为自由变量，就能指定由 变换 wi所生成的垂直尺度． 令 di为任意一个限定的实 数，则 wi其他的系数可以表示为 ai = xi － xi － 1 xN － x0 ， ei = xN xi － 1 － x0 xi xN － x0 ， ci = yi － yi － 1 xN － x0 － di ( yN － y0 ) xN － x0 ， fi = xN yi － 1 － x0 yi xN － x0 － di ( xN y0 － x0 yN ) xN － x0 ． ( 5) 可以证明此定义求取迭代函数系统总用唯一的吸 引子，且该吸引子必定是某个连续函数的图形，并同时 通过各个插值点，而这个连续函数就称为分形插值 函数［9］． 2. 2 四维空间分形插值法 分形插值基于地质数据的自相似性，具体表现为 矿床在空间上是不均匀分布的，但矿床整体品位分布 规律与其局部品位分布规律相似，即矿床在空间上呈 分形丛集分布［10］． 通过钻孔品位数据得到整体品位 分布规律，并利用其对矿床进行品位插值，能够有效地 反映局部品位变化规律，从而提高插值精度． 相对于克里金插值法，分形插值利用拼贴原理得 到迭代函数系统． 拼贴原理的实质是在矿床空间内通 过寻找合适的压缩仿射变换，使得钻孔品位数据拼贴 成一个矿床整体品位数据集合． 在拼贴过程中能够很 好地保留包含高频、局部与弱信息在内的原始信息，从 而克服克里金插值法的低通滤波性． 由于品位估算时需要利用三维空间信息，加上地 质品位就成为四维空间，而传统的分形插值最多仅能 计算三维空间［11--12］，未 充分考虑地质品位的空间信 息，在一定程度上影响了分形插值的推广能力，因此本 文提出四维空间分形插值算法． 令 I =［α，β］，J =［γ，δ］，K =［ε，ζ］，设区域 D = I × J × K = { ( x，y，z) | α≤x≤β，γ≤ y≤δ，ε≤z≤ζ} ， 以!x、!y 和!z 为步长，将 D 剖分为网格: α = x0 ＜ x1 ＜ … ＜ xN = β， γ = y0 ＜ y1 ＜ … ＜ yM = δ， ε = z0 ＜ z1 ＜ … ＜ zL = ζ { ． ( 6) ( 1) X 坐标的插值公式，其中 N 为 x 的可能取值． n ( x) = xn － 1 + ( xn － xn － 1 ) ( x － x0 ) xN － x0 ，n∈{ 1，2，…，N} ． ( 7) · 755 ·
·558· 工程科学学报，第37卷，第5期 (2)Y坐标的插值公式，其中M为y的可能取值. 误差(MSE)作为验证指标. pn(y）=ym-1+- ,（yn-ym-i）(y-yo) Yw-Yo .2-2? (12) me{1,2,…,M}. (8) n·s2 (3)Z坐标的插值公式，其中L为z的可能取值. 式中，Z为观测点估计值，Z,为观测点实际值，s2为观 -a=,le1,2,. 测点实际值的方差. 9,(z）=z1-1+- 2L-20 均方误差越小，预测值就越接近它们的真实值，插 (9) 值效果也就越好 (4)W品位的插值公式.令 本文以某钼矿为工程背景，基于MATLAB软件平 Fm（x,y,z,0）=ga,mx+hn.m,y+Pm2+ 台，对该钼矿468个钻孔钼品位样本数据进行研究，随 9my+「.mz+5.m2+【mXz+m.0, 机选取其中252个样本数据作为插值的初始值，其余 n∈{1,2，…，N},m∈{1,2，…，M,le{1,2,…,L} 216个样本数据用来和插值结果进行对照.以100m 由 水平为例，利用MATLAB软件作钼矿四维空间分形插 [Fmm1（x0y020000.0）=Z。-l,m-1.l-1y 值算法品位等值线图与钼矿克里金插值法品位等值线 F（xoy020o.l）=Z.-la- 图分别如图1和图2所示. F1（x0yn20,0o..0）=Z。-l,m1-1' 钼品位10 4.6585 40 F（xoy20o.）=Z-1' n640 4.6580 (10) 640 F(xo0o0）=Z-- 4.6575 490 90 4.6570 540590 Fa(xxyo0xo.）=Za- F.i（xyga0.0）=Zna-' 日4.6565 色4.6560 5907厂 F.m1（xxy220.Mi）=Zml 4.6555 得 4.6550 4.6545 -11=gn.m+hmo+pam+9n.moYo+ 590 4.6540 a.m.o+mmo+.0.0.0 465934815249252483524894524855248652487 a-1m18m+hmo+pa.m.L+9.m.oYo+ x/10°m a.m.o.mYozL+moYow.m.0.0.L 图1100m水平分形插值 -1.mgmhamyy+pa.m:o+9am.oy+ Fig.1 Fractal interpolation at the 100 m level 「，1X020+5mJ2g+Lnm.Xoy0+uam100.M.0y za-1.m=gn.m.o+ha.my+pa.m.+9n.m.y+ 钼品位/106 4.6585 a.m.myu.m 4.6580 640 a.m-1-1=gm.+hamYo+Pa.m1o+9m.o+ 4.6575 4.6570 「a,m,x20+snm,y020+lm.心xy020+La,m10,0.0’ 4.6565 a.m1.mha:m:o+pa.m.+9n.mo+ 兰4.6560 a.m:XNL+S.m.IYoZL+mYoL+uw.m.IWN.0.L 4.65 乙a,m1-l=gm.x+hm,yM+Pnm.10+9am式xyy+ 4.6550 4.6545 a.m.o5n.mYomyoun.m 4.6540 n.mgmha.my+pammy+ 4.6535 5.24815.24825.24835.24845.24855.24865.2487 T,m,X2L+5a,m,n2L+la,m,1心xyn2L+L。,m0X,M,上 x/105m (11) 图2100m水平克里金插值 Fig.2 Kriging interpolation at the 100 m level 3实例验证 为了评价插值方法的质量，本文将数据样本总体 对比图1和图2可以看出，四维空间分形插值算 分为两组，一组作为观测点（待估计点），另一组作为 法与克里金插值法的品位分布整体趋势相同.但在中 已知点，然后通过已知点对观测点进行预测，采用均方 间位置，分形插值中钼品位(590×106)分布呈现零星
工程科学学报，第 37 卷，第 5 期 ( 2) Y 坐标的插值公式，其中 M 为 y 的可能取值． φm ( y) = ym － 1 + ( ym － ym － 1 ) ( y － y0 ) yM － y0 ， m∈{ 1，2，…，M} ． ( 8) ( 3) Z 坐标的插值公式，其中 L 为 z 的可能取值． l ( z) = zl － 1 + ( zl － zl － 1 ) ( z － z0 ) zL － z0 ，l∈{ 1，2，…，L} ． ( 9) ( 4) W 品位的插值公式． 令 Fn，m，l ( x，y，z，w) = gn，m，lx + hn，m，ly + pn，m，lz + qn，m，lxy + rn，m，lxz + sn，m，lyz + tn，m，lxyz + un，m，lw， n∈{ 1，2，…，N} ，m∈{ 1，2，…，M} ，l∈{ 1，2，…，L} ． 由 Fn，m，l ( x0，y0，z0，w0，0，0 ) = Zn － 1，m － 1，l － 1， Fn，m，l ( x0，y0，zL，w0，0，L ) = Zn － 1，m － 1，l， Fn，m，l ( x0，yM，z0，w0，M，0 ) = Zn － 1，m，l － 1， Fn，m，l ( x0，yM，zL，w0，M，L ) = Zn － 1，m，l， Fn，m，l ( xN，y0，z0，wN，0，0 ) = Zn，m － 1，l － 1， Fn，m，l ( xN，y0，zL，wN，0，L ) = Zn，m － 1，l， Fn，m，l ( xN，yM，z0，wN，M，0 ) = Zn，m，l － 1， Fn，m，l ( xN，yM，zL，wN，M，L ) = Zn，m， l ( 10) 得 zn － 1，m － 1，l － 1 = gn，m，lx0 + hn，m，ly0 + pn，m，lz0 + qn，m，lx0 y0 + rn，m，lx0 z0 + sn，m，ly0 z0 + tn，m，lx0 y0 z0 + un，m，lw0，0，0， zn － 1，m － 1，l = gn，m，lx0 + hn，m，ly0 + pn，m，lzL + qn，m，lx0 y0 + rn，m，lx0 zL + sn，m，ly0 zL + tn，m，lx0 y0 zL + un，m，lw0，0，L， zn － 1，m，l － 1 = gn，m，lx0 + hn，m，lyM + pn，m，lz0 + qn，m，lx0 yM + rn，m，lx0 z0 + sn，m，lyM z0 + tn，m，lx0 yM z0 + un，m，lw0，M，0， zn － 1，m，l = gn，m，lx0 + hn，m，lyM + pn，m，lzL + qn，m，lx0 yM + rn，m，lx0 zL + sn，m，lyM zL + tn，m，lx0 yM zL + un，m，lw0，M，L， zn，m － 1，l － 1 = gn，m，lxN + hn，m，ly0 + pn，m，lz0 + qn，m，lxN y0 + rn，m，lxN z0 + sn，m，ly0 z0 + tn，m，lxN y0 z0 + un，m，lwN，0，0， zn，m － 1，l = gn，m，lxN + hn，m，ly0 + pn，m，lzL + qn，m，lxN y0 + rn，m，lxN zL + sn，m，ly0 zL + tn，m，lxN y0 zL + un，m，lwN，0，L， zn，m，l － 1 = gn，m，lxN + hn，m，lyM + pn，m，lz0 + qn，m，lxN yM + rn，m，lxN z0 + sn，m，lyM z0 + tn，m，lxN yM z0 + un，m，lwN，M，0， zn，m，l = gn，m，lxN + hn，m，lyM + pn，m，lzL + qn，m，lxN yM + rn，m，lxN zL + sn，m，lyM zL + tn，m，lxN yM zL + un，m，lwN，M，L ． ( 11) 3 实例验证 为了评价插值方法的质量，本文将数据样本总体 分为两组，一组作为观测点( 待估计点) ，另一组作为 已知点，然后通过已知点对观测点进行预测，采用均方 误差( MSE) ［13］作为验证指标． MSE = ∑ n k = 1 ( Z* k － Zk ) 2 n·s 2 ． ( 12) 式中，Z* k 为观测点估计值，Zk为观测点实际值，s 2 为观 测点实际值的方差． 均方误差越小，预测值就越接近它们的真实值，插 值效果也就越好． 本文以某钼矿为工程背景，基于 MATLAB 软件平 台，对该钼矿 468 个钻孔钼品位样本数据进行研究，随 机选取其中 252 个样本数据作为插值的初始值，其余 216 个样本数据用来和插值结果进行对照． 以 100 m 水平为例，利用 MATLAB 软件作钼矿四维空间分形插 值算法品位等值线图与钼矿克里金插值法品位等值线 图分别如图 1 和图 2 所示． 图 1 100 m 水平分形插值 Fig． 1 Fractal interpolation at the 100 m level 图 2 100 m 水平克里金插值 Fig． 2 Kriging interpolation at the 100 m level 对比图 1 和图 2 可以看出，四维空间分形插值算 法与克里金插值法的品位分布整体趋势相同． 但在中 间位置，分形插值中钼品位( 590 × 10 － 6 ) 分布呈现零星 · 855 ·
胡乃联等：四维空间分形插值算法在品位估算中的应用 ·559· 分布现象，而克里金插值则是表现为相对集中.这是 性而未能反映出品位分布的局部性 由于分形插值具有自相似性，能够体现相邻两插值点 用均方误差对这两种插值方法的插值精度进行对 之间的局部变化特征，而克里金插值则由于低通滤波 比.表1列出代表性观测点的插值误差对比 表1钼矿两种插值方法结果 Table 1 Results of two interpolation methods in molybdenum ore grade estimation y/m x/m m 钼品位观测值106 分形插值误差 克里金插值误差 465620 5248180 100 594.8263 0.011751521 0.012588814 465700 5248250 100 594.4836 0.011729090 0.013762319 465380 5248670 100 491.3848 -0.007394643 0.007019079 465460 5248670 100 557.1922 0.002370298 0.004807368 465540 5248670 100 576.2724 0.008065163 0.005567949 465700 5248110 200 594.733 -0.005177851 0.004436979 465780 5248110 200 595.4168 0.055252780 -0.059105252 465380 5248600 200 923.9639 0.026653377 0.019697220 465620 5248600 200 1016.716 -0.041018674 0.062061316 465540 5248670 200 995.4822 -0.010591899 0.034773408 465780 5248110 300 701.9488 -0.004437747 0.019388099 465620 5248320 300 1064.548 -0.002040820 0.030878965 465780 5248460 300 1010.303 0.035103251 -0.053138138 465540 5248530 300 1937.489 0.032239192 0.042725035 465540 5248670 300 957.8489 0.023094816 0.035682428 465540 5248110 400 587.7758 -0.000612992 -0.007780875 465780 5248110 400 660.2681 -0.004770689 0.001987888 465460 5248390 400 658.8874 -0.008764155 -0.018133457 465620 5248390 400 893.1652 0.001446660 0.032450752 465380 5248530 400 1024.34 0.006630153 0.011699876 均方误差 1.14098×10-8 2.20978×10-8 由表1可以看出，采用四维空间分形插值算法估 础上，提出了四维空间分形插值算法 算的均方误差比克里金插值法降低48.4%，说明在该 (2)该方法克服了传统分形插值算法只反映平面 钼矿品位估算中，四维空间分形插值算法明显优于克 信息的局限性和克里金插值法的低通滤波性，不仅能 里金插值法 反映相邻两插值点之间的局部变化特征，而且能充分 从插值原理上看：克里金插值法中任意两相邻插 体现地质品位变化规律复杂的特点，使插值结果更加 值点间的信息是通过光滑曲线连接的，从而掩盖了相 精确 邻两插值点间的局部变化特征，因而具有一定程度的 (3)以某钼矿为例，将该方法应用于品位估算过 光滑作用：而四维空间分形插值算法则是根据整体与 程中.与克里金插值法相比，应用四维空间分形插值 局部相似的原理，将插值数据点的变化特征映射到相 算法进行品位估算具有更高的预测精度，说明该方法 邻点之间的局部区域，从而得到两点间的局部变化特 在品位估算方面具有良好的应用前景. 征，因而用四维空间分形插值算法更加准确。结合上 述实例验证，四维空间分形插值算法在精度上优于克 参考文献 里金插值法 1]Li G Q,Hu N L.Reserve evaluation model of polymetallic depos- 结论 its based on multidimensional probability distributions.J Univ Sci Technol Beijing,2011,33(3):257 (1)本文通过对分形插值理论的研究，在充分利 (李国清，胡乃联.基于多维概率分布的多金属矿床储量估算 用地质品位的空间信息和发挥分形理论自相似性的基 模型.北京科技大学学报，2011,33(3)：257)
胡乃联等: 四维空间分形插值算法在品位估算中的应用 分布现象，而克里金插值则是表现为相对集中． 这是 由于分形插值具有自相似性，能够体现相邻两插值点 之间的局部变化特征，而克里金插值则由于低通滤波 性而未能反映出品位分布的局部性． 用均方误差对这两种插值方法的插值精度进行对 比． 表 1 列出代表性观测点的插值误差对比． 表 1 钼矿两种插值方法结果 Table 1 Ｒesults of two interpolation methods in molybdenum ore grade estimation y /m x /m z /m 钼品位观测值/10 － 6 分形插值误差 克里金插值误差 465620 5248180 100 594. 8263 0. 011751521 0. 012588814 465700 5248250 100 594. 4836 0. 011729090 0. 013762319 465380 5248670 100 491. 3848 － 0. 007394643 0. 007019079 465460 5248670 100 557. 1922 0. 002370298 0. 004807368 465540 5248670 100 576. 2724 0. 008065163 0. 005567949 465700 5248110 200 594. 733 － 0. 005177851 0. 004436979 465780 5248110 200 595. 4168 0. 055252780 － 0. 059105252 465380 5248600 200 923. 9639 0. 026653377 0. 019697220 465620 5248600 200 1016. 716 － 0. 041018674 0. 062061316 465540 5248670 200 995. 4822 － 0. 010591899 0. 034773408 465780 5248110 300 701. 9488 － 0. 004437747 0. 019388099 465620 5248320 300 1064. 548 － 0. 002040820 0. 030878965 465780 5248460 300 1010. 303 0. 035103251 － 0. 053138138 465540 5248530 300 1937. 489 0. 032239192 0. 042725035 465540 5248670 300 957. 8489 0. 023094816 0. 035682428 465540 5248110 400 587. 7758 － 0. 000612992 － 0. 007780875 465780 5248110 400 660. 2681 － 0. 004770689 0. 001987888 465460 5248390 400 658. 8874 － 0. 008764155 － 0. 018133457 465620 5248390 400 893. 1652 0. 001446660 0. 032450752 465380 5248530 400 1024. 34 0. 006630153 0. 011699876 均方误差 1. 14098 × 10 － 8 2. 20978 × 10 － 8 由表 1 可以看出，采用四维空间分形插值算法估 算的均方误差比克里金插值法降低 48. 4% ，说明在该 钼矿品位估算中，四维空间分形插值算法明显优于克 里金插值法． 从插值原理上看: 克里金插值法中任意两相邻插 值点间的信息是通过光滑曲线连接的，从而掩盖了相 邻两插值点间的局部变化特征，因而具有一定程度的 光滑作用; 而四维空间分形插值算法则是根据整体与 局部相似的原理，将插值数据点的变化特征映射到相 邻点之间的局部区域，从而得到两点间的局部变化特 征，因而用四维空间分形插值算法更加准确． 结合上 述实例验证，四维空间分形插值算法在精度上优于克 里金插值法． 4 结论 ( 1) 本文通过对分形插值理论的研究，在充分利 用地质品位的空间信息和发挥分形理论自相似性的基 础上，提出了四维空间分形插值算法． ( 2) 该方法克服了传统分形插值算法只反映平面 信息的局限性和克里金插值法的低通滤波性，不仅能 反映相邻两插值点之间的局部变化特征，而且能充分 体现地质品位变化规律复杂的特点，使插值结果更加 精确． ( 3) 以某钼矿为例，将该方法应用于品位估算过 程中． 与克里金插值法相比，应用四维空间分形插值 算法进行品位估算具有更高的预测精度，说明该方法 在品位估算方面具有良好的应用前景． 参 考 文 献 ［1］ Li G Q，Hu N L． Ｒeserve evaluation model of polymetallic deposits based on multidimensional probability distributions． J Univ Sci Technol Beijing，2011，33( 3) : 257 ( 李国清，胡乃联． 基于多维概率分布的多金属矿床储量估算 模型． 北京科技大学学报，2011，33( 3) : 257) · 955 ·
·560· 工程科学学报，第37卷，第5期 2]Li J,Li C P,Li Z X.Grade interpolation in orebody based on 机振动识别方法.北京科技大学学报，2013,35(8)：1064) support vector regression.J Unir Sci Technol Beijing,2009,31 [8]Li F,Liu H F,Zhang X J,et al.Determination of spontaneous (12):1498 combusion extent in coal seams on the basis of the fractal theory (李娟，李翠平，李仲学.基于支持向量回归机的矿体品位插 Coal Geol Explor.2013,41(3):15 值.北京科技大学学报，2009,31(12)：1498) (李峰，刘鸿福，张新军，等.基于分形理论确定地下煤层自 B]Li C P.Zheng Y X.Zhang J,et al.Ore grade interpolation model 燃火区范围.煤田地质与勘探，2013,41(3)：15) based on support vector machines optimized by geneticalgorithms. Sun H Q.Fractal Geometry and Fractal Interpolation.Beijing: JUniv Sci Technol Beijing,2013,35(7)837 Science Press,2011 (李翠平，郑瑶瑕，张佳，等.基于遗传算法优化的支持向量 (孙洪泉.分形几何与分形插值.北京：科学出版社，2011) 机品位插值模型.北京科技大学学报，2013,35(7)：837) [10]Shi J F,Wang C N.Fractal analysis of gold deposits in China: 4]Li Q M.Multifractal-Krige interpolation method.Adrances in implication for giant deposit exploration.Earth Sci China Unit Earth Science,2005,20 (2)248 Geosci,1998,23(6):616 (李庆谋.多维分形克里格方法.地球科学进展，2005,20 (施俊法，王春宁.中国金矿床分形分布及对超大型矿床的 (2):248) 勘查意义.地球科学一中国地质大学学报，1998,23(6)： 5]Fan Y H,Luan YZ,Wang Y,et al.The method study of combi- 616) nation of fractal interpolation and linear interpolation.Sci Surr [11]Chen Y,Hu Y A,Lin T.Two-dimensional improved fractal Mapp,2005,30(2):76 model of the sea surface and sea spectrum evaluation.J Com- (范玉红，栾元重，王永，等分形插值与传统插值相结合的 mn,2013,34(2):177 方法研究.测绘科学，2005,30(2)：76) (陈瑜，胡云安，林涛.二维改进分形海面模型及海谱分析. 6]Li J,Wu S C.Gao YT,et al.Rock and soil layer interface fit- 通信学报，2013,34(2)：177) ting method based on a Kriging andclosest point coupled algo- [2]Sun H Q.The study of fractal interpolation on rectangular area rithm.J Univ Sci Technol Beijing,2012,34(5)500 Acta Math Sci.2009,29A(3):773 (李健，吴顺川，高永涛，等.基于Kriging与Closest Point脸 (孙洪泉.矩形域上分形插值研究.数学物理学报，2009， 合算法的边坡岩土层界面拟合.北京科技大学学报，2012， 29A(3):773) 34(5):500) [13]Wu J H,Wang Y J.The spatial interpolation method efficiency ]Mi K F,Zhang J,Cao JG,et al.Vibration identification technol- evaluation of the coal mining isoline based on surfer.China Min ogy of tandem cold rolling mills based on wavelet and fractal analy- Mag,2007,16(1):109 sis.J Univ Sci Technol Beijing,2013,35(8):1064 (武俊红，汪云甲.基于Surfer的煤矿等值线空间插值方法 (米凯夫，张杰，曹建国，等.基于小波和小波分形的冷连轧 有效性评价.中国矿业，2007,16(1)：109)
工程科学学报，第 37 卷，第 5 期 ［2］ Li J，Li C P，Li Z X． Grade interpolation in orebody based on support vector regression． J Univ Sci Technol Beijing，2009，31 ( 12) : 1498 ( 李娟，李翠平，李仲学． 基于支持向量回归机的矿体品位插 值． 北京科技大学学报，2009，31( 12) : 1498) ［3］ Li C P，Zheng Y X，Zhang J，et al． Ore grade interpolation model based on support vector machines optimized by geneticalgorithms． J Univ Sci Technol Beijing，2013，35( 7) : 837 ( 李翠平，郑瑶瑕，张佳，等． 基于遗传算法优化的支持向量 机品位插值模型． 北京科技大学学报，2013，35( 7) : 837) ［4］ Li Q M． Multifractal-Krige interpolation method． Advances in Earth Science，2005，20( 2) : 248 ( 李庆谋． 多维分形克里格方法． 地球科学进展，2005，20 ( 2) : 248) ［5］ Fan Y H，Luan Y Z，Wang Y，et al． The method study of combination of fractal interpolation and linear interpolation． Sci Surv Mapp，2005，30( 2) : 76 ( 范玉红，栾元重，王永，等． 分形插值与传统插值相结合的 方法研究． 测绘科学，2005，30( 2) : 76) ［6］ Li J，Wu S C，Gao Y T，et al． Ｒock and soil layer interface fitting method based on a Kriging andclosest point coupled algorithm． J Univ Sci Technol Beijing，2012，34( 5) : 500 ( 李健，吴顺川，高永涛，等． 基于 Kriging 与 Closest Point 融 合算法的边坡岩土层界面拟合． 北京科技大学学报，2012， 34( 5) : 500) ［7］ Mi K F，Zhang J，Cao J G，et al． Vibration identification technology of tandem cold rolling mills based on wavelet and fractal analysis． J Univ Sci Technol Beijing，2013，35( 8) : 1064 ( 米凯夫，张杰，曹建国，等． 基于小波和小波分形的冷连轧 机振动识别方法． 北京科技大学学报，2013，35( 8) : 1064) ［8］ Li F，Liu H F，Zhang X J，et al． Determination of spontaneous combusion extent in coal seams on the basis of the fractal theory． Coal Geol Explor，2013，41( 3) : 15 ( 李峰，刘鸿福，张新军，等． 基于分形理论确定地下煤层自 燃火区范围． 煤田地质与勘探，2013，41( 3) : 15) ［9］ Sun H Q． Fractal Geometry and Fractal Interpolation． Beijing: Science Press，2011 ( 孙洪泉． 分形几何与分形插值． 北京: 科学出版社，2011) ［10］ Shi J F，Wang C N． Fractal analysis of gold deposits in China: implication for giant deposit exploration． Earth Sci J China Univ Geosci，1998，23( 6) : 616 ( 施俊法，王春宁． 中国金矿床分形分布及对超大型矿床的 勘查意义． 地球科学—中国地质大学学报，1998，23 ( 6) : 616) ［11］ Chen Y，Hu Y A，Lin T． Two-dimensional improved fractal model of the sea surface and sea spectrum evaluation． J Commun，2013，34( 2) : 177 ( 陈瑜，胡云安，林涛． 二维改进分形海面模型及海谱分析． 通信学报，2013，34( 2) : 177) ［12］ Sun H Q． The study of fractal interpolation on rectangular area． Acta Math Sci，2009，29A( 3) : 773 ( 孙洪泉． 矩形域上分形插值研究． 数学物理学报，2009， 29A( 3) : 773) ［13］ Wu J H，Wang Y J． The spatial interpolation method efficiency evaluation of the coal mining isoline based on surfer． China Min Mag，2007，16( 1) : 109 ( 武俊红，汪云甲． 基于 Surfer 的煤矿等值线空间插值方法 有效性评价． 中国矿业，2007，16( 1) : 109) · 065 ·