2. 同济大学 空间信息与可持续发展应用中心，上海 200092；
3. 中国科学院 上海天文台，上海 200030；
4. 中国地震局 地壳运动监测工程研究中心，北京，100036
2. Center for Spatial Information Science and Sustainable Development,Tongji University,Shanghai 200092,China;
3. Shanghai Astronomical Observatory,Chinese Academic of Sciences,Shanghai 200030,China;
4. National Earthquake Infrastructure Service,China Earthquake Administration,Beijing 100036,China
1 引 言
GNSS与卫星激光测距(satellite laser ranging，SLR)和甚长基线干涉测量(very long baseline interferometry，VLBI)并址站空间归心基线是形成国际地球参考框架(international terrestrial reference frame，ITRF)相关产品的基础[1, 2]。因此，全球大部分GNSS/SLR/VLBI并址站都进行了归心测量[3, 4, 5, 6, 7, 8, 9, 10]。我国共有7个并址站，分别位于上海、武汉、长春、北京、西安、昆明和乌鲁木齐。上海并址站自1995年起进行过多次归心测量[3, 5, 11, 17, 18]，法国国家地理所(Institute Geographique National，IGN)在2003年也对上海和武汉并址站进行了归心测量[11, 12]，求得的空间归心基线已用于ITRF框架计算。
并址站的空间归心基线是指在ITRF框架下GNSS天线几何中心、SLR望远镜和VLBI天线旋转中心之间的基线向量。SLR与VLBI的旋转中心是指SLR望远镜与VLBI天线跟踪目标转动的垂直轴与水平轴的交点。由于GNSS天线中心比较容易测定，SLR和VLBI的旋转中心不能直接观测，因此通过固定在SLR或VLBI上的标志间接测定其旋转中心是空间归心基线测定的关键技术。标志点坐标可以利用精密三维控制网通过极坐标测量或前方交会求得，相关技术比较成熟[5, 17, 18, 20]。目前，根据标志点坐标确定旋转中心的方法可分两类：① 将标志点的ITRF框架三维坐标变换到测站的站心坐标系，通过拟合求出绕垂直轴旋转圆心和绕水平轴旋转圆心的站心系三维坐标，各个标志的垂直轴旋转圆心的平面坐标均值和水平轴旋转圆心的高程均值就是旋转中心的站心坐标，再进行逆变换得到旋转中心的ITRF框架三维坐标[5, 11, 12, 17, 18, 20]。② 采用第一类约束条件分别求出垂直轴旋转圆心和水平轴旋转圆心的ITRF框架三维坐标，然后将各标志点的垂直轴与水平轴旋转圆心坐标变换到站心系，并分别取其水平坐标均值与高程的均值，构成旋转中心的三维站心坐标，再变换到ITRF框架。在第一类约束条件的基础之上，可以附加额外的约束条件，如两轴平行、两轴偏差条件等可以进一步提高旋转中心的解算精度[21, 22, 23, 24, 25]。
本文除了用第一类约束条件建立观测方程外，还通过引入第二类约束条件建立标志点的垂直轴与水平轴旋转圆心与旋转中心的关系作为条件方程，直接解算旋转中心的ITRF框架三维坐标及协方差矩阵，不需要引入站心坐标系。该解法理论严密、数学模型简单，方便进行程序设计。经计算与已有的分步解法差值不超过1 mm。
2 施测方案精密三维控制网的建立：VLBI站周围至少布设4个控制点，SLR站附近至少布设两个控制点。具备GNSS观测条件的各点均采用Trimble双频GNSS接收机观测两个时段，每个时段至少观测12 h。每个控制点观测所有通视点的水平方向、垂直角和距离，均采用TCA2003全站仪观测4个测回。上海佘山站VLBI和SLR站周围的控制点分别如图 1和图 2所示。
|
| 图 1 VLBI观测控制网图 Fig. 1 Co-location diagram about VLBI |
|
| 图 2 SLR观测控制点示意图 Fig. 2 Co-location diagram about SLR |
VLBI 几何旋转中心测量原理如图 3所示，VLBI天线的水平轴和竖直轴的交点就是VLBI天线的几何旋转中心。具体在VLBI天线的一侧安装观测标志，每个标志外侧一端带有红颜色小球，以便于方向观测(如图 4所示)。VLBI天线绕竖直轴每旋转15°，3台全站仪同时测定标志点水平方向和垂直角。当天线旋转1周时，测定出标志位于同一水平面内的一组离散点轨迹。同理，当VLBI天线绕水平轴旋转，测定出标志位于同一垂直面内的一组离散点轨迹，由于VLBI天线绕水平轴只能旋转180°，因此观测点间隔为10°。
|
| 图 3 几何示意图 Fig. 3 Geometric center diagram |
|
| 图 4 VLBI上安置观测标志 Fig. 4 Observation mark on the VLBI |
由于SLR望远镜上可安装棱镜或反射片(如图 5所示),因此用全站仪进行极坐标测量，即直接测定水平方向，垂直角和距离，其绕垂直轴与水平轴的旋转方法与VLBI相同。
|
| 图 5 SLR上放置棱镜 Fig. 5 Prism on the SLR |
原有的解法(分步解法)是将标志点ITRF框架三维坐标变换到测站的站心坐标系，拟合求出绕垂直轴旋转圆心和绕水平轴旋转圆心的站心系三维坐标，最后将标志点绕垂直轴旋转圆心的平面坐标均值和水平轴旋转圆心的高程均值组合成旋转中心的站心坐标，再逆变换得到ITRF框架的三维坐标，具体可详见文献[18]。本文提出的直接解法将平差后的标志点三维坐标作为相关观测量，根据标志点绕垂直轴或水平轴旋转形成的平面圆约束为第一类约束条件，平面圆心与旋转中心间的几何约束为第二类约束条件，直接求解旋转中心的三维坐标及其协方差阵。由于控制点布设在SLR和VLBI天线周围，观测站离SLR望远镜不到20 m，离VLBI天线不到100 m，垂线偏差的影响可以忽略。
3.1 第一类约束条件当固定在VLBI/SLR上的观测标志绕同一旋转轴旋转时，按一定间隔观测其三维坐标，各观测点可以表示为由一个平面与一个球面相割而形成的平面圆。因此，平面圆上的每个观测点同时满足球面和平面方程，且球心位于平面上(如图 6和图 7所示)。绕垂直轴旋转时，其球面与平面约束方程可表示为
式中，(aH,bH,cH)为平面法向量系数；(XH,YH,ZH)为旋转中心坐标；rH为球半径。同理，绕水平轴旋转时，其球面与平面约束方程为 式中，(aV,bV,cV)为平面法向量系数；(XV,YV,ZV)为旋转中心坐标；rV为球半径。若(Xi,Yi,Zi)和(Xj,Yj,Zj)分别为i和j点改正后的三维空间坐标，即 式中，(Xi,Yi,Zi)和(Xj,Yj,Zj)分别为i、j点的三维空间坐标观测值；(vXi,vYi,vZi)和(vXj,vYj,vZj)为其改正数向量。将式(3)、式(4)代入式(1)和式(2)进行线性化。若观测标志绕垂直轴与水平轴旋转时分别测定m与n个观测点，绕垂直轴与水平轴分别观测了g圈和h圈，则线性化后的观测方程可表示为
|
| 图 6 绕垂直轴旋转 Fig. 6 Rotation around the vertical axis |
|
| 图 7 绕水平轴旋转 Fig. 7 Rotation around the horizontal axis |
根据式(1)和式(2)，每旋转一圈绕垂直轴和水平轴分别需要引入7个未知参数，要求出这14个未知参数，至少分别需要测定3个标志点的三维坐标。因此x1为7(g+h)维参数向量；A为2(mg+nh)×7(g+h)阶列满秩矩阵；v为3(mg+nh)维改正数向量；B为2(mg+nh)×3(mg+nh)矩阵；y为2(mg+nh)维不符值向量。根据误差传播定律，y的协方差阵为
式中，Σ为观测点坐标的协方差阵，其维数和排列顺序与v相同，由三维控制网平差求得。 3.2 第二类约束条件
当天线绕垂直轴旋转时，其圆心与几何旋转中心位于同一铅垂线上；当天线绕水平轴旋转时，其圆心与几何旋转中心位于同一水平面上，构成的约束条件称之为第二类约束条件。
以VLBI天线或SLR望远镜旋转中心O(0,0,0)为原点，铅垂线方向为U轴，水平面北向为N轴建立局部坐标系。设O点空间直角坐标和经纬度分别为(X0,Y0,Z0)、(B0,L0)。在局部坐标系中绕垂直轴旋转圆心的坐标(NH,EH,UH)和绕水平轴旋转圆心的坐标(NV,EV,UV)分别可表示为
由于天线旋转中心与垂直轴圆心位于同一铅垂线上，与水平轴旋转圆心位于同一水平面上，因此NH=0，EH=0，UV=0，根据式(7)和式(8)整理后得
显然，分别绕垂直轴和水平轴观测1圈，可分别形成条件方程式(9)和式(10)。当分别观测g圈和h圈时，可形成2g+h个条件方程，其一般式可表示为
式中，参数x1为7(g+h)维向量；C1T为(2g+h)×7(g+h)阶矩阵，由于式(9)只包含x1中的旋转圆心坐标参数，因此C1T中的系数包含大量0；x2=(X0,Y0,Z0)T为旋转中心坐标参数；C2T为(2g+h)×3阶矩阵；w为(2g+h)维常数项向量。综合式(4)和式(11)组成 用附有约束的最小二乘准则求得公式(12)的平差值为 则x1和x2协方差矩阵为 4 成果分析
根据中国长春和上海两个并址站的GNSS与SLR或VLBI的空间归心测量数据，先利用GNSS和常规观测值三维联合平差解算出各个标志点的三维坐标，然后用本文的直接解法进行了处理。表 1给出了解算的空间归心基线及其精度；表 2给出了与ITRF2008给出的归心基线的比较结果；表 3给出了与分步解法的比较结果。其中，SHAO、CHAN分别为上海佘山和长春GNSS站的IGS(international GNSS service)站名，SLR和VLBI分别表示并址的SLR和VLBI站。
| 并址站 | 基线 | ΔX/m | ΔY/m | ΔZ/m | MΔX/mm | MΔY/mm | MΔZ/mm |
| 长春 | CHAN-SLR | 40.300 0 | 46.015 5 | -13.340 5 | 0.2 | 0.2 | 0.2 |
| 上海 | SHAO-SLR | 989.058 2 | 914.354 3 | -296.572 8 | 1.4 | 0.6 | 0.5 |
| 上海 | SHAO-VLBI | 46.346 2 | 67.642 2 | -41.815 4 | 0.8 | 0.8 | 0.9 |
| 并址站 | 基线 | ITRF框架下成果 | ITRF与直接解法较差 | ||||
| ΔX/m | ΔY/m | ΔZ/m | ΔX/mm | ΔY/mm | ΔZ/mm | ||
| 长春 | CHAN-SLR | 40.298 1 | 46.009 1 | -13.352 7 | 2.0 | 6.4 | 12.1 |
| 上海 | SHAO-SLR | 989.056 6 | 914.350 9 | -296.573 | 1.6 | -3.4 | 0.2 |
| 上海 | SHAO-VLBI | 46.351 9 | 67.638 0 | -41.816 5 | -5.7 | 4.3 | 1.1 |
| 并址站 | 基线 | ΔX/m | ΔY/m | ΔZ/m | 备注 |
| 长春 | CHAN-SLR | 40.300 0 | 46.015 5 | -13.340 5 | 直接解法 |
| CHAN-SLR | 40.299 3 | 46.016 1 | -13.339 6 | 分步解法 | |
| 差值/mm | 0.7 | -0.6 | -0.9 | ||
| 上海 | SHAO-SLR | 989.058 2 | 914.354 3 | -296.572 8 | 直接解法 |
| SHAO-SLR | 989.058 0 | 914.354 9 | -296.572 4 | 分步解法 | |
| 差值/mm | 0.2 | -0.6 | 0.4 | ||
| 上海 | SHAO-VLBI | 46.346 2 | 67.642 2 | -41.815 4 | 直接解法 |
| SHAO-VLBI | 46.346 0 | 67.642 8 | -41.815 3 | 分步解法 | |
| 差值/mm | 0.2 | -0.6 | -0.1 |
从以上给出的成果表中可以看出：
(1) 归心基线的内符合精度很高，其三维坐标的中误差最大只有1.4 mm。
(2) 分步解法和直接解法的成果相比，归心基线差值不超过1 mm，因此两种方法计算成果基本一致。
(3) 直接解法与ITRF框架下的成果值相比，存在不同程度的差值。其中，上海SLR站归心基线差值在Y轴方向为3.4 mm；上海VLBI 站差值在X轴方向为5.7 mm；长春SLR站归心基线差值更大，Z轴方向差值达到12.1 mm。
5 结 论在本次中国并址站空间归心测量中，从外业控制网的布设、观测到解算等各个环节都进行了严格的质量控制，为最终得到正确、高精度的归心结果提供了保证。
(1) 分步解法是在站心坐标系中拟合求出绕垂直轴旋转圆心平面坐标值和绕水平轴旋转圆心高程值，再逆变换得到ITRF框架的三维坐标，需要分步进行计算。
(2) 本文通过引入两类约束条件建立标志点与旋转中心之间的关系，构建观测方程和条件方程，不需要引入站心坐标系，按附有约束条件的间接平差解算旋转中心的ITRF框架三维坐标及协方差矩阵。解法理论严密、数学模型简单，方便进行程序设计。
(3) 利用上海佘山的GNSS/VLBI/SLR和长春GNSS/SLR的实测数据，计算了归心基线及其精度，该计算结果与分步解法相差不超过1 mm。与ITRF2008的基线结果相比，在长春站相差最大，达到12.1 mm，而上海SLR与VLBI站差值最大为5.7 mm。这些计算成果具有一定的参考价值。
| [1] | ALTAMIMI Z, SILLARD P, BOUCHER C. ITRF2000: A New Release of the International Terrestrial Reference Frame for Earth Science Applications[J]. Journal of Geophysical Research, 2002,107(B10): 2214. |
| [2] | RAY J, ALTAMIMI Z. Evaluation of Co-location Ties Relating the VLBI and GPS Reference Frames[J]. Journal of Geodesy, 2005, 79(4-5): 189-195. |
| [3] | XIONG Yongqing, ZHU Wenyao. The Last Measurement Results and Analysis among Shanghai Observatory VLBI, SLR and GPS Stations[J].Annals of Shanghai Observatory Academia Sinica, 1996,17: 67-72. (熊永清,朱文耀. 上海天文台VLBI,SLR和GPS站的归心结果与分析[J]. 中国科学院上海天文台年刊, 1996,17: 67-72.) |
| [4] | NOTHNAGEL A, WIRTZ C, SAUERBIER M, et al. Local Survey at the Effelsberg Radio Telescope 1997: Preliminary Results[C]//Proceedings of the 13th Working Meeting on European VLBI for Geodesy and Astrometry. Viechtach: EVGA, 1997: 25-31. |
| [5] | WANG Jiexian, JI Shanbiao, SHI Yiming. The Determination of Space Coordinate Differences among VLBI, SLR and GPS Stations of Shanghai Observatory[J]. Journal of the PLA Institute of Surveying and Mapping, 1997, 14(1): 7-10. (王解先,季善标,施一民. 上海天文台的VLBI、SLR、GPS站的空间归心测量[J]. 解放军测绘学院学报, 1997, 14(1): 7-10. |
| [6] | SARTI P, VITTUARI L, TOMASI P. GPS and Classical Survey of the VLBI Antenna in Medicina: Invariant Point Determination[C]// Proceedings of 14th Working Meeting on European VLBI for Geodesy and Astrometry. Bologna: EVGA, 2000: 67-72. |
| [7] | VITTUARI L, SARTI P, TOMASI P. 2001 GPS and Classical Survey at Medicina Observatory: Local Tie and VLBI Antenna Reference Point Determination[C]// Proceedings of the 15th Workshop Meeting on European VLBI for Geodesy and Astrometry. Barcelona: Institut dEstudis Espacials de Catalunya, 2001: 161-167. |
| [8] | JOHNSTON G, DAWSON J, TWILLEY B, et al. Accurate Survey Connections between Co-located Space Geodesy Techniques at Australian Fundamental Geodetic Observatories: Geodesy Technical Report 3[R]. Canberra: Australian Surveying and Land Information Group, 2000. |
| [9] | JOHNSTON G, DAWSON J. The 2003 Yarragadee (Moblas 5) Local Tie Survey[R]. Canberra: Geoscience Australia, 2003. |
| [10] | LONG J, MICHEL V, CORBIèRE C, et al. Hartebeesthoek Co-location Survey Reports and Results[R]. Paris: Institut Géographique National, 2005. |
| [11] | GARAYT B, KALOUSTIAN S, LONG J, et al. Sheshan Co-location Survey Report[R]. Paris: Institut Géographique National, 2005. |
| [12] | MICHEL V, KALOUSTIAN S, VERGEZ P. Wuhan Co-location Survey Report[R]. Paris: Institut Géographique National, 2005. |
| [13] | LEINEN S, BECKER M, DOW J, et al. Geodetic Determination of Radio Telescope Antenna Reference Point and Rotation Axis Parameters[J]. Journal of Surveying Engineering, 2007, 133(2): 41-51. |
| [14] | LONG J, CARPENTER T. Goddard Geophysical and Astronomic Observatory Co-location Survey Report[R]. Washington: National Aeronautics and Space Administration, 2007. |
| [15] | LANOTTE R, PIRRI M, BIANCO G. Matera Site Survey and VLBI Invariant Point Determination[C]// The 5th IVS General Meeting Proceedings. Saint Petersburg: IVS, 2008: 87-92. |
| [16] | POYARD J C. Grasse ITRF Co-location Survey Report[R]. Paris: Institut Géographique National, 2009. |
| [17] | LI Jinling, QIAO Shubo, LIU Peng, et al. Site Survey at Sheshan 25 m Radio Telescope in 2008[J]. Geomatics and Information Science of Wuhan University, 2010, 35(12): 1387-1391. (李金岭,乔书波,刘鹏,等. 2008年佘山25 m射电天线归心测量[J]. 武汉大学学报: 信息科学版, 2010, 35(12): 1387-1391.) |
| [18] | SHEN Yunzhong, CHEN Tingwu. Determination of Space Coordinate Differences of Co-location Sites in Shanghai Observatory[J]. Journal of Tongji University: Natural Science, 2006, 34(2): 217-222. (沈云中,陈廷武. 上海天文台并址站的空间归心测量[J]. 同济大学学报: 自然科学版, 2006, 34(2): 217-222.) |
| [19] | DAWSON J, SARTI P, JOHNSTON G, et al. Indirect Approach to Invariant Point Determination for SLR and VLBI Systems: An Assessment[J]. Journal of Geodesy, 2007, 81(6-8): 433-441. |
| [20] | WOODS A R. Indirect Determination of the Invariant Reference Point (IVP)of SLR and VLBI Observing Systems[C]//Proceedings of 13th FIG Symposium on Deformation Measurement and Analysis & 4th IAG Symposium on Geodesy for Geodesy for Geotechnical and Structural Engineering. Lisbon: [s. n.], 2008: 42-57. |
| [21] | LSLER M, HENNER M. An Innovative Mathematical Solution for a Time-efficient IVS Reference Point Determination[C]//Proceedings of 13th FIG Symposium on Deformation Measurement and Analysis & 4th IAG Symposium on Geodesy for Geotechnical and Structural Engineering. Lisbon: [s. n.], 2008. |
| [22] | LSLER M. New Mathematical Model for Reference Point Determination of an Azimuth-elevation Type Radio Telescope[J]. Journal of Surveying Engineering, 2009, 135(4): 131-135. |
| [23] | SARTI P, SILLARD P, VITTUARI L. Surveying Co-located Space Geodetic Instruments for ITRF Computation[J]. Journal of Geodesy, 2004, 78(3): 210-222. |