变异函数模型参数的加权总体最小二乘回归法

赵英文 王乐洋，2，3

1 东华理工大学测绘工程学院，南昌市广兰大道418号，330013

2 流域生态与地理环境监测国家测绘地理信息局重点实验室，南昌市广兰大道418号，330013

3 江西省数字国土重点实验室，南昌市广兰大道418号，330013

变异函数能同时描述区域变量的随机性和结构性［1］，在研究数据的分布特性和数据插值等方面发挥着重要的作用。变异函数模型参数估计方法中，人工拟合法采用肉眼来确定变异函数模型的参数，效率低且可靠性不高。文献［2］提出用加权回归多项式法来拟合变异函数模型的参数，但没有解决参数正负号问题。线性规划法［3－4］、最小二乘法［5－6］和加权最小二乘法［7－9］虽然提高了参数的计算效率，但只认为变异函数值含有误差，没有考虑到变异函数模型中距离值的随机误差。文献［10］首次提出总体最小二乘概念，从数值分析的角度解决系数矩阵也含有误差的平差问题。当系数矩阵含有误差时，最小二乘解是有偏的，而总体最小二乘解是无偏的［11－12］。考虑到距离值的随机误差，文献［13］提出用总体最小二乘法求解变异函数模型参数，认为变异函数值和距离值是等精度的，把变异函数值的权阵作为行尺度矩阵左乘观测向量和系数矩阵，然后用SVD 分解法进行解算。

本文将加权总体最小二乘回归法引入到变异函数模型参数估计中。以幂函数模型为例，通过协方差传播律［14］发现分组后变异函数值和距离值是不等精度的，并给出距离值的定权方法。再结合熵权法和点对数法进行参数的迭代求解，最后通过模拟数据和实测数据验证加权总体最小二乘回归法的合理性和优越性。

1 变异函数模型及参数估计

1.1 变异函数模型

离散型变异函数γ（h）可由式（1）求解：

式中，h为两点间的分隔距离，N（h）为距离为h的点对数目，Z（xi，yi）和Z（（xi，yi）＋h）分别为位置（xi，yi）和（xi，yi）＋h处的区域变化量。

根据不同的形状和结构，变异函数模型可分为球状模型、高斯模型、指数模型、幂函数模型和孔穴效应模型等［1］。幂函数模型为：

式中，M为常系数，α为幂指数。由文献［1］知，α必须小于2；若α≥2，则幂函数不再是变异函数。

1.2 变异函数模型参数估计

变异函数模型事先是未知的，一般的做法是利用已知样本点计算出的距离值h和变异函数值γ（h）拟合出相应的变异函数模型。实际中点位分布散乱，需在计算出任意点对间的距离后对距离进行分组。给定一个容许误差T，落在（h－T，h＋T）内的距离值个数为N（h），则该组距离值为所有距离值的均值，即。

以幂函数模型为例，对式（2）两边取对数：

令

式中，y∈Rm×1，A∈Rm×n，x∈Rn×1，m为分组后距离的个数，n为未知参数个数，幂函数模型n取2。

考虑到观测值的误差，幂函数模型线性化后的回归模型简记为：

式中，y表示变异函数值，A为包含距离的系数阵，x为幂函数模型参数。

式中，e为向量y的误差向量，EA为系数矩阵A的误差矩阵，为参数最佳估值，vec（·）为矩阵按列拉直运算，为单位权方差，QE为A的协因数阵，Qe为y的协因数阵。

如果向量y和系数矩阵A的元素等精度，根据总体最小二乘平差准则求得参数的总体最小二乘解：

式中，e为向量y的误差向量，EA为系数矩阵A的误差矩阵，为参数最佳估值，vec（·）为矩阵按列拉直运算，为单位权方差，QE为A的协因数阵，Qe为y的协因数阵。

如果向量y和系数矩阵A的元素等精度，根据总体最小二乘平差准则求得参数的总体最小二乘解：

总体最小二乘解法主要有SVD 分解法和迭代解法［12］。顾及到系数矩阵A中含有常数列，文献［15］提出一种顾及系数矩阵部分含有误差的加权总体最小二乘迭代解法，把算法中的协因数阵变为单位阵，得到该平差问题的总体最小二乘解：

如果顾及y和A中元素不同的精度和贡献程度，变异函数模型参数估计问题成为加权总体最小二乘平差问题。文献［16］提出一种加权总体最小二乘算法，根据加权总体最小二乘平差准则：

该算法的迭代过程为［16］：

1）计算加权最小二乘解作为迭代初始值：

3）迭代终止，得到解xi。

1.3 权的确定

1.3.1 向量y对应权阵Py的确定

Py的确定方法主要有熵权法［7］和点对数法［9］。依据文献［7］，把分组后的距离h和其对应的个数N（h）作为影响因子构造评价矩阵R，取＝［hN（h）］，得到变异函数值的权Wi。文献［9］提出的点对数法定权公式为：

1.3.2 系数矩阵A对应权阵PA的确定

假设观测值z的随机误差，令zi和zj分别等于式（1）中的Z（xi，yi）和Z（（xi，yi）＋h），由一个点对求得变异函数值γ（h）的方差D1：

顾及到式（1），由N（h）个点对求得变异函数值γ（h）的方差D（γ（h））：

分组后距离h由N（h）个符合分组条件的hij先求和再除以N（h）得到，则距离h的方差为：

由式（15）和（17），且一般情况下m个距离h对应的N（h）不完全相同，所以分组后变异函数值和距离值是不等精度的。从而，距离的权可定义为方差的倒数，即，W（h）为m维列向量。

在计 算 中，取Py＝diag（Wi），PA＝diag

上述提到的权适用于幂函数模型，而对于其它的变异函数模型，可在模型线性化后，按照文献［7］和［9］的方法确定权阵Py，依据本文的推导思路确定权阵PA。因为不同的函数模型对应的系数矩阵会含有不同的距离表现形式，所以权阵PA需依据具体的函数模型确定。

2 算例分析

2.1 模拟算例

模拟一组规则分布的坐标（x，y），并计算任意两点间的距离。假设变异函数模型为幂函数模型，通过给定的幂函数γ（h）＝0.2h1.5代入h计算得到γ（h）。由MATLAB 产生一组均值为0、标准差初始值为0.009、步长为0.005、终止值为0.049的正态分布随机误差序列，根据式（15）和（17）分别在h和γ（h）加上相应的误差，共产生9组不同标准差下的数据，每组模拟200次。采用4种方法进行参数估计，其中WLS1 和WTLS1的定权方法是熵权法，WLS2 和WTLS2 的定权方法是点对数法，TLS采用文献［15］方法，WTLS采用文献［16］方法，‖ΔX‖为参数M和α的估值与真值之差的2范数。估计结果见表1。

表1 不同方法的估计结果Tab.1 The results estimated by different methods

表1中结果为9种不同标准差下的均值。由表1可知，加权总体最小二乘法所得的参数残差范数均小于其他3种方法，采用点对数法定权的结果优于熵权法定权的结果。由图1可知，当所加误差较小时，几种方法所得结果比较接近。随着所加误差的增大，加权总体最小二乘法所得结果更接近于参数的真值，说明加权总体最小二乘法在参数估计方面具有更高的精度和合理性。

2.2 高程异常插值算例

区域一数据来自于文献［17］，其GPS控制网由17个同精度GPS水准点构成，高程异常值变化平缓。选取5个点作为已知点，剩余12个作为检核点。区域二数据来自于文献［18］，其GPS控制网由24个同精度GPS水准点构成，高程异常值变化较大。选取10 个点作为已知点，剩余14个作为检核点。点位分布如图2。

图2 点位分布图（左图为区域一，右图为区域二）Fig.2 Distribution of points（area 1on the left side，area 2on the right side）

对上述两组数据分别计算距离值和变异函数值，在距离分组后选取幂函数模型作为变异函数模型，参数估计方法同上，使用RMS 和STD 作为评价指标。

RMS为均方根预报误差：

STD 为预报残差标准差：

式中，Zi、Zi′分别为k个检核点的观测值和通过Kriging插值计算得到的检核点估计值，ri是检核点预报残差，是检核点预报残差的均值。结果如表2、表3所示。

表2 区域一不同方法的检核点预报结果Tab.2 The predicted results of check points calculated by different methods in area 1

表3 区域二不同方法的检核点预报结果Tab.3 The predicted results of check points calculated by different methods in area 2

图3 区域一检核点预报值与观测值差值曲线图Fig.3 Difference of check points between the predicted values and the observed ones in area 1

图4 区域二检核点预报值与观测值差值曲线图Fig.4 Difference of check points between the predicted values and the observed ones in area 2

由表2、表3可知，无论RMS还是STD 都表明，加权总体最小二乘法的精度最高，总体最小二乘法和加权最小二乘法次之，最小二乘法效果最差。加权总体最小二乘法使两个区域的变异函数模型参数估计精度分别提高70%和60%左右，均方根预报误差分别减少5mm 和98mm。图3和图4中，WLS和WTLS对应于表中的WLS2和WTLS2，可知加权总体最小二乘法对应残差分布曲线较其他两种方法变化更平缓，更接近横坐标轴。两个区域的数据均表明，使用点对数法对应的加权最小二乘法所得到的结果优于按熵权法对应的结果，但是两种定权方法对应的加权总体最小二乘法结果却十分接近。算例中每组数据的两种WTLS法得到的幂指数α均大于2，考虑到§1.1提到的幂函数模型适用条件，为保持模型的变异函数特性，在进行Kriging插值过程中，重新把α取值为1.999 999，这也与文献［7］的做法一致。两种WTLS法得到的幂函数模型常系数不同、幂指数相同，这就使得到的RMS 和STD差别很小。数据一结果的精度维持在mm 级，数据二结果的精度由dm 级提高到了cm 级，加权总体最小二乘法在这两种情形下都能有效地提高参数的估计精度。以上结果证明，把加权总体最小二乘法引入到变异函数领域进行参数估计是可行和有效的。

3 结 语

本文在考虑距离值误差的基础上，通过协方差传播律发现分组后的变异函数值和距离值是不等精度的，并给出距离值的定权方法，再结合熵权法和点对数法两种变异函数值的定权方法，把加权总体最小二乘回归法引入到变异函数模型参数估计中。模拟数据和实测数据证明了加权总体最小二乘回归法的可行性和有效性，相对于最小二乘法、加权最小二乘法和总体最小二乘法，加权总体最小二乘回归法能得到更高精度的变异函数模型参数估值。本文仅对幂函数模型进行了算例讨论，而对于其他变异函数模型的适用性算例验证以及距离误差对函数模型的影响机制、进一步提高加权总体最小二乘法的参数估计精度和解算效率，还有待于研究。

［1］徐建华.现代地理学中的数学方法［M］.北京：高等教育出版社，2002（Xu Jianhua.Mathmatica Methods in Contemporary Geography［M］.Beijing：Higher Education Press，2002）

［2］王仁铎，胡光道.线性地质统计学［M］.北京：地质出版社，1988（Wang Renduo，Hu Guangdao.Linear Geostatistics［M］.Beijing：Geological Press，1988）

［3］矫希国，刘超.变差函数的参数模拟［J］.物化探测技术，1996，18（2）：157－161（Jiao Xiguo，Liu Chao.Estimation of Variation Parameter［J］.Computing Techniques for Geophysical and Geochemical Exploration，1996，18（2）：157－161）

［4］李玲，何涛，张武，等.变异函数线性化的统一参数估计方法研究［J］.长江大学学报：自然科学版（理工卷），2010，7（2）：127－129（Li Ling，He Tao，Zhang Wu，et al.Study on the Unity Parameter Estimation Method of Linear Variogram［J］.Journal of Yangtze University：Nat Sci Edit，2010，7（2）：127－129）

［5］李明，高星伟，文汉江，等.Kriging方法在GPS水准拟合中的应用［J］.测绘科学，2009，34（1）：106－107（Li Ming，Gao Xingwei，Wen Hanjiang，et al.The Application of Kriging Method in GPS Leveling Fitting［J］.Science of Surveying and Mapping，2009，34（1）：106－107）

［6］郭泉河，李秀海.不同变异函数的泛Kriging法的GPS高程拟合结果［J］.黑龙江工程学院学报：自然科学版，2011，25（4）：26－28（Guo Quanhe，Li Xiuhai.Application of Universal Kriging Technology with Different Semivariation Function Models to GPS Height Anomaly Fitting［J］.Journal of Heilongjiang Institute of Technology，2011，25（4）：26－28）

［7］潘家宝，戴吾蛟，章浙涛，等.变异函数模型参数估计的信息熵加权回归法［J］.大地测量与地球动力学，2014，34（3）：125－128（Pan Jiabao，Dai Wujiao，Zhang Zhetao，et al.Parameter Estimation of Variogram Model by Using Information Entropy Weighted Regression［J］.Journal of Geodesy and Geodynamics，2014，34（3）：125－128）

［8］严华雯，吴健平.加权最小二乘法改进遗传克里金插值方法研究［J］.计算机技术与发展，2012，22（3）：92－95（Yan Huawen，Wu Jianping.Reasearch on Genetic Algorithm Kriging Optimized by Weight Least Square［J］.Computer Technology and Development，2012，22（3）：92－95）

［9］曾怀恩，黄声享.基于Kriging方法的空间数据插值研究［J］.测绘工程，2007，16（5）：5－13（Zeng Huaien，Huang Shengxiang.Research on Spatial Data Interpolation Based on Kriging Interpolation［J］.Engineering of Surveying and Mapping，2007，16（5）：5－13）

［10］Golub G H，Loan C V.An Analysis of the Total Least－Squares Problem［J］.SIAM Journal on Numerical Analysis，1980，17（6）：883－893

［11］王乐洋.总体最小二乘解性质研究［J］.大地测量与地球动力学，2012，32（5）：48－52（Wang Leyang.Research on Properties of Total Least Squares Estimation［J］.Journal of Geodesy and Geodynamics，2012，32（5）：48－52）

［12］王乐洋，许才军.总体最小二乘研究进展［J］.武汉大学学报：信息科学版，2013，38（7）：850－856（Wang Leyang，Xu Caijun.Progress in Total Least Squares［J］.Geomatics and Information Science of Wuhan University，2013，38（7）：850－856）

［13］Felus Y A，Schaffrin B.A Total Least－Squares Approach in Two Stages for Semivariogram Modeling of Aeromagnetic Data［C］.IAMG2005，Toronto，2005

［14］王乐洋，鲁铁定.总体最小二乘平差法的误差传播定律［J］.大地测量与地球动力学，2014，34（2）：55－59（Wang Leyang，Lu Tieding.Propagation Law of Errors in Total Least Squares Adjustment［J］.Journal of Geodesy and Geodynamics，2014，34（2）：55－59）

［15］Fang X.Weighted Total Least Squares Solutions for Applications in Geodesy［D］.Hanover：Leibniz University of Hanover，2011

［16］Jazaeri S，Amiri－Simkooei A R，Sharifi M A.Iterative Algorithm for Weighted Total Least Squares Adjustment［J］.Survey Review，2014，46（334）：19－27

［17］朱卫东，李全海.基于标准化动量BP神经网络的GPS高程转换［J］.大地测量与地球动力学，2010，30（1）：123－125（Zhu Weidong，Li Quanhai.Conversion of GPS Height Based on Standardization Momentum BP Neural Network［J］.Journal of Geodesy and Geodynamics，2010，30（1）：123－125）

［18］黎剑.区域GPS高程异常拟合及建模方法研究［D］.昆明：昆明理工大学，2013（Li Jian.Research on Regional GPS Height Anomaly Fitting and Modeling［D］.Kunming：Kunming University of Science and Technology，2013）