甚长基线干涉测量(very long baseline interferometry，VLBI)、全球卫星导航系统（Global Navigation Satellite System，GNSS)和卫星激光测距(satellite laser ranging，SLR)是确定地球自转参数(Earth rotation parameters, ERP)的主要技术。ERP是一个随时间变化的物理量，不同处理方法、软件可得到的ERP往往有较大差异，这种ERP的差异给用户使用带来不便。因此，需要对不同技术、不同分析中心得到的ERP进行融合处理，以得到唯一的ERP解。

1994年，Blewitt等提出了一种数据交换格式（solution independent exchange format, SINEX）。基于SINEX格式的ERP融合模型主要有两类，一类是约化其他参数仅单独融合ERP参数^[1-2]，另一类则将ERP与测站坐标等其他参数统一进行融合^[3-8]，不同融合模型解的性质略有差异^[9]。1998年，喷气推进实验室提出一种ERP融合的卡尔曼（Kalman）滤波方法，把ERP变化看成一种随机过程，利用Kalman滤波融合不同时期的ERP得到统一解^[1]。为消除高频噪声，Vondrák等^[2]提出一种平滑组合方法。Fernandez等^[10]提出一种顾及频率噪声的天尺度ERP序列的融合方法。为确保ERP与地球参考框架的一致性，2003年起学者们基于多空间大地测量技术进行测站坐标和ERP时间序列的组合研究^[3-7]，这一研究直接促成新的国际地球参考框架ITRF2005^[8]（International Terestrial Reference Frame 2005）的建立。2007年，有学者提出了基于法方程的测站坐标、ERP、对流层参数的融合方法^[11]。目前，ERP官方产品IERS C04（International Earth Rotation and Reference Systems Service Combined 04）是基于多技术ERP的融合结果^{[1, 10, 12]}。IERS C04以不同技术结果EOP(Earth orientation parameters)序列作为输入，采用三角帽法确定各技术相对权因子，通过滑动平均、高频滤波、插值等过程，综合求得EOP序列解^[13]。我国学者也开展了ERP融合研究，基于多分析中心结果组合了全球定位系统(Global Positioning System, GPS)站坐标和ERP^[14-17]。

本文从批处理角度，提出了一种多分析中心ERP结果的附加边界约束和转换参数内约束融合模型，其基本思路是：通过参数变换模型把各分析中心的结果转换到相同时刻，对相邻观测时段边界点处的ERP施加边界约束，对各分析中心的长期解增加的转换参数施加内约束，分两步得到ERP的融合解。

1 多时段ERP序列的边界约束模型构建

设某一大地测量观测技术第i个观测时段数据（其观测开始时刻为t₁ⁱ，技术时刻为t₂ⁱ，中心时刻为t₁₂ⁱ）所建立的误差方程为：

$ {\mathit{\boldsymbol{V}}^i} = \left( {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{A}}_1^i}&{\mathit{\boldsymbol{A}}_2^i}&{\mathit{\boldsymbol{A}}_3^i} \end{array}} \right)\left( {\begin{array}{*{20}{c}} {\hat x_{12}^i}\\ {\hat x_{12}^i}\\ {\hat y} \end{array}} \right) - {\mathit{\boldsymbol{L}}^\mathit{i}} $

(1)

式中，Vⁱ为误差向量；Lⁱ为该观测时段观测向量；$ \hat x_{12}^i$、$\widehat {\dot x}_{12}^i $为该观测时段t₁₂ⁱ时刻的ERP及其速率；$ {\hat y}$为其他待估参数；$A_1^i $、$A_2^i $、$A_3^i $则为其相应的设计矩阵。

约化掉其他待估参数${\hat y} $，则该观测时段i仅保留$ \hat x_{12}^i$、$\widehat {\dot x}_{12}^i$的法方程为：

$ {{\mathit{\boldsymbol{\tilde N}}}_i}\left( {\begin{array}{*{20}{c}} {\hat x_{12}^i}\\ {\widehat {\dot x}_{12}^i} \end{array}} \right) = {{\mathit{\boldsymbol{\tilde U}}}^i} $

(2)

式中，${{\mathit{\boldsymbol{\tilde N}}}_i} = \left( {\begin{array}{*{20}{l}} {\tilde N_{11}^i}&{\tilde N_{12}^i}\\ {\tilde N_{21}^i}&{\tilde N_{22}^i} \end{array}} \right);{{\mathit{\boldsymbol{\tilde U}}}^i} = \left( {\begin{array}{*{20}{l}} {\tilde U_1^i}\\ {\tilde U_2^i} \end{array}} \right)$。

时刻t₁₂ⁱ的x₁₂ⁱ、$ \widehat {\dot x}_{12}^i$可以表示为该时段i的边界值x₁ⁱ、x₂ⁱ的形式，即：

$ \left( {\begin{array}{*{20}{l}} {x_{12}^i}\\ {\dot x_{12}^i} \end{array}} \right) = {\mathit{\boldsymbol{C}}_i}\left( {\begin{array}{*{20}{l}} {x_1^i}\\ {x_2^i} \end{array}} \right) $

(3)

式中，$ {\mathit{\boldsymbol{C}}_i} = \left( {\begin{array}{*{20}{c}} {\frac{{t_2^i - t_{12}^i}}{{t_2^i - t_1^i}}}&{\frac{{t_{12}^i - t_2^i}}{{t_2^i - t_1^i}}}\\ {\frac{{ - 1}}{{t_2^i - t_1^i}}}&{\frac{1}{{t_2^i - t_1^i}}} \end{array}} \right)$。

把式（3）代入式（1），经参数变换，式（2）可转换为在新参数$ \hat X_i^{}$下的法方程，即：

$ \mathit{\boldsymbol{C}}_i^{\rm{T}}{{\mathit{\boldsymbol{\tilde N}}}_i}{\mathit{\boldsymbol{C}}_i}{{\mathit{\boldsymbol{\hat X}}}_i} = \mathit{\boldsymbol{C}}_i^{\rm{T}}{{\mathit{\boldsymbol{\tilde U}}}_i} $

(4)

式中，$ \mathit{\boldsymbol{\hat X}}_i^{} = \left( {\begin{array}{*{20}{l}} {\hat x_1^i}\\ {\hat x_2^i} \end{array}} \right)$。

观测时段i以新参数X_i为估值的法方程可由式（4）改写：

$ \left( {\begin{array}{*{20}{c}} {N_{11}^i}&{N_{12}^i}\\ {N_{21}^i}&{N_{22}^i} \end{array}} \right)\left( {\begin{array}{*{20}{c}} {\hat x_1^i}\\ {\hat x_2^i} \end{array}} \right) = \left( {\begin{array}{*{20}{c}} {U_1^i}\\ {U_2^i} \end{array}} \right) $

(5)

式中，${\mathit{\boldsymbol{N}}_i} = \left( {\begin{array}{*{20}{l}} {N_{11}^i}&{N_{12}^i}\\ {N_{21}^i}&{N_{22}^i} \end{array}} \right) = \mathit{\boldsymbol{C}}_i^{\rm{T}}{{\mathit{\boldsymbol{\tilde N}}}_i}{\mathit{\boldsymbol{C}}_i};{\mathit{\boldsymbol{U}}_i} = \left( {\begin{array}{*{20}{c}} {U_1^i}\\ {U_2^i} \end{array}} \right) = \mathit{\boldsymbol{C}}_i^{\rm{T}}{{\mathit{\boldsymbol{\tilde U}}}_i}$。

多个观测时段的法方程可进行叠加，则连续m个观测时段的总法方程为：

$ \begin{array}{*{20}{l}} {\underbrace {\left( {\begin{array}{*{20}{c}} {N_{11}^1}&{N_{12}^1}&0&0&0&0&0\\ {N_{21}^1}&{N_{22}^1}&0&0&0&0&0\\ 0&0&{N_{11}^2}&{N_{12}^2}&0&0&0\\ 0&0&{N_{21}^2}&{N_{22}^2}&0&0&0\\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 0&0&0&0&0&{N_{11}^m}&{N_{12}^m}\\ 0&0&0&0&0&{N_{21}^m}&{N_{22}^m} \end{array}} \right)}_\mathit{\boldsymbol{N}} \cdot \underbrace {\left( {\begin{array}{*{20}{c}} {\hat x_1^1}\\ {\hat x_2^1}\\ {\hat x_1^2}\\ {\hat x_2^2}\\ \vdots \\ {\hat x_1^m}\\ {\hat x_2^m} \end{array}} \right)}_{\mathit{\boldsymbol{\hat X}}} = }\\ {\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\underbrace {\left( {\begin{array}{*{20}{c}} {U_1^1}\\ {U_2^1}\\ {U_1^2}\\ {U_2^2}\\ \vdots \\ {U_1^m}\\ {U_2^m} \end{array}} \right)}_\mathit{\boldsymbol{U}}} \end{array} $

(6)

在进行多时段ERP融合时，相邻时段i、i+1的相同边界点处的ERP值应当一样，即有$ \hat x_2^i = \hat x_1^{i + 1}$，则边界约束条件^[12]为：

$ \begin{array}{l} \underbrace {\left( {\begin{array}{*{20}{c}} 0&1&{ - 1}&0&0& \cdots &0&0&0\\ 0&0&0&1&{ - 1}& \cdots &0&0&0\\ 0&0&0&0&0& \cdots &0&0&0\\ 0&0&0&0&0& \cdots &1&{ - 1}&0 \end{array}} \right)}_\mathit{\boldsymbol{B}}\underbrace {\left( {\begin{array}{*{20}{c}} {\hat x_1^1}\\ {\hat x_2^1}\\ {\hat x_1^2}\\ {\hat x_2^2}\\ \vdots \\ {\hat x_1^m}\\ {\hat x_2^m} \end{array}} \right)}_{\mathit{\boldsymbol{\hat X}}} = \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\underbrace {\left( {\begin{array}{*{20}{c}} 0\\ 0\\ \vdots \\ 0 \end{array}} \right)}_\mathit{\boldsymbol{W}} \end{array} $

(7)

联合法方程式（6）和约束条件式（7），可采用附加约束条件的参数平差进行求解，其解为：

$ \mathit{\boldsymbol{\hat X}} = {\mathit{\boldsymbol{N}}^{ - 1}}\left( {\mathit{\boldsymbol{B}}_{}^{\rm{T}}\mathit{\boldsymbol{K}} - \mathit{\boldsymbol{U}}} \right) $

(8)

式中，$ \mathit{\boldsymbol{K}} = (\mathit{\boldsymbol{BN}}_{}^{ - 1}\mathit{\boldsymbol{B}}_{}^{\rm{T}})_{}^{ - 1}\left( {\mathit{\boldsymbol{W}} - \mathit{\boldsymbol{BN}}_{}^{ - 1}\mathit{\boldsymbol{U}}} \right)$。

联合法方程式（6）和约束条件式（7），约化掉相同参数，其解为：

$ \left( {\begin{array}{*{20}{c}} {\hat x_1^1}\\ {\hat x_2^1}\\ {\hat x_2^2}\\ {\hat x_2^3}\\ \vdots \\ {\hat x_2^m} \end{array}} \right) = {\left( {\begin{array}{*{20}{c}} {N_{11}^1}&{N_{12}^1}&0&0&0&0\\ {N_{21}^1}&{N_{22}^1 + N_{11}^2}&{N_{12}^2}&0&0&0\\ 0&{N_{21}^2}&{N_{22}^2 + N_{11}^3}&{N_{12}^3}& \vdots & \vdots \\ \vdots & \vdots & \vdots & \vdots &{N_{12}^{m - 1}}&0\\ 0&0&0&{N_{21}^{m - 1}}&{N_{22}^{m - 1} + N_{11}^m}&{N_{12}^m}\\ 0&0&0&0&{N_{21}^m}&{N_{22}^m} \end{array}} \right)^{ - 1}}\left( {\begin{array}{*{20}{c}} {U_1^1}\\ {U_2^1 + U_1^2}\\ {U_2^2 + U_1^3}\\ \vdots \\ {U_2^{m - 2} + U_1^{m - 1}}\\ {U_2^{m - 1} + U_1^m}\\ {U_2^m} \end{array}} \right) $

(9)

通过式（8）、式（9）可得到一个长期解。理论上，附加约束条件式（7）后，式（8）与式（9）的解一致^[16]。

2 施加转换参数内约束的多中心解融合模型

通过式（8）或式（9）得到的长期解的精度应当均匀，但结果仍有差异。可认为不同技术、不同分析中心的ERP间的差异主要来源于参考系统间的定义差，不同系统下的ERP可相互转换，其转换模型为：

$ {\left( {\begin{array}{*{20}{l}} {x_p^i}\\ {y_p^i} \end{array}} \right)_A} = \left( {\begin{array}{*{20}{l}} {{\rm{\Delta }}x}\\ {{\rm{\Delta }}y} \end{array}} \right) + \left( {1 + m} \right)\left( {\begin{array}{*{20}{c}} 1&{ - \alpha }\\ \alpha &1 \end{array}} \right){\left( {\begin{array}{*{20}{l}} {x_p^i}\\ {y_p^i} \end{array}} \right)_C} $

(10)

式中，${x_{A, i}} = {\left( {\begin{array}{*{20}{l}} {x_p^i}\\ {y_p^i} \end{array}} \right)_A} $为转换前的ERP；$ {x_{C, i}} = {\left( {\begin{array}{*{20}{l}} {x_p^i}\\ {y_p^i} \end{array}} \right)_C}$为转换后的ERP；Δx、Δy、m、a为转换参数。

由于ERP参数不可避免存在误差，则式（10）可改写为：

$ {x_{A, i}} = {x_{C, i}} + E_i^{} \cdot \theta + \varepsilon $

(11)

式中，$ E_i^{} = \left( {\begin{array}{*{20}{c}} 1&0&{y_p^{}}&{x_p^{}}\\ 0&1&{ - x_p^{}}&{y_p^{}} \end{array}} \right)$；ε为误差；$\theta = {\left( {\begin{array}{*{20}{c}} {{\rm{\Delta }}x_p^{}}&{{\rm{\Delta }}y_p^{}}&\alpha &m \end{array}} \right)^{\rm{T}}} $。

ERP序列的转换方程为：

$ {X_A} = {X_C} + E_A^{} \cdot \theta + \mathit{\varepsilon } $

(12)

通过式（8）或式(9)可知，不同分析中心（或不同技术）结果都可得到一个长期解X_Ai，则r个分析中心的长期解可建立r个观测方程，即：

$ \left\{ {\begin{array}{*{20}{c}} {{X_{A1}} = {X_C} + {E_{A1}} \cdot {\theta _1} + {\varepsilon _1}}\\ {{X_{A2}} = {X_C} + {E_{A2}} \cdot {\theta _2} + {\varepsilon _2}}\\ \vdots \\ {{X_{Ar}} = {X_C} + {E_{Ar}} \cdot {\theta _r} + {\varepsilon _r}} \end{array}} \right. $

(13)

显然，式（13）是秩亏的。要得到唯一解X_c，必须施加一定基准，施加的基准不同，则得到的结果不同。为了保留ERP的固有物理参数信息（原点、长度与旋转）而不引入任何的外部物理信息，使用极移参数的内在约束条件^[18]，可对转换参数θ_i构造如下约束：

$ \sum\limits_{i = 1}^r {{\theta _i}} = 0 $

(14)

3 算例与分析

采用GNSS 7个分析中心（co2、em2、es2、gt2、jp2、mi2、ul2）的第二次精密处理结果为例，时间跨度从2005—2011年共6 a，精密处理结果文件格式为SINEX格式。由于IERS C04结果是目前公认的最好结果，本文将其作为参考值来分析比较本文计算结果。

图 1、图 2给出了7个分析中心结果与IERS C04结果的比较。可以看出，各个分析中心估计的ERP结果与IERS C04发布的相应历元ERP结果有一定差异，其中co2、em2、es2、gt2、jp2、mi2分析中心估计的结果从整体上看差别不大，但ul2与其他分析中心的结果间存在明显差异，在极移X、Y方向其最大值分别为1.23 mas、1.25 mas，最小值分别为-0.97 mas、-1.08 mas，而其他分析中心的结果最小值和最大值均在0.3 mas以内。这种差异产生可能主要是由于不同分析中心采用了不同的数据处理软件、不同的数据处理策略造成的。

为了验证本文提出方法的效果，按如下方案进行计算：

方案一：各分析中心的结果按其精度进行加权平均；

方案二：采用本文所述的边界约束和转换参数内约束的融合模型进行处理。

本文中也给出IGS（International GNSS Service）融合结果情况。如图 3给出了上述数据的融合结果与IERS C04结果的极移X方向的差异，图 4给出了上述数据的融合结果与IERS C04结果的极移Y分量的差异。表 1给出了上述两种方案以及IGS融合结果的统计情况。表 1中XPO表示极移X分量，YPO表示极移Y分量。方案二在获得多分析中心融合解的同时，得到了不同分析中心结果与融合解的系统转换参数。表 2给出了7个分析中心结果相对于融合解的转换参数的最佳估计值。

对上述结果进行分析可知：

1）从表 1和图 3、图 4可以看出，将多分析中心结果施加边界约束和内约束的融合模型计算的极移参数（方案二）与IERS C04结果进行比较，方案二极移参数的标准差在X、Y方向分别为0.052 6 mas、0.050 7 mas，其平均偏差分别为-0.000 9 mas、0.001 1 mas，其最大值分别为0.155 9 mas、0.193 5 mas。而多分析中心结果加权平均计算的极移参数（方案一）的标准差，在X、Y方向分别为0.057 6 mas、0.071 1 mas，其平均偏差分别为0.000 9 mas、0.001 1 mas，其最大值分别为0.231 4 mas、0.260 9 mas。

2）从表 1中最大值、平均偏差、标准差看，方案二的极移参数精度都较方案一有明显提高。极移X方向标准差从0.057 6 mas提高到0.052 6 mas，提高了8.6%，而极移Y方向的标准差从0.071 1 mas提高到0.050 7 mas，提高了28.6%。这主要是因为方案二改进了融合函数模型，对各分析中心的长期解在融合时增加了系统参数，并且对该系统参数采用内约束方式，并没有受到外部基准的干扰；而方案一仅对各分析中心的长期解结果加权平均，权的大小依据各分析中心的长期解的内部精度确定，但各中心解的协方差阵通常并未采用同一方差因子。

3）从表 1和图 3、图 4也可看出，边界约束和转换参数内约束的融合模型（方案二）确定极移参数精度较IGS的极移参数精度略有提高。极移参数X方向标准差从0.054 8 mas提高到0.052 6 mas，Y方向的标准差从0.051 3 mas提高到0.050 7 mas，但提高不明显。IGS组合结果以9个分析中心的结果为输入，融合模型采用的是与站坐标一起进行组合的相对严密的融合模型^[19]，在函数模型中考虑了七参数转换参数，采用最小二乘加权平均得到了综合解^[20]。

4）从表 2可以看出，ul2分析中心结果较其他分析中心结果存在较大的系统偏差。ul2分析中心结果相对于综合解在极移X方向存在0.063 9 mas的平移，在极移Y方向存在0.025 6 mas的平移，0.111 6 arcs的旋转和-0.123 0×10^-6的尺度差异，而其他分析中心结果的系统偏差较小，较ul2分析中心结果小一个数量级。这种差异的产生可能是由于该分析中心采用了不同的数据处理软件和不同的数据处理策略，也可能是由于数据中采用了不同物理参数。尺度参数是由数据处理中采用的物理参数和相对论改正模型来确定的，而尺度参数和平移参数具有较强的相关性^[21]。

4 结语

ERP是联系天球参考框架与地球参考框架^[22]的最重要参数，但不同技术分析中心确定的ERP往往不同，为方便用户使用、确保计算参数的可靠，常对不同分析中心结果进行综合^[12]。本文提出了一种基于多分析中心ERP结果的分布融合模型，第一步通过参数变换把各分析中心结果转换到相同时刻得到单中心结果的长期解，第二步对各中心结果的长期解增加系统参数，对系统参数施加内约束得到融合解。利用6 a的7个GNSS分析中心结果进行融合得到了0.05 mas极移精度，验证了该方法的有效性。我国对ERP的独立自主确定具有迫切的应用需求，本文提出的融合模型为提高国内ERP的融合提供了一种新的思路，是学习国际ERP融合模型的结果，但如何建立可靠高精度的ERP融合模型仍需要进一步深入研究。