一种Kriging时空统一预测模型及其构建方法和应用与流程

文档序号：16212061发布日期：2018-12-08 07:52阅读：618来源：国知局

本发明属于特别适用于特定功能的数字计算设备或数据处理设备或数据处理方法技术领域，尤其涉及一种kriging时空统一预测模型及其构建方法和应用。

背景技术

目前，业内常用的现有技术是这样的：隧道围岩位移预测预报是围岩稳定性研究中最重要的内容，尽管国内外学者在这方面已经积累了一些研究成果，但对位移预测的研究仍然是围岩变形研究领域中最薄弱的一环。隧道围岩位移预测预报对施工决策有特别重要的意义。在实际监控量测过程中，传统检测方法容易受到外界因素的干扰，再加上仪器设备系统误差和技术人员本身读数误差的不可消除性，使得量测数据很难接近围岩应力应变的真实情况。在这种情况之下，必须采用有效的数据分析方法对离散性较大的量测数据进行排查，以确保数据的真实性，以期实现较为准确的预测，避免对围岩和支护结构的稳定性产生误判，影响隧道施工安全。

综上所述，现有技术存在的问题是：当前，国内对隧道围岩变形预测的研究，多是通过位移反演进行岩体力学参数敏感性分析及数值模型优化，或是分析位移的时序变化趋势来预测其稳定收敛值，很难真正实现对围岩变形的有效预测。隧道围岩变形规律本身具有非线性特征，存在很多不确定的随机因素，以致很多理论公式的推导和解释都很难有效的达到准确预测的效果。在隧道实际施工中，往往忽略了隧道围岩的物理力学性质与现场监控量测信息之间的内在关联，之于围岩空间变异性对隧道收敛变形的影响没有引起足够的重视，而这恰恰是对后续未开挖段施工、设计变更与否等具有重要指导意义。

解决上述技术问题的难度和意义：

kriging插值法是一种运用结构分析与变异函数相关理论，在有限空间内针对区域化变量进行最优和无偏估值计算的预测方法。相比于上述其它传统预测方法，kriging插值能更好的描述描述空间的连续关联性并反映数据点的整体变化趋势，预测更精确，更优越。尽管kriging插值法在地质勘探和地质统计领域应用较为广泛，但在岩土工程领域还未引起国内外专家学者们足够的关注和研究。

技术实现要素：

针对现有技术存在的问题，本发明提供了一种kriging时空统一预测模型及其构建方法和应用。

本发明是这样实现的，一种kriging时空统一预测模型，所述kriging时空统一预测模型为：

式中：z^*(s0,t0)为空间点(s0,t0)处的估计值，z(si,tj)为时间点(si,tj)处的估计值j＝1,2,…m，λi为临近观测值z(si,tj)的加权系数，i是空间点数(i＝1,2,…n)，j是空间点数()；引入拉格朗日系数μ进得：

式中γ(ab,xy)＝γ(a-b,x-y)是时空点对(sa-sb,tx-ty)的变异函数；μ是拉格朗日乘数。

本发明的另一目的在于提供一种所述kriging时空统一预测模型的构建方法，所述kriging时空统一预测模型的构建方法包括：

步骤一，计算测点之间的距离；

步骤二，按照距离大小排序；实现距离分组；

步骤三，确定变异函数，计算变差，组成系数矩阵；

步骤四，计算权值系数，计算预测点的估计值，所有测点完成则结束；所有测点未完成则返回步骤三。

进一步，所述kriging时空统一预测模型的构建方法具体包括：

步骤一，选取kriging插值；

步骤二：时空变异函数模型创建：函数模型：

cst(hs,ht)＝k1[cs(0)-γs(hs)][ct(0)-γt(ht)]

+k2(cs(0)-γs(hs))+k3(ct(0)-γt(ht))；

步骤三：时空普通kriging预测：预测公式如下：

式中：z^*(s0,t0)为空间点处的估计值，z(si,tj)(s0,t0)为时间点(si,tj)处的估计值j＝1,2,…m，λi为临近观测值z(si,tj)的加权系数，i是空间点数(i＝1,2,…n)，j是空间点数()；

步骤四：kriging时空插值模型建模。

进一步，所述步骤四时空插值包括：

计算时空实验变异函数，估计时空变异模型参数；

确定最优的时空变异函数及相关参数；

计算出k1、k2、k3，并计算出时空变异函数模型；

计算预测值。

进一步，所述时空变异函数模型的创建方法为：空间变异函数只能描述在某一固定时刻的不同监测点的空间变异结构，与空间协方差函数cs(hs)的关系式为：

γs(hs)＝cs(0)-cs(hs)；

与时间协方差函数ct(ht)的关系式为：

γt(ht)＝ct(0)-ct(ht)；

与时空协方差函数的关系式为：

γst(hs,ht)＝cst(0,0)-cst(hs,ht)；

函数模型如下：

当h＝0，且t＝0时，由“product–sum”时空变异函数模型得：

cst(0,0)＝k1cs(0)ct(0)+k2cs(0)+k3ct(0)；

得：

γst(hs,ht)＝(k1cs(0)+k3)γt(ht)+(k1ct(0)+k2)γs(hs)-k1γs(hs)γt(ht)；

其中，cs和ct分别为空间和时间所对应的协方差，cs(0)和ct(0)分别为空间和时间所对应的变异函数基台值，cst(0,0)为时空变异函数γst的基台值；

参数k1、k2、k3的表达式如下：

k1＝[cs(0)+ct(0)-cst(0,0)]/cs(0)ct(0)；

k2＝[cst(0,0)-ct(0)]/cs(0)；

k3＝[cst(0,0)-cs(0)]/ct(0)。

本发明的另一目的在于提供一种应用所述kriging时空统一预测模型的隧道围岩位移预测预报系统。

综上所述，本发明的优点及积极效果为：检测方法不会受到外界因素的干扰；相比于上述其它传统预测方法，kriging插值能更好的描述描述空间的连续关联性并反映数据点的整体变化趋势，预测更精确，更优越。通过交叉验证和相关精度指标对时空kriging预测模型、灰色verhulst预测模型两种模型进行对比评定，验证结果表明：相对于灰色verhulst预测模型，时空kriging模型预测精度更高，能更准确地预测拱顶下沉位移的发展变化趋势，在实际隧道工程的围岩位移预测中是有效、可行的，同时也为隧道工程监测预报技术提供了一种全新的研究手段。(具体见下表13～表15)

附图说明

图1是本发明实施例提供的kriging时空统一预测模型的构建方法流程图；

图2是本发明实施例提供的kriging时空统一预测模型的构建方法实现流程图；

图3是本发明实施例提供的kriging时空统一预测模型及构建方法周边位移及拱顶下沉测点布置图；

图4是本发明实施例提供的kriging时空统一预测模型及构建方法yk31+872拱顶沉降实测值与灰色预测值对比图；

图5是本发明实施例提供的kriging时空统一预测模型及构建方法yk31+875拱顶沉降实测值与灰色预测值对比图；

图6是本发明实施例提供的kriging时空统一预测模型及构建方法yk31+878拱顶沉降实测值与灰色预测值对比图；

图7是本发明实施例提供的kriging时空统一预测模型及构建方法yk31+872拱顶沉降实测值与kriging预测值对比图；

图8是本发明实施例提供的kriging时空统一预测模型及构建方法yk31+875拱顶沉降实测值与kriging预测值对比图；

图9是本发明实施例提供的kriging时空统一预测模型及构建方法yk31+878拱顶沉降实测值与kriging预测值对比图；

图10是本发明实施例提供的kriging时空统一预测模型及构建方法yk31+872拱顶沉降kriging预测值日沉降值误差图；

图11是本发明实施例提供的kriging时空统一预测模型及构建方法yk31+872拱顶沉降kriging预测值累积沉降值误差图；

图12是本发明实施例提供的kriging时空统一预测模型及构建方法yk31+875拱顶沉降kriging预测值日沉降值误图；

图13是本发明实施例提供的kriging时空统一预测模型及构建方法yk31+875拱顶沉降kriging预测值累积沉降值误差图；

图14是本发明实施例提供的kriging时空统一预测模型及构建方法yk31+878拱顶沉降kriging预测值日沉降值误差图；

图15是本发明实施例提供的kriging时空统一预测模型及构建方法yk31+878拱顶沉降kriging预测值累积沉降值误差图；

图16是本发明实施例提供的kriging时空统一预测模型及构建方法yk31+878拱顶沉降实测值与kriging预测值对比图；

图17是本发明实施例提供的kriging时空统一预测模型及构建方法yk31+878拱顶沉降kriging预测值日沉降值误差图；

图18是本发明实施例提供的kriging时空统一预测模型及构建方法yk31+878拱顶沉降kriging预测值累积沉降值误差图；

图19是本发明实施例提供的kriging时空统一预测模型及构建方法yk31+881拱顶沉降实测值与kriging预测值对比图；

图20是本发明实施例提供的kriging时空统一预测模型及构建方法yk31+881拱顶沉降kriging预测值日沉降值误差图；

图21是本发明实施例提供的kriging时空统一预测模型及构建方法yk31+881拱顶沉降kriging预测值累积沉降值误差图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

本发明将经典的kriging空间插值进行时空扩展，应用时空变异函数建立了隧道围岩位移预测模型；分析了参数与变形的之间的关系，给出了时空插值的计算程序和方法。

下面结合附图对本发明的应用原理作详细的描述。

本发明实施例提供的kriging时空统一预测模型为：

式中：z^*(s0,t0)为空间点(s0,t0)处的估计值，z(si,tj)为时间点(si,tj)处的估计值j＝1,2,…m，λi为临近观测值z(si,tj)的加权系数，i是空间点数(i＝1,2,…n)，j是空间点数()；引入拉格朗日系数μ进行推导可得：

式中γ(ab,xy)＝γ(a-b,x-y)是时空点对(sa-sb,tx-ty)的变异函数；μ是拉格朗日乘数，可得研究区域内任意点的插值估计。

如图1所示，本发明实施例提供的kriging时空统一预测模型的构建方法包括以下步骤：

s101：计算测点之间的距离；

s102：按照距离大小排序；实现距离分组；

s103：确定变异函数，计算变差，组成系数矩阵；

s104：计算权值系数，计算预测点的估计值，所有测点完成则结束；所有测点未完成则返回s103。

下面结合附图对本发明的应用原理作进一步的描述。

如图3所示，本发明实施例提供的kriging时空统一预测模型构建方法包括以下步骤：

步骤一：选取kriging插值方法：建模时主要关注时空自相关结构，以确保所选择的模型是有效的，模型是足够灵活，允许拟合数据并估计模型参数，普通kriging法相关假设在隧道监测数据分析预测中具有较好的适用性。

步骤二：时空变异函数模型创建：不可分离型中的“product–sum”时空变异函数模型较为符合地质学实际情况，得其函数模型如下：

cst(hs,ht)＝k1[cs(0)-γs(hs)][ct(0)-γt(ht)]

+k2(cs(0)-γs(hs))+k3(ct(0)-γt(ht))；

步骤三：时空普通kriging预测：普通kriging方法实现数据的时空插值，其预测公式如下：

步骤四：kriging时空插值模型建模思想：为了简化时空kriging插值预测模型计算过程，假设：(1)隧道变形监测数据是隧道围岩形变规律的宏观体现，其变形分析的预测值仅对围岩形变规律进行深度数据挖掘；(2)隧道围岩在动态施工过程中最终产生的应力应变是受地质体内外界条件与各种因素的相互影响、综合作用的结果；(3)隧道围岩各影响因素间相互作用、相互影响且呈现非线性规律，但能够组建隧道变形分析预测模型对其进行拟合处理。

时空插值的基本步骤如下：

(5)计算时空实验变异函数，估计时空变异模型参数；

(6)确定最优的时空变异函数及相关参数；

(7)计算出k1、k2、k3，并计算出时空变异函数模型；

(8)计算预测值。

本发明实施例提供的所述时空变异函数模型的创建方法为：空间变异函数只能描述在某一固定时刻的不同监测点的空间变异结构，其与空间协方差函数cs(hs)的关系式为：

γs(hs)＝cs(0)-cs(hs)(1)

与空间变异函数比较，时间变异函数则只能表示某一固定的监测点在不同时刻或者阶段的时间变异结构，其与时间协方差函数ct(ht)的关系式为：

γt(ht)＝ct(0)-ct(ht)(2)

时空变异函数用来表征时空点的变异结构，其与时空协方差函数的关系式为：

γst(hs,ht)＝cst(0,0)-cst(hs,ht)(3)

为了使上述函数转化为可计算形式，引入“product–sum”时空变异函数模型。经过大量研究发现，不可分离型中的“product–sum”时空变异函数模型较为符合地质学实际情况，联立式(1)、式(2)可得其函数模型如下：

当h＝0，且t＝0时，由“product–sum”时空变异函数模型的定义可得：

cst(0,0)＝k1cs(0)ct(0)+k2cs(0)+k3ct(0)(5)

联立方程式(3)可得：

γst(hs,ht)＝(k1cs(0)+k3)γt(ht)+(k1ct(0)+k2)γs(hs)-k1γs(hs)γt(ht)(6)

其中，cs和ct分别为空间和时间所对应的协方差，同时，cs(0)和ct(0)分别为空间和时间所对应的变异函数基台值，cst(0,0)为时空变异函数γst的基台值。a、c0、c1等参数并不是一成不变的，他们的具体值在建模计算过程中都是可以进行调节的，每个参数都具有其最优取值。以上诸多参数值的限值定义以确保时空变差函数的有效性。

若γst(0,0)＝γs(0)＝γt(0)＝0，则由上式式(6)可知：

γst(hs,0)＝[k2+k1ct(0)]γs(hs)(7)

γst(0,ht)＝[k3+k1cs(0)]γt(ht)(8)

又假设：

k2+k1ct(0)＝1；k3+k1cs(0)＝1(9)

由式(5)和式(9)可知参数k1、k2、k3的表达式如下：

k1＝[cs(0)+ct(0)-cst(0,0)]/cs(0)ct(0)(10)

k2＝[cst(0,0)-ct(0)]/cs(0)(11)

k3＝[cst(0,0)-cs(0)]/ct(0)(12)

本发明实施例提供的所述kriging方法预测公式中引入拉格朗日系数μ进行推导可得：

式中γ(ab,xy)＝γ(a-b,x-y)是时空点对(sa-sb,tx-ty)的变异函数；μ是拉格朗日乘数，通过预测公式可得研究区域内任意点的插值估计。

下面结合试验对本发明作进一步的叙述。

1现场监测试验

1.1、现场监测试验方案

试验区是娄衡高速公路笋安山隧道右线典型监测断面，具体试验监测区域为yk31+860-880段。该段监测区域为进口浅埋段，围岩主要为黏性土夹碎石呈松散破碎状，基本无自稳能力。由于该段围岩拱顶下沉沉降量较大，大大超出设计预留沉降量且拱顶喷射混凝土局部有开裂现象，故拱顶沉降试验监测断面间距均加密设置为3米。

(1)监测内容

a)隧道掌子面地质状况观察

其目的是及时了解和掌握隧道开挖掌子面的岩型特征、成分、颜色、结构构造、节理产状、裂隙发育程度、地下水状况等，以便对围岩级别及时做出合理判断。

(b)隧道周边收敛位移量测

其目的是为判断隧道空间稳定提供可靠信息；根据变形速度判断围岩稳定程度，为二次衬砌提供合理支护时机；指导现场设计与施工。

c)隧道拱顶下沉位移量测

其目的是根据拱顶下沉速度判断围岩稳定程度，为二次衬砌提供合理支护时机；指导现场设计与施工。

(2)测点布置

周边收敛及拱顶下沉量测根据围岩类别、隧道埋深、开挖方法等，沿隧道纵向在隧道墙中布设测点(测点与基线的布置将视具体施工方案的变化进行修改和调整，一般在围岩级别变化地段应进行布点量测)，测点间距为3m，测点布置见图2。

(3)量测频率

因试验需要，按照每天一次的量测频率进行监测。

(4)仪器设备

本监测量测试验实施过程中主要使用仪器详见下表1。

表1本试验投入仪器设备表

2、现场试验数据分析

取yk31+863、yk31+866、yk31+869、yk31+872、yk31+875、yk31+878六个断面持续观测34d的拱顶沉降数据为基础实验数据(实测数据如表2、3、4所示)。在灰色预测中，gm(1,1)模型是灰色理论的基础模型，也是应用最为广泛的灰色理论模型。由于gm(1,1)模型波动较小，只能描述简单的数据变化过程，且建模序列必须满足等时距的要求。而在隧道工程监测过程中，往往现场监测数据都是非等时距的，数据的波动存在较大的情况，从而产生较大的误差。本发明从实际出发考虑到以上因素，应用灰色verhulst模型进行基本数据分析。该模型克服了传统的gm(1,1)模型预测精度不高和预测稳定性等问题，对提高围岩位移的长时间预测精度具有非常重要的理论价值和实践意义。因此，对以上试验断面进行分组建模并采用灰verhulst预测模型预测结果与实测数据结果进行对比分析。

表2yk31+863、yk31+866断面拱顶沉降实测变形值

表3yk31+869、yk31+872断面拱顶沉降实测变形值

表4yk31+875、yk31+878断面拱顶沉降实测变形值

定义非负数列x⁽⁰⁾＝{x⁽⁰⁾(1)，x⁽⁰⁾(2)，…，x⁽⁰⁾(n)}为原始沉降数据数列，x⁽¹⁾＝{x⁽¹⁾(1)，x⁽¹⁾(2)，…，x⁽¹⁾(n)}为x⁽⁰⁾的一阶累加生成数列，z⁽¹⁾＝{z⁽¹⁾(1)，z⁽¹⁾(2)，…，z⁽¹⁾(n)}为x⁽¹⁾的紧邻均值生成数列，又那么有：

x⁽⁰⁾(k)+az⁽¹⁾(k)＝b(z⁽¹⁾(k))^α(15)

上式称为gm(1,1)幂模型，当α＝2时：

x⁽⁰⁾(k)+az⁽¹⁾(k)＝b(z⁽¹⁾(k))²(16)

上式(16)即为灰色verhulst模型。

灰色verhulst模型的白化方程如下：

式中：a为发展系数；b为灰作用量。对式(17)的微分方程进行求解可得：

对上式参数a、b应用最小二乘法进行估计，将还原得若且：

则灰色verhulst模型的时间响应函数表达式为：

由上式(15)～(19)建立灰色verhulst预测模型，可计算得到相应断面的的时间响应函数如下：

1)yk31+872拱顶沉降时间响应函数为：

2)yk31+875拱顶沉降时间响应函数为：

3)yk31+878拱顶沉降时间响应函数为：

由式(20)～(22)可计算得出预测变形值(见表5)；灰色verhulst预测变形值与实测值的位移-时间曲线对比见图4～图6。

表5yk31+872、yk31+875、yk31+878断面拱顶沉降灰色verhulst预测值

3、kriging时空统一模型预测分析

1、计算方案

(1)第一种方案：由三个已知监测点数据为基础预测前方相邻的第四个监测点的变形(三测四)。将上述现场试验数据分成三组分步导入kriging时空模型进行计算，第一组：以yk31+863、yk31+866、yk31+869三个断面34d的监测点数据预测前方yk31+872断面监测点34d内的拱顶沉降数据；第二组：以yk31+866、yk31+869、yk31+872三个断面34d的监测点数据为基础预测前方yk31+875断面监测点34d内的拱顶沉降数据；第三组：以yk31+869、yk31+872、yk31+875三个断面34d的监测点数据为基础预测前方yk31+878断面监测点34d内的拱顶沉降数据。理论分析和实践表明：预测误差与块金值(c0)的大小呈正比例关系，变程对预测误差的影响较大，而基台值(c0+c)对预测结果的影响很有限，每组计算估计的变异函数相关参数如表6、表7所示：

表6空间变异函数参数

表7时空变异函数参数

(2)第二种方案：由五个已知监测点数据为基础预测前方相邻的第六个监测点的变形(五测六)。以上述现场试验数据即yk31+863、yk31+866、yk31+869、yk31+872、yk31+875五个断面34d的监测点数据预测前方相邻第六各监测点的yk31+878断面监测点34d内的拱顶沉降数据。为了与第一种方案进行对比分析，本组计算采用的变异函数相关参数如表8。

表8空间变异函数参数

(3)第三种方案：保持第一种方案(三测四)的参数不变，试验断面的间距调整为每6米一个监测点。以上述现场试验数据即yk31+863、yk31+869、yk31+875三个断面34d的监测点数据预测前方相邻yk31+881断面监测点34d内的拱顶沉降数据。该组计算采用的半变异函数相关参数如表6、表7所示。

4、预测结果及对比分析

(1)采用时空kriging预测理论对第一种方案的三组拱顶沉降数据分别进行计算，计算结果如表8所示，对比如图7～图9所示，matlab计算程序如下：

clearall；closeall；

clc；

c0＝0.01；

c＝0.20；

a＝20；

a＝xlsread('data.xls','sheet2','e4:f37')；

b＝xlsread('data.xls','sheet2','h4:i37')；

c＝xlsread('data.xls','sheet2','b4:c37')；

d＝xlsread('data.xls','sheet2','l4:m37')；

％e＝xlsread('data.xls','sheet2','o4:p37')；

％f＝xlsread('data.xls','sheet2','r4:s37')

i＝xlsread('data.xls','sheet2','l4:m37')；

z＝[a(:,1)b(:,1)c(:,1)d(:,1)]；

％z＝[a(:,1)b(:,1)c(:,1)d(:,1)e(:,1)f(:,1)]；

g＝[a(:,2)b(:,2)c(:,2)d(:,2)]；

％g＝[a(:,2)b(:,2)c(:,2)d(:,2)e(:,2)f(:,2)]；

[d,b]＝size(z)；

fori＝1:d

z1＝z(i,1)；

z2＝z(i,2)；

z3＝z(i,3)；

z4＝z(i,4)；

％z5＝z(i,5)；

％z6＝z(i,6)；

position＝[1z1；2z2；3z3；4z4；]；

％position＝[1z1；2z2；3z3；4z4；5z5；6z6]；

c_mean＝covar_m(position,c0,c,a)；

covar_gk＝covarivance(position,c0,c,a)；

lambda_nu＝lambda(covar_gk,c_mean)；

gk＝general_k(lambda_nu,position)；

cj(i)＝gk；

wc_cj(i)＝abs(cj(i)-i(i,1))/i(i,1)；

end

fori＝1:d

z1＝g(i,1)；

z2＝g(i,2)；

z3＝g(i,3)；

z4＝g(i,4)；

％z5＝g(i,5)；

％z6＝g(i,6)；

position＝[1z1；2z2；3z3；4z4；]；

％position＝[1z1；2z2；3z3；4z4；5z5；6z6]；

c_mean＝covar_m(position,c0,c,a)；

covar_gk＝covarivance(position,c0,c,a)；

lambda_nu＝lambda(covar_gk,c_mean)；

gk＝general_k(lambda_nu,position)；

wc_ljcj(i)＝abs(ljcj(i)-i(i,2))/i(i,2)；

end

x＝1:1:34；

y＝cj；

figure(1)

plot(x,y,'k-')；

xlabel('日期')；

ylabel('沉降值')；

holdon

gridon

％随天数的增加，第k+1个桩沉降估计值与实际沉降值的误差

y1＝wc_cj；

figure(2)

plot(x,y1,'k-')；

xlabel('日期')；

ylabel('沉降值误差')；

holdon

表9yk31+872、yk31+875、yk31+878断面拱顶沉降kriging预测值

计算结果及成果图如图7～图9所示：

(2)采用时空kriging预测理论对第二种方案的拱顶沉降数据分别进行计算，计算结果如表10所示，对比如图16～图18所示，

表10yk31+878断面拱顶沉降实测值与kriging预测值

(3)采用时空kriging预测理论对第三种方案的拱顶沉降数据分别进行计算，变异函数参数见表11所示，计算结果如表11所示，对比如图19～图21所示：

表11空间变异函数参数

部分matlab程序如下：

function[covar_gk]＝covarivance(position,c0,c,a)

％该函数计算普通克里金插值的变异函数矩阵

％第k+1个测点与其他桩测点的距离

％由距离计算普通克里金插值的变异函数矩阵

表12yk31+881断面拱顶沉降实测值与kriging预测值

5、精度评定

预测结果是否符合实际，应该进行交叉验证，并对本预测模型进行检验评价。本发明根据式(23)～式(27)，分别计算出预测误差(rmse)、平均误差(re)及相关系数。由数理统计学可知：预测误差值越小则说明预测结果越准确，预测值与真实值越接近；平均误差越接近于“0”，估计值越准确，精度越高；相关系数r²值越接近“1”，表明预测曲线与实测曲线的线性相关程度越高。

均方根误差：

z(θ)＝{z(h,t)，z(h),z(t)}(24)

式中：z^*(θ)为估计值，z(θ)为在相同点的观测值；n为估计点数量。

平均误差：