总体最小二乘法将系数矩阵中随机量误差和观测向量中随机量误差一同考虑，是基于EIV^[1](errors-in-variables)模型的处理方法。实际中，很多情况下系数矩阵由随机元素和固定元素共同组成，使用Partial EIV^[2](partial errors-in-variables)模型进行未知参数的求解具有更高的效率。当参数估计具有先验信息时，需要将Partial EIV模型转换为附有约束的Partial EIV模型进行求解。Zhang等^[3]用列举法处理不等式约束的总体最小二乘问题，并穷举出所有最优解; Moor^[4]将不等式约束的平差问题转化为线性互补问题求最优参数解; 王乐洋等^[5]讨论了附有不等式约束的大地测量反演理论的技术条件; 汪奇生等^[6]给出一种迭代算法求解附有不等式约束时的参数估值，但只考虑权阵为单位阵，未考虑观测数据的相关性; 曾文宪等^[7]用罚函数法计算不等式约束加权总体最小二乘问题，但未对海森矩阵进行修正，在求解过程中会出现海森矩阵奇异的问题; Zeng等^[8]提出未知参数与系数矩阵都含有不等式约束的总体最小二乘解法。以上文献大多只强调参数估值的计算方法，很少关注参数估值的精度问题^[9]。王乐洋等^[10]使用拟牛顿修正法解算附有不等式约束的EIV模型，也未对参数估值进行精度评定。针对这个问题，本文将附有不等式约束的EIV模型转化为更具有一般性的附有不等式约束的Partial EIV模型，使用拟牛顿修正解法进行参数求解，并引入SUT^[11-12]法对参数估值进行精度评定。

1 附有不等式约束的Partial EIV模型的序列二次规划模型 1.1 附有不等式约束的Partial EIV模型

附有不等式约束的Partial EIV模型的数学模型表达式为：

$ \left\{ \begin{array}{l} \mathit{\boldsymbol{y}} = \left({{\mathit{\boldsymbol{\beta }}^{\rm{T}}} \otimes {\mathit{\boldsymbol{I}}_n}} \right)(\mathit{\boldsymbol{h}} + \mathit{\boldsymbol{B\bar a}}) + {\mathit{\boldsymbol{e}}_y}\\ \mathit{\boldsymbol{a}} = \mathit{\boldsymbol{\bar a}} + {\mathit{\boldsymbol{e}}_a}\\ {\mathop{\rm vec}\nolimits}(\mathit{\boldsymbol{\overline A}}) = \mathit{\boldsymbol{h}} + \mathit{\boldsymbol{B}}\mathit{\boldsymbol{\overline \alpha }} \\ \mathit{\boldsymbol{G\beta }}-\mathit{\boldsymbol{w}} \ge 0 \end{array} \right. $

(1)

$ \mathit{\boldsymbol{e}} = \left[{\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{e}}_y}}\\ {{\mathit{\boldsymbol{e}}_a}} \end{array}} \right]\sim \left({\left[{\begin{array}{*{20}{l}} 0\\ 0 \end{array}} \right], \sigma _0^2{\mathit{\boldsymbol{Q}}_e} = \sigma _0^2{\mathit{\boldsymbol{P}}^{-1}}} \right) $

(2)

式中，y∈R^n×1为观测向量，e_y∈R^n×1为其含有的随机误差; A∈R^n×m为系数矩阵的真值; a∈R^t×1为系数矩阵中随机元素构成的列向量，e_a∈R^t×1为其含有的随机误差，a∈R^t×1为系数矩阵中随机元素的真值; β∈R^m×1为未知参数; h∈R^m×1为系数矩阵中非随机元素构成的列向量; B∈R^nm×t为构造矩阵; ⊗表示克罗内克积，vec(·)表示将矩阵拉直为列向量; σ₀²为单位权方差; 矩阵Q_e和P分别表示误差向量e的协因数阵和权阵; G∈R^U×m、w∈R^U×1分别为不等式组系数矩阵和常数向量。

对式(1)进行改写，得到函数模型^[13]：

$ \left\{ {\begin{array}{*{20}{l}} {\left[ {\begin{array}{*{20}{l}} \mathit{\boldsymbol{y}}\\ \mathit{\boldsymbol{a}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\left({{\mathit{\boldsymbol{\beta }}^{\rm{T}}} \otimes {\mathit{\boldsymbol{I}}_n}} \right)(\mathit{\boldsymbol{h}} + \mathit{\boldsymbol{B\bar a}})}\\ {\mathit{\boldsymbol{\bar a}}} \end{array}} \right]-\mathit{\boldsymbol{e}}}\\ {\mathit{\boldsymbol{G\beta }}-\mathit{\boldsymbol{w}} \ge 0} \end{array}} \right. $

(3)

1.2 附有不等式约束Partial EIV模型的序列二次规划模型

根据观测数据的误差平方和最小原则，得到最小二乘准则：

$ \left\{ {\begin{array}{*{20}{l}} {\min {\mathit{\boldsymbol{e}}^{\rm{T}}}\mathit{\boldsymbol{Pe}}}\\ {{\rm{ s}}{\rm{.t}}{\rm{. }}\mathit{\boldsymbol{y}}-\left({{\mathit{\boldsymbol{\beta }}^{\rm{T}}} \otimes {\mathit{\boldsymbol{I}}_n}} \right)(\mathit{\boldsymbol{h}} + \mathit{\boldsymbol{Ba}}) + \mathit{\boldsymbol{Ce}} = 0} \end{array}} \right. $

(4)

根据式(4)构造拉格朗日目标函数：

$ \begin{array}{l} \mathit{\Phi }(\mathit{\boldsymbol{e}}, \mathit{\boldsymbol{\beta }}, \mathit{\boldsymbol{\lambda }}) = {\mathit{\boldsymbol{e}}^{\rm{T}}}\mathit{\boldsymbol{Pe}} + 2{\mathit{\boldsymbol{\mu }}^{\rm{T}}}[\mathit{\boldsymbol{y}}-\\ \;\;\;\left({{\mathit{\boldsymbol{\beta }}^{\rm{T}}} \otimes {\mathit{\boldsymbol{I}}_n}} \right)(\mathit{\boldsymbol{h}} + \mathit{\boldsymbol{Ba}}) + \mathit{\boldsymbol{Ce}}] \end{array} $

(5)

式中，$\mathit{\boldsymbol{C}} = \left[{{\mathit{\boldsymbol{I}}_n}, -\left({{\mathit{\boldsymbol{\beta }}^{\rm{T}}} \otimes {\mathit{\boldsymbol{I}}_n}} \right)\mathit{\boldsymbol{B}}} \right]$，μ为拉格朗日乘子。

对目标函数(5)分别求e、β、μ的偏导数并令其为0，得到：

$ \frac{{\partial \mathit{\Phi }}}{{\partial \mathit{\boldsymbol{e}}}} = 2\left({\mathit{\boldsymbol{Pe}} + {\mathit{\boldsymbol{C}}^{\rm{T}}}\mathit{\boldsymbol{\mu }}} \right) = 0 $

(6)

$ \frac{{\partial \mathit{\Phi }}}{{\partial \mathit{\boldsymbol{\beta }}}} = -2{\mathit{\boldsymbol{\overline A}} ^{\rm{T}}}\mathit{\boldsymbol{\mu }} = 0 $

(7)

$ \frac{{\partial \mathit{\Phi }}}{{\partial \mathit{\boldsymbol{\mu }}}} = 2\left[ {\mathit{\boldsymbol{y}}-\left({{\mathit{\boldsymbol{\beta }}^{\rm{T}}} \otimes {\mathit{\boldsymbol{I}}_n}} \right)(\mathit{\boldsymbol{h}} + \mathit{\boldsymbol{Ba}}) + \mathit{\boldsymbol{Ce}}} \right] = 0 $

(8)

由式(6)可得：

$ \mathit{\boldsymbol{e}} = -{(\mathit{\boldsymbol{P}})^{-1}}{\mathit{\boldsymbol{C}}^{\rm{T}}}\mathit{\boldsymbol{\mu }} $

(9)

将式(9)代入式(8)得到：

$ \mathit{\boldsymbol{\mu }} = {\left({\mathit{\boldsymbol{C}}{{(\mathit{\boldsymbol{P}})}^{-1}}{\mathit{\boldsymbol{C}}^{\rm{T}}}} \right)^{-1}}\left[ {\mathit{\boldsymbol{y}}-\left({{\mathit{\boldsymbol{\beta }}^{\rm{T}}} \otimes {\mathit{\boldsymbol{I}}_n}} \right)(\mathit{\boldsymbol{h}} + \mathit{\boldsymbol{Ba}})} \right] $

(10)

联立式(9)和式(10)，构造二次函数形式：

$ \left\{ {\begin{array}{*{20}{l}} {{{\mathit{\boldsymbol{\hat e}}}^{\rm{T}}}\mathit{\boldsymbol{P\hat e}} = \left[ {\mathit{\boldsymbol{y}} - \left( {{\mathit{\boldsymbol{\beta }}^{\rm{T}}} \otimes {\mathit{\boldsymbol{I}}_n}} \right)(\mathit{\boldsymbol{h}} + \mathit{\boldsymbol{Ba}})} \right]{{\left( {\mathit{\boldsymbol{C}}{\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}}{\mathit{\boldsymbol{C}}^{\rm{T}}}} \right)}^{ - 1}} \bullet }\\ {\;\;\;\;\;\;\left[ {\mathit{\boldsymbol{y}} - \left( {{\mathit{\boldsymbol{\beta }}^{\rm{T}}} \otimes {\mathit{\boldsymbol{I}}_n}} \right)} \right.\left. {(\mathit{\boldsymbol{h}} + \mathit{\boldsymbol{Ba}})} \right]}\\ {\mathit{\boldsymbol{f}}\left( \mathit{\boldsymbol{\beta }} \right) = {{\mathit{\boldsymbol{\hat e}}}^{\rm{T}}}\mathit{\boldsymbol{P\hat e}}} \end{array}} \right. $

(11)

此时，附有不等式约束的Partial EIV模型可以转化为一般的附有不等式约束的最优化问题：

$ \left\{ \begin{array}{l} \min f(\mathit{\boldsymbol{\beta }}) = \left[ {\mathit{\boldsymbol{y}}-\left({{\mathit{\boldsymbol{\beta }}^{\rm{T}}} \otimes {\mathit{\boldsymbol{I}}_n}} \right)(\mathit{\boldsymbol{h}} + \mathit{\boldsymbol{Ba}})} \right] \bullet \\ \;\;\;\;\;{\left({\mathit{\boldsymbol{C}}{\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{e}}}{\mathit{\boldsymbol{C}}^{\rm{T}}}} \right)^{-1}}\left[ {\mathit{\boldsymbol{y}}-\left({{\mathit{\boldsymbol{\beta }}^{\rm{T}}} \otimes {\mathit{\boldsymbol{I}}_n}} \right)(\mathit{\boldsymbol{h}} + \mathit{\boldsymbol{Ba}})} \right]\\ \mathit{\boldsymbol{G}}\mathit{\boldsymbol{\widehat \beta }}-\mathit{\boldsymbol{w}} \ge 0 \end{array} \right. $

(12)

利用拟牛顿修正解法求解附有不等式约束的Partial EIV模型参数，可以有效提高王乐洋等^[10]使用EIV模型求解参数的计算效率，并扩展其适用范围。针对上述最优化问题，与曾文宪等^[7]的方法不同，本文用含有拟牛顿修正的序列二次规划法求解附有不等式约束的Partial EIV模型，利用拟牛顿修正得到海森矩阵的代替矩阵，避免了每次迭代计算时海森矩阵出现非正定的情况，收敛速度更快。

2 基于拟牛顿修正的SQP(序列二次规划)方法

将式(12)构造成拉格朗日函数：

$ L(\mathit{\boldsymbol{\beta }}, \mathit{\boldsymbol{\lambda }}) = f(\mathit{\boldsymbol{\beta }})-{\mathit{\boldsymbol{\lambda }}^{\rm{T}}}(\mathit{\boldsymbol{G\beta }}-\mathit{\boldsymbol{w}}) $

(13)

根据文献[14]，式(13)存在互补条件或松弛条件，有λ_j≥0且λ_j(Gβ-w)≥0，可以保证拉格朗日乘子向量λ的非负性，其中λ_j为乘子向量λ的第j个元素。

在给定未知参数β_j和拉格朗日乘子向量λ处，对式(13)进行二次多项式近似，并把不等式约束函数Gβ-w≥0写为G(β)≥0线性化，便可以得到序列二次规划的子问题^[15]：

$ \left\{ \begin{array}{l} \min \frac{1}{2}\mathit{\boldsymbol{d}}_j^{\rm{T}}\mathit{\boldsymbol{H}}\left( {{\mathit{\boldsymbol{\beta }}_j}} \right){\mathit{\boldsymbol{d}}_j} + \nabla f{\left( {{\mathit{\boldsymbol{\beta }}_j}} \right)^{\rm{T}}}{\mathit{\boldsymbol{d}}_j}\\ {\rm{s}}.{\rm{t}}.\;G\left( {{\mathit{\boldsymbol{\beta }}_j}} \right) + \nabla G\left( {{\mathit{\boldsymbol{\beta }}_j}} \right){\mathit{\boldsymbol{d}}_j} \ge 0 \end{array} \right. $

(14)

式中，不等式约束函数G(β)≥0是将不等式约束条件$\mathit{\boldsymbol{G}}\mathit{\boldsymbol{\widehat \beta }} - \mathit{\boldsymbol{w}}$≥0改写为关于$\mathit{\boldsymbol{\widehat \beta }}$的函数所得; G(β_j)+$\nabla G{\left({{\mathit{\boldsymbol{\beta }}_j}} \right)^{\rm{T}}}{d_j}$≥0是对G($\mathit{\boldsymbol{\hat \beta }}$)≥0在β_j处进行线性化得到; $\nabla f\left({{\mathit{\boldsymbol{\beta }}_j}} \right)$是目标函数f(β)的一阶梯度; H(β_j)是f(β)的海森矩阵; d_j表示步长，其表达式为d_j=β_j+1-β_j。

拟牛顿法的基本思想是用某个矩阵B_j代替海森矩阵，所得的矩阵B_j对称正定，可以确保d_j的方向为目标函数的下降函数。本文采用与文献[10]中BFGS校正公式不同的WHP方法。为了保证序列二次规划法的全局收敛性，一般采用价值函数Ψ(β)来取得步长^[15]：

$ \mathit{\Psi }(\mathit{\boldsymbol{\beta }}) = f(\mathit{\boldsymbol{\beta }}) + \frac{{{{\left\| {{g_i}(\mathit{\boldsymbol{\beta }})} \right\|}_1}}}{\sigma } $

(15)

式中，σ代表价值函数的参数(σ>0);g_i(β)=max{G_iβ-w_i, 0}(i=1, 2, …L); 使用线性搜索法来取得步长，得到优化改正后的参数$\mathit{\boldsymbol{\hat \beta }}$_j+1和拉格朗日乘子$\hat \lambda $_j+1; 如果满足精度要求，则参数$\mathit{\boldsymbol{\hat \beta }}$_j+1为附有不等式约束的Partial EIV模型的参数最优解，否则通过WHP法校正矩阵B_j，直到求得最优解。

基于不等式约束Partial EIV模型的拟牛顿修正SQP解法计算步骤如下。

1) 给定初始参数解β₀、初始拉格朗日乘子向量λ₀、初始矩阵B₀、不等式方程的雅戈比矩阵。β₀为最小二乘解，λ₀为元素全为0的向量，B₀为单位阵，约束方程的雅戈比矩阵A₀=$\nabla \mathit{\boldsymbol{C}}{\left({{\mathit{\boldsymbol{\beta }}_0}} \right)^{\rm{T}}}$。

2) 求解序列二次规划子问题(14)，得到d_j。

3) 对于价值函数Ψ(β, σ)，选取价值函数的参数σ_j，使d_j是该函数在B_j处的下降方向。

4) Armijo搜索，令m_j是使下列不等式成立的最小非负整数m：

$ \mathit{\Psi }\left( {{\mathit{\boldsymbol{B}}_j} + {p^m}{\mathit{\boldsymbol{d}}_j}, {\sigma _j}} \right) - \mathit{\Psi }\left( {{\mathit{\boldsymbol{B}}_j}, {\sigma _j}} \right) \le {\eta ^m}{\mathit{\Psi }^\prime }\left( {{\mathit{\boldsymbol{B}}_j}, {\sigma _j};{\mathit{\boldsymbol{d}}_j}} \right) $

(16)

令${a_j} = {p^m}^{_j}, {\mathit{\boldsymbol{\beta }}_{j + 1}} = {\mathit{\boldsymbol{\beta }}_j} + {a_j}{\mathit{\boldsymbol{d}}_j}, \eta \in \left({0, \frac{1}{2}} \right), {\mathit{\Psi }^\prime }\left({{\mathit{\boldsymbol{B}}_j}, {\sigma _j}; {\mathit{\boldsymbol{d}}_j}} \right)$为Ψ(B_j, σ_j)的方向导数。

5) 计算A_j+1=$\nabla \mathit{\boldsymbol{C}}{\left({{\mathit{\boldsymbol{\beta }}_{j + 1}}} \right)^{\rm{T}}}$和拉格朗日乘子向量${\mathit{\boldsymbol{\lambda }}_{j + 1}} = {\left[{{\mathit{\boldsymbol{A}}_{j + 1}}\mathit{\boldsymbol{A}}_{j + 1}^{\rm{T}}} \right]^{ - 1}}{\mathit{\boldsymbol{A}}_{j + 1}}\nabla {\mathit{\boldsymbol{f}}_{j + 1}}$。

6) 校正矩阵B_j为B_j+1。令

$ {\mathit{\boldsymbol{h}}_j} = {\mathit{\boldsymbol{\widehat \beta }}_{j + 1}} - {\mathit{\boldsymbol{\widehat \beta }}_j} $

(17)

$ {\mathit{\boldsymbol{l}}_j} = {\nabla _\beta }\mathit{\boldsymbol{L}}\left( {{{\mathit{\boldsymbol{\widehat \beta }}}_{j + 1}}, {{\mathit{\boldsymbol{\hat \lambda }}}_{j + 1}}} \right) - {\nabla _\beta }\mathit{\boldsymbol{L}}\left( {{{\mathit{\boldsymbol{\widehat \beta }}}_j}, {{\mathit{\boldsymbol{\hat \lambda }}}_j}} \right) $

(18)

WHP拟牛顿校正要求向量h_j和I_j满足曲率条件，即$\mathit{\boldsymbol{h}}_j^{\rm{T}}{\mathit{\boldsymbol{I}}_j}$大于0，但是由式(17)、式(18)确定的向量h_j和I_j可能不满足这一条件。为此，需要用如下步骤对h_j进行修正^[14]：

$ {\mathit{\boldsymbol{Z}}_j} = {\theta _j}{\mathit{\boldsymbol{h}}_j} + {\mathit{\boldsymbol{B}}_j}{\mathit{\boldsymbol{l}}_j} - {\theta _j}{\mathit{\boldsymbol{B}}_j}{\mathit{\boldsymbol{l}}_j} $

(19)

参数θ_j的定义为：

$ {\theta _j} = \left\{ \begin{array}{l} 1, \mathit{\boldsymbol{I}}_j^{\rm{T}}\mathit{\boldsymbol{h}} \ge 0.2\mathit{\boldsymbol{I}}_j^{\rm{T}}{\mathit{\boldsymbol{B}}_j}{\mathit{\boldsymbol{I}}_j}\\ \frac{{0.8\mathit{\boldsymbol{I}}_j^{\rm{T}}{\mathit{\boldsymbol{B}}_j}{\mathit{\boldsymbol{I}}_j}}}{{\mathit{\boldsymbol{I}}_j^{\rm{T}}{\mathit{\boldsymbol{B}}_j}{\mathit{\boldsymbol{I}}_j} - \mathit{\boldsymbol{I}}_j^{\rm{T}}{\mathit{\boldsymbol{h}}_j}}}, \mathit{\boldsymbol{I}}_j^{\rm{T}}\mathit{\boldsymbol{h}} < 0.2\mathit{\boldsymbol{I}}_j^{\rm{T}}{\mathit{\boldsymbol{B}}_j}{\mathit{\boldsymbol{I}}_j} \end{array} \right. $

(20)

此时可以得到B_j的校正公式：

$ {\mathit{\boldsymbol{B}}_{j + 1}} = {\mathit{\boldsymbol{B}}_j} - \frac{{{\mathit{\boldsymbol{B}}_j}{\mathit{\boldsymbol{I}}_j}\mathit{\boldsymbol{I}}_j^{\rm{T}}{\mathit{\boldsymbol{B}}_j}}}{{\mathit{\boldsymbol{I}}_j^{\rm{T}}{\mathit{\boldsymbol{B}}_j}\mathit{\boldsymbol{I}}_j^{\rm{T}}}} + \frac{{{\mathit{\boldsymbol{Z}}_j}\mathit{\boldsymbol{Z}}_j^{\rm{T}}}}{{\mathit{\boldsymbol{I}}_j^{\rm{T}}{\mathit{\boldsymbol{Z}}_j}}} $

(21)

7) 替代矩阵B_j更新为B_j+1，令j=j+1，回到步骤2)，当最优解满足条件时停止循环。

3 不等式约束相关观测Partial EIV模型精度评定

由上文的拟牛顿修正法计算步骤可以看出，整个迭代过程的参数估计量是观测向量y和系数矩阵A的函数。将观测向量y和系数矩阵A表示为向量l=[y; vec(A)]，参数估值与向量I之间的非线性关系可以表示为：

$ \mathit{\boldsymbol{\widehat \beta }} = f(\mathit{\boldsymbol{l}}) $

(22)

针对文献[10]中未对参数估值结果进行精度评定，且已有附有不等式约束的精度评定方法需要进行复杂的求导运算来获取参数估值的协方差阵，本文引入精度评定的SUT法，其精度评定步骤如下^[11]。

1) 计算不含约束条件时的参数估值$\mathit{\boldsymbol{\widehat \beta }}$、改正数$\mathit{\boldsymbol{\hat V}}$、单位权方差$\hat \sigma _0^2$。

2) 根据输入量I=vec(y, A)、$\mathit{\boldsymbol{\hat V}}$、Q_l和单位权方差估值$\hat \sigma _0^2$构建观测值向量I的均值和协方差阵$\mathit{\boldsymbol{\bar l}} = \mathit{\boldsymbol{l}} + \mathit{\boldsymbol{\hat V}}, {\mathit{\boldsymbol{D}}_l} = \hat \sigma _0^2{\mathit{\boldsymbol{Q}}_l}$。

3) 通过比例对称采样策略构建2t+1个sigma点列向量I_i：

$ \left\{ \begin{array}{l} {\mathit{\boldsymbol{I}}_0} = \mathit{\boldsymbol{\bar l}}\\ {\mathit{\boldsymbol{I}}_i} = \mathit{\boldsymbol{\bar l}} + {\left( {\sqrt {t + \gamma } \sqrt {{\mathit{\boldsymbol{D}}_l}} } \right)_i}, i = 1, 2, \cdots , t\\ {\mathit{\boldsymbol{I}}_i} = \mathit{\boldsymbol{\bar l}} - {\left( {\sqrt {t + \gamma } \sqrt {{\mathit{\boldsymbol{D}}_l}} } \right)_i}, i = t + 1, t + 2, \cdots , 2t \end{array} \right. $

(23)

式中，$\sqrt {{\mathit{\boldsymbol{D}}_l}} $为对D_l进行卡洛斯基分解得到的下三角阵; ${\left({\sqrt {t + \gamma } \sqrt {{\mathit{\boldsymbol{D}}_l}} } \right)_i}$为矩阵$\left({\sqrt {t + \gamma } \sqrt {{\mathit{\boldsymbol{D}}_l}} } \right)$的第i列。γ=a²(t+k)-t为比例参数，常数a推荐值为0.001，k为一个比例常数，通常取0。

4) 计算sigma点列向量经过非线性变换后的样本：

$ \mathit{\boldsymbol{v}}_i^{\hat \beta } = f\left( {{\mathit{\boldsymbol{l}}_i}} \right), i = 0, 1, \cdots , 2t $

(24)

5) 确定样本$\mathit{\boldsymbol{v}}_i^{\hat \beta }$的权值：

$ \left\{ \begin{array}{l} W_0^m = \gamma /(t + \gamma )\\ W_0^c = \gamma /(t + \gamma ) + 1 - {a^2} + b\\ W_i^m = W_i^c = 1/(2(t + \gamma )), i = 1, 2, \cdots , 2t \end{array} \right. $

(25)

式中，b=2。

6) 加权计算参数估值的均值E($\mathit{\boldsymbol{\hat \beta }}$)、偏差${\mathit{\boldsymbol{b}}_{\hat \beta }}$、偏差改正后的参数估值$\mathit{\boldsymbol{\tilde \beta }}$、参数估值近似协方差阵${\mathit{\boldsymbol{D}}_{\hat \beta \hat \beta }}$：

$ E(\mathit{\boldsymbol{\widehat \beta }}) = \sum\limits_{i = 0}^{2t} {W_i^m} \mathit{\boldsymbol{v}}_i^{\hat \beta } $

(26)

$ {\mathit{\boldsymbol{b}}_{\hat \beta }} = E(\mathit{\boldsymbol{\widehat \beta }}) - \mathit{\boldsymbol{\widehat \beta }} $

(27)

$ \mathit{\boldsymbol{\widetilde \beta }} = \mathit{\boldsymbol{\widehat \beta }} - {\mathit{\boldsymbol{b}}_{\hat \beta }} $

(28)

$ {\mathit{\boldsymbol{D}}_{\hat \beta \hat \beta }} = \sum\limits_{i = 0}^{2t} {W_i^c} \left( {\mathit{\boldsymbol{v}}_i^{\hat \beta } - \mathit{\boldsymbol{\widetilde \beta }}} \right){\left( {\mathit{\boldsymbol{v}}_i^{\hat \beta } - \mathit{\boldsymbol{\widetilde \beta }}} \right)^{\rm{T}}} $

(29)

4 算例及分析

本文共仿真模拟了2套算例，并与文献[10]的计算结果进行比较。算例1主要讨论系数矩阵元素和观测量不等精度、相关的情况，算例2讨论不等精度不相关的情况。

4.1 算例1(不等式约束的平面拟合模型)

算例1使用平面拟合模型，系数矩阵由1列固定量和2列随机量组成。系数矩阵元素与观测量元素的互协因数阵可以使用MATLAB软件中的函数randn模拟生成，单位权方差σ₀²=0.1。再使用MATLAB函数mvnrnd生成随机误差进行计算。

算例1可以表示为：

$ \begin{array}{l} \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left[ {\begin{array}{*{20}{c}} {{y_1}}\\ {{y_2}}\\ \vdots \\ {{y_n}} \end{array}} \right] - \left[ {\begin{array}{*{20}{c}} {{e_{{y_1}}}}\\ {{e_{{y_2}}}}\\ \vdots \\ {{e_{{y_n}}}} \end{array}} \right] = \\ \left( {\left[ {\begin{array}{*{20}{c}} 1&{{c_1}}&{{b_1}}\\ 1&{{c_2}}&{{b_2}}\\ \vdots & \vdots & \vdots \\ 1&{{c_n}}&{{b_n}} \end{array}} \right] - \left[ {\begin{array}{*{20}{c}} 0&{{e_{{c_1}}}}&{{e_{{b_1}}}}\\ 0&{{e_{{c_2}}}}&{{e_{{b_2}}}}\\ \vdots & \vdots & \vdots \\ 0&{{e_{{c_n}}}}&{{e_{{b_n}}}} \end{array}} \right]} \right)\left[ {\begin{array}{*{20}{c}} {{\beta _1}}\\ {{\beta _2}}\\ {{\beta _3}} \end{array}} \right] \end{array} $

(30)

将式(30)表示为Partial EIV模型式(1)的形式，式中的矩阵和向量可以表示为：

$ \mathit{\boldsymbol{\bar a}} = \left[ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{c}}_{n \times 1}}}\\ {{\mathit{\boldsymbol{b}}_{n \times 1}}} \end{array}} \right], \mathit{\boldsymbol{h}} = \left[ {\begin{array}{*{20}{c}} {{{\bf{1}}_{n \times 1}}}\\ {{{\bf{0}}_{2n \times 1}}} \end{array}} \right], \mathit{\boldsymbol{B}} = \left[ {\begin{array}{*{20}{c}} {{{\bf{0}}_{n \times 2n}}}\\ {{\mathit{\boldsymbol{I}}_{2n \times 2n}}} \end{array}} \right] $

(31)

算例1的相应数据见表 1，表 1中的A为系数矩阵，y为观测向量，G为约束条件的系数矩阵，w为约束条件的常数向量。分别用不含约束条件时的总体最小二乘、附有不等式约束的Partial EIV模型与附有不等式约束的EIV模型求参数估值的结果$\mathit{\boldsymbol{\widehat \beta }}$_TLS、$\mathit{\boldsymbol{\widehat \beta }}$_ICPEIV、$\mathit{\boldsymbol{\widehat \beta }}$_ICEIV、迭代次数k以及参数真值$\mathit{\boldsymbol{\tilde \beta }}$，见表 2，精度评定结果见表 3。

由表 2可以看出，本文算法与附有不等式约束的EIV模型参数估值计算结果几乎完全一致。在运算效率上，本文算法中求解子问题的迭代次数为53次，母问题次数为7次，总次数为60次，而使用文献[10]中的附有不等式约束的EIV模型求解的迭代总次数为106次，本文算法大大减少了迭代次数，收敛速度有一定的提升。

由表 2可以看出，当参数处理具有约束条件时，利用不等式约束的Partial EIV模型进行参数估计，可以得到比不考虑约束条件时总体最小二乘算法偏差更小的参数估值。由表 3的近似协方差计算结果可以看到，当考虑约束条件时，得到的参数估值具有更高的精度。

4.2 算例2(不等式约束的坐标转换模型)

算例2采用比平面拟合更为复杂的四参数坐标转换模型，系数矩阵含有2列固定量和2列随机量，可以表示为：

$ \begin{array}{l} \left[ {\begin{array}{*{20}{c}} {{X_1}}\\ {{Y_1}}\\ \vdots \\ {{X_n}}\\ {{X_n}} \end{array}} \right] - \left[ {\begin{array}{*{20}{c}} {{e_{{X_1}}}}\\ {{e_{{Y_1}}}}\\ \vdots \\ {{e_{{X_n}}}}\\ {{e_{{Y_n}}}} \end{array}} \right] = \left( {\left[ {\begin{array}{*{20}{c}} {{x_1}}&{ - {y_1}}&1&0\\ {{y_1}}&{{x_1}}&0&1\\ \vdots & \vdots & \vdots & \vdots \\ {{x_n}}&{ - {y_n}}&1&0\\ {{y_n}}&{{x_n}}&0&1 \end{array}} \right] - } \right.\\ \;\;\;\;\;\;\;\;\;\left. {\left[ {\begin{array}{*{20}{c}} {{e_{{x_1}}}}&{ - {y_1}}&0&0\\ {{e_{{y_1}}}}&{{x_1}}&0&0\\ \vdots & \vdots & \vdots & \vdots \\ {{e_{{x_n}}}}&{ - {e_{{y_n}}}}&0&0\\ {{e_{{y_n}}}}&{{e_{{x_n}}}}&0&0 \end{array}} \right]} \right)\left[ {\begin{array}{*{20}{l}} {{\beta _1}}\\ {{\beta _2}}\\ {{\beta _3}}\\ {{\beta _4}} \end{array}} \right] \end{array} $

(32)

本算例包括10个起始坐标点和10个目标坐标点，通过给定4个参数真值来联系2套坐标系，然后在坐标真值上添加随机误差进行平差解算。2套坐标对应的权阵分别为：

$ \begin{array}{*{20}{l}} {\;\;\;\;\;\;\;\;{\mathit{\boldsymbol{P}}_{XY}} = {\rm{diag}}([1, 2, 3, 1, 5, }\\ {4, 2, 7, 2, 1, 1, 2, 3, 1, 5, 4, 2, 7, 2, 1{]^{ - 1}})} \end{array} $

(33)

$ \begin{array}{*{20}{l}} {\;\;\;\;\;\;\;\;\;{\mathit{\boldsymbol{P}}_{xy}} = {\rm{diag}}([1, 3, 6, 1, 1, }\\ {8, 4, 3, 6, 5, 1, 3, 6, 1, 1, 8, 4, 3, 6, 5{]^{ - 1}})} \end{array} $

(34)

单位权方差为σ₀²=0.01，所添加随机误差的协方差阵为：

$ {\mathit{\boldsymbol{D}}_{XY}} = \sigma _0^2\mathit{\boldsymbol{P}}_{XY}^{ - 1}, {\mathit{\boldsymbol{D}}_{xy}} = \sigma _0^2\mathit{\boldsymbol{P}}_{xy}^{ - 1} $

(35)

坐标转换的模型数据见表 4，分别用不含约束条件时的总体最小二乘、附有不等式约束的Partial EIV模型与附有不等式约束的EIV模型求参数估值的结果$\mathit{\boldsymbol{\hat \beta }}$_TLS、$\mathit{\boldsymbol{\hat \beta }}$_ICPEIV、$\mathit{\boldsymbol{\hat \beta }}$_ICEIV、迭代次数k以及参数真值$\mathit{\boldsymbol{\tilde \beta }}$，见表 5，精度评定结果见表 6。

从表 5可见，用不含约束条件的总体最小二乘算法求得的估计值$\hat \beta $₂、$\hat \beta $₃、$\hat \beta $₄不满足不等式的约束条件，而本文算法可以很好地将参数估值控制在约束范围内。本文算法中，求解子问题的迭代次数为83次，母问题次数为6次，总次数为89次，而使用文献[10]中的ICEIV模型求解的迭代总次数为109次，本文算法在一定程度上减少了迭代次数。

由表 6可见，对$\hat \beta $₂、$\hat \beta $₃、$\hat \beta $₄进行一定的约束时，得到的参数估值比不考虑约束条件时总体最小二乘算法得到的精度更高。

5 结语

本文推导了求解不等式约束的Partial EIV模型的序列二次规划方法，并用WHP拟牛顿法求解最优估值，能够采用统一的方法计算结构性系数矩阵，有效减少迭代次数; 引入精度评定的SUT法，得到估值参数的近似协方差阵，其原理简单、适用性强，进一步完善了附有不等式约束的总体最小二乘问题的理论研究，是拟牛顿修正求解附有不等式约束总体最小二乘方法的一个重要补充，具有一定的理论意义。本文方法只适用于系数矩阵为列满秩的情况，对于秩亏的附有不等式约束的平差问题需要进一步研究。