病态可分离非线性最小二乘问题的一种解算方法及其应用

本发明公开了病态可分离非线性最小二乘问题的一种解算方法及其应用，属于测量平差技术领域

背景技术：

病态问题常存在于在大地网平差、变形监测网平差以及gps定位等测量数据处理中，影响平差成果的精度，其危害性十分严重。分析测量平差中的病态矩阵、减弱病态性对平差成果影响是测量数据处理面临的一个重要课题。病态矩阵往往具有较大的条件数，以其为系数阵构成的方程组的解对微小扰动十分敏感，使得方程组的解数值稳定性较差，对方程求解带来很大的困难。其中，基于奇异值分解的变量投影法能够在一定程度上克服方程组系数阵的病态性，若矩阵病态程度较为严重，通常可采用奇异值截断或修正的方法。然而，无论是单一地对奇异值截断或是先修正后截断方法，都是一种有偏估计方法，而对全部奇异值进行统一修正的方法在降低方差的同时也引入了偏差。

技术实现要素：

本发明公开了病态可分离非线性最小二乘问题的一种解算方法及其应用，以解决现有技术中，奇异值截断方法引入偏差降低平差结果精度的问题。

病态可分离非线性最小二乘问题的一种解算方法，使用变量投影法分离非线性函数中的线性参数与非线性参数，结合矩阵的奇异值分解法对两类参数分别求解；

所述奇异值分解法中，对分解后的u、s、v三个矩阵中的s矩阵进行修正。

优选地，所述变量投影法分离非线性参数的步骤为：

可分离非线性最小二乘问题的模型是非线性函数的线性组合形式，其观测方程用矩阵表示为：

y＝φβ+r(1)；

式中，为非线性函数，为φ的每个列向量，包含了模型变量t和非线性参数α＝(α1,α2,…,αn)^t；β＝(β1,β2,…,βn)^t是线性参数；n为非线性参数或线性参数的个数；y为t对应的已知观测数据，观测值个数为m；r为残差向量；

式(1)的求解式为：

将矩阵φ通过奇异值分解法分解为u、s、v三个矩阵相乘的形式：φ＝usv^t，令u＝[u1u2]，其中u1、u2分别为m×r、m×(m-r)阶矩阵，r为φ的秩；

令式(2)中的β计算按下式：β＝φ⁺y，结合φ的奇异值分解式，得变量投影函数如下：

式中，φ⁺是φ的广义逆矩阵，rsvd表示基于奇异值分解的残差向量表达形式；利用lm算法迭代，求解使得残差平方和最小的非线性参数α。

优选地，所述变量投影法分离线性参数的步骤为：

对于良态方程组，求解线性参数β如下式：β＝φ⁺y(4)；

对于病态方程组，求解线性参数β如下式：β＝φ⁺y＝vs^-1u^ty(5)。

优选地，所述对s矩阵进行修正的具体方法为：

令式(5)中的φ⁺按下式计算：

式中，表示φ的每个奇异值，m(k)为奇异值的修正项，k为修正参数。

优选地，式(5)更换为：

式中，表示φ的第i个奇异值，表示最小奇异值，表示最大奇异值，s为截断参数，s<r，k为修正参数。

优选地，所述截断参数s的确定方法为：

利用l曲线法确定截断参数s，以为横坐标，为纵坐标，假定函数模型y＝f(x)，以为样本点进行拟合，求解模型参数并绘制曲线，由下式计算曲线的曲率ρ：曲率最大的点对应的参数s即为截断参数。

优选地，所述修正参数k的确定方法为：

利用h-k公式确定修正参数k，可分离非线性最小二乘问题的观测方程为y＝φβ+r，设(ω1,…,ωn)为φ^tφ特征根(λ1,…,λn)对应的特征向量，由下式计算修正参数k：

式中，表示中第i个元素，ω＝(ω1，…，ωn)，λ＝diag(λ1，…，λn)。

优选地，一种局部地区坐标转换参数的解算方法，包括：

s1.选取一组点位坐标作为原点，将坐标转换的实验数据代入局部地区坐标转换模型：

式中，(x10，t10，z10)，(x20，y20，z20)为原点分别在两空间直角坐标系下的坐标值，(x1，y1，z1)，(x2，y2，z2)表示点位分别在两空间直角坐标系下的坐标值，c为缩放尺度，即线性参数，(α，β，γ)旋转参数，即非线性参数，模型对应的残差函数表达式的矩阵形式如式(1)所示；

s2.将非线性矩阵φ进行svd分解，并结合式(3)得到其对应的变量投影函数表达式；

s3.在最小二乘原则下对基于svd分解的变量投影函数进行迭代解算，求得旋转参数，即非线性参数；

s4.利用式(4)解得缩放参数，即线性参数。

优选地，一种高程异常的计算方法，包括：

(1)将实验数据代入对数函数拟合模型形成可分离非线性方程组，模型对应的残差函数表达式的矩阵形式如式(1)所示；

(2)将矩阵φ进行svd分解，φ＝usv^t，得到其对应的变量投影函数表达式；

(3)在最小二乘原则下对基于svd分解的变量投影函数进行迭代解算，求得非线性参数；

(4)按照式(7)处理矩阵s中的奇异值，得到其对应逆阵按式(6)计算得到线性参数。

优选地，所述对数函数拟合模型为：

ξ(α，β；x，y)＝β0+β1lnα1x+β2lnα2y+β3ln²α3x+β4ln²α4y+β5lnα5xlnα6y；

式中，x，y分别表示控制点的坐标值，待求参数为：

非线性参数α＝(α1，α2，α3，α4，α5，α6)^t和线性参数β＝(β0，β1，β2，β3，β4，β5)^t。

与现有技术相比，本发明探究了计算可分离模型的变量投影算法及其基于矩阵分解的改进算法，并将基于奇异值分解的变量投影法应用于局部地区坐标转换参数解算中，使病态性较弱的方程组保持较好的解算稳定性，提高了坐标转换问题的转换精度；

在此基础上进一步公开了病态可分离问题中求解线性参数的奇异值截断法和修正法，在经典的奇异值截断和修正的基础上，结合计算岭参数的h-k公式提出了一种改进的变量投影算法以克服非线性方程组解算遇到的严重病态问题，提高了高程异常的拟合精度。

附图说明

图1为本发明的s矩阵修正的技术流程图；

图2和图3为八种模型的高程异常拟合成果图，具体模型为：

图2(a)：vptsvd；图2(b)：vpmsvd+m1；图2(c)：vpmsvd+m2；图2(d)：vpmtsvd；

图3(a)：qsls；图3(b)：csls；图3(c)：hyperboloidls；图3(d)：rhls。

具体实施方式

下面结合具体实施方式对本发明作进一步详细说明：

测量平差中常用的guass-markov模型为：

l+v＝ax(1)

式中，l为m×1阶观测值向量，v为m×1阶残差向量，a为m×n(m≥n)阶系数阵，x为n×1阶待估参数向量。依据最小二乘准则，可导出式(1)的法方程、最小二乘解及其协方差阵为：

若g-m模型中的矩阵a和观测向量l有微小扰动，则会引起式(2)中矩阵n和w分别产生误差δn和δw，相应地x产生误差δx，即：

(n+δn)(x+δx)＝w+δw(3)

顾及nx＝w，将式(3)整理可得：

δx＝-n^-1δn(x+δx)+n^-1δw(4)

将式(4)右端取范数，依据范数次可加性及矩阵与向量范数相容性，有如下不等式：

令c(n)＝‖n‖‖n^-1‖，将式(5)整理可得：

其中，c(n)称为矩阵n的条件数，可作为矩阵病态程度的判断依据。从式(6)可以看出，条件数c(n)表示矩阵n和w的相对扰动对解的相对扰动的影响程度：若扰动较小，而引起解的变动较大，此时方程组病态程度较为严重，即不等式右边数值较小，左边数值较大，则需要较大c值以保证不等式成立；反之，c值较小。条件数衡量的是矩阵病态性的相对程度，当条件数大于1000时，可认为矩阵n病态程度较为严重，此时方程组解抗扰动的数值稳定性很差。在这种情况下，须对矩阵进行调整。

可分离非线性最小二乘问题的函数模型为非线性函数的线性组合形式，一般可表示为：

式中，表示非线性函数部分，yi(i＝1,2,…,m)为ti对应的观测数据，β＝(β1,β2,…,βn)^t和α＝(α1，α2，…，αn)^t分别是待估的线性和非线性参数。在给定观测数据(ti，yi)(i＝1，2，…，m)的情况下，依据最小二乘原则，未知参数的求解即求使残差函数平方和最小的和用矩阵表示为：

式中β＝(β1，β2，…βj)^t(j＝1,2,…,n)，‖·‖2表示欧式范数。线性参数β可通过式(9)求得：

β＝φ⁺y(9)

其中，φ⁺是φ的广义逆矩阵。将式(9)代入式(8)，可得：

式(10)称为变量投影函数，仅与非线性参数α有关，非线性参数可通过迭代计算求得。为提高求解的运算效率，通常将m×n的矩阵φ通过奇异值分解法(singularvaluedecomposition,svd)分解为u、s、v三个矩阵相乘的形式。令u＝[u1u2]，φ的奇异值分解式与基于奇异值分解的变量投影表达式为：

式中，矩阵u＝(u1u2…um)、v＝(v1v2…vn)分别为m×m、n×n阶正交矩阵，s为m×n阶矩阵，是由矩阵φ的奇异值由大到小排列而成，进而，可利用levenberg-marquardt(lm)算法迭代求解使得残差平方和最小的参数估值其算法步骤如下：

(1)给定迭代初值α0，设置收敛准则：迭代终止阈值ε，最大迭代次数kmax；

(2)计算第k次迭代时残差函数rsvd的雅可比矩阵j，令阻尼因子μ＝τ·max{diag(j^tj)}，计算hlm＝-(j^tj+μi)j^tr，αk+1＝αk+hlm；

(3)若hlm使得r(αk+1)<r(αk)，则接受αk+1并减小阻尼因子μ，进入下一步，否则，增大阻尼因子μ并重新计算hlm，直至使得r(αk+1)<r(αk)；

(4)若r(αk)-r(αk+1)<ε或k>kmax，则终止迭代，否则，执行(2)-(3)。

在解算方程组的过程中，我们常常遇到病态程度严重较为严重(即条件数c(φ)很大)的系数矩阵，此时单纯采用svd方法求解方程组并不能有效地削弱和克服病态性的影响。将φ的奇异值分解式分别结合式(9)和式(2)可得：

当矩阵φ严重病态时，其奇异值单调趋向于零，由式(12)可知，分母位置上较小的奇异值放大了观测误差对线性参数估值的影响，导致其解不稳定。为此，有学者用截断或修正奇异值的方法对其进行改进。

截断奇异值(truncatedsingularvaluedecomposition,tsvd)是将矩阵s中较小的奇异值舍弃，即令s1中较小的为0。若令进而用s×s、m×s、n×s阶矩阵sts、uts、vts来近似描述原矩阵，则：

式中，该方法通过截去较小的奇异值减小了估值方差，使得参数估值更加稳定，可在一定程度上克服矩阵的病态性。在采用截断奇异值法时，按照何种原则进行截断是应用该方法的关键。修正奇异值(modifiedsingularvaluedecomposition,msvd)是对每个奇异值的数值大小进行调整而非舍弃，即每个奇异值加入一修正项m(k)，奇异值个数不会减少。对奇异值进行修正后φ⁺的表达形式为：

式中，

在参数分离及解算过程中，若利用截断奇异值方法舍弃较小的奇异值，则改变了非线性矩阵φ的结构，是一种有偏估计方法，会使得线性参数解的精度受到一定的影响。本发明在原有奇异值截断与修正方法的基础上提出了一种改进算法，新算法旨在修正较小奇异值以降低方差，同时相应地减小非线性矩阵的条件数。在欧几里得范数下，矩阵的条件数可表达为奇异值最大值与最小值之比由式(12)知，较小奇异值出现在分母位置上，显著放大了参数估值的方差，也导致了解的不稳定性。因此在较小奇异值中加入修正项依据具体数值调整奇异值大小，即矩阵病态程度越严重，奇异值修正幅度越大；同时在较大奇异值中加入修正项对其作微小调整，缩小奇异值最大值与最小值差距，可适当地减小矩阵条件数，进而削弱其病态程度。因此，将非线性矩阵φ进行奇异值分解后，取s(s<r)为截断参数，矩阵s及其逆阵为：

则线性参数为：

l曲线(l-curve)核心是认为曲线上曲率最大点对应的参数即为最优参数。在可分离非线性最小二乘问题中，利用l曲线法确定截断奇异值算法的截断参数s，是以为横坐标，为纵坐标，并假定函数模型y＝f(x)，以为样本点进行拟合，求解模型参数并绘制曲线，再由式(17)计算曲线上曲率最大的点：

曲率最大的点对应的参数s即为截断参数。

改进算法中的修正项参数k可采用岭估计中岭参数的选取方法。本发明利用h-k公式选取修正参数k。在可分离非线性最小二乘问题中，观测方程为y＝φβ+r，设(ω1,…,ωt)为φ^tφ特征根(λ1,…,λt)对应的特征向量，记：

ω＝(ω1,…,ωt)，λ＝diag(λ1,…,λt)(18)

则

式中，

本发明包括两种实施例，实施例一为变量投影法分离非线性函数中的线性参数与非线性参数，在测量学科的空间直角坐标转换中，两个不同的空间直角坐标系进行转换需要用七个参数来描述，其中包括三个坐标轴的平移参数、三个旋转参数和一个坐标轴缩放参数，若已知两空间直角坐标系下的三个及以上的公共点坐标值，那么这七个参数可以唯一确定。在进行坐标转换时，常用的七参数模型的矩阵形式为：

当点位分布集中在一个范围较小的局部地区时，应用七参数法求得的转换参数，尤其是平移参数精度是不高的，公共点坐标很小的变化会引起参数较大的变化。因此为了保证参数解算的准确性，局部地区的坐标转换参数的计算按照如下方案求解。选取一组点作为“原点”并分别求出各点对原点的坐标差值。设原点坐标在两空间直角坐标系下的坐标值分别为(x10,y10,z10)，(x20,y20,z20)，对于原点由式(20)可得：

式(20)与式(21)作差可得：

式中含有线性参数c和非线性参数(α,β,γ)，可视为非线性函数的线性组合形式。解算非线性模型待估参数的常规算法是将模型线性化，再将所有参数统一进行迭代计算，而本实施例将应用变量投影算法将参数分离开来逐步进行解算。实施例采用独立空间直角坐标系和cgcs2000坐标系下的9组公共点进行实施例，经纬度范围分别为116°39’33.5334”e～116°59’23.1124”e和34°38’17.4987”n～34°54’06.6574”n。9组坐标值如表1所示。

表1.独立空间直角坐标系和cgcs2000坐标系下9组公共点坐标值

将数据代入模型(22)，迭代计算时旋转参数的初始值在弧度制下取α0＝0.1,β0＝0.1,γ0＝0.1，缩放尺度参数为c0＝0.1。利用公式c＝‖φ‖‖φ^-1‖计算非线性矩阵条件数得c(φ)＝1.401，其病态程度极弱，矩阵呈现良态，对φ进行奇异值分解即可使得求解过程稳定。实施例将采用基于奇异值分解的变量投影法(vpsvd)求解模型中的两类参数。实施例过程如下：

(1)选取第2组点位坐标作为“原点”，将表1实施例数据代入模型(22)形成可分离非线性方程组；

(2)通过式(11)将非线性矩阵进行svd分解φ(α,β,γ)＝usv^t，得到其对应的变量投影函数表达式

(3)在最小二乘原则下，对基于svd分解的变量投影函数进行迭代解算，求得旋转参数

(4)利用式(9)解得缩放参数

实施例将vpsvd算法与不分离参数算法(nons)、经典变量投影算法(vp)和基于qr分解的变量投影法(vpqr)所得结果进行了对比。

实施例结果从迭代次数，原始坐标与估值坐标的残差平方和、均方根误差，拟合优度和运算时间五个方面进行了比较，如表2所示，其中拟合优度用决定系数r²表示。四种算法对应的参数解算结果如表3所示。由表2决定系数可知，表3所示的四种算法解算结果都是可靠的。参数不分离算法的残差平方和相对较大，相比之下，参数分离算法对应残差平方和较小，由此判断模型参数解算结果更为精确。由于变量投影算法降低了待估参数的维度，使得迭代计算的复杂度有所减少，因此三种参数分离算法对应的迭代次数和运算时间相比参数不分离算法都有所减少，计算效率有所提升，其中vpsvd算法最为显著。

表2.迭代次数、残差平方和(ssr)、均方根误差(rmse)和决定系数(r²)和运算时间的对比

表3.参数解算结果

实施例二为基于本发明所述的s矩阵修正的高程异常拟合，高程异常拟合模型通常采用二次曲面模型，模型自变量通常为中心化、同量级化后的坐标。本节实施例采用的对数函数拟合模型以平面坐标为自变量、高程异常值为因变量，其表达式具有线性与非线性参数可分离的特性。拟合模型表达式为：

式中，待求参数为α＝(α1,α2,α3,α4,α5,α6)^t，β＝(β0,β1,β2,β3,β4,β5)^t。本实施例将借助已知控制点的平面坐标与高程异常值进行拟合。已知点平面坐标、大地高、正常高、高程异常值如表4所示。

表4已知点平面坐标(x,y)、大地高(h)、正常高(h正常)、高程异常值(ξ)

实施例利用d01～d08数据进行拟合，d09、d10验证参数的可靠性与模型的预测能力。将数据代入模型(25)，设置非线性参数迭代初值均为100，求得非线性矩阵条件数c(φ)＝2.159×10¹¹，矩阵病态程度较为严重。因此实施例将采用模型(25)和本发明提出的变量投影改进算法(vpmtsvd)求解模型参数。过程如下：

(1)将表3实施例数据代入模型(25)形成可分离非线性方程组；

(2)通过式(11)将非线性矩阵进行svd分解φ(α)＝usv^t，得到其对应的变量投影函数表达式

(3)在最小二乘原则下，对基于svd分解的变量投影函数进行迭代解算，求得非线性参数

(4)按照式(15)处理矩阵s中的奇异值，得到其对应逆阵最后利用式(16)计算得到线性参数

实施例将vpmtsvd算法与基于奇异值截断的变量投影法(vptsvd)、基于奇异值修正式的变量投影法(vpmsvd+m1)和基于奇异值修正式的变量投影法(vpmsvd+m2)所得结果，以及二次曲面模型(quadraticsurface,qs)、三次曲面模型(cubicsurface,cs)、以双曲面(hyperboloid)和倒曲面(reciprocalhyperboloid,rh)为核函数的多面函数所得最小二乘解(ls)进行对比。

经迭代计算后得到非线性参数估值，进而计算得到线性参数。各算法对应高程异常函数拟合图像如图2和图3所示。表5列出了各种算法对应拟合与预测的残差平方和、均方根误差以及模型决定系数。由于实施例数据的限制，三次曲面模型只存在拟合结果。对比表5中的数据ssr和rmse可知，所有模型与算法都具有较高的拟合精度，其中，高程异常拟合常用模型qs、cs、hyperboloid和rh精度极高，然而由ssrf和rmsef可知常规模型预测结果可靠性较差，而本发明介绍的拟合模型和变量投影算法在高程异常预测中表现出了更好的性能。在vptsvd、vpmsvd+m1、vpmsvd+m2和vpmtsvd算法中，vptsvd算法精度相对较低，其拟合与预测精度数量级分别为10^-3和10^-2，相比之下，单纯的修正方法vpmsvd+m1、vpmsvd+m2和本发明提出的vpmtsvd的拟合与预测精度数量级为10^-4，算法精度更高。观察图2和图3可得出一致的结论。

表5各算法对应残差平方和(ssr、ssrf)、均方根误差(rmse、rmsef)和决定系数(r²、)

本发明探究了计算可分离模型的变量投影算法及其基于矩阵分解的改进算法，并将基于奇异值分解的变量投影法应用于局部地区坐标转换参数解算中，使病态性较弱的方程组保持较好的解算稳定性。在此基础上进一步研究了病态可分离问题中求解线性参数的奇异值截断法和修正法，在经典的奇异值截断和修正的基础上，结合计算岭参数的h-k公式提出了一种改进的变量投影算法以克服非线性方程组解算遇到的严重病态问题，并通过高程异常拟合实施例对各类算法的可行性与适用性进行研究，得出如下结论：

(1)在局部地区坐标转换参数解算实施例中，常规的参数不分离算法精度稍低，相比之下，参数分离算法求解结果的残差平方和较小，旋转参数和缩放尺度解算结果更为精确，且参数分离算法对应的迭代次数和运算时间相比不分离算法都有所减少，其中基于奇异值分解的变量投影法vpsvd效果最为显著；

(2)在高程异常拟合实施例中，常用的拟合模型qs、cs、hyperboloid和rh拟合精度极高，然而其预测结果可靠性较差，而本发明介绍的对数形式拟合模型和改进的变量投影算法具有较高的预测精度。在四种改进的变量投影算法中，vptsvd算法不论是高程异常值的拟合或是预测，其精度都相对较低，而单纯的修正方法和本发明提出的vpmtsvd算法具有更高的拟合和预测精度。

当然，上述说明并非是对本发明的限制，本发明也并不仅限于上述举例，本技术领域的技术人员在本发明的实质范围内所做出的变化、改型、添加或替换，也应属于本发明的保护范围。