1 引 言

天文定位是以太阳、月球、行星和恒星等自然天体作为信标，以天体的地平坐标(方位或高度)作为观测量，进而确定测量点地理位置的技术和方法^{[1, 2]}。因其是一种绝对定位手段，天文点在大地控制网中常用于控制测量误差的累积^{[3, 4, 5, 6]}，又因其具有隐蔽性好、不易受干扰和独立性强等优点，所以可作为一种独立的自主导航方式^{[7, 8, 9, 10]}。

若仪器已经精密整平，则利用测量仪器可以直接观测得到自然天体的地平方位角和高度角，由于自然天体总是按其固有规律运动，某个时刻它们在赤道天球坐标系中的位置矢量是可以精确计算得到的，通过对天体地平坐标观测值和在赤道天球坐标系中的理论位置进行处理，就可获测站的天文经纬度，可见天文定位需要精密的水平基准^[11]。传统的天文定位在观测之前都需要将仪器精密整平，如专门用于天文测量的T4经纬仪上带有5种水准器，需要在观测前进行一系列复杂的整平操作，这种手动整平方式不仅耗时长，而且不容易实现^[12]。因此，研究在不需精确整平的条件下进行天文定位的方法具有重要意义。

经过标校的倾角传感器可以直接测得其水平倾角值，若将两个倾角传感器相互垂直放置，则可以同时获得二维倾角值^{[13, 14]}。将二维倾角传感器与测量仪器固定在同一平台之上，就可以通过二维倾角传感器的读数变化来检测测量仪器的姿态变化。更进一步，对上述二维倾角传感器和测量仪器的相对位置进行标定，使传感器读数为零时测量仪器垂直轴恰好与当地铅垂线重合，则传感器的倾斜值即为测量仪器的倾斜值。在此基础上就可以根据传感器的读数对未经整平的天体位置观测值进行二维倾斜改正，得到地平坐标中的天体水平角和高度角。为了实现无需精密整平的天文测量，本文推导了根据二维倾角传感器数据改正天体位置观测量的公式；为了减小系统误差的影响，在天文定位模型中加入了天顶距改正量；为了限制测量粗差的影响，引入抗差估计的算法，建立了抗差天文定位模型；最后通过两个算例分析了本方法的可行性和精度。无需精密整平的天文定位流程如图 1所示。

2 观测值的二维倾斜改正

将二维倾角传感器和天文测量仪器放置于同一测量平台之上，并经过标校使两者的坐标轴相互平行，如图 2所示。

图 2中，O_iX_i、O_iY_i、O_iZ_i分别为天文测量仪器的3个坐标轴，构成右手系；O_gX_g、O_gY_g、O_gZ_g分别为二维倾角传感器的3个坐标轴，分别与测量仪器的三轴相互平行。设X_g轴和Y_g轴的倾斜读数分别为θ₁和θ₂，而测量平台与水平面的二面角为θ。当X_g轴和Y_g轴正向偏高时，对应的输出倾斜值为正，则如图 3所示。

图 3中，OX_g、OY_g、OZ_g分别为二维倾角传感器的3个相互垂直的坐标轴，O-X_gY_gZ_g构成测量坐标系；OZ_h为当地铅垂线，OX_h、OY_h分别为OX_g、OY_g轴在水平面上的投影；∠OX_hX_g=θ₁、∠OY_hY_g=θ₂ 。需要注意的是，OX_g轴与OY_g轴相互垂直，但是OX_h轴和OY_h轴不一定垂直，它们之间的夹角与测量平面倾角θ有关。OT_g在测量平面内，并且与测量平面和水平面的交线相互垂直。OT_h轴为OT_g轴在水平面的投影，它们之间的夹角即为测量平面的倾角θ。记∠X_gOT_g和∠X_hOY_h分别为φ_g和φ_h，则根据二面角的相关定理，有

根据式(1)可求得测量平面的倾角，根据式(2)可求得测量平面内倾斜方向与二维倾角传感器OX_g轴的夹角，根据式(3)可求得水平面上倾角方向与OX_h轴的夹角。上述公式中只推导了θ₁、θ₂均为正的情况，实际使用时需要根据θ₁、θ₂的正负判断倾斜方向。

假设测量仪器在倾斜状态下获得的观测量为水平角L_g和高度角h_g。为了改正测量平面倾斜的影响，下面探讨将观测值改正为真实水平角L_h和真实高度角h_h的方法，如图 4所示。

图 4中，Z_g为二维倾角传感器的垂直轴；Z_h为铅垂线；θ为测量平面的倾角；σ为某一观测目标；h_g为测量仪器倾斜时的观测高度角；h_h为真实的高度角。在测量平面上，以OT_g方向为X轴，以OR的反方向为Y轴，以OZ_g方向为Z轴构成右手直角坐标系。由于测量仪器的水平零轴与二维倾向传感器的OX_g轴重合，则根据式(2)，测量目标与X轴在测量平面内的夹角为L_g－φ_g，则将观测目标在测量平面内的直角坐标可转换为

以OR为轴，将测量坐标系按图中所指方向旋转θ角，从而测量平面与水平面重合，根据右手法则，可得转换关系式为 从而求得了相应的真实地平直角坐标，再化算为观测目标在地平坐标中的水平角L_h和高度角h_h，如式(6)所示

结合式(4)、式(5)、式(6)可得L_h和高度角h_h与测量平面倾角θ的关系式

式中，θ、φ_g、φ_h可根据式(1)、式(2)、式(3)求得。

3 抗差天文定位模型 3.1 天文定位模型

经过上述步骤计算得到了天体的水平角和高度角，结合观测时刻和天体的理论位置，即可解算测站的天文经纬度。在天球上，北天极、天顶和恒星可构成定位三角形，如图 5所示。

图 5中，P为北天极；Z为天顶；W为西点；WMQ为地平圈；WNQ′为天球赤道；σ为天球上的一颗恒星。由球面三角公式中的边余弦公式，可以得到天顶距z、测站天文纬度φ、经度λ、恒星的赤经α、赤纬δ与时角t之间的关系式^{[15, 16]}

式中，天顶距z由经过二维倾斜改正的天体高度角得到 t为天体的时角 式中，S为观测瞬间的格林尼治真恒星时，可以由观测瞬间的UTC时刻T化算得到。首先根据时间频率公报将UTC时刻T加上UT₁与UTC的差值化算为UT₁时刻T′，然后将其化为平恒星时，最后加上赤经章动Δu(Δu=Δψcosε，Δψ为黄经章动，可由公式计算得到，ε为黄赤交角)将其化为真恒星时S。恒星的赤经α，赤纬δ可由视位置计算程序从星表中加恒星自行改正得到，时刻T(或T′)由观测同时记录。

因此，式(8)中只含有未知数φ、λ，只需观测2颗恒星即可求得测站的天文经纬度。但由于天文测量中大气折射改正不彻底，求得的天顶距z含有系统误差，故引入一个微小量Δz，即观测的天顶距应为z+Δz，因此，至少需观测3颗恒星才能解算这3个参数^[15]。设Δφ为测站纬度φ与其初值φ₀之差，即

Δλ为测站经度λ与其初值λ₀之差，即
因此，方程中只需解3个参数Δφ、Δλ、Δz。若对n个天体进行观测，并以φ₀、λ₀为经纬度初始值，对式(8)进行线性化^[17]，可得到

式中，V为n维残差向量；A是n×t阶系数矩阵；是t维参数估计向量；L是n维观测值向量
式中

对式(13)进行解算即可得到测站的天文经度λ和天文纬度φ。

3.2 抗差M估计

式(13)的解法通常采用最小二乘平差，由于测量条件影响和天体识别错误等问题，部分观测量中不可避免的含有较大误差，通常也存在少数粗差，导致天文经纬度的最小二乘估值往往与真值偏差较大。本文引入抗差M估计的方法，充分利用高精度信息，限制利用低精度信息，排除粗差的影响，获得较为可靠、有效、且具有实际意义的天文经纬度估值。

抗差估计也称稳健估计，源于统计学中的稳健性概念，它是针对最小二乘估计不具有抗干扰性这一缺陷提出来的。它是一种既能减免粗差影响，或抗拒异常值的干扰，又具有较高效率的估计^{[18, 19, 20, 21, 22, 23]}。

天文测量时对各天体进行分别测量，因此观测值相互独立，先验权矩阵P为对角阵。为了消除或减弱粗差和较大偏差对估计值的影响，引入抗差M估计的原理，用等加权矩阵p代替先验权矩阵P，则可解得未知参数的抗差M估值为

参数的验后协方差阵为

上式的解法一般采用迭代法，第k+1次的迭代解为

抗差估计的关键是选择合适参数初值和有效的等价权函数。针对天文测量的特点，本文首先以具有强淘汰性质的L₁范数最小为准则求解参数初值^[24]，设初始等价权

式中，p_i为各观测量的先验权；v_i为对应的最小二乘残差。将该等价权函数代入式(19)即可求得基本不受粗差影响的参数初值，根据式(13)即可解得各观测量的初始残差，然后将残差标准化，即取 式中，v_i是观测值残差；m_vi是v_i的中误差，m_vi由下式计算 式中，σ₀为单位权中误差，可采用理论值或经验值；p_i是观测值的初始权；A_i是矩阵A的第i行，即第i个误差方程的系数向量^[25]。

根据式(23)解得的标准化残差，选择IGG3方案来确定观测值的等价权^{[26, 27]}

式中，k₀、k₁的取值一般为1.5和3.0。

将等价权函数代入式(21)迭代解算即可求得参数的最终抗差估值，利用式(11)、(12)即得测站经度λ和纬度φ。

4 算例与分析 4.1 算例1

为验证本文所提出的无需精密整平的抗差天文定位方法的有效性和精度，根据目前的倾角传感器精度和天文测量仪器测角精度模拟天文观测数据，并利用本文所提出的方法进行天文定位解算。

参照大地天文测量规范，制定一个时段内的观测步骤为^[28]：

(1) 概略整平测量平台。将测量平台稳定放置，使二维倾角传感器的读数在其量程范围内，并实时记录其读数。

(2) 生成观测星表。根据测站位置和观测时间，遴选待测天体，给出待测天体的概略水平角和高度角，要求所有待测天体高度角互差不超过±1°，各象限数量最大互差不超过观测天体总数的10%。

(3) 观测天体。利用电子经纬仪照准某一个天体，并记录其水平角和高度角观测值，观测的同时记录观测时刻，直至完成所有天体观测。

完成后，即可按照本文所提出方法进行天文定位解算。目前二维倾角传感器的双轴测角精度最高可达±0.2″(1σ)，常用的Leica NIVEL200系列双轴倾角传感器的量程为±220″，其测角中误差m_θ1和m_θ2均为±1″^[29]。天文测量常用的仪器为电子经纬仪，如Leica TCA 2003和Leica TM 5100A等，其测角中误差均为±0.5″，故设定以下仿真条件：

(1) 测站的真实天文坐标经纬度分别为东经120°和北纬30°，解算时所有观测值均归算到同一时刻，UTC时间2013年1月1日20时整。

(2) 二维倾角传感器与测量仪器固连在统一测量平台上，并且其相对位置已经经过精密标校，其坐标轴相互平行。观测时测量平台概略整平，二维倾角θ₁和θ₂分别为+150.0″和+200.0″。

(3) 二维倾角传感器的测角中误差m_θ1和m_θ2均为±1″，天文测量仪器的测角中误差±0.5″。

首先根据仿真条件(1)、(2)计算得到各颗恒星在测量平面内的理论水平角和高度角；然后按照仿真条件(3)加入相应的二维倾角测量误差和测量仪器测角误差，再对各观测值以5%的概率加入5″的粗差、以5%的概率加入-5″的粗差，以25%的概率加入1″的较大误差，以25%的概率加入-1″的较大误差；最后利用本文算法解算得到测站经纬度及其误差。根据以上条件，选取观测18颗恒星，其水平角和高度角如图 6所示。

从图 6中可以看出，恒星高度角均在60°左右，方位角分布比较均匀。以这18颗恒星为观测对象，按照上述条件和步骤进行100次重复试验，得到的经纬度估值与测站真值之差如图 7所示。

从图 7中可以看出，利用本文算法得到的经度和纬度误差基本都在±1″范围内，经度和纬度解算值的中误差分别为±0.46″和±0.42″，误差均值分别为-0.06″和-0.05″。说明本文提出的方法不仅能够正确地将恒星水平角和高度角观测值从测量平面改正到水平面，而且能够以较高的精度解算测站的经纬度。

4.2 算例2

为进一步验证本文所提出的抗差天文定位方法具有较好的定位精度，分别利用本文算法和传统算法分别对2011年11月16日在华北某地的实测数据进行解算分析。

该次试验一个时段内共观测了18颗恒星，对每颗恒星进行大约10次采样测量，其观测高度角误差随方位角的变化如图 8所示。

由于各颗恒星的方位角并不相同，致使每颗恒星的高度角误差在图 8中呈现为列状，图中所有高度角误差的均值为+5.64″，中误差为±2.45″，均值存在较大偏差主要是由于大气折射改正不彻底，导致改正恒星高度角时存在系统差。

从第3节的分析可知，本文在天文定位模型中加入了天顶距改正数Δz，可以有效地减弱系统差的影响；在参数解算时采用抗差估计的方法，能够有效地抑制粗差的影响，限制利用误差偏大的观测值信息。而传统算法首先利用最小二乘法解算得到参数初值和观测值残差，然后以2倍中误差为限剔除粗差，最后再利用最小二乘法求解参数估值。分别利用本文算法和传统算法对上述实验数据进行解算，其结果对比如表 1所示。

从表 1可以看出，无论是定位经纬度误差，还是其中误差，本文算法均优于传统算法，能够有效地抑制粗差的影响，减弱系统误差的影响，充分利用高精度的观测信息，得到较高精度的测站经纬度估值。

5 结 论

(1) 按照本文所提出的方法，能够实现无需精密整平的天文导航定位。天文测量时，只需将测量仪器概略整平，然后根据本文算法对天体位置观测值进行二维倾斜改正，最后进行定位解算。

(2) 本文算法的定位精度优于传统算法。本文算法引入了天顶距改正数Δz和抗差估计的方法，能够减弱系统误差的影响，抑制粗差的影响，在解算经纬度估值的误差和中误差两方面，本文算法均优于传统算法。

(3) 若利用星敏感器取代经纬仪进行天文导航定位，它能够同时获得大量的天体观测数据，但是不易实现精密整平，其中不可避免地包含误差偏大的观测值和粗差，采用本文算法将更具优势。