基于Kriging与Closest Point融合算法的边坡岩土层界面拟合

D0L:10.13374/.issn1001-053x.2012.05.006 第34卷第5期 北京科技大学学报 Vol.34 No.5 2012年5月 Journal of University of Science and Technology Beijing May 2012 基于Kriging与Closest Point融合算法的边坡岩土层 界面拟合 李健”吴顺川”s 高永涛” 周 喻”邓增兵) 1)北京科技大学土木与环境工程学院，北京1000832)中煤平朔煤业有限责任公司，朔州036006 通信作者，E-mail:wushunchuan(@usth.edu.cm 摘要以先后发生过两次土质滑坡的平朔公司东露天矿边坡为研究背景，针对目前高精度岩土层界面拟合时有效数据不 足和人工评价误差较大的现状，充分利用岩土层界面的地质调查点数据，有效扩充钻孔数据范围，为高精度的岩土层界面拟 合提供有效数据，并利用基于Closest Point算法的结构分析模型进行Kriging插值计算，实现Kriging算法与Closest Point算法 的融合，完成对露天矿边坡各岩土层的三维拟合.现场选取测量点的比对结果表明，基于Kriging与Closest Point融合算法的 界面拟合，与传统的人工经验评价方法相比，能大幅提高拟合精度，且简单易行. 关键词露天开采：边坡稳定：界面：算法：计算机模拟 分类号TD164.1:TP391.9 Rock and soil layer interface fitting method based on a Kriging and Closest Point coupled algorithm LI Jian,WU Shun-chuan,GAO Yong-tao",ZHOU Yu,DENG Zeng-bing? 1)School of Civil and Environmental Engineering,University of Science and Technology Beijing,Beijing 100083,China 2)Pingshuo Coal Industry Co.,Shuozhou 036006,China Corresponding author,E-mail:wushunchuan@ustb.edu.cn ABSTRACT Taking the slope at the East Open Pit Mine of Pingshuo Coal Industry Co.where two soil landslides occurred succes- sively as the research background,in consideration of the current situation that valid data lacks and relatively large error exists in expe- rience evaluation,data of geological survey points in interfaces between rock and soil layers were fully used and the drilling data range was effectively expanded to provide valid data for high-precision interface fitting.A Closest Point algorithm-based structural analysis model was applied to conduct Kriging interpolation calculations,which made the two algorithms coupled.Thus a three-dimensional fit- ting model of all rock and soil layers can be completed.By comparing the data of measuring points selected from the field,the accuracy of the interface fitting method is much higher than the traditional method of human experience KEY WORDS open pit mining:slope stability:interfaces:algorithms:computer simulation 在露天矿山开采设计中，需要对矿山各边帮基 Krige在20世纪50年代提出的Kriging插值算法 础地质条件进行评价，其中最重要的资料是各岩土 应用最广泛.美国Schuenemeyer等回以及国内李晓 层的分布情况.国内传统矿山多依赖有限的钻孔数 军和靳国栋等B-l针对Kriging算法的应用研究都 据人工经验评价，在此基础上绘制的各岩土层顶底 表明，Kriging插值算法在数据表达和地层分布状况 板等高线及等厚线图往往与实际情况相差甚远. 拟合等方面优于普通插值算法.然而，由于露天矿 近年来，随着科学计算可视化研究的不断深入， 各边帮钻孔数据分布离散性强，Kriging插值算法在 图形仿真技术及其插值算法被引入到矿山三维边坡 露天矿的三维边坡岩土层拟合中鲜有成功的实例. 岩土层拟合领域，国内外学者针对不同插值算法展 随着算法综合应用技术的不断发展，两种甚至多种 开了广泛的探索性应用研究，尤其以南非工程师 算法的综合应用在国内外已展开了相关研究5-) 收稿日期：201108-28 基金项目：国家高技术研究发展计划资助项目(2009AA11Z105):国家自然科学基金资助项目(51074014)：长江学者和创新团队发展计划资助 项目(IRT0950)

第 34 卷 第 5 期 2012 年 5 月 北京科技大学学报 Journal of University of Science and Technology Beijing Vol． 34 No． 5 May 2012 基于 Kriging 与 Closest Point 融合算法的边坡岩土层 界面拟合 李 健1) 吴顺川1) 高永涛1) 周 喻1) 邓增兵2) 1) 北京科技大学土木与环境工程学院，北京 100083 2) 中煤平朔煤业有限责任公司，朔州 036006 通信作者，E-mail: wushunchuan@ ustb． edu． cn 摘 要 以先后发生过两次土质滑坡的平朔公司东露天矿边坡为研究背景，针对目前高精度岩土层界面拟合时有效数据不 足和人工评价误差较大的现状，充分利用岩土层界面的地质调查点数据，有效扩充钻孔数据范围，为高精度的岩土层界面拟 合提供有效数据，并利用基于 Closest Point 算法的结构分析模型进行 Kriging 插值计算，实现 Kriging 算法与 Closest Point 算法 的融合，完成对露天矿边坡各岩土层的三维拟合． 现场选取测量点的比对结果表明，基于 Kriging 与 Closest Point 融合算法的 界面拟合，与传统的人工经验评价方法相比，能大幅提高拟合精度，且简单易行． 关键词 露天开采; 边坡稳定; 界面; 算法; 计算机模拟 分类号 TD164 + . 1; TP391. 9 Rock and soil layer interface fitting method based on a Kriging and Closest Point coupled algorithm LI Jian1) ，WU Shun-chuan1) ，GAO Yong-tao 1) ，ZHOU Yu1) ，DENG Zeng-bing2) 1) School of Civil and Environmental Engineering，University of Science and Technology Beijing，Beijing 100083，China 2) Pingshuo Coal Industry Co． ，Shuozhou 036006，China Corresponding author，E-mail: wushunchuan@ ustb． edu． cn ABSTRACT Taking the slope at the East Open Pit Mine of Pingshuo Coal Industry Co． where two soil landslides occurred succes￾sively as the research background，in consideration of the current situation that valid data lacks and relatively large error exists in expe￾rience evaluation，data of geological survey points in interfaces between rock and soil layers were fully used and the drilling data range was effectively expanded to provide valid data for high-precision interface fitting． A Closest Point algorithm-based structural analysis model was applied to conduct Kriging interpolation calculations，which made the two algorithms coupled． Thus a three-dimensional fit￾ting model of all rock and soil layers can be completed． By comparing the data of measuring points selected from the field，the accuracy of the interface fitting method is much higher than the traditional method of human experience． KEY WORDS open pit mining; slope stability; interfaces; algorithms; computer simulation 收稿日期: 2011--08--28 基金项目: 国家高技术研究发展计划资助项目( 2009AA11Z105) ; 国家自然科学基金资助项目( 51074014) ; 长江学者和创新团队发展计划资助 项目( IRT0950) 在露天矿山开采设计中，需要对矿山各边帮基 础地质条件进行评价，其中最重要的资料是各岩土 层的分布情况． 国内传统矿山多依赖有限的钻孔数 据人工经验评价，在此基础上绘制的各岩土层顶底 板等高线及等厚线图往往与实际情况相差甚远． 近年来，随着科学计算可视化研究的不断深入， 图形仿真技术及其插值算法被引入到矿山三维边坡 岩土层拟合领域，国内外学者针对不同插值算法展 开了广泛的探索性应用研究，尤其以南非工程师 Krige ［1］在 20 世纪 50 年代提出的 Kriging 插值算法 应用最广泛． 美国 Schuenemeyer 等［2］以及国内李晓 军和靳国栋等［3--4］针对 Kriging 算法的应用研究都 表明，Kriging 插值算法在数据表达和地层分布状况 拟合等方面优于普通插值算法． 然而，由于露天矿 各边帮钻孔数据分布离散性强，Kriging 插值算法在 露天矿的三维边坡岩土层拟合中鲜有成功的实例． 随着算法综合应用技术的不断发展，两种甚至多种 算法的综合应用在国内外已展开了相关研究［5--9］． DOI:10.13374/j.issn1001-053x.2012.05.006

第5期 李健等：基于Kriging与Closest Point融合算法的边坡岩土层界面拟合 ·501· 该技术可以充分发挥各算法的优势，有效克服各算 受到主观影响，往往误差较大，究竞采取何种方法才 法的不足.基于算法融合技术，在Kriging算法插值 能最大程度地开发和利用确定性数据，一直是国内 计算离散数据过程中引入Closest Point算法，可以 外研究的热点.加拿大的Houlding较早地提出了 很好地解决露天矿边坡数据离散的问题，实现露天 基于钻孔数据进行三维地质建模的概念；Lemon 矿的三维边坡岩土层拟合 等田基于钻孔数据建模后又采用自定义剖面的方 本文以平朔公司东露天矿边坡为研究背景，通 法进行动态修正：国内学者刘振平等回基于钻孔数 过对钻孔数据和地质点调查数据的充分挖掘利用， 据对三维可视化快速建模技术进行了应用研究.这 使用Kriging与Closest Point融合算法，实现该露天 些研究大多针对数据离散性不大的工程数据进行插 矿边坡各岩土层的拟合 值拟合，对于露天矿边坡这种广域离散型数据的插 值拟合，一直缺乏相关的研究探索. 1 地质数据利用与工程背景 针对露天矿边坡钻孔数据离散性大的特点，将 露天矿地质数据种类繁多，根据数据来源的可 地质调查点数据转化成微型钻孔数据，实现钻孔数 靠性可将其分为三类：确定性数据（如工程测量数 据的广域分布，最大程度开发和利用已有的钻孔和 据、钻孔数据)、不确定性数据（如各岩土层厚度的 地质点调查等确定性数据.在本次应用研究中，原 空间变化特征，主要依靠主观解释、插值等手段获 始钻孔共97个，基本均匀分布于露天矿的境界周边 得)以及知识推理型数据（如煤层底板等高线图、剖 (图1(a)),将针对发生过两次滑坡的东北帮附近 面地质图，是地质工程师根据己有数据结合专家经 进行的704个地质调查点数据转化成704个虚拟微 验知识推理推断得出) 型钻孔，完成钻孔共计801个（如图1(b)所示，图中 在露天煤矿地质数据中，确定性数据有效性高， 东北帮附近小圆点即为虚拟微型钻孔)，有效扩充 是最值得信赖的数据来源，但其采集成本高，数量有 了钻孔数据范围，为后期高精度的拟合计算提供了 限，而不确定性数据和知识推理型数据很大程度上 数据来源 北 492000493600494000 495000 496000497000 8491000 492000 496000 497000 (b) 是 490 4940 495000 497000 49300 4940049300496000 南 图1确定性数据分布.(a)原始钻孔分布：(b)加工处理后钻孔分布 Fig.I Deterministic data distribution:(a)primitive drilling distribution:(b)processed drilling distribution 2边坡岩土层拟合模型与方法 函数，记为y(x,h): 2.1区域化变量 y,h)=2yarZ)-Zx+月]. (1) 设露天矿边坡区域为A,则以边坡范围内空间 2.3平稳假设和内蕴假设 点x的三个直角坐标(x,x,x）为自变量的随机场 对于露天矿边坡岩土层拟合问题，为进行下一 Z(x。,x,x)=Z(x)即为该边坡范围内的区域化变 步的插值计算，提出二阶平稳假设和内蕴假设 量，其中，Z(x)∈A (1)平稳假设.当区域化变量Z(x)满足以下 2.2变异函数 两个条件时，称边坡范围内区域化变量Z(x)满足二 当边坡范围内空间点x在一维x轴上变化时， 阶平稳假设 把区域化变量在x与x+h处的值Z(x)与Z(x+h) ①区域化变量Z(x)的数学期望存在且等于常 的差的方差之半定义为Z(x)在x轴方向上的变异 数，即

第 5 期 李 健等: 基于 Kriging 与 Closest Point 融合算法的边坡岩土层界面拟合 该技术可以充分发挥各算法的优势，有效克服各算 法的不足． 基于算法融合技术，在 Kriging 算法插值 计算离散数据过程中引入 Closest Point 算法，可以 很好地解决露天矿边坡数据离散的问题，实现露天 矿的三维边坡岩土层拟合． 本文以平朔公司东露天矿边坡为研究背景，通 过对钻孔数据和地质点调查数据的充分挖掘利用， 使用 Kriging 与 Closest Point 融合算法，实现该露天 矿边坡各岩土层的拟合． 1 地质数据利用与工程背景 露天矿地质数据种类繁多，根据数据来源的可 靠性可将其分为三类: 确定性数据( 如工程测量数 据、钻孔数据) 、不确定性数据( 如各岩土层厚度的 空间变化特征，主要依靠主观解释、插值等手段获 得) 以及知识推理型数据( 如煤层底板等高线图、剖 面地质图，是地质工程师根据已有数据结合专家经 验知识推理推断得出) ． 在露天煤矿地质数据中，确定性数据有效性高， 是最值得信赖的数据来源，但其采集成本高，数量有 限，而不确定性数据和知识推理型数据很大程度上 受到主观影响，往往误差较大，究竟采取何种方法才 能最大程度地开发和利用确定性数据，一直是国内 外研究的热点． 加拿大的 Houlding ［10］较早地提出了 基于钻孔数据进行三维地质建模的概念; Lemon 等［11］基于钻孔数据建模后又采用自定义剖面的方 法进行动态修正; 国内学者刘振平等［12］基于钻孔数 据对三维可视化快速建模技术进行了应用研究． 这 些研究大多针对数据离散性不大的工程数据进行插 值拟合，对于露天矿边坡这种广域离散型数据的插 值拟合，一直缺乏相关的研究探索． 针对露天矿边坡钻孔数据离散性大的特点，将 地质调查点数据转化成微型钻孔数据，实现钻孔数 据的广域分布，最大程度开发和利用已有的钻孔和 地质点调查等确定性数据． 在本次应用研究中，原 始钻孔共 97 个，基本均匀分布于露天矿的境界周边 ( 图 1( a) ) ，将针对发生过两次滑坡的东北帮附近 进行的 704 个地质调查点数据转化成 704 个虚拟微 型钻孔，完成钻孔共计 801 个( 如图 1( b) 所示，图中 东北帮附近小圆点即为虚拟微型钻孔) ，有效扩充 了钻孔数据范围，为后期高精度的拟合计算提供了 数据来源． 图 1 确定性数据分布． ( a) 原始钻孔分布; ( b) 加工处理后钻孔分布 Fig． 1 Deterministic data distribution: ( a) primitive drilling distribution; ( b) processed drilling distribution 2 边坡岩土层拟合模型与方法 2. 1 区域化变量 设露天矿边坡区域为 A，则以边坡范围内空间 点 x 的三个直角坐标( xu，xv，xw ) 为自变量的随机场 Z( xu，xv，xw ) = Z( x) 即为该边坡范围内的区域化变 量，其中，Z( x) ∈A． 2. 2 变异函数 当边坡范围内空间点 x 在一维 x 轴上变化时， 把区域化变量在 x 与 x + h 处的值 Z( x) 与 Z( x + h) 的差的方差之半定义为 Z( x) 在 x 轴方向上的变异 函数，记为 γ( x，h) : γ( x，h) = 1 2 Var［Z( x) － Z( x + h) ］． ( 1) 2. 3 平稳假设和内蕴假设 对于露天矿边坡岩土层拟合问题，为进行下一 步的插值计算，提出二阶平稳假设和内蕴假设． ( 1) 平稳假设． 当区域化变量 Z( x) 满足以下 两个条件时，称边坡范围内区域化变量 Z( x) 满足二 阶平稳假设． ① 区域化变量 Z( x) 的数学期望存在且等于常 数，即 ·501·

·502 北京科技大学学报 第34卷 E[☑(]=m(m为常数) (2) 0, h=0 ②区域化变量Z(x)的空间协方差函数存在且 y(h)= Co +C 0a. C(h),Yx,h. (3) (10) 在二阶平稳假设下， 式中，C。为块金常数，C为基台值，a为变程 E[☑(x+h)]=m=E[☑(x)],Hh,Hx.(4) (2)高斯模型.通式为 变异函数可改写为 0, h=0: y(h)= y(s.)=Var [z()-Z(s+]-E(Z()- C。+c(1-e5),h>0. (11) Zx+月]2-EZ(国]-EZx+)]2= 式中a不是变程，由于当h=5a,l-e兰-1-e3≈ 0.95≈1，即当h=3a时，y(h)≈C。+C,所以该模 E()-7 (5) 型的变程为3a,其他符号含义同式(10) (2)内蕴假设.当区域化变量Z(x)的增量 (3)指数模型.一般公式为 Z(x)-Z(x+h)满足以下两个条件时，称边坡范围 rO, h=0: y(h)= (12) 内区域化变量Z(x)满足内蕴假设. C。+c(1-e),h>0. ①区域化变量Z(x)的增量Z(x)-Z(x+h)的 式中，当h=3a时，1-e是=1-e=0.95=1,即 数学期望为0，即 当h=3a时，y(h)≈C。+C,所以该模型变程为3a, E [Z (x)-Z(x+h)]=0,h,Vx. (6) 其他符号含义同式(10). ②对于所有矢量的增量Z(x)-Z(x+h)的方 2.6基于Closest Point算法的结构分析 差函数存在且平稳，即 (I)Closest Point算法.Closest Point算法是在 Var☑(x)-Z(x+)]=E[☑(x)-Z(x+h)]2= 最小生成树算法Pim算法基础上改造而成的，其基 2y(x,h)=2y(h),Yx,Yh. (7) 本原理为：选取区域内任一点，为初始环，找出离环最 在内蕴假设下，E[☑(x)-Z(x+)]2仅依赖 近的点u,将u插到v后面得到一个新的环，重复以上 于分割他们的距离IhI和方向a,与点x的位置无 的步骤直到设定的点间距离或所有点都包在环内. 关，因此变异函数可改写为 (2)结构分析.在露天矿边坡岩土层拟合过程 yh,a)=EZ()-zx+0]2. (8) 中，由于区域化变量Z(x)的变化性很复杂，在不同 方向上有不同的变化性，或者在一个方向上包含着 2.4试验变异函数 不同尺度的多层次变化性，因而无法用一个理论模 在露天矿边坡岩土层拟合过程中，钻孔得到的 型拟合。为全面了解区域化变量的变异性，必须进 数据点是有限的，把所有已知点数据构制的变异函 行结构化分析，即构造一个变异函数模型对于全部 数称之为试验变异函数，记为y(h): 有效结构信息作定量化的概括，以表征区域化变量 1 N(A) V=2N三)-z6+0]3o) 的主要特征.结构分析的主要方法是套合结构，即 把分别出现在不同距离Ih1上和（或）不同方向a上 式中，N(h)为被矢量h相分隔的试验数据对的 同时起作用的点数据变异性组合起来.套合结构可 数目. 以表示为多个变异函数之和，每一个变异函数代表 2.5理论模型 一种特定尺度上的变异性，套合结构的表达式为 变异函数构制完成后，可以解决边坡范围内区 y(h)=yo(h)+y1(h)+…+y:(h)+….(13) 域化变量Z(x)的变化特征和结构性状的问题，要完 传统的套合结构在露天矿边坡拟合过程中没有 成整个露天矿边坡岩土层的拟合，就必须给试验变 考虑其数据离散的特点，将离散性很大的两侧边坡 异函数配以相应的理论模型，直接参与到Kriging插 钻孔数据组合到一起进行拟合，拟合结果与实际相 值计算中.最常用的理论模型有球状模型、高斯模 差较大. 型及指数模型闭 本文基于Closest Point算法，设定露天矿边坡 (1)球状模型.一般公式为 实际宽度H为终止的点间距离，每一个岩土层设置

北 京 科 技 大 学 学 报 第 34 卷 E［Z( x) ］= m ( m 为常数) ． ( 2) ② 区域化变量 Z( x) 的空间协方差函数存在且 平稳，即 Cov［Z( x) ，Z( x + h) ］= E［Z( x) Z( x + h) ］－ m2 = C( h) ，x，h． ( 3) 在二阶平稳假设下， E［Z( x + h) ］= m = E［Z( x) ］，h，x． ( 4) 变异函数可改写为 γ( x，h) = 1 2 Var［Z( x) － Z( x + h) ］= 1 2 E［Z( x) － Z( x + h) ］2 － 1 2 { E［Z( x) ］－ E［Z( x + h) ］} 2 = 1 2 E［Z( x) － Z( x + h) ］2 ． ( 5) ( 2) 内蕴假设． 当区域化变量 Z ( x) 的增量 Z( x) %Z( x + h) 满足以下两个条件时，称边坡范围 内区域化变量 Z( x) 满足内蕴假设． ① 区域化变量 Z( x) 的增量 Z( x) %Z( x + h) 的 数学期望为 0，即 E［Z( x) － Z( x + h) ］= 0，h，x． ( 6) ② 对于所有矢量的增量 Z( x) %Z( x + h) 的方 差函数存在且平稳，即 Var［Z( x) － Z( x + h) ］= E［Z( x) － Z( x + h) ］2 = 2γ( x，h) = 2γ( h) ，x，h． ( 7) 在内蕴假设下，E ［Z( x) － Z( x + h) ］2 仅依赖 于分割他们的距离 | h | 和方向 α，与点 x 的位置无 关，因此变异函数可改写为 γ( h，α) = 1 2 E［Z( x) － Z( x + h) ］2 ． ( 8) 2. 4 试验变异函数 在露天矿边坡岩土层拟合过程中，钻孔得到的 数据点是有限的，把所有已知点数据构制的变异函 数称之为试验变异函数，记为 γ* ( h) : γ* ( h) = 1 2N( h) ∑ N( h) i = 1 ［Z( xi ) － Z( xi + h) ］2 ． ( 9) 式中，N( h) 为被矢量 h 相分隔的试验数据对的 数目． 2. 5 理论模型 变异函数构制完成后，可以解决边坡范围内区 域化变量 Z( x) 的变化特征和结构性状的问题，要完 成整个露天矿边坡岩土层的拟合，就必须给试验变 异函数配以相应的理论模型，直接参与到 Kriging 插 值计算中． 最常用的理论模型有球状模型、高斯模 型及指数模型［13］． ( 1) 球状模型． 一般公式为 γ( h) = 0， h = 0; C0 + C ( 3 2 × h a － 1 2 × h3 a3 ) ， 0 ＜ h≤a; C0 + C， h ＞      a． ( 10) 式中，C0为块金常数，C 为基台值，a 为变程． ( 2) 高斯模型． 通式为 γ( h) = 0， h = 0; C0 + C( 1 － e － h2 a2 { ) ， h ＞ 0． ( 11) 式中 a 不是变程，由于当 h =槡3a，1 － e － h2 a2 = 1 － e －3 ≈ 0. 95≈1，即当 h = 槡3a 时，γ( h) ≈C0 + C，所以该模 型的变程为槡3a，其他符号含义同式( 10) ． ( 3) 指数模型． 一般公式为 γ( h) = 0， h = 0; C0 + C( 1 － e － h a { ) ， h ＞ 0． ( 12) 式中，当 h = 3a 时，1 － e － h2 a2 = 1 － e － 3 ≈0. 95≈1，即 当 h = 3a 时，γ( h) ≈C0 + C，所以该模型变程为 3a， 其他符号含义同式( 10) ． 2. 6 基于 Closest Point 算法的结构分析 ( 1) Closest Point 算法． Closest Point 算法是在 最小生成树算法 Prim 算法基础上改造而成的［14］，其基 本原理为: 选取区域内任一点 v 为初始环，找出离环最 近的点 u，将 u 插到 v 后面得到一个新的环，重复以上 的步骤直到设定的点间距离或所有点都包在环内． ( 2) 结构分析． 在露天矿边坡岩土层拟合过程 中，由于区域化变量 Z( x) 的变化性很复杂，在不同 方向上有不同的变化性，或者在一个方向上包含着 不同尺度的多层次变化性，因而无法用一个理论模 型拟合． 为全面了解区域化变量的变异性，必须进 行结构化分析，即构造一个变异函数模型对于全部 有效结构信息作定量化的概括，以表征区域化变量 的主要特征． 结构分析的主要方法是套合结构，即 把分别出现在不同距离 | h | 上和( 或) 不同方向 α 上 同时起作用的点数据变异性组合起来． 套合结构可 以表示为多个变异函数之和，每一个变异函数代表 一种特定尺度上的变异性，套合结构的表达式为 γ( h) = γ0 ( h) + γ1 ( h) + … + γi ( h) + …． ( 13) 传统的套合结构在露天矿边坡拟合过程中没有 考虑其数据离散的特点，将离散性很大的两侧边坡 钻孔数据组合到一起进行拟合，拟合结果与实际相 差较大． 本文基于 Closest Point 算法，设定露天矿边坡 实际宽度 H 为终止的点间距离，每一个岩土层设置 ·502·

第5期 李健等：基于Kriging与Closest Point融合算法的边坡岩土层界面拟合 ·503· 一个变异函数y:(h)以表征其变异性，如下式所示： 推导可得估计方差的计算公式为 「Y:(h),0<Ihl≤H: σ=C(A,A)+C(v,)-2C(A,),( 23) y:(h)= (14) 0, Ihl H. 用变异函数表示为 在构造套合结构的过程中，可以把结构模型 o2=y(A,v）-y(A,A)-2y(u,u). (24) Y(h)看成由N个各向同性结构{y:(Ih1),i=1,2, 式中，C(4,v)及y(A,v)分别为当矢量的两个端点 …,}套合而成，即 各自独立扫过域A和信息域v的协方差函数平均值 之Y.hD 及变异函数平均值，C(v,v)及y(v,v)分别为当矢 y(h)= (15) 量的两个端点各自独立扫过任两个信息域v的协方 式中各个组成结构y:(Ih:I)的特有的各向异性是用 差函数平均值及变异函数平均值，C(A,A)及y(A, 线性变换矩阵A]来表示的，它将矢量h变换成h, A)分别为当矢量两个端点各自独立的在域A扫过 Ch]=[A 0h,]. (16) 时的协方差函数平均值及变异函数平均值.当估计 从而使结构变为各向同性。构成代替露天矿边坡范 量为加权平均值时，估计方差的公式可表示为 围内区域化变量Z(x)变异特征的结构模型. 2=C(A,A)-2∑，C(A,v)+ 2.7 Kriging插值算法 ∑∑，c(u) (25) 如前文所述，露天矿边坡区域为A,所有露天矿 边坡已知点x:(i=1,2,…,n)的三个直角坐标 要求在无偏条件下，使估计方差σ达到极小 (xmi,xa,xa)为自变量对应的属性值Z(x)(i=1, 的诸权系数入：(i=1,2,…,n)是个求条件极值的 2,…,n),Z(x,)∈A,未知点对应的属性真值Z(x。) 问题，为了便于求解，将E{☑(x)-Z(x)]}化 的估计值Z(x。)为己知点属性值的加权组合，即 为无约束的拉格朗日乘数法求极值的问题，即将 Z(x)=∑A,Z(x) (17) ∑入：=1引入目标函数之中，令 = 其中，入，(i=1,2,…,n)为待求权系数，Kriging插 (26) 值计算目的是求出权系数入，使Z'(xo)是Z(xo)的 F=2-24(-1） 无偏估计量，且估计方差最小. 式中，F是n个权系数入：(i=1,2,…,n)和u的 (1)无偏条件.要使Z(x。)是Z(x)的无偏 (n+1)元函数，σ为估计方差，-24是拉格朗日乘 估计量，即要求 数.求出F对入，和u的偏导数并令其为0，最后得 E{Z(x)-Z(x)}=0, (18) 到Kriging方程组，表达式为 在二阶平稳条件下， E{Z(x)}=E{Z(x)}=m, 含Ado-gaW. (19) (i=1,2,…,n). 而 Eiz()Z (27) 基于内蕴假设，用y(h)表示如下 AE(Z())m (20) Ajr(v;v)+u=y(v;,A), 要使 (i=1,2,…,n). EZ1=EZ,1,则三4m=m,即 名A=引 (28) (21) 这是一个n+1个未知数(n个入：和u),n+1个方 式(21)即无偏条件. 程的方程组，求出所有未知数即可完成插值计算. (2)Kriging方程组求解.在边坡范围内要使 Z(xo)对Z(x)的估计方差σ最小，其表达式为 3工程应用 oi=Var[☑(xo)-Z(xo）]→0， (22) 本文选取的工程对象为平朔公司东露天矿边 而Z(x)是一个二阶平稳的随机函数，其期望为m, 坡.该矿东北帮附近的土质边坡在当地季风气候雨 协方差为C(h),变异函数为y(h)且只依赖于h. 季的三个月时间内，先后发生过两次规模较大的滑

第 5 期 李 健等: 基于 Kriging 与 Closest Point 融合算法的边坡岩土层界面拟合 一个变异函数 γi ( h) 以表征其变异性，如下式所示: γi ( h) = γi ( h) ， 0 ＜ | h |≤H; {0， | h | ＞ H． ( 14) 在构造套合结构的过程中，可以把结构模型 γ( h) 看成由 N 个各向同性结构{ γi ( | h | ) ，i = 1，2， …，N} 套合而成，即 γ( h) = ∑ N i = 1 γi ( | hi | ) ． ( 15) 式中各个组成结构 γi ( | hi | ) 的特有的各向异性是用 线性变换矩阵［Ai ］来表示的，它将矢量 h 变换成 h'， ［h' i ］ = ［Ai ］［hi ］． ( 16) 从而使结构变为各向同性． 构成代替露天矿边坡范 围内区域化变量 Z( x) 变异特征的结构模型． 2. 7 Kriging 插值算法 如前文所述，露天矿边坡区域为 A，所有露天矿 边坡已知点 xi ( i = 1，2，…，n) 的三个直角坐标 ( xui，xvi，xwi ) 为自变量对应的属性值 Z( xi ) ( i = 1， 2，…，n) ，Z( xi ) ∈A，未知点对应的属性真值 Z( x0 ) 的估计值 Z* ( x0 ) 为已知点属性值的加权组合，即 Z* ( x0 ) = ∑ n i = 1 λiZ( xi ) ． ( 17) 其中，λi ( i = 1，2，…，n) 为待求权系数，Kriging 插 值计算目的是求出权系数 λi，使 Z* ( x0 ) 是 Z( x0 ) 的 无偏估计量，且估计方差最小． ( 1) 无偏条件． 要使 Z* ( x0 ) 是 Z( x0 ) 的无偏 估计量，即要求 E{ Z( xi ) － Z* ( x0 ) } = 0， ( 18) 在二阶平稳条件下， E{ Z( xi ) } = E{ Z* ( x0 ) } = m， ( 19) 而 E{ Z* ( x0 ) } = E { ∑ n i = 1 λiZ( xi ) } = ∑ n i = 1 λi ·E{ Z( xi ) } = ∑ n i = 1 λi ·m， ( 20) 要使 E{ Z( xi ) } = E{ Z* ( x0 ) } ，则 ∑ n i = 1 λi ·m = m，即 ∑ n i = 1 λi = 1． ( 21) 式( 21) 即无偏条件． ( 2) Kriging 方程组求解． 在边坡范围内要使 Z* ( x0 ) 对 Z( x0 ) 的估计方差 σ2 E 最小，其表达式为 σ2 E = Var［Z( x0 ) － Z* ( x0) ］0， ( 22) 而 Z( x) 是一个二阶平稳的随机函数，其期望为 m， 协方差为 C( h) ，变异函数为 γ( h) 且只依赖于 h． 推导可得估计方差的计算公式为 σ2 E = C( A，A) + C( υ，υ) － 2C( A，υ) ， ( 23) 用变异函数表示为 σ2 E = γ( A，υ) － γ( A，A) － 2γ( υ，υ) ． ( 24) 式中，C( A，υ) 及 γ( A，υ) 分别为当矢量的两个端点 各自独立扫过域 A 和信息域 υ 的协方差函数平均值 及变异函数平均值，C( υ，υ) 及 γ( υ，υ) 分别为当矢 量的两个端点各自独立扫过任两个信息域 υ 的协方 差函数平均值及变异函数平均值，C( A，A) 及 γ( A， A) 分别为当矢量两个端点各自独立的在域 A 扫过 时的协方差函数平均值及变异函数平均值． 当估计 量为加权平均值时，估计方差的公式可表示为 σ2 E = C( A，A) － 2 ∑i λiC( A，υi ) + ∑i ∑ j λiλj C( υi，υj ) ． ( 25) 要求在无偏条件下，使估计方差 σ2 E 达到极小 的诸权系数 λi ( i = 1，2，…，n) 是个求条件极值的 问题，为了便于求解，将 E{ ［Z( x0 ) － Z* ( x0) ］2 } 化 为无约束的拉格朗日乘数法求极值的问题，即将 ∑ n i λi = 1 引入目标函数之中，令 F = σ2 E － 2μ ( ∑ n i = 1 λi － ) 1 ． ( 26) 式中，F 是 n 个权系数 λi ( i = 1，2，…，n) 和 μ 的 ( n + 1) 元函数，σ2 E 为估计方差，－ 2μ 是拉格朗日乘 数． 求出 F 对 λi和 μ 的偏导数并令其为 0，最后得 到 Kriging 方程组，表达式为 ∑ n j = 1 λj C( υi，υj ) － μ = C( υi，A) ， ∑ n i = 1 λi = 1 { ， ( i = 1，2，…，n) ． ( 27) 基于内蕴假设，用 γ( h) 表示如下 ∑ n j = 1 λjγ( υi，υj ) + μ = γ( υi，A) ， ∑ n j = 1 λj = 1 { ， ( i = 1，2，…，n) ． ( 28) 这是一个 n + 1 个未知数( n 个 λi 和 μ) ，n + 1 个方 程的方程组，求出所有未知数即可完成插值计算． 3 工程应用 本文选取的工程对象为平朔公司东露天矿边 坡． 该矿东北帮附近的土质边坡在当地季风气候雨 季的三个月时间内，先后发生过两次规模较大的滑 ·503·

·504 北京科技大学学报 第34卷 坡，塌方总量达到1.78×10m3,严重影响到矿山的 Kriging与Closest Point融合算法，可实现对露天矿 正常生产建设.考虑到当地土层的工程特殊性，为 边坡的岩土层拟合。由于之前发生的两次滑坡均发 尽快展开下一步的滑坡处理和待开挖边坡的稳定性 生在土层尤其是粉质黏土层，本次岩土层拟合重点 研究工作，亟需掌握矿区内各岩土层分布情况. 为粉质黏土层及其上覆粉土层两个关键土层，其三 通过前期对确定性数据的开发和利用，基于 维模型见图2. (a) 粉土层 粉质黏土层 图2边坡关键层三维模型.()粉土层三维模型：(b)粉质黏土层三维模型 Fig.2 3D models of the key strata in the slope:(a)3D model of the silty soil:(b)3D model of the silty clay layer 为验证岩土层拟合效果，随着边坡的不断开挖， 经验评价相比（准确率一般为50%~70%左右），有 本文在三个月时间内先后随机选取了200个测量点 了大幅度的提高，可以作为矿山判定各边坡岩土层 进行地质点调查，并将调查结果与拟合模型进行对 分布情况的重要依据. 比，发现有179个调查点岩土层性质与模型预测相 由三维模型可知，区域内粉质黏土层厚度及分 符，即本次界面拟合整体准确率达到89.5%.在有 布变化十分明显，需对其进行重点分析.基于三维 限的确定性数据情况下，本文基于Kriging与Closest 模型数据，可生成粉质黏土层顶底板等高线图和等 Point融合算法对边坡各岩土层的拟合与传统人工 厚线图，见图3和图4. 8491000 492000 49g000 49500 496000 470 84910 492000 49300 4970 8 (b) 491000 400049400 4950 496000 497000 49时00049400049560 49600 497000 南 南 图3粉质黏土层等高线图.()粉质黏土层顶板等高线图：(b)粉质黏土层底板等高线图 Fig.3 Contour maps of the silty clay layer:(a)contour map of the silty clay layer roof:(b)contour map of the silty clay layer floor 由图3和图4可知：该矿边坡范围内，粉质黏土 沿矿坑边界方向不足200m的区域内，粉质黏土厚 层厚度分布呈现北帮最厚、东北帮次之、南端帮和西 度介于20~60m之间，厚度差达40m左右，而此处 端帮较薄的特点：从粉质黏土层厚度变化来看，呈现 正是发生第二次滑坡的区域. 东北帮和北帮变化最为明显，南端帮次之，西端帮变 基于两次滑坡实例分析，对该矿其他未开挖边 化较小的特点.以东北帮与北帮结合部为例，在沿 帮的土层分布情况做出初步预测，划定两个区域为 矿坑边界方向不足500m的区域内，粉质黏土厚度 重点关注区.目前，工程己开挖至南侧重点关注区， 介于30~100m,厚度变化达70m左右，而此处正是 经现场生产实践表明，该区域遇到了与滑坡区类似 发生第一次滑坡的区域：东北帮东段粉质黏土层，在 的工程难题，区域内土层厚度及分布变化明显，导致

北 京 科 技 大 学 学 报 第 34 卷 坡，塌方总量达到 1. 78 × 105 m3 ，严重影响到矿山的 正常生产建设． 考虑到当地土层的工程特殊性，为 尽快展开下一步的滑坡处理和待开挖边坡的稳定性 研究工作，亟需掌握矿区内各岩土层分布情况． 通过前期对确定性数据的开发和利用，基于 Kriging 与 Closest Point 融合算法，可实现对露天矿 边坡的岩土层拟合． 由于之前发生的两次滑坡均发 生在土层尤其是粉质黏土层，本次岩土层拟合重点 为粉质黏土层及其上覆粉土层两个关键土层，其三 维模型见图 2． 图 2 边坡关键层三维模型． ( a) 粉土层三维模型; ( b) 粉质黏土层三维模型 Fig． 2 3D models of the key strata in the slope: ( a) 3D model of the silty soil; ( b) 3D model of the silty clay layer 为验证岩土层拟合效果，随着边坡的不断开挖， 本文在三个月时间内先后随机选取了 200 个测量点 进行地质点调查，并将调查结果与拟合模型进行对 比，发现有 179 个调查点岩土层性质与模型预测相 符，即本次界面拟合整体准确率达到 89. 5% ． 在有 限的确定性数据情况下，本文基于 Kriging 与 Closest Point 融合算法对边坡各岩土层的拟合与传统人工 经验评价相比( 准确率一般为 50% ～ 70% 左右) ，有 了大幅度的提高，可以作为矿山判定各边坡岩土层 分布情况的重要依据． 由三维模型可知，区域内粉质黏土层厚度及分 布变化十分明显，需对其进行重点分析． 基于三维 模型数据，可生成粉质黏土层顶底板等高线图和等 厚线图，见图 3 和图 4． 图 3 粉质黏土层等高线图． ( a) 粉质黏土层顶板等高线图; ( b) 粉质黏土层底板等高线图 Fig． 3 Contour maps of the silty clay layer: ( a) contour map of the silty clay layer roof; ( b) contour map of the silty clay layer floor 由图 3 和图 4 可知: 该矿边坡范围内，粉质黏土 层厚度分布呈现北帮最厚、东北帮次之、南端帮和西 端帮较薄的特点; 从粉质黏土层厚度变化来看，呈现 东北帮和北帮变化最为明显，南端帮次之，西端帮变 化较小的特点． 以东北帮与北帮结合部为例，在沿 矿坑边界方向不足 500 m 的区域内，粉质黏土厚度 介于 30 ～ 100 m，厚度变化达 70 m 左右，而此处正是 发生第一次滑坡的区域; 东北帮东段粉质黏土层，在 沿矿坑边界方向不足 200 m 的区域内，粉质黏土厚 度介于 20 ～ 60 m 之间，厚度差达 40 m 左右，而此处 正是发生第二次滑坡的区域． 基于两次滑坡实例分析，对该矿其他未开挖边 帮的土层分布情况做出初步预测，划定两个区域为 重点关注区． 目前，工程已开挖至南侧重点关注区． 经现场生产实践表明，该区域遇到了与滑坡区类似 的工程难题，区域内土层厚度及分布变化明显，导致 ·504·

第5期 李健等：基于Kriging与Closest Point融合算法的边坡岩土层界面拟合 ·505· B]Schuenemeyer J H,Power HC.Uncertainty estimation for resource 491 assessment:an application to coal.Math Geol,2000,32(5):521 188 第一次滑坡 Li X J,Hu J H,Zhu HH,et al.The estimation of coal thickness based on Kriging technique and 3D coal seam modeling.China Coal Soc,2008,33(7):765 30505 (李晓军，胡金虎，朱合华，等.基于Kriging方法的煤层厚度估 计及三维煤层建模.煤炭学报，2008,33(7)：765) [4] Jin C D,Liu YC,Niu WJ.Comparison between inverse distance 第二次滑坡 weighting method and Kriging.Changchun Uni Technol,2003. 2E 409 49 49 496497 24(3):53 (新国栋，刘衍聪，牛文杰.距离加权反比插值法和克里金插值 图4粉质黏土层等厚线图 法的比较.长春工业大学学报，2003,24(3)：53) Fig.4 Isopachous map of the silty clay layer [5] Zou G D.Chen Y P.Coupling algorithm of simulated annealing al- gorithm and random search method for slope stability analysis.Chin 工程现场出现小规模塌方，见图5 J Rock Mech Eng,2004,23 (12),2032 (邹广电，陈永平.滑坡和边坡稳定性分析的模拟退火随机搜 索耦合算法.岩石力学与工程学报，2004,23(12)：2032) [6]Shi L,Li X C.Ren W,et al.Hybrid of ant colony algorithm and ge- netic algorithm and its application to searching critical slope slip surface.Rock Soil Mech,2009.30(11):3486 (石露，李小春，任伟，等.蚁群算法与遗传算法融合及其在边 坡临界滑动面搜索中的应用.岩土力学，2009,30(11)：3486) 7]Liu Y,Li Y,Xu J T.Application of a multi-algorithm fusion real 图5南侧重点关注区现场图 time filter in FOG.Acta Photonica Sin,2010,39(6):1116 Fig.5 Picture of the souther area which we pay close attention to (刘颖，李言，徐金涛.一种多算法融合的实时滤波在光纤陀螺 与两次滑坡不同的是，开挖本区域时，当地的季 中的应用.光子学报，2010,39(6)：1116) 8] Hu J J,Tang C J,Li C,et al.Efficient multi-ayer clustering algo- 风性气候已进入旱季，因而水对边坡影响较小，加之 rithm based on nearest neighbors first.I Sichuan Unir Eng Sci, 施工过程中，生产单位高度重视，对本区域的施工情 2004,36(6):93 况重点关注，及时采取了相应安全措施，因而未发生 (胡建军，唐常杰，李川，等.基于最近邻优先的高效聚类算法 大规模滑坡. 四川大学学报：工程科学版，2004,36(6)：93) Wu J H,Zhang J,Liu Z H.Solution of TSP problem based on the 4结论 combination of ant colony algorithm and immune algorithm./Hu- nan Unir Nat Sci,2009,36(10):81 (1)充分利用岩土层界面的地质调查点数据， (吴建辉，章兢，刘朝华.基于蚁群算法和免疫算法融合的TSP 可有效扩充钻孔数据范围，为高精度的岩土层界面 问题求解.湖南大学学报：自然科学版，2009,36(10)：81) 拟合提供有效数据. [10]Houlding S W.3D Geosciences Modeling:Computer Technique for (2)利用基于Closest Point算法的结构分析模 Geological Characterization.Berlin:Springer-Verlag,1994 型进行了Kriging插值计算，实现了Kriging算法和 [11]Lemon A M,Jones N L.Building solid models from boreholes and user defined eross section.Comput Geosci,2003,29(5):547 Closest Point算法的融合，完成了对东露天矿边坡各 02] Liu Z P,He H J,Zhu F H.Study of technology of fast 3D model- 岩土层的三维拟合.经测量验证，与实际情况基本 ing and visualization based on borehole data.Rock Soil Mech, 相符，可为矿区后续生产及边坡防护设计提供有效 2009,30(Suppl1):260 的数据支撑 (刘振平，贺怀建，朱发华.基于钻孔数据的三维可视化快速 (3)基于Kriging与Closest Point融合算法的界 建模技术的研究.岩土力学，2009,30（增刊1）：260) [13]Hou J R,Yin Z N,Li W M,et al.Practical Geological Statistics 面拟合方法具有较高的精度，且简单易行，可为类似 Beijing:Geological Publishing House,1998 矿山的岩土层拟合提供参考. (侯景儒，尹镇南，立维明，等.实用地质统计学.北京：地质 出版社，1998) 参考文献 [14]Liu Y J,Xie X G,Chen S D.Improved approximation algorithm [1]Krige D G.A statistical approach to some basic mine valuation for TSP.Comput Eng Appl,2006,42(33):71 problems on the Witwatersrand.J Chem Metall Min Soc South Afr, (刘艳娟，谢晓钢，陈胜达.一种改进的求解TSP问题的近似 1951,52(6):119 算法.计算机工程与应用，2006,42(33)：71)

第 5 期 李 健等: 基于 Kriging 与 Closest Point 融合算法的边坡岩土层界面拟合 图 4 粉质黏土层等厚线图 Fig． 4 Isopachous map of the silty clay layer 工程现场出现小规模塌方，见图 5． 图 5 南侧重点关注区现场图 Fig． 5 Picture of the southern area which we pay close attention to 与两次滑坡不同的是，开挖本区域时，当地的季 风性气候已进入旱季，因而水对边坡影响较小，加之 施工过程中，生产单位高度重视，对本区域的施工情 况重点关注，及时采取了相应安全措施，因而未发生 大规模滑坡． 4 结论 ( 1) 充分利用岩土层界面的地质调查点数据， 可有效扩充钻孔数据范围，为高精度的岩土层界面 拟合提供有效数据． ( 2) 利用基于 Closest Point 算法的结构分析模 型进行了 Kriging 插值计算，实现了 Kriging 算法和 Closest Point 算法的融合，完成了对东露天矿边坡各 岩土层的三维拟合． 经测量验证，与实际情况基本 相符，可为矿区后续生产及边坡防护设计提供有效 的数据支撑． ( 3) 基于 Kriging 与 Closest Point 融合算法的界 面拟合方法具有较高的精度，且简单易行，可为类似 矿山的岩土层拟合提供参考． 参 考 文 献 ［1］ Krige D G． A statistical approach to some basic mine valuation problems on the Witwatersrand． J Chem Metall Min Soc South Afr， 1951，52( 6) : 119 ［2］ Schuenemeyer J H，Power H C． Uncertainty estimation for resource assessment: an application to coal． Math Geol，2000，32( 5) : 521 ［3］ Li X J，Hu J H，Zhu H H，et al． The estimation of coal thickness based on Kriging technique and 3D coal seam modeling． J China Coal Soc，2008，33( 7) : 765 ( 李晓军，胡金虎，朱合华，等． 基于 Kriging 方法的煤层厚度估 计及三维煤层建模． 煤炭学报，2008，33( 7) : 765) ［4］ Jin G D，Liu Y C，Niu W J． Comparison between inverse distance weighting method and Kriging． J Changchun Univ Technol，2003， 24( 3) : 53 ( 靳国栋，刘衍聪，牛文杰． 距离加权反比插值法和克里金插值 法的比较． 长春工业大学学报，2003，24( 3) : 53) ［5］ Zou G D，Chen Y P． Coupling algorithm of simulated annealing al￾gorithm and random search method for slope stability analysis． Chin J Rock Mech Eng，2004，23( 12) ，2032 ( 邹广电，陈永平． 滑坡和边坡稳定性分析的模拟退火--随机搜 索耦合算法． 岩石力学与工程学报，2004，23( 12) : 2032) ［6］ Shi L，Li X C，Ren W，et al． Hybrid of ant colony algorithm and ge￾netic algorithm and its application to searching critical slope slip surface． Rock Soil Mech，2009，30( 11) : 3486 ( 石露，李小春，任伟，等． 蚁群算法与遗传算法融合及其在边 坡临界滑动面搜索中的应用． 岩土力学，2009，30( 11) : 3486) ［7］ Liu Y，Li Y，Xu J T． Application of a multi-algorithm fusion real time filter in FOG． Acta Photonica Sin，2010，39( 6) : 1116 ( 刘颖，李言，徐金涛． 一种多算法融合的实时滤波在光纤陀螺 中的应用． 光子学报，2010，39( 6) : 1116) ［8］ Hu J J，Tang C J，Li C，et al． Efficient multi-layer clustering algo￾rithm based on nearest neighbors first． J Sichuan Univ Eng Sci， 2004，36( 6) : 93 ( 胡建军，唐常杰，李川，等． 基于最近邻优先的高效聚类算法． 四川大学学报: 工程科学版，2004，36( 6) : 93) ［9］ Wu J H，Zhang J，Liu Z H． Solution of TSP problem based on the combination of ant colony algorithm and immune algorithm． J Hu￾nan Univ Nat Sci，2009，36( 10) : 81 ( 吴建辉，章兢，刘朝华． 基于蚁群算法和免疫算法融合的 TSP 问题求解． 湖南大学学报: 自然科学版，2009，36( 10) : 81) ［10］ Houlding S W． 3D Geosciences Modeling: Computer Technique for Geological Characterization． Berlin: Springer-Verlag，1994 ［11］ Lemon A M，Jones N L． Building solid models from boreholes and user defined cross section． Comput Geosci，2003，29( 5) : 547 ［12］ Liu Z P，He H J，Zhu F H． Study of technology of fast 3D model￾ing and visualization based on borehole data． Rock Soil Mech， 2009，30( Suppl 1) : 260 ( 刘振平，贺怀建，朱发华． 基于钻孔数据的三维可视化快速 建模技术的研究． 岩土力学，2009，30( 增刊 1) : 260) ［13］ Hou J R，Yin Z N，Li W M，et al． Practical Geological Statistics． Beijing: Geological Publishing House，1998 ( 侯景儒，尹镇南，立维明，等． 实用地质统计学． 北京: 地质 出版社，1998) ［14］ Liu Y J，Xie X G，Chen S D． Improved approximation algorithm for TSP． Comput Eng Appl，2006，42( 33) : 71 ( 刘艳娟，谢晓钢，陈胜达． 一种改进的求解 TSP 问题的近似 算法． 计算机工程与应用，2006，42( 33) : 71) ·505·