设有如下的线性系统：

式中，$\mathit{\boldsymbol{\tilde x}}\mathit{\boldsymbol{、\tilde y}} $分别表示q×1、m×1阶的实数输入和输出向量的真值；$\mathit{\boldsymbol{A}}\left( {\mathit{\boldsymbol{\tilde x}}} \right) $表示m×p阶映射矩阵，即其元素为常数或向量${\mathit{\boldsymbol{\tilde x}}} $的函数；ξ为p×1阶的参数向量。在实际问题中，往往只能得到含误差的输入输出变量，设观测值满足以下关系：

式中，x、y分别表示输入和输出向量的观测值；e_x、e_y分别表示相应的误差向量，并满足如下的随机模型：

式中，N[μ, D]表示数学期望为μ、方差为D的正态分布。很显然，A(x)的元素可以是观测值的线性或非线性形式，本文仅考虑A(x)为观测值向量的线性组合的情况。目前TLS方法应用中的直线拟合、仿射变换、Bursa变换与自回归模型参数估计等A矩阵中含误差的元素均为观测值的线性函数。在测量平差问题中，由于x、y向量均含有误差的模型属于EIV模型，宜采用WTLS方法求解，现有的WTLS平差方法可表示为求解以下EIV问题：

式中，A表示含有随机误差E_A的m×p阶设计矩阵；Q_B⁺表示协因数矩阵Q_B的广义逆矩阵(由于B矩阵常包含常数或重复元素，因此Q_B多为奇异矩阵)。A中含有常数的EIV模型也被称为部分EIV模型^[20]。在大地测量中，含有误差的输入输出向量可等价为观测值，其先验误差可以近似或精确给定，根据误差传播原理，当设计矩阵中的随机元素是观测向量的线性组合时，可得到如下的等价随机模型：

Q_B表示增广设计矩阵的协因数阵，可以通过以下方法计算^[8]：

式中，$\mathit{\boldsymbol{J = }}{\left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{J}}_{11}^{\rm{T}}}& \cdots &{\mathit{\boldsymbol{J}}_{m1}^{\rm{T}}}& \cdots &{\mathit{\boldsymbol{J}}_{1p}^{\rm{T}}}& \cdots &{\mathit{\boldsymbol{J}}_{mp}^{\rm{T}}} \end{array}} \right]^{\rm{T}}} \in {{\bf{R}}^{mp \times q}} $，q表示观测值向量的维数；J_ij表示矩阵元素A_ij对应于向量x的Jacobi矩阵，即：

显然，当A的随机元素可表示为x的线性形式时，J为常数矩阵。若向量x和y相关，即Q_xy≠0，A和y也相关。

在测量平差问题中，很多情况下，参数向量ξ需满足一定的限制条件，限制条件用来提供观测系统的先验误差信息，或表示参数间的物理或几何约束，或用来减弱观测系统的病态，通常可以表示为等式条件或不等式约束形式^[13]。本文仅考虑等式约束，设约束条件的一般形式可表示为:

式中, c(ξ)表示s×1的函数向量，表示约束条件由s个独立函数组成。得到附限制条件的WTLS平差的通用模型：

根据拉格朗日乘数法原理，构造拉格朗日函数：

式中，λ、μ分别表示m×1、s×1阶拉格朗日乘数向量。为方便表达，令$ \mathit{\boldsymbol{\eta = }}{\left( {{\bf{ \pmb{\mathit{ ξ}} }} - {\bf{1}}} \right)^{\rm{T}}} \in {{\bf{R}}^{\mathit{p + }{\rm{1}}}}$，不考虑限制条件的WTLS问题等价于以下无约束条件的极值问题^[8]：

其中，⊗表示Kronecker积，得到考虑限制条件的目标函数为：

式(12)是一个无约束条件的非线性函数的极值问题，可以用常用的数值解算方法，如梯度法、拟牛顿法、牛顿法、Levenberg-Marquardt法求解^[21]。由于牛顿法具有二阶收敛速度，本文采用牛顿迭代法，为此需要对式(12)求二阶导数。为方便表达，分别对f(ξ)和c(ξ)求导数，得到：

式中，J_f∈R^p×1、H_f∈R^p×p分别表示函数f(ξ)的Jacobi向量和Hessian矩阵，其解析表达式在文献[8]中已给出。限制条件是一个由若干函数组成的向量，利用矩阵求导公式：

式中，R、S分别表示含有未知参数Γ的矩阵；Γ表示s×t的参数矩阵。由此得到限制条件的一阶和二阶导数分别为：

式中，J_c∈R^p×s、H_c∈R^p×ps分别表示函数c(ξ)的一阶、二阶导数。若用J和H分别表示式(12)中的函数l(ξ, μ)相对于参数[ξ^T    μ^T]^T向量的梯度向量和Hessian矩阵，则：

根据牛顿迭代法的基本原理，设广义参数向量为$\mathit{\boldsymbol{\overline \xi }} = {\left[ {{\mathit{\boldsymbol{\xi }}^{\rm{T}}}{\mathit{\boldsymbol{\mu }}^{\rm{T}}}} \right]^{\rm{T}}} $，在已知参数当前值的情况下，得到参数的改正数为：

以上给出了算法的解析表达式，需要经数值迭代实现。附有限制条件的EIV模型的WTLS平差方法归纳如下：

1) 计算参数的初值${\mathit{\boldsymbol{\hat \xi }}^{\left( 0 \right)}}、{\mathit{\boldsymbol{\hat \mu }}^{\left( 0 \right)}} $，由式(6)计算增广设计矩阵的协方差阵Q_B, 并设初始化迭代次数i=0。

2) 按以下顺序计算：$\mathit{\boldsymbol{\eta = }}{\left[ {{\mathit{\boldsymbol{\xi }}^{\left( i \right)}} - {\bf{1}}} \right]^{\rm{T}}} $，$\mathit{\boldsymbol{R = }}\left( {{{\mathit{\boldsymbol{\hat \eta }}}^{\rm{T}}}{\mathit{\boldsymbol{I}}_m}} \right){\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{B}}}{\left( {\mathit{\boldsymbol{\hat \eta }} \otimes {\mathit{\boldsymbol{I}}_\mathit{m}}} \right)^{ - 1}} $，$\mathit{\boldsymbol{u = RB\hat \eta }} $，${\mathit{\boldsymbol{N}}_\mathit{i}} = \left( {\mathit{\boldsymbol{e}}_i^{\rm{T}}{\mathit{\boldsymbol{I}}_\mathit{m}}} \right)\mathit{\boldsymbol{QB}}\left( {\mathit{\boldsymbol{\hat \eta }} \otimes {\mathit{\boldsymbol{I}}_\mathit{m}}} \right) + \left( {{{\mathit{\boldsymbol{\hat \eta }}}^{\rm{T}}}{\mathit{\boldsymbol{I}}_\mathit{m}}} \right)\mathit{\boldsymbol{QB}}\left( {{\mathit{\boldsymbol{e}}_i} \otimes {\mathit{\boldsymbol{I}}_\mathit{m}}} \right)$(e_i为p+1的单位向量，除第i个元素为1外，其他的元素均为零)。

式中，i=1…m，j=1…p；[J_f]_i表示向量J_f的第i行元素；[H_f]_ij表示矩阵H_f的第i行、第j列元素；由式(15)计算限制条件的J_c、H_c，并进一步计算参数的改正数：

3) 判定$\frac{{\left\| {\delta \mathit{\boldsymbol{\hat {\overline \xi} }} } \right\|}}{{\left\| {\mathit{\boldsymbol{\hat {\overline \xi} }} } \right\|}} \le \varepsilon $(ε为给定的阈值)是否成立，若成立则结束计算；否则更新参数${\mathit{\boldsymbol{\hat \xi }}^{\left( {i + 1} \right)}} = {\mathit{\boldsymbol{\hat \xi }}^{\left( i \right)}} + \mathit{\delta }\mathit{\boldsymbol{\hat \xi }} $，${\mathit{\boldsymbol{\hat \mu }}^{\left( {i + 1} \right)}} = {\mathit{\boldsymbol{\hat \mu }}^{\left( i \right)}} + \mathit{\delta }\mathit{\boldsymbol{\hat \mu }} $，并令迭代次数i：=i+1，继续执行步骤2)。

上文给出了附有限制条件的EIV模型的WTLS平差方法，由于需要求目标函数的二阶导数，导致算法的解析表达式较为复杂，但算法清晰地给出了一阶、二阶偏导的解析表达式，易于编程实现，而且牛顿算法具有较高的收敛速度。本文算法可解算与设计矩阵和观测向量相关的问题，适用于任何形式的限制条件或无限制条件的EIV模型的平差问题，是一种通用的WTLS平差方法。

相关研究已证明总体最小二乘目标函数是一个非凸函数，初值的选取将影响收敛速度，甚至决定能否准确收敛。由于测量平差问题属于超定问题而且误差通常较小，因此可以利用部分观测值确定参数$\mathit{\boldsymbol{\hat \xi }} $初值，然后通过将附有限制条件的EIV模型线性化为附有限制条件的条件问题确定联系数$\mathit{\boldsymbol{\hat \mu }} $的初值，具体的方法在有关测量平差的教材中已有详细推导，本文不再赘述。

精度评定是测量平差的必要步骤，诸多文献已证明WTLS平差方法是有偏的^{[4, 7, 20]}。在小误差情况下，参数的估值偏差极小可以忽略。文献[8]将EIV模型线性化为附有参数的条件平差法问题来进行近似的精度评定。由于本文采用的牛顿迭代算法在计算过程中已经得到了目标函数的一阶、二阶导数，因此可以在精度评定中加以利用。本文提出一种新的近似精度评定方法，适用于小误差情况。首先计算单位权方差的估值：

式中，m-p+s表示多余观测数，也称为平差系统的自由度。在EIV模型中，参数的误差是由[Ay]的误差引起的，若不考虑参数的误差，式(20)中只有J_f受观测值误差的影响：

式中，ΔJ_f表示ΔB对J_f产生的影响量。将式(22)代入式(20)中，并运用误差传播定律得到：

其中，

本文提出的附有限制条件的平差方法适用于设计矩阵A是观测值的线性组合、A和y相关的EIV模型，而且限制条件可以是任意非线性函数的情况。函数模型、随机假设、限制条件均可以拓展或简化以满足特定应用。如在房屋重建问题中，限制条件通常只包含一个二次型限制条件，文献[14]中附有线性限制条件和二次型条件也可以简化为二次型限制条件，设二次型限制条件的通式为：

式中，M是一个p×p的方阵，按照本文方法可以得到以下的迭代公式：

类似地可以得到相应的数值计算方法，通过模型和算法的简化，可以拓展EIV模型和WTLS方法的应用。

本文算例模拟某矩形房屋的重建问题，设房屋长10 m、宽5 m, AB边与Y轴正向的夹角为π/6，A点的坐标为(10, 10)，沿长度方向测量10个点，沿宽度方向测量5个点，其坐标的理论值分别沿长度和宽度方向等距离选取。为模拟观测值的不等精度问题，AB、BC、CD、DA边的x、y坐标的误差分别为±0.2、±0.4、±0.6、±0.8 m，同一点的x、y坐标的相关系数在[-0.5，0.5]之间随机产生，用于模拟设计矩阵和观测向量相关的情况。将模拟得到的测量误差和坐标理论值相加得到观测值。设直线AB和BC的斜率分别为k_u和k_v，4条边的截距分别为b₁、b₂、b₃、b₄。为保证房屋的两条邻边互相垂直，斜率需满足k_uk_v=-1的限制条件。由此得到附有二次型限制条件的EIV模型，其中各矩阵可表示为：

式中，(x_ABⁱ, y_ABⁱ)表示线段AB上第i个点的坐标值，其他边上的测量点的表达方法类似。至此，二维建筑物重建问题已经被表示为附有限制条件的EIV模型，和文献[18]的传统平差方法相比，本文方法构建简单，计算时无需对非线性的条件方程线性化。表 1给出了一组模拟数据及其误差，将其作为观测值组成附有限制条件的EIV模型。很显然，上述问题属于附有线性和二次型限制条件的WTLS方法的特殊形式，将其用文献[14]方法求解，不能准确收敛。由此可见，文献[14]方法在某些情况下不能适用，而用本文方法在经过5次迭代后能准确收敛，从而说明本文方法是一种有效且通用的方法。

考虑二次型限制条件的WTLS平差方法与不考虑限制条件，采用WTLS方法单独拟合4条边界的结果见表 2，其中$\mathit{\hat a}、\mathit{\hat b} $分别表示各边界的斜率和截距，Δa、Δb表示估计参数与相应的真值间的差值。然后通过求出直线交点得到房屋角点坐标，上述两种方法的重建效果见图 1。从表 2的结果可以看出，附二次型限制条件的估计结果与真值的最大差值在0.4 m之内，而不考虑正交限制条件的单独拟合方法与真值的最大差值达7.2 m，说明前者的重建精度明显高于后者。增加模拟次数，比较参数的均方根误差也可以得出相同的结论，其主要原因是：附限制条件的WTLS方法顾及了矩形参数间的几何约束条件，并将所有的点整体参与平差，从而增加了整个平差系统的多余观测数，因而提高了重建精度。从图 1中可以看出，附限制条件的重建方法能保持矩形的形状，而无限制条件的WTLS方法无法保持，说明了附有限制条件的WTLS方法的应用价值。

理论研究方面，本文针对基于遥感数据的二维建筑物重建问题，考虑观测向量和设计矩阵相关的情况，推导了附限制条件的EIV模型的WTLS平差方法，并将二次型限制条件作为限制条件的特例给出了相应的算法和近似精度评定方法。本文采用的牛顿迭代算法具有二阶收敛速度，研究成果对于完善WTLS平差方法、拓展总体最小二乘方法的应用有一定的理论价值。本文仅考虑设计矩阵是观测值的线性组合的情况，而且采用的牛顿算法需要计算目标函数的二阶导数，从而增加了解析表达式和相应算法的复杂程度，所采用的精度评定方法是近似方法，只适用于小误差的情况，以上问题尚需进一步研究和完善。

在应用方面，本文以基于LiDAR或高分遥感数据的二维直角建筑物高精度重建为应用背景，研究了EIV模型构建及相应的WTLS平差算法。和现有平差方法相比，本文方法具有建模简洁、解算精度高的优势，可进一步拓展至复杂建筑的三维重建问题。本文研究成果在测绘技术从以确定点的坐标和高程为主的单点测量模式逐步转换为以测量二维或三维几何模型为主的信息测绘模式、遥感技术大规模应用的技术背景下，具有重要的应用价值。

图  1  有限制条件和无限制条件的WTLS平差方法的矩形房屋重建效果

Figure 1.  Performance of Constrained and Unconstrained WTLS Algorithms for Building Reconstruction