测绘学报 | 王路遥：基于矩阵分解的可分离非线性最小二乘问题求解方法及其应用
本文内容来源于《测绘学报》2022年第11期（审图号GS京（2022）1229号）
基于矩阵分解的可分离非线性最小二乘问题求解方法及其应用
王路遥
山东科技大学测绘与空间信息学院, 山东 青岛 266590
基金项目：国家自然科学基金(42074009);山东省自然科学基金(ZR2020MD043;ZR2020MD044)
摘要：针对可分离非线性函数模型的特殊结构, 本文使用变量投影法(VP)将线性参数与非线性参数分离开来, 并分别与矩阵的满秩分解、QR分解、奇异值分解和施密特正交化相结合, 对两类参数分别求解, 缩短了计算机解算方程组的运算时间, 使算法更加高效, 同时也使得具有一定病态程度的方程组在解算过程中保持相对较好的稳定性。本文利用Mackey-Glass时间序列拟合试验和空间直角坐标转换参数解算试验对比分析了基于不同矩阵分解方法的算法优劣性。试验结果表明, 基于矩阵分解的改进变量投影法具有高效的运算效率与稳定的解算过程, 也适用于解算空间直角坐标转换参数问题。
关键词：可分离非线性最小二乘 变量投影 矩阵分解 Mackey-Glass时间序列 空间直角坐标转换
WANG Luyao, LIU Guolin, WANG Fengyun, et al. The method and application for solving separable nonlinear least squares based on matrix decomposition[J]. Acta Geodaetica et Cartographica Sinica, 2022, 51(11): 2317-2327. DOI: 10.11947/j.AGCS.2022.20200502
引 言
最小二乘估计在线性与非线性模型的参数估计理论与实际应用中都占有非常重要的地位[1]，而在数据处理中大多数实际问题都涉及非线性函数模型[2]。非线性最小二乘问题是最优化理论的重要分支，也可视为无约束条件下的二次规划问题。若要在非线性数据处理中广泛使用非线性模型参数估计理论，必须进行大量深入细致的工作[3]。文献[4]推导了不等式约束整体最小二乘(inequality constrained WTLS, ICWTLS)最优解的一阶必要条件和二阶充分条件，采用有效集(active set, AS)法和序列二次规划(sequential quadratic programming, SQP)法求ICWTLS解。序列二次规划法和有效集法等算法常常被应用于求解等式或不等式约束，以及混合约束的非线性规划问题。而针对可分离非线性最小二乘问题，Golub和Pereyra提出变量投影算法，此后众多学者对变量分离后再进行解算的方法进行了深入研究。在可分离问题中，应用变量投影算法并结合Gauss-Newton方法对模型参数进行求解，函数的计算效率会有一定程度的提高[5]。文献[6]证明了变量投影和联合优化是求解可分离非线性最小二乘问题的两种截然不同的方法，其最重要的区别在于前者的不平衡信赖域假设。文献[7]通过蒙特卡洛模拟试验验证了Ruano提出的简化雅可比矩阵JR=－DFF－y不适用于变量投影算法，这可能使算法难以收敛。此外，基于矩阵的梯形分解[8]、约束方程的广义LU(lower-upper, LU)分解[9]和奇异值分解[10]等分解方法，多种变量投影的改进算法也被相继提出。
经国内外学者长期的潜心钻研，变量投影算法在数值分析、化学工程、神经网络、电子通信与工程、环境科学和时间序列分析等领域得到了广泛的发展和应用。文献[11]将变量投影法用于有理函数或指数模型的非线性超定方程组求解的试验，并与参数未分离的方法进行对比，在利用公式J=－PF⊥DFF－y化简雅可比矩阵后，采用Gauss-Newton(GN)算法迭代计算非线性参数，结果表明在结果渐近收敛性质相同的情况下，简化雅可比矩阵方法的计算效率更高。文献[12]采用施密特正交分解方法对变量投影过程中的非线性参数矩阵进行求解，简化了伪逆的求解过程，并在RBF-AR模型的参数求解中得到验证。此外，变量投影方法也被应用在LiDAR测深波形参数确定[13]、未知源函数全波形的反演问题[14]、大规模不适定问题中解决图像处理中的失真问题[15]及GPS定位技术中估算由多路径效应引起的误差[16]等。学者们在不同科研领域对变量投影法进行了多种尝试，并且取得了较为丰富的研究成果。在空间直角坐标转换中，七参数模型符合可分离非线性模型一般形式，而许多坐标转换参数解算方法为避免求解复杂的非线性模型，通常是将非线性函数在近似值处泰勒展开到一次项进行线性化处理[17]，并采用迭代法以实现提高解算精度和简化计算过程[18]。常用的Gauss-Newton迭代法利用了目标函数的一阶导数近似二阶导数信息，表现出近似二次收敛的特性[19]，且能够以指定精度收敛[20]，在求解坐标转换参数中得到了很好的应用。测绘领域中存在较多的可分离非线性模型，但变量投影算法的应用相对较少，需进一步拓宽其应用范围。
在变量投影算法与矩阵分解相结合的研究文献中，学者们在经典算法的基础上优化算法来提高求解可分离问题的效率，但其中多数只针对某种单一的矩阵分解方法与变量投影算法结合的实用性进行了探究，缺少多种方法的对比分析，且常用的满秩分解方法未与变量投影算法相结合。在空间直角坐标转换参数解算中，七参数模型具有可分离非线性模型的一般形式，而变量投影算法并未应用其中。因此，本文在经典变量投影算法的基础上进行改进，对参数分离后的非线性矩阵进行满秩分解、QR分解、奇异值分解和施密特正交化方法分解处理，导出4种变量投影公式。在数值试验部分，通过Mackey-Glass时间序列模拟试验，验证变量投影算法与基于不同矩阵分解的改进算法的优缺点，并将变量投影算法及其改进算法应用于空间直角坐标转换参数解算中，采用基于Gauss-Newton法改进的Levenberg-Marquardt(LM)法迭代计算旋转参数。通过独立坐标系和CGCS2000坐标系转换试验，验证改进的变量投影算法在此类实际问题中的可行性，并对比不同改进算法解决此类问题的优劣性，为坐标转换模型参数解算提供一种新方法。
1 可分离非线性最小二乘及其变量投影算法
可分离非线性最小二乘问题的函数模型为非线性函数的线性组合形式，一般可表示为
2 基于矩阵分解的变量投影算法
2.1 满秩分解满秩分解是将矩阵表示成一个列满秩矩阵C与一个行满秩矩阵B的乘积形式。设矩阵A∈ℂrm×p(r>0)，则存在C∈ℂrm×r(r>0)，B∈ℂrr×p(r>0)，使得A=C×B，此为A的满秩分解。在变量投影算法中，将m×p阶矩阵Φ(α)通过满秩分解方法分解为3 数值试验
为了验证基于矩阵分解的变量投影算法在可分离非线性最小二乘问题求解的有效性及在测绘领域中的应用，分别采用Mackey-Glass时间序列模拟试验和独立坐标系与CGCS2000坐标系之间参数求解试验，对参数不分离方法(Nons)、经典变量投影法(VP)、基于满秩分解的变量投影法(VPCB)、基于QR分解的变量投影法(VPQR)、基于奇异值分解的变量投影法(VPSVD)和基于施密特正交化的变量投影法(VPGSO)进行了对比分析，试验环境为Matlab R2016a，1.80 GHz PC，Windows7系统。3.1 Mackey-Glass时间序列试验利用指数模型拟合混沌型Mackey-Glass时间序列，此模型在非线性时间序列预测中具有优良的性能，且已从理论上被证明可以任意精度逼近任何函数[24]。设函数模型为图 1 混沌型Mackey-Glass时间序列 Fig. 1 Chaotic Mackey-Glass time series 图选项
令指数模型中的p=5，则模型中有9个非线性参数和5个线性参数。利用区间[24, 123]内的数据进行训练，即将数据代入模型(23)形成非线性方程组，给定相同的迭代初值，使用LM迭代算法，限制最大迭代次数为500次，自变量误差限||αk+1－αk|| < 10－6，以及残差平方和误差限ssr(αk+1)－ssr(αk) < 10－6，对模型参数进行迭代求解。6种算法结合LM算法训练及预测数据的曲线拟合程度如图 2所示，其中，区间[24, 123]内的图像为训练曲线，区间[124, 200]内的图像为预测曲线，红色散点和黄色散点是从时间序列图像中提取的数据点，从图中可以判断VPCB+LM、VPSVD+LM对应曲线与数据拟合程度较低，其余算法拟合效果较好。图 2 曲线拟合 Fig. 2 Curve fitting 图选项
表 1列出了所有方法对应的迭代次数、函数计算次数、均方根误差及运算时间。图 3给出了迭代计算时，参数分离的5种算法对应的非线性参数的收敛过程。
表 1 迭代次数、函数计算次数、均方根误差(RMSE)、运算时间对比Tab. 1 Comparison of iteration times, function calculation times, RMSE and operation time
表选项
对比表 1中的数据可知，模型在短期预测中所得均方根误差比训练时要小，这与文献[24]的试验结果一致，说明模型能够从训练数据中获取较好的预测能力。在解算结果精确度方面，所有方法都具有足够小的均方根误差，其结果都有较高的可信度，其中VPQR+LM算法对应的训练和预测的均方根误差数值最小，说明此种方法在拟合数据和预测趋势的过程中，其结果最为精准，但其计算效率不占据优势；Nons+LM方法迭代计算次数最多，所需时间最长，而经过分离参数的变量投影算法在这两方面有了较大幅度的减少，本文提出的VPCB+LM具有最少的迭代次数；在迭代次数与函数计算次数水平相近的情况下，VPSVD+LM和VPGSO+LM方法在运算时间上有较为显著的优势。文献[25]分别利用200、400、600组模拟试验数据，在解算结果精度相近的条件下对比验证了基于施密特正交化的变量投影法的运算效率[25]，与本文所得结论相同。给定相同的迭代初值，参数分离的5种方法对应的非线性参数迭代过程如图 3所示。图 3 非线性参数迭代过程 Fig. 3 Iteration process of nonlinear parameters 图选项
由图 3可知，在给定相同初值的条件下，所有方法都有明显的收敛趋势，且VP+LM、VPCB+LM和VPGSO+LM的收敛速度较快。在区间[400, 500]内进行了多步预测试验，从Mackey-Glass时间序列中的第400个点开始，每10个点一组，利用5种分离参数算法求得的模型对时间序列进一步预测至第500个点，绘制出该区间内曲线拟合效果图，并求得每一组的均方根误差RMSE，进而绘制成折线图，如图 4所示。图 4 预测曲线拟合图及分组多步预测的RMSE Fig. 4 Curve fitting of prediction and RMSE of multi-step prediction 图选项
由图 4可知，VPCB+LM和VPSVD+LM算法随着预测次数的增加，其均方根误差越来越大且变化速率较快，表明其预测精度越来越低，预测能力较弱，而其他算法对应折线相对平缓，预测能力更加稳定，其结果具有较高的精度，其中VPQR+LM算法最稳定，预测效果最佳。3.2 空间直角坐标转换模型参数求解试验在测量学科的空间直角坐标转换中，两个不同的空间直角坐标系进行转换需要用7个参数来描述，其中包括3个坐标轴的平移参数、3个旋转参数和1个坐标轴缩放参数，如果已知两空间直角坐标系下的3个及以上的公共点坐标值，那么这7个参数可以唯一确定。下面对七参数模型进行介绍。如图 5所示，对于两个不同的空间直角坐标系O1X1Y1Z1和O2X2Y2Z2，它们的原点不一致，对应的坐标轴相互不平行。两个坐标系之间存在平移参数(ΔX, ΔY, ΔZ)、旋转参数(α, β, γ)和缩放参数s。图 5 坐标转换 Fig. 5 Coordinate transformation 图选项
在进行坐标转换时，常用的七参数模型的一般形式为
表 2 独立空间直角坐标系和CGCS2000坐标系下9组公共点坐标值Tab. 2 Nine sets of common point coordinate values in independent rectangular coordinate system and CGCS2000 coordinate system m
表选项
利用参数不分离算法(Nons+LM)、经典变量投影法(VP+LM)、基于满秩分解的变量投影法(VPCB+LM)、基于QR分解的变量投影法(VPQR+LM)、基于奇异值分解的变量投影法(VPSVD+LM)和基于施密特正交化的变量投影法(VPGSO+LM)求解模型中的旋转和缩放参数。由大量公共点求得参数的最优值为
表 3 参数解算结果的对比Tab. 3 Comparison of parameters calculation results
表选项
由表 3的结果可以看出，每种方法得到的参数值几乎相等；当计算每次迭代的残差平方和时，原始模型的残差实际上包含两种参数类型的误差，这也是在收敛分析时图 6(a)中Nons+LM方法在最开始时残差平方和较大的原因。参数不分离算法的参数估计值与最优值的残差平方和最小，其他方法也都具有较小的残差，得到的结果都是足够精确的。图 6 非线性参数的迭代过程及残差平方和变化过程 Fig. 6 Iteration process of nonlinear parameters and variation process of residual sum of squares 图选项
3.2.2 计算过程分析除了比较数值解算结果外，试验过程中还比较了每种方法的迭代次数、函数计数次数、残差平方和及运算时间，结果见表 4。
表 4 迭代次数、函数计算次数、残差平方和及运算时间的对比Tab. 4 Comparison of iteration times, function calculation times, residual sum of squares and operation time
表选项
由表 4可知，每种方法的迭代次数均相等，而经过参数分离的算法对应的函数计算次数均小于Nons+LM算法。Nons+LM对应的残差稍大而参数分离的方法相对较小，这表明相比于不分离而直接进行迭代的方法，变量投影算法在估计坐标转换参数时，能够减小观测值与转换模型之间的残差，使得参数计算结果更加可靠。由计算时间可以看出，参数分离的方法对应的时间相比分离前都有所减少，在解算结果精度相同的条件下，改进算法的计算速度有不同程度的提升，其中VPGSO+LM最为显著。试验利用了9组公共点求解转换参数，展开式对应形成了含有27个方程的方程组，经变量投影分离参数后，非线性函数构成27×3阶矩阵，当VPGSO+LM对目标函数进行分解时，施密特正交化中计算的矩阵G为27×3阶矩阵，比用其他方法得到的矩阵(如QR分解得到矩阵Q为27×27阶矩阵)更加简单一些，因此在计算机存储或调用此矩阵并代入函数进行计算时，时间会有所减少。在迭代次数相同、函数计算次数和残差平方和相近的情况下，VPGSO+LM方法相比于Nons+LM方法降低了估计参数维度，相比于其他改进的变量投影法又减少了计算时间。3.2.3 收敛过程的对比分析本节以试验中的旋转角α为例，分析了6种算法在迭代过程中的收敛性。由于所有方法都是基于LM算法进行计算的，因此参数分离的算法在迭代步长、方向和循环过程是基本一致的。在加入参数不分离算法进行对比时，只在图中显示了Nons+LM和VPGSO+LM对应的计算过程。Nons+LM算法的迭代过程中也包含线性参数，其迭代方向和迭代步长与参数分离的算法在数值上存在差异。参数的具体迭代过程及残差平方和的变化如图 6所示。由图 6可以看出，Nons+LM方法在初次计算时的残差要远大于VPGSO+LM，这是由于在迭代过程中线性与非线性参数统一进行计算所导致的。在参数α的迭代过程中，VPGSO+LM方法在第2次迭代过程中存在振荡，第3次迭代接近最优值，而Nons+LM方法相对稳定一些；图 6(c)显示了分离变量后的5种方法迭代过程，其过程比较一致，都是在第2次迭代出现震荡之后开始接近最优值，由于迭代过程差异较小，因此将第2次迭代结果放大至图 6(d)，可以看出VPSVD+LM算法结果较差，而VPCB+LM相比之下最稳定。在参数迭代过程中，LM迭代法计算的是所有非线性参数每次迭代后的总体变化，因此，迭代过程中某个参数的变化出现抬升或下降并不意味着算法不好。一般采用迭代法求解非线性方程组往往只是局部收敛，这使得迭代初值很难确定[28]。当随机给定不同的初值时，分别用分离参数的5种算法估计非线性转换参数，将程序独立运行30次，计算得到每次结果的均方根误差RMSE并绘制箱形图，如图 7所示。由图 7可知，VPSVD+LM对应的RMSE数值的中位数最小，对应的箱形图所表示的数据整体分布区域更接近横轴，由此可见，VPSVD+LM算法受初值影响较小，程序运行得更加稳定一些。图 7 不同算法对应均方根误差分布状况 Fig. 7 Distribution of RMSE corresponding to different algorithms 图选项
4 结论
本文结合混沌型Mackey-Glass时间序列模拟试验，全面地对比分析了变量投影法及其基于矩阵分解的4种改进算法的优缺点。参数不分离方法在解算模型时，迭代计算次数较多，运算时间最长，而分离变量的VP算法在保证解算结果精度的前提下，可有效削弱这一弊端；当非线性矩阵存在一定程度的病态时，对矩阵进行分解可使得算法迭代过程稳定，以保障解算的有效性；改进的算法在迭代次数、解算精度、计算效率、预测能力及稳定性等方面也有不同程度的提升。在解算可分离问题的模型参数时，若参数最优值所在区间不可知，在不追求运算效率的情况下可选用VPQR+LM，若考虑运算效率而对精度要求较低，则建议使用VPCB+LM或VPSVD+LM，若两者都须顾及则可选用VPGSO+LM；若已知其最优值所在区间，可使用VPGSO+LM算法并设置迭代初值在最优值附近，以追求参数求解的计算效率；在解得模型参数后，可使用VPQR+LM进行数据预测。在基于矩阵分解推导变量投影公式的过程中发现，当非线性矩阵Φ(α)秩亏时满秩分解具有较大优势，其余3种矩阵分解法最终都将非线性矩阵及其逆阵的乘积PΦ(α)导出为正交矩阵与含有单位阵矩阵相乘，其形式具有统一性，而正交矩阵的逆矩阵与其转置相同这一性质使得算法在求解模型参数时简化了计算过程，节约了运算时间。本文将变量投影法及改进算法应用于空间直角坐标转换模型参数求解中，并对算法的实用性进行了探究。经试验发现，该算法及改进算法适用于解算坐标转换模型参数，且4种改进算法在解决此类问题的运算效率方面有所提升。在转换参数求解中应用改进的变量投影算法，可在保证求解结果精度的同时降低待估计参数的维数，简化非线性函数矩阵大小，提高计算机的变量存储、调用和计算的速度，使运算时间最小化，为坐标转换模型参数的求解提供了新的方法，对丰富可分离非线性最小二乘理论在测绘领域中的应用具有实际意义。
引 言
最小二乘估计在线性与非线性模型的参数估计理论与实际应用中都占有非常重要的地位[1]，而在数据处理中大多数实际问题都涉及非线性函数模型[2]。非线性最小二乘问题是最优化理论的重要分支，也可视为无约束条件下的二次规划问题。若要在非线性数据处理中广泛使用非线性模型参数估计理论，必须进行大量深入细致的工作[3]。文献[4]推导了不等式约束整体最小二乘(inequality constrained WTLS, ICWTLS)最优解的一阶必要条件和二阶充分条件，采用有效集(active set, AS)法和序列二次规划(sequential quadratic programming, SQP)法求ICWTLS解。序列二次规划法和有效集法等算法常常被应用于求解等式或不等式约束，以及混合约束的非线性规划问题。而针对可分离非线性最小二乘问题，Golub和Pereyra提出变量投影算法，此后众多学者对变量分离后再进行解算的方法进行了深入研究。在可分离问题中，应用变量投影算法并结合Gauss-Newton方法对模型参数进行求解，函数的计算效率会有一定程度的提高[5]。文献[6]证明了变量投影和联合优化是求解可分离非线性最小二乘问题的两种截然不同的方法，其最重要的区别在于前者的不平衡信赖域假设。文献[7]通过蒙特卡洛模拟试验验证了Ruano提出的简化雅可比矩阵JR=－DFF－y不适用于变量投影算法，这可能使算法难以收敛。此外，基于矩阵的梯形分解[8]、约束方程的广义LU(lower-upper, LU)分解[9]和奇异值分解[10]等分解方法，多种变量投影的改进算法也被相继提出。
经国内外学者长期的潜心钻研，变量投影算法在数值分析、化学工程、神经网络、电子通信与工程、环境科学和时间序列分析等领域得到了广泛的发展和应用。文献[11]将变量投影法用于有理函数或指数模型的非线性超定方程组求解的试验，并与参数未分离的方法进行对比，在利用公式J=－PF⊥DFF－y化简雅可比矩阵后，采用Gauss-Newton(GN)算法迭代计算非线性参数，结果表明在结果渐近收敛性质相同的情况下，简化雅可比矩阵方法的计算效率更高。文献[12]采用施密特正交分解方法对变量投影过程中的非线性参数矩阵进行求解，简化了伪逆的求解过程，并在RBF-AR模型的参数求解中得到验证。此外，变量投影方法也被应用在LiDAR测深波形参数确定[13]、未知源函数全波形的反演问题[14]、大规模不适定问题中解决图像处理中的失真问题[15]及GPS定位技术中估算由多路径效应引起的误差[16]等。学者们在不同科研领域对变量投影法进行了多种尝试，并且取得了较为丰富的研究成果。在空间直角坐标转换中，七参数模型符合可分离非线性模型一般形式，而许多坐标转换参数解算方法为避免求解复杂的非线性模型，通常是将非线性函数在近似值处泰勒展开到一次项进行线性化处理[17]，并采用迭代法以实现提高解算精度和简化计算过程[18]。常用的Gauss-Newton迭代法利用了目标函数的一阶导数近似二阶导数信息，表现出近似二次收敛的特性[19]，且能够以指定精度收敛[20]，在求解坐标转换参数中得到了很好的应用。测绘领域中存在较多的可分离非线性模型，但变量投影算法的应用相对较少，需进一步拓宽其应用范围。
在变量投影算法与矩阵分解相结合的研究文献中，学者们在经典算法的基础上优化算法来提高求解可分离问题的效率，但其中多数只针对某种单一的矩阵分解方法与变量投影算法结合的实用性进行了探究，缺少多种方法的对比分析，且常用的满秩分解方法未与变量投影算法相结合。在空间直角坐标转换参数解算中，七参数模型具有可分离非线性模型的一般形式，而变量投影算法并未应用其中。因此，本文在经典变量投影算法的基础上进行改进，对参数分离后的非线性矩阵进行满秩分解、QR分解、奇异值分解和施密特正交化方法分解处理，导出4种变量投影公式。在数值试验部分，通过Mackey-Glass时间序列模拟试验，验证变量投影算法与基于不同矩阵分解的改进算法的优缺点，并将变量投影算法及其改进算法应用于空间直角坐标转换参数解算中，采用基于Gauss-Newton法改进的Levenberg-Marquardt(LM)法迭代计算旋转参数。通过独立坐标系和CGCS2000坐标系转换试验，验证改进的变量投影算法在此类实际问题中的可行性，并对比不同改进算法解决此类问题的优劣性，为坐标转换模型参数解算提供一种新方法。
1 可分离非线性最小二乘及其变量投影算法
可分离非线性最小二乘问题的函数模型为非线性函数的线性组合形式，一般可表示为1 可分离非线性最小二乘及其变量投影算法
2 基于矩阵分解的变量投影算法
2.1 满秩分解满秩分解是将矩阵表示成一个列满秩矩阵C与一个行满秩矩阵B的乘积形式。设矩阵A∈ℂrm×p(r>0)，则存在C∈ℂrm×r(r>0)，B∈ℂrr×p(r>0)，使得A=C×B，此为A的满秩分解。在变量投影算法中，将m×p阶矩阵Φ(α)通过满秩分解方法分解为2 基于矩阵分解的变量投影算法
3 数值试验
为了验证基于矩阵分解的变量投影算法在可分离非线性最小二乘问题求解的有效性及在测绘领域中的应用，分别采用Mackey-Glass时间序列模拟试验和独立坐标系与CGCS2000坐标系之间参数求解试验，对参数不分离方法(Nons)、经典变量投影法(VP)、基于满秩分解的变量投影法(VPCB)、基于QR分解的变量投影法(VPQR)、基于奇异值分解的变量投影法(VPSVD)和基于施密特正交化的变量投影法(VPGSO)进行了对比分析，试验环境为Matlab R2016a，1.80 GHz PC，Windows7系统。3.1 Mackey-Glass时间序列试验利用指数模型拟合混沌型Mackey-Glass时间序列，此模型在非线性时间序列预测中具有优良的性能，且已从理论上被证明可以任意精度逼近任何函数[24]。设函数模型为令指数模型中的p=5，则模型中有9个非线性参数和5个线性参数。利用区间[24, 123]内的数据进行训练，即将数据代入模型(23)形成非线性方程组，给定相同的迭代初值，使用LM迭代算法，限制最大迭代次数为500次，自变量误差限||αk+1－αk|| < 10－6，以及残差平方和误差限ssr(αk+1)－ssr(αk) < 10－6，对模型参数进行迭代求解。6种算法结合LM算法训练及预测数据的曲线拟合程度如图 2所示，其中，区间[24, 123]内的图像为训练曲线，区间[124, 200]内的图像为预测曲线，红色散点和黄色散点是从时间序列图像中提取的数据点，从图中可以判断VPCB+LM、VPSVD+LM对应曲线与数据拟合程度较低，其余算法拟合效果较好。
表 1列出了所有方法对应的迭代次数、函数计算次数、均方根误差及运算时间。图 3给出了迭代计算时，参数分离的5种算法对应的非线性参数的收敛过程。
对比表 1中的数据可知，模型在短期预测中所得均方根误差比训练时要小，这与文献[24]的试验结果一致，说明模型能够从训练数据中获取较好的预测能力。在解算结果精确度方面，所有方法都具有足够小的均方根误差，其结果都有较高的可信度，其中VPQR+LM算法对应的训练和预测的均方根误差数值最小，说明此种方法在拟合数据和预测趋势的过程中，其结果最为精准，但其计算效率不占据优势；Nons+LM方法迭代计算次数最多，所需时间最长，而经过分离参数的变量投影算法在这两方面有了较大幅度的减少，本文提出的VPCB+LM具有最少的迭代次数；在迭代次数与函数计算次数水平相近的情况下，VPSVD+LM和VPGSO+LM方法在运算时间上有较为显著的优势。文献[25]分别利用200、400、600组模拟试验数据，在解算结果精度相近的条件下对比验证了基于施密特正交化的变量投影法的运算效率[25]，与本文所得结论相同。给定相同的迭代初值，参数分离的5种方法对应的非线性参数迭代过程如图 3所示。
由图 3可知，在给定相同初值的条件下，所有方法都有明显的收敛趋势，且VP+LM、VPCB+LM和VPGSO+LM的收敛速度较快。在区间[400, 500]内进行了多步预测试验，从Mackey-Glass时间序列中的第400个点开始，每10个点一组，利用5种分离参数算法求得的模型对时间序列进一步预测至第500个点，绘制出该区间内曲线拟合效果图，并求得每一组的均方根误差RMSE，进而绘制成折线图，如图 4所示。
由图 4可知，VPCB+LM和VPSVD+LM算法随着预测次数的增加，其均方根误差越来越大且变化速率较快，表明其预测精度越来越低，预测能力较弱，而其他算法对应折线相对平缓，预测能力更加稳定，其结果具有较高的精度，其中VPQR+LM算法最稳定，预测效果最佳。3.2 空间直角坐标转换模型参数求解试验在测量学科的空间直角坐标转换中，两个不同的空间直角坐标系进行转换需要用7个参数来描述，其中包括3个坐标轴的平移参数、3个旋转参数和1个坐标轴缩放参数，如果已知两空间直角坐标系下的3个及以上的公共点坐标值，那么这7个参数可以唯一确定。下面对七参数模型进行介绍。如图 5所示，对于两个不同的空间直角坐标系O1X1Y1Z1和O2X2Y2Z2，它们的原点不一致，对应的坐标轴相互不平行。两个坐标系之间存在平移参数(ΔX, ΔY, ΔZ)、旋转参数(α, β, γ)和缩放参数s。
在进行坐标转换时，常用的七参数模型的一般形式为
利用参数不分离算法(Nons+LM)、经典变量投影法(VP+LM)、基于满秩分解的变量投影法(VPCB+LM)、基于QR分解的变量投影法(VPQR+LM)、基于奇异值分解的变量投影法(VPSVD+LM)和基于施密特正交化的变量投影法(VPGSO+LM)求解模型中的旋转和缩放参数。由大量公共点求得参数的最优值为
由表 3的结果可以看出，每种方法得到的参数值几乎相等；当计算每次迭代的残差平方和时，原始模型的残差实际上包含两种参数类型的误差，这也是在收敛分析时图 6(a)中Nons+LM方法在最开始时残差平方和较大的原因。参数不分离算法的参数估计值与最优值的残差平方和最小，其他方法也都具有较小的残差，得到的结果都是足够精确的。
表选项
由表 4可知，每种方法的迭代次数均相等，而经过参数分离的算法对应的函数计算次数均小于Nons+LM算法。Nons+LM对应的残差稍大而参数分离的方法相对较小，这表明相比于不分离而直接进行迭代的方法，变量投影算法在估计坐标转换参数时，能够减小观测值与转换模型之间的残差，使得参数计算结果更加可靠。由计算时间可以看出，参数分离的方法对应的时间相比分离前都有所减少，在解算结果精度相同的条件下，改进算法的计算速度有不同程度的提升，其中VPGSO+LM最为显著。试验利用了9组公共点求解转换参数，展开式对应形成了含有27个方程的方程组，经变量投影分离参数后，非线性函数构成27×3阶矩阵，当VPGSO+LM对目标函数进行分解时，施密特正交化中计算的矩阵G为27×3阶矩阵，比用其他方法得到的矩阵(如QR分解得到矩阵Q为27×27阶矩阵)更加简单一些，因此在计算机存储或调用此矩阵并代入函数进行计算时，时间会有所减少。在迭代次数相同、函数计算次数和残差平方和相近的情况下，VPGSO+LM方法相比于Nons+LM方法降低了估计参数维度，相比于其他改进的变量投影法又减少了计算时间。3.2.3 收敛过程的对比分析本节以试验中的旋转角α为例，分析了6种算法在迭代过程中的收敛性。由于所有方法都是基于LM算法进行计算的，因此参数分离的算法在迭代步长、方向和循环过程是基本一致的。在加入参数不分离算法进行对比时，只在图中显示了Nons+LM和VPGSO+LM对应的计算过程。Nons+LM算法的迭代过程中也包含线性参数，其迭代方向和迭代步长与参数分离的算法在数值上存在差异。参数的具体迭代过程及残差平方和的变化如图 6所示。由图 6可以看出，Nons+LM方法在初次计算时的残差要远大于VPGSO+LM，这是由于在迭代过程中线性与非线性参数统一进行计算所导致的。在参数α的迭代过程中，VPGSO+LM方法在第2次迭代过程中存在振荡，第3次迭代接近最优值，而Nons+LM方法相对稳定一些；图 6(c)显示了分离变量后的5种方法迭代过程，其过程比较一致，都是在第2次迭代出现震荡之后开始接近最优值，由于迭代过程差异较小，因此将第2次迭代结果放大至图 6(d)，可以看出VPSVD+LM算法结果较差，而VPCB+LM相比之下最稳定。在参数迭代过程中，LM迭代法计算的是所有非线性参数每次迭代后的总体变化，因此，迭代过程中某个参数的变化出现抬升或下降并不意味着算法不好。一般采用迭代法求解非线性方程组往往只是局部收敛，这使得迭代初值很难确定[28]。当随机给定不同的初值时，分别用分离参数的5种算法估计非线性转换参数，将程序独立运行30次，计算得到每次结果的均方根误差RMSE并绘制箱形图，如图 7所示。由图 7可知，VPSVD+LM对应的RMSE数值的中位数最小，对应的箱形图所表示的数据整体分布区域更接近横轴，由此可见，VPSVD+LM算法受初值影响较小，程序运行得更加稳定一些。
4 结论
本文结合混沌型Mackey-Glass时间序列模拟试验，全面地对比分析了变量投影法及其基于矩阵分解的4种改进算法的优缺点。参数不分离方法在解算模型时，迭代计算次数较多，运算时间最长，而分离变量的VP算法在保证解算结果精度的前提下，可有效削弱这一弊端；当非线性矩阵存在一定程度的病态时，对矩阵进行分解可使得算法迭代过程稳定，以保障解算的有效性；改进的算法在迭代次数、解算精度、计算效率、预测能力及稳定性等方面也有不同程度的提升。在解算可分离问题的模型参数时，若参数最优值所在区间不可知，在不追求运算效率的情况下可选用VPQR+LM，若考虑运算效率而对精度要求较低，则建议使用VPCB+LM或VPSVD+LM，若两者都须顾及则可选用VPGSO+LM；若已知其最优值所在区间，可使用VPGSO+LM算法并设置迭代初值在最优值附近，以追求参数求解的计算效率；在解得模型参数后，可使用VPQR+LM进行数据预测。在基于矩阵分解推导变量投影公式的过程中发现，当非线性矩阵Φ(α)秩亏时满秩分解具有较大优势，其余3种矩阵分解法最终都将非线性矩阵及其逆阵的乘积PΦ(α)导出为正交矩阵与含有单位阵矩阵相乘，其形式具有统一性，而正交矩阵的逆矩阵与其转置相同这一性质使得算法在求解模型参数时简化了计算过程，节约了运算时间。本文将变量投影法及改进算法应用于空间直角坐标转换模型参数求解中，并对算法的实用性进行了探究。经试验发现，该算法及改进算法适用于解算坐标转换模型参数，且4种改进算法在解决此类问题的运算效率方面有所提升。在转换参数求解中应用改进的变量投影算法，可在保证求解结果精度的同时降低待估计参数的维数，简化非线性函数矩阵大小，提高计算机的变量存储、调用和计算的速度，使运算时间最小化，为坐标转换模型参数的求解提供了新的方法，对丰富可分离非线性最小二乘理论在测绘领域中的应用具有实际意义。4 结论
作者简介
第一作者简介：王路遥(1994—), 男, 博士生, 研究领域为现代测量数据处理。E-mail: 541355164@qq.com
通信作者：刘国林, E-mail：gliu@sdust.edu.cn
终审：金 君
往期推荐
资讯
○ 测绘学报 | 王乐洋：加乘性混合误差模型精度评定的SUT法
○ 湖北大学2023年“沙湖论坛”暨全球人才云聘会重磅开启！