变异函数自动拟合初探

D0I:10.13374/j.issn1001053x.1984.02.012 北京钢铁学院学报 1984年 第2期 变异函数自动拟合初探 数学教研室 戚国安 地质教研室 吴雨沛 理化系 陈廷琨 摘 要 克立格法是估计矿石储量的一种较新方法，目前使用该方法时，求实验变异函数及作克立格估值均可在电 子计算机上自动进行，唯由实验变异函数拟合理论模型时还要依靠手工，它坊碍了由原始数据到储量估算整个 过程的自动化。本文从变异函数的统计性质出发，探讨了变异函数的自动拟合问题，得到了一些初步的结果。 克立格法是目前世界上开始广泛使用的一种计算矿石储量的较新方法。这个方法的大 致步骤是：根据地质资料求出实验变异函数，然后用一定的理论变异函数模型去拟合它， 再利用理论模型对矿石储量作克立格估计。求实验变异函数及进行克立格估计均已有电算 程序，可以在电子计算机上自动进行。而用理论模型拟合实验变异函数时，还主要是靠手工 进行，由于手工拟合缺乏统一的标准，结果往往因人而异，有时差别还比较大。而且拟合时需 要反复试配，比较耗工费时。更主要的是手工拟合不利于在电子计算机上自动进行，妨碍了 从原始数据到克立格估计整个过程的自动化。针对这个问题，本文探讨了区域化变量的一些 统计性质，利用这些性质对理论相关系数作区间估计，从而确定了变程的范围。再利用实验 变异函数的区间估计对理论模型进行筛选，最后用剩余标准离差确定理论变异函数。整个过 程可以在电子计算机上进行，这样便有可能实现从原始数据到储量计算全过程的自动化。 本文中对区域化变量Z(x)作如下假定： (1)平稳性假设。因而均值与方差都是常数，E[Z(x)]=m及DCZ(x)]=σ2.而协 方差C(h)=Cov[Z(x),Z(x+h)]是h的函数。 (2)内蕴假设。根据平稳性假设及关系式：Y(h)=g2-C(h)可知：变异函数Y(h) =合E([Z(x+h)-Z(x)]:}只是h的函数。 (3)正态性假设。即设Z(x)与Z(x+h)服从二维正态分布。 (4)襁立增量性假设。即设增量Z(x+h)-Z(x)是独立的。 一、变程a的区间估计 在本文的假设下，区域化变量Z(x)的方差σ2，协方差c(h)与变异函数Y(h)之间有下 列关系： c(h)=g2-Y(h). 125

北 京 钢 铁 学 院 学 报 啤 第 · 期 ‘ 变异函数 自动拟合初探 数学教研室 戚国安 地质教研 室 吴雨沛 理 化 系 陈廷混 、 夕 、 、 摘 要 克立 格法是估计矿 石储盆的一 种较新方法 ， 目前 使用该方法时 ，求实验变异函数及作克立格估值均可在电 子计算机上 自动进行 ， 唯由实验变异函数拟合理论棋型时还要依靠手工 ，它妨碍了 由原始数据到储 估算盆个 过程的 自动化 。 木文从变异函数的统计性质出发 ， 探讨了变异函致的 自动拟合间题 ， 得到 了一些初步的结果 。 克立格法是 目前世界上开始广泛使用的一种计算矿石储量的较新方法 。 这个方法的大 致步骤是 根据地质资料求出实验变异 函数 ， 然后用一定的理论变异函数模型去拟合它 ， 再利用理论模型对矿石储量作克立格估计 。 求实验变异 函数及进行克立格估计均 已有电算 程序 ， 可 以在 电子计算机上 自动进行 。 而用理 论模型拟 合实验变异 函数时 ，还主要是靠手工 进行 ， 由于手工拟 合缺乏统一 的标准 ， 结果往往因人而异 ， 有时差别还 比较大 。 而且拟合 时需 要反复试配 ， 比较耗工费时 。 更主要 的是手工拟合不利于在 电子计算机上 自动进行 ， 妨碍了 从 原始数据 到克立 格估计整个过程 的 自动 化 。 针对这个间题 ， 本文探讨了 区域化变量的一些 统计性质 ， 利用这 些性质对理论相关系数作 区 间估计 ， 从而确定了变程 的范 围 。 再 利用 实验 变异 函数的区 间估计对理论模型进行筛选 ， 最 后用 剩余标准离差 确定理 论变异 函数 。 整个过 程可 以在 电子计算机上进行 ， 这样便有可能实现 从原始数据到储量计算全 过程的 自动 化 。 本文 中对 区域化变量 作如下假定 平稳性假设 。 因而均值与方差都是常数 〔 〕 及 〔 ， 而协 方差 【 ， 」是 的 函数 。 幻 内蕴假设 。 根据平稳性假设及关系式 丫 二 一 ‘ 可知 变异函数 一 合 ‘ 仁 ‘ ，一 ， 只 是 的函数 。 正态性假设 。 即设 与 服从二维正态分 布 。 七 ‘ 独立增量 性假设 。 即设增量 一 是独立 的 。 一 、 变程 的区间估计 在本文 的假设下 ， 区域化变量 的方差 ， 协方差 与变异 函数 之 间有下 列关系 一 DOI ：10．13374／j ．issn1001－053x．1984．02．012

从而理论相关系数为 e(h)=c(h)=1-Y(h) 02 由数理统计知道，在正态性假设下，样本相关系数 是EZ(x1+h)-Zx+五)][Z(x,)-Zx] Y(h)= ial [Z(x1+h)-Z(x+h五)]2√2[Z(x)-Z(x) i=1 i=1 的概率密度函数是 0-1 0-4 f0y)=2a-31-p2)2g1-y)7 (n-3)1 (）(-2+）. a=0 其中p是理论相关系数，而 Zx+h)=上2Z(x1+h), n i=1 Zx)=122(x. n i=1 而当n较大时， μ=(th-1Y-th-lp)√n-3 渐近地服从标准正态分布N(01)。于是对于给定的a,设对应的双侧分位数为u1:,则 理论相关系数P的置信区间为： h(--是<h仙y+是. (1) 因为p还是变程a的函数，例如对于球状模型 p1-费+ (0≤h≤a). 将样本相关系数Y及理论相关系数p的表达式代入(1)式，对于不同的h,解不等式(1)， 便可对变程a作区间估计，从而求出a的范围：a1≤a≤az。 二、理论变异函数的确定 现在以球状模型为例，说明确定理论变异函数的方法。众所周知，球状模型的理论变 异函数为 Y)-C+C(驶-六)，0≤h<a. (2) 令 bg=co,b1=,b-2点. 126

从而理论相关系数为 口 ， 五 一 一 一 尸万一 由数理统计知道 ， 在正态性假设下 ， 样本相关系数 三〔 一 〕 一 ， 〕 刃 ， 一 。 〕 忿 了刃 ， 一 ， 〕 的概率密度 函数是 一 一 丫 一 一 一 丫 留 一 兀 “ 二 、 。 。 一 、 气、 花二一 夕 ‘ 气一一一石一一 认 ， 、 ‘ 其中 是理 论相关系 数 ， 而 二 ， 、 冬 ， ， 、 乙 气 十 － ‘ 乙 气盆 十 声， 二 布 ， 、 乙 口 一 艺 ， 。 而 当 较大时 ， 卜 一 ‘ 丫一 一 ‘ 召云二厄 渐近 地服从标准正态分布 。 · 。 于是对于给定的 ， 设对应的双侧分位数 为 ， 则 理论相关系数 的置信区 间为 一 一 一 亿五二落 一 ，丫 谬里斗一 石 因为 还是变程 的 函数 ， 例如对于球状棋童 ， 一 、 一 一二一 十 二一万 气 飞 气 ， 乙 乙 ， 将样本相关系数丫及理论相关系数 的表达式代入 式 ， 对 于不 同的 ， 解不等式 ， 便可对变程 作 区间估计 ， 从而求出 的范 围 ，簇 《 。 二 、 理论变异函数的确定 现在 以球状模型为例 ， 说明确定理论变异 函数的方法 。 众所周知 ， 球状模型的理论变 异函数为 ， 。 、 李乙 卫 、 歹 ， 一 《 镇 。 一 不下一 ‘ 典

或当b1与b:异号时，有 co=b,a=√-30.，c=gV-61 3V3b. (3) 于是(2)式变成多项式： Y(h)=b。+b1h+b2h3. (4) 用最小二乘法求出(4)之后，再由(3)求出参数c。,a及c。 设实验变异函数(h)的计算结果为： h1 h2 hk ha (h1)(hz)… f(h)…Y(h) 若h在变程a的置信区间内：a1≤hk≤a2,且k≥4，则用最小二乘法求(4)时，只要粮取 下列数据，即 h h2…hk f(h)(h2)…(h) 就可以了。这里要求k≥4是因为这样的多项式回归要求的数据对不得少于4。但是满足上 述关系的k往往不止一个，因此求得的理论变异函数也有若干个。这时再按下面的方法选 出最优的理论变异函数。 三、实验变异函数的统计性质 定理设 1N(h) ↑(h)=2N(h)1=1 [Z(x1+h)-Z(x)]2 为实验变异函数，而Y(h)是理论变异函数，则 N(h)分(h)~x2[N(h)门。 Y(h) 证由本文所作的关于区域化变量的假设及正态分布的性质知Z(x+h)-Z(x)服从一 维正态分布。且其均值与方差分别为 E[Z(x+h)-Z(x)]=m-m=0. D[Z(x+h)-Z(x)]=E{[Z(x+h)-Z(x)]2} =E{([Z(x+h)-m]-[Z(x)-m])} =E{CZ(x+h)-m]2}+E{[Z(x)-m]2}-2E{E(x+h)-m] .[Z(x)-m]} =2[g2-c(h)] =2Y(h)。 于是 Z(x+h)-Z(x2N(0.1)。 V2Y(h) 又由独立性假设知道 127

或当 与 异号时 ， 有 了一飞了 ， ‘ 。 。 。 一 产万下万 舟二 勺 一 丝竺 ，一 劣 一 一 “ 。 ’ “ 一 丫 ’ “ 一飞厂 于是 式变成多项式 丫 。 。 用最小二乘法求出 之 后 ， 再 由 求 出参数 。 ， 及 设实验变异 函数今 的计算结果为 … 、 … 宁 宁 … 宁 、 … 宁 若 ，在变程 的置信区 间内 《 、 镇 ， 且 》 ， 则用最小二乘法求 下列数据 ， 即 一 … ， 令 宁 … 今 ， 《 时 ， 只 要截取 卜， 就可 以 了 。 这里要求 是 因为这样 的多项式 回归要求的数据对不得少 于 。 但是满足上 述关系 的 往往不止一 个 ， 因此求得 的理论变异 函数也有若千个 。 这 时再按下 面的方法选 出最优的理论变异 函数 。 三 、 实验变异函数的统计性质 定理 设 了、 ， 、 ，屯， 尸 ， ， 。 ， 、 尸了 ， 、 ， ‘ ” ’ 气厄瓦 乃 雪 “ 乙 ‘ 不 ‘ 十 “ ’ 一 ‘ “ ” “ 为实验变异 函数 ， 而 是理论变异函数 ， 则 兮 二 ， 〕 。 证 由本文所作 的关于 区域化变量的假设及正态分 布的性质知 一 服从一 维正态分布 。 且其均值与方差分别为 〔 一 一 〔 一 〕 一 〕 盆 〔 一 一 〔 一 〕 吕 〔 一 〕 名 〔 一 〕 一 一 〕 〔 一 〕 〔 一 〕 。 于是 一 一万芬而硒 闷 。 。 叉由独立性假设知道

黑T Nrz(x+h)-Z(xx(N(h)). 根据这个定理，对于给定的a(0<a<1),可求得对应的x2分布的双侧分位数 Xt-a2(N(h)及X2r-a/2(N(h)。于是得到实验变异函数的置信区间为： (N()()<X (N()) N(h) (6) 将前面求得的各个理论变异函数及实验变异函数代入(6)式逐个进行检验。将其中 不满足(6)：的点数最少的保留下来，然后求出其中利余标准离差 S=√12[Y(h)-(hD] (7) n-31=1 最小的理论变异函数，即为最优理论变异函数。其中n为使h:≤a≤h:+1成立的正整数。 四、实 例 下面表1中列出某矿体的样品间距h,品位Z(),实验变异函数？(h),样本相关系数 Y(h),样本相关系数Y(h)的下限Y(h)与上限Y*(h): ,)=轴[-y)-2g], Y*(h)=tht iv(h) 以及变程a的下限a1及上限a2。由表1可见变程a的范围是1.74≤a≤9.27。因此在h的下列 范围内截取数据作多项式回归，即[1,4幻，[1,5门，，[1,9]。共得到6个回归 方程。经过电算检验，实验变异函数均在由它们所确定的置信区间即(6)式之内。 表2中列出了这6个理论模型的参数a,c,及si1和剩余标准离差s。其中s最小的各参 数与手工拟合的参数对比如下： 参 数 a Co sill 手工拟合 6.5 0.003 0.040924 自动拟合 6.97066 0.005009 0.041339 由此可见结果是比较接近的，最后用电子计算机绘出了实验变异函数，手工拟合的理论变 异函数及自动拟合的理论变异函数图形，从图形看出，结果还是令人满意的。 我们虽然只列举了球状模型一级套合的例子，但所阐述的基本理论也适用于球状模型 的多级套合、或其它类型模型的自动拟合。 陈希康副教授对本文提供了宝贵意见，并得到曹诺农同志的帮助。在电算过程中，得 到许刚，尚小卫同志的大力协助，在此表示衷心的感谢。 我们仅对变异函数自动拟合问题进行了初步探讨，由于作者水平有限，不当之处恳请 128

， ， 、 ， ， 、 一 ， 望墨粉 三 丝聋号端箭 生丛 一” ‘“ ‘“ ” · 根据这个定理 ， 对 于给定的 ， 可求得对 应 的 名 分 布 的 双 侧 分 位 数 渺 卜 及 ’ 卜 。 于是得到实验变异 函数的置信区 间为 一 “ 蹲噢之 宁 、 气 少 。 ， 犀上、华气 少 将前面求得的各个理论变异函数及实脸变异 函数代入 式逐个进行检脸 。 将 其 中 不满足 的 点数最少的保留下来 ， 然后求 出其中剩余标准离差 了击 套 〔“ ，，一 今‘ ，，” 最小 的理论变异 函数 ， 即为最优理论变异 函数 。 其 中 为使 《 五 不 成立的正整数 。 四 、 实 例 下 面表 中列 出某矿体 的样品间距 ， 品位 ， 实验变异 函数净 ， 样本相关系数 ， 样本相关系数 的下 限 与上限丫 丫· ， 【 ， 一 ‘ 丫 丫· 【 一 丫 岁景清 ， 步器了 〕 · 以及变程 的下 限 及上限 。 由表 可见变程 的范 围是 《 《 。 。 因此在 的下列 范围内截取数据作多项式 回归 ， 即 〔 ， 〕 ， ， 〕 ， … ， ， 〕 。 共得到 个回归 方程 。 经过 电算检验 ， 实验变异 函数均在 由它们所确定的置信区 间即 式之内 。 表 中列 出了这 个理论模型 的参数 ， 。 及 和剩余标准离差 。 其中 最小的各参 数与手工拟合的参数对比如下 手 工 拟 合 自 动 拟 合 。 。 。 。 。 。 由此可见结果是 比较接近 的 ， 最后用 电子计算机绘出了实验变异 函数 ， 手工拟合的理论变 异 函数及 自动拟合的理论变异 函数图形 ， 从 图形看出 ， 结果还是令人满意的 。 我们 虽然只 列举 了球状模型一级套合的例子 ， 但所 阐述 的基本理论也适用于球状模型 的多级套合 、 或其它类型模型 的 自动拟合 。 陈希廉副教授对本文提供 了宝 贵意见 ， 并得到曹诺农 同志 的带助 。 在 电算过程 中 ， 得 到许刚 ， 尚小卫 同志 的大力协助 ， 在此表示衷心 的感谢 。 我们 仅对变异 函 数 自动 拟 合 问题进行 了 初步探讨 由于作者水平有 限 ， 不 当之处恳请 吕

批评指正。 Z(h) (h) r(h) r.(h) r◆(h) &s 0.2000 0.012051 0.628879 0.235080 0.845330 19 1.749.66 2 0.2664 0.021462 0.305584 0.198077 0.681583 18 09.27 3 0.5320 0.026159 0.160264 -0.356495 0.601938 17 0 11 0.5320 0.034126 0.066835 -0.452441 0.552265 16 0 12 5 0.5320 0.036802 0.113007 -0.433224 0.598523 15 0 18 6 0.5320 0.039119 0.139847 -0.431919 0.631431 14 0 >19 7 0.5320 0.042159 0.211712 -0.394816 0.689718 13 0 >19 8 0.5320 0.038988 0.031439 -0.561634 0.603171 12 0 >19 9 0.5320 0.036871 0.245850 -0.426913 0.743427 11 0 >19 10 0.3728 0.042204 0.356219 -0.365639 0.810498 10 0 >19 11 0.1340 0.044314 0.297337 -0.469868 0.808642 9 0 >19 12 0.1340 0.039953 0.347439 -0.486831 0.850223 8 0 >19 13 0.1340 0.039605 0.169774 -0.679705 0.824731 7 0 >19 14 0.1340 0.043304 0.190272 -0.745204 0.873419 6 0 >19 15 0.2726 0.049180 0.367238 -0.773502 0.946748 5 0 >19 16 0.3320 0.067132 0.800822 -0.715859 0.995957 0 >19 17 0.3320 0.057046 0.496726 18 0.4925 0.016814 19 0.0005 0.004930 20 0.1007 表 2 a c0 sill 8.07892 0.004599 0.047419 0.007149 5.98365 0.003739 0.038344 0.006658 6.18853 0.003966 0.039055 0.006094 6.97066 0.005009 0.041339 0.005903 6.82183 0.004762 0.040987 0.005911 6.97950 0.005156 0.041218 0.005916 129

批评指正 裹 一八 一工 一一计峨 ‘ 人孟土上上，，上，上， ” ， “ “ ， ·‘ ， 二丝 曰内 月九 白吕一 。 厅‘ ，月 二甘上几‘， 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 一 。 一 。 一 。 一 。 一 。 一 。 一 。 一 。 一 。 一 。 一 。 一 。 一 。 一 。 一 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 勺自 丹︸甘‘ ﹄， 上几二，，怪勺甘口，︸丹几︸‘月扩 自几，，几，占，几上，‘， ‘目尸 衰 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。