四维空间分形插值算法在品位估算中的应用

工程科学学报，第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 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 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 depos￾its 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 combi￾nation 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 fit￾ting method based on a Kriging andclosest point coupled algo￾rithm． 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 technol￾ogy of tandem cold rolling mills based on wavelet and fractal analy￾sis． 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 Com￾mun，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 ·