2018, Vol. 38
表1(Table 1)
测量平差建模所得函数模型系数阵大多为稀疏矩阵，如仿射变换[1-2]、坐标转换[3-4]、IGS/CORS站坐标时序分析[5-6]、三维相机检校[7-8]等。事实上，上述平差模型均属多维观测建模问题，常规建模方法一般将全部待估参数表示为向量形式(简称向量参数法)，使得系数阵必须填充大量的零元素，导致其稠密度仅为1/d(d为观测维数)，并降低平差计算效率[9]。一些学者尝试研究矩阵分块降维[10-11]、复数平差建模[12]等快速建模新思路，但本质上仍属向量参数法。
对同一个测量平差问题，采用稀疏矩阵不及稠密矩阵有利于提高平差计算效率。已有学者探讨了多维观测变量的建模与最小二乘方法(简称矩阵参数法)[13]以及多维AR模型、多维灰色模型最小二乘[14]和仿射变换模型的总体最小二乘方法。这些研究主要针对某个特殊多变量观测模型开展，没有指出多变量观测的矩阵参数建模条件，也没有考虑多维观测之间的相关性。本文分析指出多维观测模型不同观测维度参数独立且个数相等特点(即独立同构特征)，满足该特征的矩阵参数建模方法可将系数矩阵化为完全稠密矩阵，从而利用Kronecker积运算性质[13]和加权最小二乘平差原理推导简便实用的矩阵参数及协因数阵最小二乘估计。与此同时，通过仿射变换模型、空间直线模型和IGS站坐标时序模型说明矩阵参数法建模过程，分析其在增大系数阵稠密度、降低法矩阵求逆阶数及提升解估计效率等方面的特点，并用算例验证新方法的有效性。
1 矩阵参数法平差模型及方法 1.1 矩阵参数法建模 1.1.1 仿射变换模型二维仿射变换模型为：
| $ \left\{ \begin{array}{l} {u_i} = {\beta _{11}} + {\beta _{12}}{x_i} + {\beta _{13}}{y_i}\\ {v_i} = {\beta _{21}} + {\beta _{22}}{x_i} + {\beta _{23}}{y_i} \end{array} \right. $ | (1) |
式中，(xi, yi)、(ui, vi)为变换前后坐标，β =[β11, β12, β13, β21, β22, β23]T为仿射变换参数。
按向量参数建模方法将式(1)写成方程组形式(待求参数用一维向量表示)，则：
| $ \left[ {\begin{array}{*{20}{c}} u\\ {{v_i}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} 1&{{x_i}}&{{y_i}}&0&0&0\\ 0&0&0&1&{{x_i}}&{{y_i}} \end{array}} \right]\mathit{\boldsymbol{\beta }} $ | (2) |
由式(2)可看出，利用参数向量法建模时，系数矩阵需补充一半零元素，导致其稠密度低，大大降低了建模效率和计算效率。仿射变换模型中不同维度的待估参数独立、个数相等，且系数矩阵各行向量的非零元素相同(不同维参数独立同构)。顾及该同构特征，将待估参数写成矩阵形式，则单点观测的仿射变换模型为：
| $ \left( {\begin{array}{*{20}{c}} {{u_i}}&{{v_i}} \end{array}} \right) = \left( {\begin{array}{*{20}{c}} 1&{{x_i}}&{{y_i}} \end{array}} \right)\left[ {\begin{array}{*{20}{c}} {{\beta _{11}}}&{{\beta _{21}}}\\ {{\beta _{12}}}&{{\beta _{22}}}\\ {{\beta _{13}}}&{{\beta _{23}}} \end{array}} \right] $ | (3) |
对比式(2)、式(3)，利用参数向量法进行二维仿射变换建模时，系数矩阵的稠密度仅为1/2，而矩阵参数法的系数矩阵不含零元素，为完全稠密矩阵，表明利用不同维之间的独立同构特征进行矩阵参数法建模，可以大大提高矩阵稠密度，从而提高计算效率。将式(3)拓展至多点观测，可得二维仿射变换模型的一般形式[13]:
| $ \left[ {\begin{array}{*{20}{c}} {{u_1}}&{{v_1}}\\ \vdots&\vdots \\ {{u_n}}&{{v_n}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} 1&{{x_1}}&{{y_1}}\\ \vdots&\vdots&\vdots \\ 1&{{x_n}}&{{y_n}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\beta _{11}}}&{{\beta _{21}}}\\ {{\beta _{12}}}&{{\beta _{22}}}\\ {{\beta _{13}}}&{{\beta _{23}}} \end{array}} \right] $ | (4) |
假设空间直线模型为：
| $ {x_i} = a + b{z_i}, {y_i} = c + d{z_i} $ | (5) |
式中，xi、yi、zi分别表示空间直线的3个方向，a、c、b、d为4个未知参数。按常规建模方法将式(5)写成方程组形式(待求参数用一维向量表示)，则空间直线的单点观测模型为：
| $ \left[ {\begin{array}{*{20}{c}} {{x_i}}\\ {{y_i}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} 1&{{z_i}}&0&0\\ 0&0&1&{{z_i}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} a\\ b\\ c\\ d \end{array}} \right] $ | (6) |
可以看出，向量参数法所得系数矩阵中零元素与非零元素个数各占一半，矩阵稠密度为1/2，表明大量无效的间接运算降低了空间直线模型的建模效率。综合式(5)、式(6)，该空间直线模型xi、yi分量之间的待估参数独立、个数相等，且系数矩阵各行向量的非零元素相同(不同维参数独立同构)。顾及这种同构特征，将待估参数重构为矩阵形式，则该空间直线的单点观测模型为：
| $ \left[ {\begin{array}{*{20}{c}} {{x_i}}&{{y_i}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} 1&{{z_i}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} a&c\\ b&d \end{array}} \right] $ | (7) |
对比式(6)、(7)，新方法有效利用不同维参数独立同构特征，完全排除无效间接运算，提高了空间直线模型的建模效率。将式(7)拓展至n个点，得空间直线模型矩阵参数法的一般形式：
| $ \left[ {\begin{array}{*{20}{c}} {{x_1}}&{{y_1}}\\ \vdots&\vdots \\ {{x_n}}&{{y_n}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} 1&{{z_1}}\\ \vdots&\vdots \\ 1&{{z_n}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} a&c\\ b&d \end{array}} \right] $ | (8) |
IGS站坐标各分量时序分析模型可以表示为[6]:
| $ \left\{ \begin{array}{l} x\left( {{t_i}} \right) = {a_x} + {b_x}{t_i} + \sum\limits_{j = 1}^m {\left( {{s_{x, j}}\sin \frac{{2\pi {t_i}}}{{{T_j}}} + {c_{x, j}}\cos \frac{{2\pi {t_i}}}{{{T_j}}}} \right)} + \sum\limits_{k = 1}^{{n_{x, k}}} {{g_{x, k}}H\left( {{t_i} - {t_{x, k}}} \right)} + {v_{x, {t_i}}}\\ y\left( {{t_i}} \right) = {a_y} + {b_y}{t_i} + \sum\limits_{j = 1}^m {\left( {{s_{y, j}}\sin \frac{{2\pi {t_i}}}{{{T_j}}} + {c_{y, j}}\cos \frac{{2\pi {t_i}}}{{{T_j}}}} \right)} + \sum\limits_{k = 1}^{{n_{y, k}}} {{g_{y, k}}H\left( {{t_i} - {t_{y, k}}} \right)} + {v_{y, {t_i}}}\\ z\left( {{t_i}} \right) = {a_z} + {b_z}{t_i} + \sum\limits_{j = 1}^m {\left( {{s_{z, j}}\sin \frac{{2\pi {t_i}}}{{{T_j}}} {c_{z, j}}\cos \frac{{2\pi {t_i}}}{{{T_j}}}} \right)} + \sum\limits_{k = 1}^{{n_{z, k}}} {{g_{z, k}}H\left( {{t_i} - {t_{z, k}}} \right)} + {v_{z, {t_i}}} \end{array} \right. $ | (9) |
式中，ti为时间，a*、b*为各分量线初始位置和运动速率，m为周期项个数，s*, j、c*, j为周期项系数，g*, k为阶跃项，n*, k表示阶跃项个数，v*, ti为误差项。下标表示x、y、z3个空间分量。
为叙述方便，假设3个分量坐标序列中均已进行阶跃探测、缺值插补和粗差剔除等预处理。按向量参数法将式(9)写成单历元观测模型：
| $ \left[ {\begin{array}{*{20}{c}} {x\left( {{t_i}} \right)}\\ {y\left( {{t_i}} \right)}\\ {z\left( {{t_i}} \right)} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{T}}_i}}&{\mathop 0\limits_{1 \times 7} }&{\mathop 0\limits_{1 \times 7} }\\ {\mathop 0\limits_{1 \times 7} }&{{\mathit{\boldsymbol{T}}_i}}&{\mathop 0\limits_{1 \times 7} }\\ {\mathop 0\limits_{1 \times 7} }&{\mathop 0\limits_{1 \times 7} }&{{\mathit{\boldsymbol{T}}_i}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{\alpha }}_x}}\\ {{\mathit{\boldsymbol{\alpha }}_y}}\\ {{\mathit{\boldsymbol{\alpha }}_z}} \end{array}} \right] $ | (10) |
其中，
由式(10)可得出与式(6)一致的结论，即利用向量参数法建模时，系数矩阵零元素个数远大于非零元素，矩阵稠密度为1/3，表明大量无效间接运算降低了站点三维时序建模效率。进一步分析可知，坐标时序模型中三维分量之间的待估参数独立且个数相等，系数矩阵各行向量的非零元素相同(不同维参数结构相同)。与式(7)同理，若顾及该同构特征，将待估参数重构为矩阵形式，则该坐标时序的单历元观测模型为：
| $ \left[ {\begin{array}{*{20}{c}} {x\left( {{t_i}} \right)}&{y\left( {{t_i}} \right)}&{z\left( {{t_i}} \right)} \end{array}} \right] = {\mathit{\boldsymbol{T}}_i}\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{\alpha }}_x}}&{{\mathit{\boldsymbol{\alpha }}_y}}&{{\mathit{\boldsymbol{\alpha }}_z}} \end{array}} \right] $ | (11) |
比较式(10)与(11)设计矩阵，前者系数矩阵的稠密度为1/3，后者系数矩阵TiT为完全稠密矩阵(稠密度为1)，表明矩阵参数法有效利用不同维参数独立同构特征, 完全排除无效间接运算，提高了空间直线模型建模效率。将式(11)拓展至n个历元，则矩阵参数法建立的站坐标时序模型为：
| $ \left[ {\begin{array}{*{20}{c}} {x\left( {{t_1}} \right)}&{y\left( {{t_1}} \right)}&{z\left( {{t_1}} \right)}\\ \vdots & \vdots & \vdots \\ {x\left( {{t_n}} \right)}&{x\left( {{t_n}} \right)}&{z\left( {{t_n}} \right)} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{T}}_1}}\\ \vdots \\ {{\mathit{\boldsymbol{T}}_n}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{\alpha }}_x}}&{{\mathit{\boldsymbol{\alpha }}_y}}&{{\mathit{\boldsymbol{\alpha }}_z}} \end{array}} \right] $ | (12) |
综合式(8)和式(12)，利用矩阵参数法可将d维单点或单历元观测模型表示为：
| $ \mathop {{{\mathit{\boldsymbol{V}}}_k}}\limits_{1 \times d} = \mathop {{{\mathit{\boldsymbol{B}}}_k}}\limits_{1 \times u} \mathop {\mathit{\boldsymbol{X}}}\limits_{u \times d} - \mathop {{{\mathit{\boldsymbol{L}}}_k}}\limits_{1 \times d} $ | (13) |
式中，X为待估参数，Bk为系数矩阵，Lk为观测值，Vk为观测值残差。下标k表示第k次观测。
将上式拓展到n个点或n个历元，并顾及d维观测模型不同维参数独立同构性，则d维观测模型的矩阵参数法平差函数模型为：
| $ \mathop {\mathit{\boldsymbol{V}}}\limits_{n \times d} = \mathop {\mathit{\boldsymbol{B}}}\limits_{n \times u} \mathop {\mathit{\boldsymbol{X}}}\limits_{u \times d} - \mathop {\mathit{\boldsymbol{L}}}\limits_{n \times d} $ | (14) |
式中，
分析式(14)，由于Bk稠密度为1，则系数矩阵B稠密度也为1，表明矩阵参数法所建模型的系数阵不含零元素，排除了无效存储和间接运算的可能。为进一步阐明矩阵参数法与向量参数法建模效率，按常规向量参数建模思路将式(14)进行误差方程的等价变换：
| $ \mathop {\mathit{\boldsymbol{\bar V}}}\limits_{\left( {d \cdot n} \right) \times 1} = \mathop {\mathit{\boldsymbol{\bar B}}}\limits_{\left( {d \cdot n} \right) \times \left( {d \cdot u} \right)} \mathop {\mathit{\boldsymbol{\bar X}}}\limits_{\left( {d \cdot u} \right) \times 1} - \mathop {\mathit{\boldsymbol{\bar L}}}\limits_{\left( {d \cdot n} \right) \times 1} $ | (15) |
式中，V=vec(V), L=vec(L), B=Id⊗B, X=vec(X)。⊗表示Kronecker积算子，vec(·)表示矩阵按列拉直算子。
比较矩阵参数法所建模型式(14)与向量参数法所建模型式(15)，前者待估参数以矩阵形式表示，后者待估参数以向量表示，这一差异直接导致前者系数阵B为稠密矩阵，后者系数阵B为稀疏矩阵且稠密度为1/d。因此，矩阵参数法可直接从建模层次提高平差计算效率。但是，常规最小二乘平差方法仅适用于式(15)，需探讨适用于式(14)的多维观测最小二乘平差方法。
设各维观测向量的协因数阵均为Ql，不同维观测之间协因数阵为Qd，则L的协因数阵为：
| $ {\mathit{\boldsymbol{Q}}_{\bar L}} = {\mathit{\boldsymbol{Q}}_d} \otimes {\mathit{\boldsymbol{Q}}_l} $ | (16) |
结合式(15)、(16)，根据VTPLV =min可得常规向量参数法的待估参数、协因数阵及单位权方差估计式为：
| $ \widehat {\mathit{\boldsymbol{\bar X}}} = {\mathit{\boldsymbol{Q}}_{\bar X}}{\mathit{\boldsymbol{\bar B}}^{\rm{T}}}{\mathit{\boldsymbol{P}}_{\bar L}}\mathit{\boldsymbol{\bar L}} $ | (17) |
| $ {\mathit{\boldsymbol{Q}}_{\bar X}} = {\left( {{{\mathit{\boldsymbol{\bar B}}}^{\rm{T}}}{\mathit{\boldsymbol{P}}_{\bar L}}\mathit{\boldsymbol{\bar B}}} \right)^{ - 1}} $ | (18) |
| $ \hat \sigma _0^2 = \frac{{{{\widehat {\mathit{\boldsymbol{\bar V}}}}^{\rm{T}}}{\mathit{\boldsymbol{P}}_{\bar L}}\widehat {\mathit{\boldsymbol{\bar V}}}}}{{d\left( {n - u} \right)}} $ | (19) |
其中，PL=(Qd⊗Ql)-1。
对于矩阵
| $ {\left( {\mathit{\boldsymbol{\tilde A}} \otimes \mathit{\boldsymbol{\tilde B}}} \right)^{ - 1}} = {\mathit{\boldsymbol{\tilde A}}^{ - 1}} \otimes {\mathit{\boldsymbol{\tilde B}}^{ - 1}} $ | (20) |
| $ {\left( {\mathit{\boldsymbol{\tilde A}} \otimes \mathit{\boldsymbol{\tilde B}}} \right)^T} = {\mathit{\boldsymbol{\tilde A}}^T} \otimes {\mathit{\boldsymbol{\tilde B}}^T} $ | (21) |
| $ \left( {\mathit{\boldsymbol{\tilde A}} \otimes \mathit{\boldsymbol{\tilde B}}} \right)\left( {\mathit{\boldsymbol{\tilde C}} \otimes \mathit{\boldsymbol{\tilde D}}} \right) = \left( {\mathit{\boldsymbol{\tilde A\tilde C}}} \right) \otimes \left( {\mathit{\boldsymbol{\tilde B\tilde D}}} \right) $ | (22) |
利用以上Kronecker积运算性质对式(17)~(19)进行化简，可得待估参数、协因数阵及单位权方差估计新公式：
| $ \begin{array}{l} \;\;\;\;\;\;\;\widehat {\mathit{\boldsymbol{\bar X}}} = {\mathit{\boldsymbol{Q}}_{\bar X}}{\left( {{\mathit{\boldsymbol{I}}_d} \otimes \mathit{\boldsymbol{B}}} \right)^{\rm{T}}}\left( {\mathit{\boldsymbol{Q}}_d^{ - 1} \otimes \mathit{\boldsymbol{Q}}_l^{ - 1}} \right)\mathit{\boldsymbol{\bar L}} = \\ \left[ {{\mathit{\boldsymbol{Q}}_d} \otimes {{\left( {{\mathit{\boldsymbol{B}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_l^{ - 1}\mathit{\boldsymbol{B}}} \right)}^{ - 1}}} \right]\left( {{\mathit{\boldsymbol{I}}_d} \otimes {\mathit{\boldsymbol{B}}^{\rm{T}}}} \right)\left( {\mathit{\boldsymbol{Q}}_d^{ - 1} \otimes \mathit{\boldsymbol{Q}}_l^{ - 1}} \right)\mathit{\boldsymbol{\bar L}} = \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\left\{ {{\mathit{\boldsymbol{I}}_d} \otimes \left[ {{{\left( {{\mathit{\boldsymbol{B}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_l^{ - 1}\mathit{\boldsymbol{B}}} \right)}^{ - 1}}{\mathit{\boldsymbol{B}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_l^{ - 1}} \right]} \right\}\mathit{\boldsymbol{\bar L}} \end{array} $ | (23) |
| $ \begin{array}{l} {\mathit{\boldsymbol{Q}}_{\bar X}} = {\left[ {{{\left( {{I_d} \otimes \mathit{\boldsymbol{B}}} \right)}^{\rm{T}}}\left( {\mathit{\boldsymbol{Q}}_d^{ - 1} \otimes \mathit{\boldsymbol{Q}}_l^{ - 1}} \right)\left( {{\mathit{\boldsymbol{I}}_d} \otimes \mathit{\boldsymbol{B}}} \right)} \right]^{ - 1}} = \\ \;\;\;\;{\left[ {\left( {{\mathit{\boldsymbol{I}}_d} \otimes {\mathit{\boldsymbol{B}}^{\rm{T}}}} \right)\left( {\mathit{\boldsymbol{Q}}_d^{ - 1} \otimes \mathit{\boldsymbol{Q}}_l^{ - 1}} \right)\left( {{\mathit{\boldsymbol{I}}_d} \otimes \mathit{\boldsymbol{B}}} \right)} \right]^{ - 1}} = \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\mathit{\boldsymbol{Q}}_d} \otimes {\left( {{\mathit{\boldsymbol{B}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_l^{ - 1}\mathit{\boldsymbol{B}}} \right)^{ - 1}} \end{array} $ | (24) |
| $ \begin{array}{l} \hat \sigma _0^2 = \frac{{{\rm{vec}}{{\left( \mathit{\boldsymbol{V}} \right)}^{\rm{T}}}\left( {\mathit{\boldsymbol{Q}}_d^{ - 1} \otimes \mathit{\boldsymbol{Q}}_l^{ - 1}} \right){\rm{vec}}\left( \mathit{\boldsymbol{V}} \right)}}{{d\left( {n - u} \right)}} = \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\frac{{{\rm{tr}}\left( {\mathit{\boldsymbol{Q}}_d^{ - 1}{\mathit{\boldsymbol{V}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_i^{ - 1}\mathit{\boldsymbol{V}}} \right)}}{{d\left( {n - u} \right)}} \end{array} $ | (25) |
式中，tr(·)表示求迹算子。顾及向量参数X与矩阵参数X、观测向量L与观测矩阵L的分块矩阵关系，可将式(23)进一步简化，得到矩阵参数平差公式：
| $ \mathit{\boldsymbol{\hat X}} = {\left( {{\mathit{\boldsymbol{B}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_l^{ - 1}\mathit{\boldsymbol{B}}} \right)^{ - 1}}{\mathit{\boldsymbol{B}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_l^{ - 1}\mathit{\boldsymbol{L}} $ | (26) |
式(24)~(26)即为矩阵参数加权最小二乘平差方法，简称矩阵参数法。比较其与参数向量法，参数估值与多维观测互相关性Qd无关，参数协因数阵估值与Qd有关。其次，矩阵参数
算例1 平差函数模型选用§1.1.2的空间直线模型，其中a、c、b、d分别为2、1、1、-1，第三维分量z取[1: 1: 200]。利用式(5)仿真空间直线数据，分别加入中误差为10的白噪声(相关系数为0)和相关噪声(相关系数为0.5)，则协因数阵Qd分别为单位阵(白噪声)和[1, 0.5;0.5, 1](相关噪声)，协因数阵Ql均为单位阵。比较方案如下：方案A，利用矩阵参数法进行平差，其中Qd和Ql如前所述；方案B，利用向量参数法进行平差，其中QL=Qd⊗Ql。
为使结果具有统计意义，仿真1 000次，利用中值法统计参数估值、参数中误差、单位权中误差、耗时比(矩阵参数法和向量参数法计算时间的比值)，结果见表 1。
|
表 1 算例1不同方案统计结果 Tab. 1 Statistical results of two schemes in example 1 |
由表 1不难看出，无论所加噪声为白噪声还是相关噪声，向量参数法与矩阵参数法均具有相同的参数估值和中误差估值，两种算法等价。同时，相同噪声条件下不同维的参数中误差完全相同，从而验证了多维观测模型参数的独立同构特征。此外，矩阵参数法提高了计算效率，其计算时间仅为向量参数法的44.23%和46.47%。
算例2 IGS站坐标时序选用SOPAC(http://sopac-ftp.ucsd.edu/pub/timeseries/measures/)提供的ITRF2008框架下GUAO站和ONSA站2002.5y~2009.5y原始坐标序列(样本长度2 558，文件格式RawNEU，单位cm)。
基于预处理坐标序列和模型(9)，分别采用向量参数法和矩阵参数法进行平差计算(取年周期和半年周期)，并统计参数估值、参数中误差、单位权中误差、耗时比。两种计算方案为：方案C，采用矩阵参数法进行平差计算，Ql取高程方向的先验方差构成对角阵，Qd取各方向先验方差与高程方向先验方差比值构成的对角阵；方案D，采用向量参数法进行平差计算，其中QL=Qd⊗Ql。为节省篇幅，仅显示GUAO站的预处理坐标序列、残差序列(图 1)和平差结果(表 2)。
|
图 1 GUAO站坐标时序和残差时序 Fig. 1 Time series and residual time series of station GUAO |
|
|
表 2 算例2不同方案统计结果 Tab. 2 Statistical results of two schemes in example 2 |
由表 2可以看出，对于IGS坐标时序模型平差，矩阵参数法和矩阵向量法得出相同的参数估值和参数中误差，进一步验证了矩阵参数法的正确性。此外，与表 1相比，由于坐标时序历元数远高于仿真数据长度，使系数矩阵、法矩阵阶数更高，故矩阵参数法对计算效率的提高更加明显，其计算时间仅为向量参数法的8.97%和9.25%。
综合算例1和算例2可得，矩阵参数法与向量参数法平差结果相同，两种方法等价。同时，由于算例1数据量远少于算例2，故矩阵参数法效率提高更加显著。
3 结语本文提出顾及多维观测模型参数独立同构的矩阵参数建模与加权最小二乘平差方法。该方法将仿射变换模型、空间直线模型和IGS站坐标模型进行多维观测建模，系数阵转化为完全稠密矩阵，建模简单高效。
| [1] |
李博峰. 无缝仿射基准转换模型的方差分量估计[J]. 测绘学报, 2016, 45(1): 30-35 (Li Bofeng. Variance Component Estimation in the Seamless Affine Transformation Model[J]. Acta Geodaetica et Cartographica Sinica, 2016, 45(1): 30-35)
(  1) |
| [2] |
刘志平, 余前勇, 査剑锋. 空间直角坐标至两类常用坐标的快速变换[J]. 测绘科学, 2015, 40(3): 8-11 (Liu Zhiping, Yu Qianyong, Zha Jianfeng. Fast Coordinate Transformations for Both XYZ-BLH and XYZ-RhA[J]. Science of Surveying and Mapping, 2015, 40(3): 8-11)
( 1) |
| [3] |
Amiri-Simkooei A R. Non-Negative Least-Squares Variance Component Estimation with Application to GPS Time Series[J]. Journal of Geodesy, 2016, 90(5): 1-16
( 1) |
| [4] |
曹良中, 杨辽, 李欣, 等. 基于附加参数光束法平差的非量测型数码相机检校研究[J]. 遥感信息, 2014, 29(3): 14-18 (Cao Liangzhong, Yang Liao, Li Xin, et al. Calibration of a Nonmetric Camera Based on Bundle Adjustment with Additional Parameters[J]. Remote Sensing Information, 2014, 29(3): 14-18)
( 1) |
| [5] |
周凌焱, 刘成龙, 张强, 等. 稀疏矩阵压缩技术在CPⅢ网平差中的应用研究[J]. 铁道科学与工程学报, 2014, 11(6): 142-147 (Zhou Lingyan, Liu Chenglong, Zhang Qiang, et al. Research on Application of Sparse Matrix to the Compression Technology in the Adjustment of CPⅢ Network[J]. Journal of Railway Science and Engineering, 2014, 11(6): 142-147)
( 1) |
| [6] |
宋力杰, 欧阳桂崇. 超大规模大地网分区平差快速解算方法[J]. 测绘学报, 2003, 32(3): 204-207 (Song Lijie, Ouyang Guichong. A Fast Method of Solving Partitioned Adjustment for Super Large-scale Geodetic Network[J]. Acta Geodaetica et Cartographica Sinica, 2003, 32(3): 204-207)
( 2) |
| [7] |
刘志平, 李思达. 复数域与实数域最小二乘平差的等价性研究[J]. 大地测量与地球动力学, 2016, 36(8): 741-744 (Liu Zhiping, Li Sida. Research on Equivalence between Complex-Valued and Real-Valued Least Squares Adjustment Method[J]. Journal of Geodesy and Geodynamics, 2016, 36(8): 741-744)
( 1) |
| [8] |
Koch K R. Parameter Estimation and Hypothesis Testing in Linear Models[M]. Springer-Verlag, 1999
( 1) |
| [9] |
张朝玉. 多维AR序列的最小二乘建模方法[J]. 武汉大学学报:信息科学版, 2002, 27(4): 377-381 (Zhang Chaoyu. Multi-Dimensional AR Series Modeled by Least Square Criterion[J]. Geomatics and Information Science of Wuhan University, 2002, 27(4): 377-381)
( 1) |
| [10] |
刘志平, 何秀凤. 无偏扩展灰色模型及其在高边坡形变预测中的应用[J]. 大地测量与地球动力学, 2007, 28(1): 41-44 (Liu Zhiping, He Xiufeng. Unbiased Estimation of Extended Grey Model and Its Application in Deformation Prediction of Steep Slopes[J]. Journal of Geodesy and Geodynamic, 2007, 28(1): 41-44)
( 1) |
| [11] |
王乐洋, 赵英文, 陈晓勇, 等. 多元总体最小二乘问题的牛顿解法[J]. 测绘学报, 2016, 45(4): 411-417 (Wang Leyang, Zhao Yingwen, Chen Xiaoyong, et al. A Newton Algorithm for Multivariate Total Least Squares Problems[J]. Acta Geodaetica et Cartographica Sinica, 2016, 45(4): 411-417 DOI:10.11947/j.AGCS.2016.20150246)
( 1) |
| [12] |
马超, 李斐, 张胜凯, 等. 顾及共性误差的南极半岛地区连续GPS站坐标时间序列分析[J]. 地球物理学报, 2016, 59(8): 2783-2795 (Ma Chao, Li Fei, Zhang Shengkai, et al. The Coordinate Time Series Analysis of Continuous GPS Stations in the Antarctic Peninsula with Consideration of Common Mode Error[J]. Chinese Journal Geophysics, 2016, 59(8): 2783-2795 DOI:10.6038/cjg20160806)
( 1) |
| [13] |
蒋志浩, 张鹏, 秘金钟, 等. 顾及有色噪声影响的CGCS2000下我国CORS站速度估计[J]. 测绘学报, 2010, 39(4): 355-363 (Jiang Zhihao, Zhang Peng, Bei Jinzhong, et al. Velocity Estimation on the Colored Noise Properties of CORS Network in China Based on the CGCS 2000 Frame[J]. Acta Geodaetica et Cartographica Sinica, 2010, 39(4): 355-363)
( 3) |
| [14] |
武艳强, 黄立人. 时间序列处理的新插值方法[J]. 大地测量与地球动力学, 2004, 24(4): 43-47 (Wu Yanqiang, Huang Liren. A New Interpolation Method in Time Series Analyzing[J]. Journal of Geodesy and Geodynamics, 2004, 24(4): 43-47)
( 1) |