-
在测量数据处理领域，根据观测对象的几何或物理特性，往往能建立起参数间的不等式约束关系。合理地利用这些先验信息，能提高平差结果的精度和可靠性[1]。不等式约束最小二乘平差在GPS数据处理[2-3]、变形分析[4]领域得到成功应用。不等式约束最小二乘平差算法主要包括规划类算法、椭圆约束法[5]、最小距离法[6]、贝叶斯方法[3]、凝聚函数法[1]等。参数精度评定的方法主要有贝叶斯方法[3]、将有效约束视为等式约束[7]、凝聚函数法[1]和Monte-Carlo模拟法[8]。
由于能考虑系数阵元素含误差的情况，变量含误差(errors-in-variables, EIV)模型成为近年来的研究热点[9-12]。当参数间含有不等式约束时，就形成了不等式约束EIV(inequality cvonstrained EIV, ICEIV)模型。de Moor给出了不等式约束TLS问题的KKT条件，并转化为广义线性补问题采用组合式方法求解[13]。Zhang等采用穷举搜索求解不等式约束EIV模型，其计算量随约束数的增加呈指数增长，且未考虑不等权和相关观测量[14]; Fang提出了不等式约束加权整体最小二乘(weighted total least squares with inequality constraints, ICWTLS)问题，给出了最优解的一阶必要和二阶充分条件，将ICWTLS转换为标准形式的约束最优化问题，采用起作用集方法和SQP方法求解[15]; Zeng等给出了不等式约束部分EIV(inequality constrained partial EIV, ICPEIV)模型，转化为非线性约束最优化问题采用罚函数法求解[16, 17], 当约束同时包含待求参数和系数阵元素时，将ICPEIV模型转化为线性补(linear complementarity problem, LCP)问题求解，并得到了解的近似精度。
Fang的算法基于最优化方法，未考虑解的精度，而Zeng等的方法没有考虑系数阵元素与右端项相关的情况。本文将所有的不等式约束转换成一个等价的凝聚约束函数，将不等式约束加权整体最小二乘平差转换成等式约束非线性最优化问题，引入乘子罚函数策略进一步表达成无约束非线性规划问题，采用BFGS算法求解。推导了不等式约束加权整体最小二乘解与观测向量的一阶线性近似函数关系，采用协方差传播律得到了解的近似方差。并用数值实例对该方法进行了验证。
HTML
-
不等式约束加权整体最小二乘平差模型的函数模型为[17]：
式中，y和ey分别表示n维观测向量和误差向量；A和EA表示n×u维的设计矩阵和误差矩阵；x表示u维参数向量；G表示s×u的约束矩阵；W表示相应的s维常数向量。将约束矩阵按行表达成G=[g1 g2…gs]T，其中，giT=[gi1 gi2…giu]是G的第i行(i=1, 2…s)。W=[w1 w2…ws]T表示约束右端项。式(2)可以写成gi(x)=giTx-wi。
相应的随机模型为：
式中，QAA和Qyy表示nu×nu系数阵元素的协因数阵和n×n观测值协因数阵；QAy=QyAT表示nu×n的系数阵元素和观测值向量的互协因数阵；Qe是所有随机误差向量的对称正定协因数矩阵；σ02为单位权方差。式(1)-(4)的最优解应使所有随机误差的平方和最小，且满足不等式约束条件，即
式(5)中极小化函数的变量是eA和ey，而约束条件中的变量是待求参数x。为了使最小化函数和约束条件中的变量一致，首先不考虑不等式约束条件，函数模型式(1)和随机模型在整体最小二乘准则下，组成Lagrange目标函数为：
分别对x、e、λ求偏导，并令其值为0，最优解的Euler-Lagrange必要条件为[15]：
式中，${\mathit{\boldsymbol{\hat x}}}$、${\mathit{\boldsymbol{\hat \lambda }}}$分别表示参数和Lagrange乘子的估值；${\mathit{\boldsymbol{\tilde e}}}$表示观测误差的预测值。并且令$\mathit{\boldsymbol{\hat B = }}\left[ {{{\mathit{\boldsymbol{\hat x}}}^{\rm{T}}} \otimes {\mathit{\boldsymbol{I}}_n}, - {\mathit{\boldsymbol{I}}_n}} \right]$，由式(8)可以得到：
代入式(9)，则有：
将式(11)重新代入式(10)得：
根据式(12)，可以得到整体最小二乘准则函数的另一种形式如下：
因此，由式(13)可以得到和式(5)等价的约束非线性最优化问题：
式(14)中极小化目标函数和约束条件仅含待估参数，是一个典型的非线性最优化问题。可以采用最优化理论与方法中的积极约束法、序列二次规划算法、牛顿法、拟牛顿法等计算。但这些方法都基于迭代策略，不能得到解与观测值之间的函数关系，难以获得解的统计性质。因此，采用凝聚函数法来解式(14)并推导解的统计性质。
-
在最优化理论与方法中，处理不等式约束非线性最优化问题的一种策略是将其转换为相应的等式约束最优化问题求解。凝聚函数法是将多个不等式约束合并成一个等价的凝聚约束函数，然后用拉格朗日乘子迭代法求解。根据文献[1]的结论，式(2)可以转换为如下的凝聚约束函数：
式中，p是趋向于无穷大的正数，在计算中根据计算精度和计算效率的要求，一般取p=106。联合式(14)、式(15)，不等式约束的加权整体最小二乘平差模型等价于求以下等式约束的加权整体最小二乘平差模型：
式(16)是一个附有等式约束的非线性最优化问题，它的解可以通过极小化如下的拉格朗日目标函数得到：
其中，α是等式约束对应的拉格朗日乘子，是一个标量。当式(17)不收敛的情况下，可以在此目标函数的基础上加入一个二次惩罚项，得到一个新的最小化目标函数:
式中，c是一个充分大的正常数，称为罚因子，可取c=106。其作用是：当gp(x)≠0时施以惩罚。在迭代计算过程中，α通过下式更新：
实际运算中，初值可以取α0=0，p=c=106，可以取相应的不等式约束最小二乘解作为迭代初值，当‖xk+1-xk‖≤ε(ε是迭代终止阈值，可取ε=10-6)时，停止迭代。然后计算▽Φ(x, α)，若▽Φ(x, α)较大，根据式(19)计算新的α值，迭代直到‖▽Φ(x, α)‖也满足限差，最终得到不等式约束加权整体最小二乘解。由于拟牛顿算法解非线性规划问题具有较快的收敛速度，且无需求解目标函数的Hesse矩阵，近似的二阶校正公式更新规则简单，从而得到广泛应用。本文采用目前最流行有效的拟牛顿算法BFGS方法来求解无约束非线性最优化问题。
-
式(16)是关于参数的非线性等式约束平差模型，为了得到解的近似精度，需要将目标函数做线性近似处理。设参数x可以表示成它的近似值和改正数之和，即x=x0+Δx，注意，经过非线性平差有限次迭代计算后x0无限趋近于不等式约束加权整体最小二乘解${{\hat x}_{{\rm{IC}}}}$，Δx趋向于0。且x和Δx具有相同的方差，则式(1)可以表示为[18]：
将式(20)展开，并忽略二阶微小量EAΔx，则式(20)可以重新表示为：
令：
e′=EAx0-ey，y′=y-Ax0，则式(21)可以写成：
由于EAx0-ey=Ke，根据协方差传播律有：
将凝聚约束函数式(15)在x0处线性化，有：
即
联立式(22)与式(25)，组成只含有一个等式约束的线性最小二乘平差模型，参数Δx的协因数矩阵可以表示为：
式中，Naa=ATQe-1′A，Ncc=CTNaa-1C是一个常数；矩阵K的元素用最终迭代解${{\mathit{\boldsymbol{\hat x}}}_{{\rm{IC}}}}$代替。单位权方差的有偏估值近似写成：
协因数矩阵式(26)与文献[1]相比较完全一致。由于加权整体最小二乘模型线性化过程中忽略了二阶项的影响，且凝聚函数只有当p趋向无穷大时才与不等式约束等价，所以，式(26)是不等式约束加权整体最小二乘的近似协因数估值。对于单位权方差式(27)，文献[1]考虑了偏差项，而本文是把不等式约束加权整体最小二乘平差模型近似成等式约束的非线性平差模型，然后一阶线性近似成附有等式约束的间接平差模型来评定参数的精度。
2.1. 计算方法
2.2. 精度评定
-
算例引自文献[1]，观测方程(1)和约束条件(2)中的A、y、G1和w1如表 1所示。区间约束条件为-0.1≤xi≤2(i=1, 2, 3, 4)。设系数阵元素和右端向量的随机误差都是独立同分布，即Qe=I25为25阶单位矩阵。首先，采用二次规划方法计算相应的不等式约束最小二乘解：
A y 0.950 1 0.762 0 0.615 3 0.405 7 0.057 8 0.231 1 0.456 4 0.791 9 0.935 4 0.352 8 0.606 8 0.018 5 0.921 8 0.916 9 0.813 1 0.485 9 0.821 4 0.738 2 0.410 2 0.009 8 0.891 2 0.444 7 0.176 2 0.893 6 0.138 8 G1 w1 0.202 7 0.272 1 0.746 7 0.465 9 0.525 1 0.198 7 0.198 8 0.445 0 0.418 6 0.202 6 0.603 7 0.015 2 0.931 8 0.846 2 0.672 1 Table 1. A, y, G1, w1 of Observation Vector(1), Constraints Condition(2)
然后将上述不等式约束最小二乘解作为迭代初始值，计算相应的不等式约束整体最小二乘解。比较文献[17]提出的有效集方法、序列二次规划算法以及文献[16]的穷举搜索法和本文提出的凝聚函数法，数值计算结果如表 2所示。
有效集方法 序列二次规划 穷举搜索法 凝聚函数法 ${{\hat x}_1}$ -0.100 000 -0.100 000 -0.100 000 -0.099 991 ${{\hat x}_2}$ -0.100 000 -0.100 000 -0.100 000 -0.100 003 ${{\hat x}_3}$ 0.168 547 0.168 547 0.168 500 0.168 561 ${{\hat x}_4}$ 0.399 777 0.399 777 0.399 800 0.399 765 ${{\hat e}^{\rm{T}}}\mathit{\boldsymbol{P}}\hat e$ 0.139 737 0.139 737 0.139 742 0.139 741 Table 2. The Solutions of Different Algorithms and the Residuals
采用有效集方法外部迭代3次，内部迭代平均8次；序列二次规划算法外部迭代4次；内部迭代平均4次，而穷举搜索法外部迭代2 048次，内部迭代平均23次。采用凝聚函数方法，取p=c=105，经过9次迭代就能得到上述解，相应的目标函数梯度为：
‖▽Φ(x, α)‖满足阈值，停止迭代，可见本文的凝聚约束函数目标函数采用BFGS算法具有较快的收敛速度。
计算凝聚函数在最优解处对参数x的梯度值为：
根据式(27)得到单位权方差的有偏估计量为：σ02=0.069 9。根据式(26)可以得到不等式约束整体最小二乘解的近似协因数矩阵为：
在运用凝聚函数法进行计算时，如果p取值过大，采用普通不等式约束最小二乘解作为初值进行计算时，容易出现不收敛的情况，并且随着p值的增大计算量会显著增加，可以根据精度和计算效率选取合适的p值。根据计算结果可以看出，解算值并不是严格满足所有不等式约束，这是由于p值是有限值时，凝聚函数并不是完全等价于最大的不等式约束。前面两种方法的结果在小数点后前6位都相同，凝聚函数法也能得到接近一致的结果。这说明凝聚函数法求解不等式约束的加权整体最小二乘问题是可行的。
-
本文提出了一般权矩阵情形下的不等式约束加权整体最小二乘模型的计算方法和近似精度评定公式。对于不等式约束加权整体最小二乘平差模型，且权矩阵为一般对称正定矩阵的情形，将残差加权平方和最小准则先转化成仅含有参数的非线性目标函数极小化问题，并将所有不等式约束表达成一个近似的凝聚等式约束，构造出无约束拉格朗日目标函数，采用非线性最优化BFGS算法得到不等式约束加权整体最小二乘平差的解。数值实例结果表明，该方法能收敛到最优解，且具有较快的收敛速度。
DownLoad: