留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

非线性参数估计的自适应松弛正则化算法

曲国庆 孙振 苏晓庆 杜存鹏

曲国庆, 孙振, 苏晓庆, 杜存鹏. 非线性参数估计的自适应松弛正则化算法[J]. 武汉大学学报 ● 信息科学版, 2019, 44(10): 1491-1497. doi: 10.13203/j.whugis20170425
引用本文: 曲国庆, 孙振, 苏晓庆, 杜存鹏. 非线性参数估计的自适应松弛正则化算法[J]. 武汉大学学报 ● 信息科学版, 2019, 44(10): 1491-1497. doi: 10.13203/j.whugis20170425
QU Guoqing, SUN Zhen, SU Xiaoqing, DU Cunpeng. Adaptive Relaxation Regularization Algorithm for Nonlinear Parameter Estimation[J]. Geomatics and Information Science of Wuhan University, 2019, 44(10): 1491-1497. doi: 10.13203/j.whugis20170425
Citation: QU Guoqing, SUN Zhen, SU Xiaoqing, DU Cunpeng. Adaptive Relaxation Regularization Algorithm for Nonlinear Parameter Estimation[J]. Geomatics and Information Science of Wuhan University, 2019, 44(10): 1491-1497. doi: 10.13203/j.whugis20170425

非线性参数估计的自适应松弛正则化算法

doi: 10.13203/j.whugis20170425
基金项目: 

国家重点研发计划 2016YFB051700

国家自然科学基金 41674014

国家自然科学基金 41704003

详细信息
    作者简介:

    曲国庆, 教授, 主要从事大地测量数据处理研究。qgq@sdut.edu.cn

    通讯作者: 孙振, 硕士。986254102@qq.com
  • 中图分类号: P228.4

Adaptive Relaxation Regularization Algorithm for Nonlinear Parameter Estimation

Funds: 

The National High-tech Research and Development Program of China 2016YFB051700

the National Natural Science Foundation of China 41674014

the National Natural Science Foundation of China 41704003

More Information
    Author Bio:

    QU Guoqing, professor, specializes in the geodesy data processing, qgq@sdut.edu.cn

    Corresponding author: SUN Zhen, postgraduate. E-mail:986254102@qq.com
图(4) / 表(3)
计量
  • 文章访问数:  1309
  • HTML全文浏览量:  147
  • PDF下载量:  110
  • 被引次数: 0
出版历程
  • 收稿日期:  2018-12-17
  • 刊出日期:  2019-10-05

非线性参数估计的自适应松弛正则化算法

doi: 10.13203/j.whugis20170425
    基金项目:

    国家重点研发计划 2016YFB051700

    国家自然科学基金 41674014

    国家自然科学基金 41704003

    作者简介:

    曲国庆, 教授, 主要从事大地测量数据处理研究。qgq@sdut.edu.cn

    通讯作者: 孙振, 硕士。986254102@qq.com
  • 中图分类号: P228.4

摘要: 非线性方程参数估计存在的弊端在于非线性观测方程存在不适定问题时,以线性化平差估计和高斯牛顿为代表的经典数值算法会产生较强的不稳定特征。因此,针对传统非线性最小二乘求解不稳定且可靠性低的特点,基于稳定泛函极小准则最优化思想,提出了一种自适应松弛正则化数值算法。该算法采用正则化参数几何递增计算方法和残差最小步长准则,实现了正则参数和迭代步长计算的完全自适应,提高了非线性迭代收敛效率。以病态仿真数据和水下实测数据为例,验证了该方法的数值收敛解优于线性平差估计解,收敛效率优于迭代Tikhonov正则化方法。

English Abstract

曲国庆, 孙振, 苏晓庆, 杜存鹏. 非线性参数估计的自适应松弛正则化算法[J]. 武汉大学学报 ● 信息科学版, 2019, 44(10): 1491-1497. doi: 10.13203/j.whugis20170425
引用本文: 曲国庆, 孙振, 苏晓庆, 杜存鹏. 非线性参数估计的自适应松弛正则化算法[J]. 武汉大学学报 ● 信息科学版, 2019, 44(10): 1491-1497. doi: 10.13203/j.whugis20170425
QU Guoqing, SUN Zhen, SU Xiaoqing, DU Cunpeng. Adaptive Relaxation Regularization Algorithm for Nonlinear Parameter Estimation[J]. Geomatics and Information Science of Wuhan University, 2019, 44(10): 1491-1497. doi: 10.13203/j.whugis20170425
Citation: QU Guoqing, SUN Zhen, SU Xiaoqing, DU Cunpeng. Adaptive Relaxation Regularization Algorithm for Nonlinear Parameter Estimation[J]. Geomatics and Information Science of Wuhan University, 2019, 44(10): 1491-1497. doi: 10.13203/j.whugis20170425
  • 大地测量观测模型大部分是非线性模型[1],其模型参数解算一般采用线性化最小二乘法。线性化最小二乘的数学基础为非线性函数模型的一阶泰勒函数逼近。但当参数初值误差较大或模型非线性强度较大时,非线性模型的线性近似有可能会导致模型误差过大而致解算结果无效,此时,若同时伴随出现不适定问题,会进一步降低参数估计的稳定性和可靠性[2-3]。因此,根据问题的秩亏或病态程度选择合适的非线性平差计算模型就显得尤为重要。对于非线性平差模型,常用固有曲率和参数效应来刻画非线性模型的非线性强度,进而评估参数估计的统计推断效果[4]

    非线性参数估计一般需要进行迭代计算,这就要求数值迭代算法快速、稳定地收敛到预期解,许多学者通过扰动分析进行了相关研究,指出非线性扰动主要来源于线性近似时系数矩阵的扰动、附加的截断误差及正交过程[5]。当非线性模型强度较高时,初始值的选取及近似正交都可能导致线性近似引起的系数矩阵扰动和截断误差超出容许范围[6],特别是在非线性模型病态时,迭代计算的不稳定将会加剧使平差计算结果失真[7]。为解决非线性最小二乘不适定问题,许多学者对其数值算法进行了大量风险评估,发展了一系列数值迭代法,如阻尼最小二乘法、非线性最小二乘同伦法和正则化牛顿法等,目的是加速非线性最小二乘迭代速度,避免迭代矩阵求逆或者改善迭代计算稳定性[8]

    无约束化最小二乘解易受到线性化初值和模型病态的影响,解不稳定且可靠性低;高斯牛顿迭代公式一定程度上降低了线性化残差的影响,但解的不适用性依然存在。为此,本文基于稳定泛函极小准则最优化思想,提出了一种自适应松弛正则化数值方法。并以仿真数据和实测数据验证了该方法的有效性和实用性。

    • 非线性模型广泛存在于大地测量数据处理和地球物理参数反演中。有效处理非线性模型是现代测量数据处理研究的难点和热点之一[9-10]。设观测方程的向量表达式为:

      $$ {\mathit{\boldsymbol{L = f}}}({\mathit{\boldsymbol{x}}}){\mathit{\boldsymbol{ + \varepsilon }}} $$ (1)

      式中,f(x):RmRnx的非线性映射函数,且f(x)=[f1(x), f2(x)…fn(x)]Txm×1未知参数向量;ε=[ε1, ε2εn]T为观测误差;L=[L1, L2Ln]T为观测数据,且互不相关。

      当模型非线性强度较低时,可将非线性模型线性化近似。对于较低精度的初值,还需要进行迭代计算,常用迭代算法包括牛顿法、高斯牛顿法、信赖域法、最速下降法、改进的高斯牛顿法[11-12]。上述非线性模型的线性化近似可表示为:

      $$ {\mathit{\boldsymbol{f}}}({\mathit{\boldsymbol{x}}}) \approx {\mathit{\boldsymbol{f}}}({{\mathit{\boldsymbol{x}}}_0}){\mathit{\boldsymbol{ + J}}}({{\mathit{\boldsymbol{x}}}_0}){\text{d}}{{\mathit{\boldsymbol{x}}}_0} $$ (2)

      式中,x0为线性化初值;dx0为参数改正向量;

      $$ \mathit{\boldsymbol{J}}(\mathit{\boldsymbol{x}}) = (\frac{{\partial \mathit{\boldsymbol{f}}(\mathit{\boldsymbol{x}})}}{{\partial {\mathit{\boldsymbol{x}}_1}}}, \frac{{\partial \mathit{\boldsymbol{f}}(\mathit{\boldsymbol{x}})}}{{\partial {\mathit{\boldsymbol{x}}_2}}}, \frac{{\partial \mathit{\boldsymbol{f}}(\mathit{\boldsymbol{x}})}}{{\partial {\mathit{\boldsymbol{x}}_3}}}, \cdots , \frac{{\partial \mathit{\boldsymbol{f}}(\mathit{\boldsymbol{x}})}}{{\partial {\mathit{\boldsymbol{x}}_m}}}) $$ (3)

      式(3)称为Jacobian矩阵。由式(1)和式(2)联立可得线性误差方程:

      $$ {\mathit{\boldsymbol{V}}}({{\mathit{\boldsymbol{x}}}_0}) = {\mathit{\boldsymbol{J}}}({{\mathit{\boldsymbol{x}}}_0}){\text{d}}{{\mathit{\boldsymbol{x}}}_0} - {\mathit{\boldsymbol{l}}} $$ (4)

      式中,l=L-f(x0)。上述方程的最小二乘解使残差在平方意义上取极小值,即

      $$ \varphi ({\text{d}}{{\mathit{\boldsymbol{x}}}_0}) = {{\mathit{\boldsymbol{V}}}^T}({{\mathit{\boldsymbol{x}}}_0}){\mathit{\boldsymbol{V}}}({{\mathit{\boldsymbol{x}}}_0}) = \min $$ (5)

      式中,φ:RnR称为最小二乘目标函数。

      运用无约束优化解法,最小二乘解满足下述正交条件:

      $$ {\mathit{\boldsymbol{h}}}({\text{d}}{{\mathit{\boldsymbol{x}}}_0}) = \varphi '({\text{d}}{{\mathit{\boldsymbol{x}}}_0}) = 2{{\mathit{\boldsymbol{J}}}^T}({{\mathit{\boldsymbol{x}}}_0}){\mathit{\boldsymbol{V}}}({{\mathit{\boldsymbol{x}}}_0}) = 0 $$ (6)

      可采用高斯牛顿法求解正交条件方程式(6),其迭代公式为:

      $$ {{\mathit{\boldsymbol{x}}}_{k + 1}} = {{\mathit{\boldsymbol{x}}}_k} + {\text{d}}{{\mathit{\boldsymbol{\hat x}}}_k} = {\mathit{\boldsymbol{x}}}{}_k + {{\mathit{\boldsymbol{N}}}^{ - 1}}({{\mathit{\boldsymbol{x}}}_k}){{\mathit{\boldsymbol{J}}}^T}({{\mathit{\boldsymbol{x}}}_k}){{\mathit{\boldsymbol{l}}}_k} $$ (7)

      式中,N(xk)=JT(xk)J(xk)称为设计矩阵; lk=L-f(xk)。

      高斯牛顿法的收敛条件相对于牛顿法更加严格。对于一般非线性最小二乘问题,Teunissen利用几何性质证明高斯牛顿法适应条件取决于残差的大小和非线性强度[13]。当残差较大或者设计矩阵存在病态时,此时高斯牛顿法易陷入死循环或者得到的解并非有效解。

    • 非线性最小二乘参数平差实质是求泛函极小值的最优化问题[14]。当利用无约束优化解算时,其优化问题可描述为:

      $$ \left\{ \begin{array}{l} \varphi ({\text{d}}{{\mathit{\boldsymbol{x}}}_0}) = {{\mathit{\boldsymbol{V}}}^T}({{\mathit{\boldsymbol{x}}}_0}){\mathit{\boldsymbol{V}}}({{\mathit{\boldsymbol{x}}}_0}) = \min \\ {\text{s}}{\text{.t}}{\text{. }}{\mathit{\boldsymbol{V}}}({x_0}) = {\mathit{\boldsymbol{J}}}({{\mathit{\boldsymbol{x}}}_0}){\text{d}}{{\mathit{\boldsymbol{x}}}_0} - {\mathit{\boldsymbol{l}}} \\ \end{array} \right. $$ (8)

      非线性方程组依据泰勒级数线性化后,直接采用无约束化思想解算易受到线性化误差和方程组态性的影响。若方程组Jacobian矩阵J(x0)秩亏或病态严重,将导致经典非线性模型线性估计方法解算结果不唯一、不稳定甚至求解失败。为了改善Jacobian矩阵秩亏或病态对参数估值的影响,可引入稳定泛函极小准则,将上述最优化问题转化为:

      $$ \left\{ \begin{array}{l} s({\rm{d}}{{\mathit{\boldsymbol{x}}}_0}) = {{\mathit{\boldsymbol{V}}}^T}({{\mathit{\boldsymbol{x}}}_0}){\mathit{\boldsymbol{V}}}({{\mathit{\boldsymbol{x}}}_0}) + \alpha {\rm{d}}{\mathit{\boldsymbol{x}}}_0^{\rm{T}}{\mathit{\boldsymbol{R}}}({{\mathit{\boldsymbol{x}}}_0}){\rm{d}}{{\mathit{\boldsymbol{x}}}_0} = \min \\ {\rm{s}}{\rm{.t}}{\rm{. }}{\mathit{\boldsymbol{V}}}({{\mathit{\boldsymbol{x}}}_0}) = {\mathit{\boldsymbol{J}}}({{\mathit{\boldsymbol{x}}}_0}){\rm{d}}{{\mathit{\boldsymbol{x}}}_0} - {\mathit{\boldsymbol{l}}} \end{array} \right. $$ (9)

      式中,α为正则化参数,用来平衡解的不稳定性和光滑性;R (x0)为正则化矩阵,代表着未知参数事先已知的先验信息,主要用来抑制伪目标,降低小特征值对应的奇异分量对数值解的影响[15-16]。若正则化矩阵R(x0)=I,则

      $$ s({\rm{d}}{{\mathit{\boldsymbol{x}}}_0}) = {{\mathit{\boldsymbol{V}}}^T}({{\mathit{\boldsymbol{x}}}_0}){\mathit{\boldsymbol{V}}}({{\mathit{\boldsymbol{x}}}_0}) + \alpha {\rm{d}}{\mathit{\boldsymbol{x}}}_0^{\rm{T}}{\rm{d}}{{\mathit{\boldsymbol{x}}}_0} = \min $$ (10)

      称为标准稳定泛函极小准则。将标准稳定泛函目标函数对待估计量dx0求一阶偏导,并令其等于0,即

      $$ \frac{{\partial s({\rm{d}}{{\mathit{\boldsymbol{x}}}_0})}}{{\partial {\rm{d}}{{\mathit{\boldsymbol{x}}}_0}}} = 2{{\mathit{\boldsymbol{J}}}^T}({{\mathit{\boldsymbol{x}}}_0}){\mathit{\boldsymbol{V}}}({{\mathit{\boldsymbol{x}}}_0}) + 2\alpha {\rm{d}}{\mathit{\boldsymbol{x}}}_0^{\rm{T}} = 0 $$ (11)

      联合式(4)对式(11)化简为:

      $$ ({\mathit{\boldsymbol{N}}}({{\mathit{\boldsymbol{x}}}_0}) + \alpha {\mathit{\boldsymbol{I}}}{\rm{)d}}{\mathit{\boldsymbol{\hat x}}} - {{\mathit{\boldsymbol{J}}}^T}({{\mathit{\boldsymbol{x}}}_0}){\mathit{\boldsymbol{l}}} = 0 $$ (12)

      式(12)的Tikhonov正则化解为:

      $$ {\rm{d}}{\mathit{\boldsymbol{\hat x}}} = {({\mathit{\boldsymbol{N}}}({{\mathit{\boldsymbol{x}}}_0}) + \alpha {\mathit{\boldsymbol{I}}})^{ - 1}}{{\mathit{\boldsymbol{J}}}^T}({{\mathit{\boldsymbol{x}}}_0}){\mathit{\boldsymbol{l}}} $$ (13)

      考虑到线性化误差和初值精度,根据式(13)可得迭代Tikhonov正则化公式:

      $$ {{\mathit{\boldsymbol{x}}}_{k + 1}} = {{\mathit{\boldsymbol{x}}}_k} + {({\mathit{\boldsymbol{N}}}({{\mathit{\boldsymbol{x}}}_k}) + \alpha {\mathit{\boldsymbol{I}}})^{ - 1}}{{\mathit{\boldsymbol{J}}}^T}({{\mathit{\boldsymbol{x}}}_k}){{\mathit{\boldsymbol{l}}}_k} $$ (14)

      由式(14)可知,迭代过程中各分量变化而正则化参数不变,由于数据扰动影响,会降低正则化参数的作用效果,从而导致结果改善有限。为此,本文对式(14)进行如下改进:

      $$ {{\mathit{\boldsymbol{x}}}_{k + 1}} = {{\mathit{\boldsymbol{x}}}_k} + {\beta _k}{({\mathit{\boldsymbol{N}}}({{\mathit{\boldsymbol{x}}}_k}) + \alpha {\mathit{\boldsymbol{I}}})^{ - 1}}{{\mathit{\boldsymbol{J}}}^T}({{\mathit{\boldsymbol{x}}}_k}){{\mathit{\boldsymbol{l}}}_k} $$ (15)

      式(15)称为自适应松弛正则化迭代算法,其中β为线性搜索步长,主要基于病态问题自适应调整迭代步长,提高数值解的整体效果;αk为第k次迭代的正则化参数。

    • 文献[17]采用正则化参数衰减的几何级数方法伴随迭代不断更新正则化参数,这就必须要求初始正则化参数要足够大,然而正则化参数的最大界限难以确定。测量数据处理中,L曲线法是确定初始正则化参数最优的方法且精度最高,而且正则化参数数值取值一般比较小。若迭代过程中,数据产生扰动加剧了解的不适定性,导致迭代序列不稳定,由此可见,迭代法在计算正则化参数时,应该让其在一定范围内不断增加,提高迭代整体的收敛特性。为此,本文建议采用几何递增级数方法获取正则化参数序列:

      $$ {\alpha _k} = {\alpha _0}{q^{k - 1}} $$ (16)

      式中,α0为正则化参数初始设定值,且α0>0,准确确定正则化参数初始值是伴随迭代不断更新正则化参数的关键,本文采用L曲线法确定正则化参数初值[18]q为增长因子[19]

      自适应松弛正则迭代公式中βk的自适应更新选取显然也非常重要,过大会使迭代步长过大而导致迭代发散,过小而影响其收敛速度,增加了计算时间。文献[20-21]提出了残差最小准则确定步长的思想,验证了方法的计算效率。为此,本文采用了残差最小准则来确定迭代步长。残差最小准则方法充分采用目标函数的导数信息和函数值信息,相比于传统线搜索方法如黄金分割法和抛物线法具有一定优势。具体推导如下:设第k+1次的拟合残差为步长的函数,即

      $$ \begin{array}{c} \min R(\beta ) = {\left\| {{\mathit{\boldsymbol{f}}}({{\mathit{\boldsymbol{x}}}_{k + 1}}) - {\mathit{\boldsymbol{L}}}} \right\|^2} = \\ {\left\| {{\mathit{\boldsymbol{f}}}({{\mathit{\boldsymbol{x}}}_k} + \beta {\rm{d}}{{\mathit{\boldsymbol{x}}}_k}) - {\mathit{\boldsymbol{L}}}} \right\|^2} \end{array} $$ (17)

      式中,f(xk+1)为xk+1处的函数计算值,可将其在xk依据泰勒级数展开得:

      $$ {\mathit{\boldsymbol{f}}}({{\mathit{\boldsymbol{x}}}_{k + 1}}) \approx {\mathit{\boldsymbol{f}}}({{\mathit{\boldsymbol{x}}}_k}) + {\mathit{\boldsymbol{J}}}({{\mathit{\boldsymbol{x}}}_k})({{\mathit{\boldsymbol{x}}}_{k + 1}} - {{\mathit{\boldsymbol{x}}}_k}) $$ (18)

      考虑到xk+1=xk+βdxk,联合式(17)和式(18),可将步长函数化简为:

      $$ \min R(\beta ) = {\left\| {{\mathit{\boldsymbol{f}}}({{\mathit{\boldsymbol{x}}}_k}) + \beta {\mathit{\boldsymbol{J}}}({{\mathit{\boldsymbol{x}}}_k}){\rm{d}}{{\mathit{\boldsymbol{x}}}_k} - {\mathit{\boldsymbol{L}}}} \right\|^2} $$ (19)

      为简要计算,本文定义lk=L-f(xk),dβ=βdxk,则上述方程可化简为:

      $$ \min R(\beta ) = {\left\| {{\mathit{\boldsymbol{J}}}({{\mathit{\boldsymbol{x}}}_k}){{\mathit{\boldsymbol{d}}}_\beta } - {{\mathit{\boldsymbol{l}}}_k}} \right\|^2} $$ (20)

      将式(20)关于β求其偏导为:

      $$ \begin{array}{c} \frac{{\partial R(\beta )}}{{\partial \beta }} = {({\mathit{\boldsymbol{J}}}({{\mathit{\boldsymbol{x}}}_k}){{\mathit{\boldsymbol{d}}}_\beta } - {{\mathit{\boldsymbol{l}}}_k})^T}{\mathit{\boldsymbol{J}}}({{\mathit{\boldsymbol{x}}}_k})\frac{{\partial {\mathit{\boldsymbol{d}}_\beta }}}{{\partial \beta }} = \\ {({\mathit{\boldsymbol{J}}}({{\mathit{\boldsymbol{x}}}_k}){{\mathit{\boldsymbol{d}}}_\beta } - {{\mathit{\boldsymbol{l}}}_k})^T}{\mathit{\boldsymbol{J}}}({{\mathit{\boldsymbol{x}}}_k}){\rm{d}}{{\mathit{\boldsymbol{x}}}_k} \end{array} $$ (21)

      若要k+1次残差最小,需满足∂R(β)/∂β=0,即:

      $$ \frac{{\partial R(\beta )}}{{\partial \beta }} = {({\mathit{\boldsymbol{J}}}({{\mathit{\boldsymbol{x}}}_k}){{\mathit{\boldsymbol{d}}}_\beta } - {{\mathit{\boldsymbol{l}}}_k})^T}{\mathit{\boldsymbol{J}}}({{\mathit{\boldsymbol{x}}}_k}){\rm{d}}{{\mathit{\boldsymbol{x}}}_k} = 0 $$ (22)

      采用正则化数值方法处理非线性方程,其搜索方向dxk表示为:

      $$ {\rm{d}}{{\mathit{\boldsymbol{x}}}_k} = {({{\mathit{\boldsymbol{J}}}^T}({{\mathit{\boldsymbol{x}}}_k}){\mathit{\boldsymbol{J}}}({{\mathit{\boldsymbol{x}}}_k}) + {\alpha _k}{\mathit{\boldsymbol{I}}})^{ - 1}}{{\mathit{\boldsymbol{J}}}^T}({{\mathit{\boldsymbol{x}}}_k}){{\mathit{\boldsymbol{l}}}_k} $$ (23)

      此时,dβ=βdxk可化简为:

      $$ {{\mathit{\boldsymbol{d}}}_\beta } = \beta {({{\mathit{\boldsymbol{J}}}^T}({{\mathit{\boldsymbol{x}}}_k}){\mathit{\boldsymbol{J}}}({{\mathit{\boldsymbol{x}}}_k}) + {\alpha _k}I)^{ - 1}}{{\mathit{\boldsymbol{J}}}^T}({{\mathit{\boldsymbol{x}}}_k}){{\mathit{\boldsymbol{l}}}_k} $$ (24)

      将式(23)和式(24)代入式(22)可得步长计算公式:

      $$ \beta = \frac{{\left\| {{\mathit{\boldsymbol{l}}}_k^T({\mathit{\boldsymbol{J}}}({{\mathit{\boldsymbol{x}}}_k}){{({{\mathit{\boldsymbol{J}}}^T}({{\mathit{\boldsymbol{x}}}_k}){\mathit{\boldsymbol{J}}}({{\mathit{\boldsymbol{x}}}_k}) + {\alpha _k}I)}^{ - 1}}{{\mathit{\boldsymbol{J}}}^T}({{\mathit{\boldsymbol{x}}}_k}){{\mathit{\boldsymbol{l}}}_k})} \right\|}}{{{{\left\| {{\mathit{\boldsymbol{J}}}({{\mathit{\boldsymbol{x}}}_k}){{({{\mathit{\boldsymbol{J}}}^T}({{\mathit{\boldsymbol{x}}}_k}){\mathit{\boldsymbol{J}}}({{\mathit{\boldsymbol{x}}}_k}) + {\alpha _k}I)}^{ - 1}}{{\mathit{\boldsymbol{J}}}^T}({{\mathit{\boldsymbol{x}}}_k}){{\mathit{\boldsymbol{l}}}_k}} \right\|}^2}}} $$

      考虑到迭代过程,可将步长记为如下形式:

      $$ {\beta _k} = \frac{{\left\| {{\mathit{\boldsymbol{l}}}_k^T({\mathit{\boldsymbol{J}}}({{\mathit{\boldsymbol{x}}}_k}){{({{\mathit{\boldsymbol{J}}}^T}({{\mathit{\boldsymbol{x}}}_k}){\mathit{\boldsymbol{J}}}({{\mathit{\boldsymbol{x}}}_k}) + {\alpha _k}I)}^{ - 1}}{{\mathit{\boldsymbol{J}}}^T}({{\mathit{\boldsymbol{x}}}_k}){{\mathit{\boldsymbol{l}}}_k})} \right\|}}{{{{\left\| {{\mathit{\boldsymbol{J}}}({{\mathit{\boldsymbol{x}}}_k}){{({{\mathit{\boldsymbol{J}}}^T}({{\mathit{\boldsymbol{x}}}_k}){\mathit{\boldsymbol{J}}}({{\mathit{\boldsymbol{x}}}_k}) + {\alpha _k}I)}^{ - 1}}{{\mathit{\boldsymbol{J}}}^T}({{\mathit{\boldsymbol{x}}}_k}){{\mathit{\boldsymbol{l}}}_k}} \right\|}^2}}} $$ (25)

      综上所述,本文设计得到了自适应松弛正则化迭代公式的解算步骤如下:

      1) 给定迭代初值x0,将非线性模型在x0处线性化,按一般间接平差组成法方程。基于先验信息进行搜索或采用无需矩阵逆运算的梯度法确定初值。对于一些问题,可采用部分观测数据解析获取问题的初值[22]

      2) 基于条件数法判断法矩阵的病态性,基本公式为:

      $$ {\rm{cond}}({\mathit{\boldsymbol{N}}}({\mathit{\boldsymbol{x}}})) = \left\| {{\mathit{\boldsymbol{N}}}({\mathit{\boldsymbol{x}}})} \right\|\left\| {{{\mathit{\boldsymbol{N}}}^{ - 1}}({\mathit{\boldsymbol{x}}})} \right\| $$ (26)

      条件数反映了系统抗干扰性,刻画了参数解对原始数据变化的灵敏程度,即方程组法矩阵的病态程度。法方程态性良好,选择高斯牛顿迭代公式计算,并直接跳到5);法方程病态,采用自适应松弛正则化迭代公式进行计算,并进行下一步;

      3) 根据L曲线方法和残差最小准则分别确定正则化参数初值和步长βk

      4) 根据实际观测数据进行试算确定增长系数,兼顾数值解算精度和稳定性取值q

      5) 采用数值迭代方法进行求解。若满足正交条件方程‖h(dx)‖≤δ,则迭代结束,输出解算结果;否则,继续迭代,直至满足收敛条件。

    • 案例1:采用文献[18]模拟的病态测边网进行试验。

      P1P2P9为9个已知点,且已知点到待测点P10(模拟值为(0, 0, 0)m)的观测距离如表 1所示,距离为等精度观测,观测图形如图 1所示,要求根据9个观测距离来确定待测点的坐标。假设初始值为(0.1, -0.1, 0.1), 迭代终止条件为$ \left\|\boldsymbol{h}\left(\hat{\boldsymbol{x}}_{k}\right)\right\| \leqslant 1 \times 10^{-8}$。非线性代数方程构建的非线性方程组依据泰勒公式线性化后,由于已知点和未知点近似共面,方程组系数矩阵具有严重的病态性,条件数为cond(N(x0))=1.815 3×106

      表 1  测边网控制点坐标和已知观测距离

      Table 1.  Trilateration Network Known Coordinates and Observation Distance

      点号 控制点坐标 观测距离/m
      X/m Y/m Z/m
      P1 23.000 10.000 0.010 25.078 6
      P2 -10.000 9.990 0.000 14.134 5
      P3 35.000 10.010 -0.010 36.415 9
      P4 100.000 19.990 0.005 101.479 4
      P5 -36.000 10.005 -0.000 37.364 2
      P6 0.000 10.010 -0.005 10.010 0
      P7 56.000 9.995 0.010 56.996 1
      P8 -15.000 10.015 -0.010 18.035 9
      P9 -1.700 10.008 0.015 10.150 6

      图  1  病态测边网

      Figure 1.  Morbid Trilateration Network

      图 1表 1数据图形点分布初步分析可知,构建的非线性方程组线性化后,方程组病态性主要集中在Z方向,Y方向只存在微弱的病态性,X方向态性良好。分别采用线性化平差估计、高斯牛顿法、迭代Tikhonov方法和自适应松弛正则化算法计算,结果详见表 2

      表 2  不同方法的解算结果(算例1)

      Table 2.  Calculation Results of Different Methods(Example 1)

      方法 数值收敛解/m 迭代次数k 时间/s 残差‖Δ$\hat x $‖/m
      线性化平差 0.034 0 1.109 7 112.778 3 112.783 7
      高斯牛顿法 不收敛 无解
      迭代Tikhonov -0.027 5 5.370 1 8.883 1 1 097 0.146 617 10.380 2
      自适应松弛正则化算法 -0.027 5 5.370 1 8.883 1 98 0.019 011 10.380 2

      表 2可知,非线性函数模型直接采用线性化平差估计,由于模型病态和初始值精度的影响,其估计解与真解偏差较大,可信度低;高斯牛顿法由于初始值精度和方程组病态的影响而无法收敛,解算失效,为此,本文建议病态方程组解算一般不采用高斯牛顿法(案例2不再考虑该方法)。相对于高斯牛顿法,迭代Tikhonov方法和自适应松弛正则化算法采用了稳定泛函约束,改善了非线性方程的不适定性,可靠性更高。迭代Tikhonov方法和自适应松弛正则化算法均明显提高了线性平差估计解的精度,其解算结果一致,原因是正则化参数的自动更新并没有改变正则化参数的作用机理,而自适应调整步长是为了降低相邻迭代值过大对数值解的影响,并提高数值迭代法的收敛效率。

      图 2给出了迭代Tikhonov方法和自适应松弛正则化算法的点位迭代序列,k为迭代次数。由图 2可知,自适应松弛正则化算法在提高迭代Tikhonov算法收敛速度的同时,并不会降低数值算法的稳定性。结合表 2数值收敛解所耗时间发现,自适应松弛正则化算法明显提高了迭代Tikhonov方法的收敛速度,表明自适应松弛正则化算法具有更好的局部收敛性质。

      图  2  点位坐标迭代解序列(算例1)

      Figure 2.  Coordinate Iteration Sequences(Example 1)

      案例2:短程水下定位实测数据验证。

      为进一步研究自适应松弛正则化算法的实用性,本文采用水下定位实测数据验证,数据主要是测量船围绕应答器航行获取。在1圈数据选取15个相邻数据进行验证,数据坐标如图 3所示。在1圈数据中,假设观测距离为等精度观测,要求根据15个已知点坐标和观测距离求待测点的坐标。以1圈数据的单点定位结果作为真值,即(2 438 925.079 3, 491 997.435 4, 2 010.693 8) m。假定初值为(2 438 900, 492 000, 1 000) m。迭代终止条件为‖ h(xk)‖≤1×10-8

      图  3  15个相邻数据的已知点坐标

      Figure 3.  Known Point Coordinates for 15 Adjacent Data

      水下定位数据采集过程中,测量船在水平面航行,非线性代数方程构建的非线性方程组依据泰勒公式线性化后,Z方向上由于已知点共面,具有较强的病态性。此外,若采集数据相隔时间较短,X方向和Y方向也会产生微弱的病态性。该案例所构建非线性方程组依据泰勒公式转为线性方程组后,方程组的系数矩阵具有严重的病态性,条件数为cond(N(x0))=1.613 5×103。分别采用线性化平差估计、迭代Tikhonov方法和自适应松弛正则化算法进行计算,比较其数值收敛解、迭代次数k、所耗时间和待估参数的残差。

      表 3给出了线性化平差估计、迭代Tikhonov方法和自适应松弛正则化算法的解算结果。由于水下测距定位方程系数矩阵具有较强的病态性和线性化初值精度的影响,线性化平差估计解与真解相差较大,可信度低;迭代Tikhonov方法和自适应松弛正则化算法采用稳定泛函准则约束,解决了水下短程测距定位方程的不适定问题,数值收敛解明显提高了线性化平差估计解。自适应松弛正则化算法采用几何递增方法更新正则化参数,保持了正则化参数的作用特性,而自适应采用残差最小准则来自适应调整步长,降低了迭代序列之间残差过大对数值解的影响,提高了迭代Tikhonov方法的收敛速度。

      表 3  不同方法的解算结果(算例2)

      Table 3.  Calculation Results of Different Methods(Example 2)

      方法 数值收敛解/m k 时间/s 残差‖Δ$\hat x$‖/m
      线性化平差 2 553 993.560 3 -145 334.776 5 371 795 350.815 7 371 793 904.189
      迭代Tikhonov 2 438 949.021 9 492 011.361 4 2 142.926 3 512 0.087 446 135.102
      自适应松弛正则化算法 2 438 949.021 9 492 011.361 4 2 142.926 3 44 0.009 753 135.102

      图 4给出了迭代Tikhonov方法和自适应松弛正则化算法的点位迭代序列。自适应松弛正则化算法在提高迭代Tikhonov方法收敛速度的同时,并不会降低数值算法的稳定性。此外,结合表 3中迭代序列所消耗时间,能够看出自适应松弛正则化算法明显提高了迭代Tikhonov方法的收敛速度,表明自适应松弛正则化算法具有更好的局部收敛性质。

      图  4  点位坐标迭代解序列(算例2)

      Figure 4.  Coordinate Iteration Sequences(Example 2)

    • 若有非线性代数方程构成的线性化方程组直接采用无约束非线性最小二乘求解,可靠性低且解算误差大,原因是线性平差估计解易受到线性化初值精度和模型态性的影响。高斯牛顿法处理短程测距定位问题时,常由于受到初始值精度和方程组病态的影响而无法收敛。本文基于递增几何级数自适应选择正则化参数和残差最小步长准则自适应更新迭代步长,提出一种自适应松弛正则化数值算法,并以病态仿真数据和水下定位实测数据为例进行试验。结果表明,线性平差估计解可靠性低且误差大;自适应松弛正则化算法数值收敛解优于线性平差估计解,且收敛特性优于迭代Tikhonov正则化方法。当然,该方法受到增长因子经验取值限制,其合理确定有待进一步研究。

参考文献 (22)

目录

    /

    返回文章
    返回