翁烨,邵德盛
(1.云南省地震局,云南 昆明 650200; 2.昆明理工大学国土资源与工程学院,云南 昆明 650093)
在经典线性平差模型中,采用最小二乘求解未知参数称为最优线性无偏估计,而在实际中因|BTPB|≈0导致法方程的参数解变得不稳定,此时的法方程称为病态方程[1]。最小二乘估计虽然是无偏估计,但是在自变量相关性较强或者某些自变量变化范围较大时的估计结果差强人意,基于对法方程系数矩阵的病态性考量,Hoerl和Kenard通过对法方程系数矩阵对角线上添加一个参数k得到了岭估计[2,3],岭估计目的在于按照正则化法则取目标函数的最优化选择,达到减小均方误差的同时提升参数的估计可信度[4]。岭估计中主要是关于岭参数的求解,因此岭参数的选取至关重要,目前关于岭参数的选取方法主要有岭迹法、L-曲线法[5,6]、U-曲线法[7,8]、广义交叉检核法[9]。文献[10]总结了8种求取岭参数的方法、对含有较多自变量的大型回归问题的解法进行研究;文献[11]提出了一种可以利用典则参数重复计算的方法,使得均方误差逐步减小,该方法改进了Hoerl和Kennard的相应方法;文献[12]提出了一种广义岭估计的直接解法,不需要计算岭参数就可以直接求得具有最小均方误差的解。文献[13]阐述了在病态情形下对于病态矩阵的两步解法。由于广义岭估计算法确定岭参数太过繁杂,计算量庞大,而岭估计对病态问题的参数改正考虑得不够全面,缺乏针对性。本文在广义岭估计的基础上将参数进行分块拟定,根据特征值的贡献率确定分块大小,按照法方程矩阵的病态性影响大小分类计算,将不同的对角线元素分块添加岭参数。在均方误差下验证了新的岭估计优于狭义岭估计和LS估计,推导出新的岭估计岭参数选取方法,通过算例验证了该估计的可行性和适用性。
在线性模型[14]
(1)
中,L为观测向量,B为系数矩阵,X为未知参数,△为观测向量的随机误差,Q已知正定矩阵且满足Q=P-1,P为权矩阵。设G为正交方阵,使得GT(BTPB)G=Λ=diag(λ1,λ2,…,λt),λi(i=1,2,…,t)为法方程系数阵的t个特征根,引入典则参数α=GTX,并记A=BG。
典则参数α在模型(1)中的最小二乘估计和岭估计分别为[14]:
(2)
(3)
广义岭估计的参数估计为:
(4)
广义岭估计的均方误差为[14]:
(5)
(6)
本文在线性模型(1)上重新构造岭估计,相对于式(3)做出改进:
(7)
(8)
(9)
(10)
(11)
=f1(k1)+f2(k2)+f3(k3)+f4(k1)+f5(k2)+f6(k3)=F(k)
(12)
(BTPB+GKGT)-1BTB(BTPB+GKGT)-1
(13)
(14)
(15)
(16)
(17)
(18)
图1 岭回归参数估计
观测值真值以及添加微小误差后的观测值Li 表1
模型参数X的计算结果表 表2
参数估计的协因数矩阵 表3
采用文献[19]和[20]中的数值算例,含有5个未知数,病态设计矩阵和观测向量真值如下所示:
法矩阵N的条件数为cond(N)=2.084×104,严重病态性。岭迹法得出的岭参数k=0.10,如图2所示。表4中显示了求未知参数的LS估计、岭估计以及本文方法的计算结果,新的岭估计参数值的协因数矩阵如表5所示。
图2 岭迹法确定岭参数
模型参数X的计算结果表 表4
参数估计的协因数矩阵 表5
通过算例1~2可以得出以下结论:
新的岭估计可以有效地改善LS估计和狭义岭估计,提升参数估计的可信度。就算例一而言,新的岭估计相比较LS估计,残差值和均方误差值分别为:
均方误差减小明显,残差增加也在可控范围内。相比较狭义岭估计,参数估计可信度也有一定的提升。在算例二中,新的岭估计相比较LS估计,残差值和均方误差值分别为:
因此本文新的岭估计优于LS估计和狭义岭估计。
新的岭估计相比较于LS估计和普通岭估计,对参数的估计更加贴近于真值,相比较LS估计、狭义岭估计更加具有可信度;新的岭估计相比较LS估计和岭估计在均方误差上有所改善,同时残差值增加也在限制范围内。本文的广义平差方法针对贡献率较小的多数特征值有很大的改进作用,能在一定要求范围内快速求解参数估计值。可以在大型病态矩阵中减少计算量,提高计算效率。