抗差估计在多星近似等高法中的应用.pdf

第 37 卷 第 4 期2023 年 4 月 北京测绘Beijing Surveying and MappingVol.37 No.4April 2023引文格式:栗志昊,詹银虎,王成龙.抗差估计在多星近似等高法中的应用J.北京测绘,2023,37(4):475-480.Reference format:LI Zhihao,ZHAN Yinhu,WANG Chenglong.The application of anti-difference estimation in the multi-star approximate contour methodJ.Beijing Surveying and Mapping,2023,37(4):475-480.DOI:10.19580/ki.1007-3000.2023.04.002收稿日期 2022-08-15基金项目 国家自然科学基金(42074013)作者简介 栗志昊(1999),男,河北石家庄人,本科在读,研究方向为天文导航。E-mail:531662918 通信作者 詹银虎,E-mail:oscardad 抗差估计在多星近似等高法中的应用栗志昊 詹银虎 王成龙(战略支援部队信息工程大学,河南 郑州 450001)摘 要 传统天文测量中,利用多星近似等高法测定天文经纬度时默认平差解算过程中的权矩阵为单位阵,在一定程度上影响了观测定位精度。为了综合考虑已知的误差项和尚未确定的影响因素,利用抗差估计调整观测数据的权值来降低异常观测量对定位结果的影响。为此,提出将 IGG抗差估计方案引入平差解算,根据前一次解算结果的中误差大小动态调整方案阈值系数,实现权矩阵结构的动态迭代更新,并推导出详细的解算公式。通过实地观测并结合仿真试验验证了理论可行性,能够在高精度全站仪观测条件下将天文定位精度提升 10%以上,并且证明抗差方案可以有效削弱不良选星空间结构对于解算结果的影响。关键词 天文定位测量;多星近似等高法;抗差估计中图分类号 P227 文献标识码 A 文章编号 1007-3000(2023)04-0475-060 引言光学天文大地测量技术是大地测量技术的重要组成部分,在大地水准面模型验证、高精度垂线偏差测量、地震等自然灾害预报、地球自转参数监测与预报以及深空探测等领域有着广泛的应用1-2。目前,配备图像全站仪的视频测量机器人已经实现了高精度、集成化,并且在向全自动化无人干预、智能化阶段推进3-4。而在天文定位算法方面,目前广泛采用的是利用多星近似等高法同时测定天文经纬度5,该理论将大气垂直折光模型误差作为系统参数同天文经纬度一同解算6,有效削弱了大气垂直折射模型误差对天文定位的影响并大大提高了观测效率。但是现阶段成熟的解算系统中,均认为观测数据质量差异较小,平差解算时默认权矩阵为单位阵,在一定程度上影响了解算精度。为了区分不同质量的天文观测数据对于解算结果的贡献度,本文将 IGG抗差估计方案7引入平差解算程序中,利用前一次迭代的解算结果中各时段中误差大小来确定下一次解算过程中 IGG权函数的阈值系数 k0和 k18-9,实现权矩阵结构的动态调整,兼顾算法效率的同时提升抗差能力。通过处理实地观测数据得到结论,与传统多星近似等高法相比,该算法一定程度上提升了解算精度,并且选星空间结构越差时算法应用效果越好。1 基于抗差估计的多星近似等高法1.1 多星近似等高法2根据球面几何中的边余弦公式10,可以得到天顶距 z(或高度角 h),待测点天文经纬度、,观测恒星的赤经、赤纬 和时角 t 之间的关系,公式为cos z=sin sin +cos cos cos t(1)z=2-h(2)北京测绘第37 卷 第4 期式中,时角 t 与观测恒星赤经、待测点天文经度 和观测瞬间格林尼治真恒星时 S 有关,具体关系式为t=S-+(3)观测天体的赤经、赤纬 可通过查询星表并进行恒星视位置计算得到,观测天体的天顶距z(或高度角 h)则通过仪器直接测量得到,观测瞬间格林尼治真恒星时 S 可以由观测瞬间的 UTC时刻加改正换算解得11。由于大气折射干扰对于天文观测结果影响较大,在正式解算前需要对高度数据进行大气折射改正。因此,式(1)只含有待测点天文经纬度、两个未知数,理论上只需要观测 2 颗恒星即可解算。在实际天文测量中,考虑到大气折射改正不彻底和人仪差的影响,天顶距 z 精度较差,无法用于解算,为此引入天顶距误差 z2。天顶距误差主要与大气折射残留差 和由仪器指标差引起的量测高度差 zi有关,即 z=+zi。由于观测时高度相差很小(1以内),可认为 和 zi是常数。所以 z 可以作为系统误差在解算中求出,即观测的天顶距应为 z+z。因此,在解算式(1)时至少需观测 3 颗恒星构成 3 个方程,但在实际测量工作中为了提高观测结果的精度往往观测 12 24 颗星(合为一组)4,并且要求所选测星应该均匀分布在等高圈上12。1.2 IGG抗差方案的多星近似等高法在实际天文作业中,不管是外业测量还是内业数据处理,都无法避免的存在粗差。目前天文定位多采用自动化、集成化的基于视频测量的自动天文测量系统,当利用测量机器人星敏感器图像传感器(charge-coupled device,CCD)对目标进行连续多次随机拍照、并对所拍摄的星图进行图形处理的过程中,由于算法问题,星点目标识别与提取、星点质心定位、像素坐标差值与度盘差值之间的转换等模型都不可避免地带来误差,产生精度损失13。而现阶段的选星算法存在缺乏综合考虑星表编制的限制条件、恒星的空间几何分布及可测性等问题14,导致测星的几何分布结构难以满足多星近似等高法中所要求的均匀对称分布,也可以视为观测误差,带来无法量化的精度损失。1.1 节中采用经典最小二乘法用以天文定位解算,当观测样本服从正态分布时,参数在最小二乘估计中可以具有无偏性、一致性和最优性15。但是,实际天文观测量不会严格服从正态分布规律。当观测环境恶劣,观测数据中存在粗差时,采用最小二乘法进行解算会对精度产生影响。根据最小二乘准则,ni=1piv2i是残差的平方函数,若天文观测值中掺入异常粗差值,单个观测值的偏差就会导致定位结果产生较大偏差,而为了得到ni=1piv2i的极值,必然会迁就异常值,所以传统最小二乘法不具有抗干扰性16-17。基于此,本文将 IGG抗差估计方案应用于多星近似等高法,目的是检测并剔除异常观测值以获得可靠的定位结果。方案可以为离群值分配零权重,通过降低权重承载误差值较大的观测量,并充分利用高精度观测信息,得到最优的解算结果。1.2.1 解算公式推导为提高天文观测定位精度,应观测多颗星并对每颗星进行多次测量,取平差结果。在式(1)的基础上方程式可以写为cos(zi+z+vi)=sin(0+)sin i+cos(0+)cos cos(Si-i+0+)(4)式中,为测站纬度 与测站概略纬度 0之差,即:=-0;为测站经度 与测站概略经度 0之差,即:=-0;其中测站概略纬度0与测站概略经度 0从卫星计时器的导航数据中提取得到(精度要求不高);赤经 i、赤纬 i、观测瞬间格林尼治真恒星时 Si和天顶距 zi获取方式已经在 1.1 节中介绍,vi表示残差。因此,方程中只需解、和 z 三个参数,具体过程如下。设观测值数量为 n,设置 0、0作为经度和纬度的初值以便利用级数展开使得方差线性化,得到V=AX+L(5)式中,V 表示 n 维残差向量;A 表示 nt 阶系数矩阵;X为 3 维参数向量;L 为 n 维观测向量,具体构成为674第37 卷 第4 期栗志昊,詹银虎,王成龙.抗差估计在多星近似等高法中的应用V=v1v2vn n1(6)X=z31(7)A=a1b11a2b21anbn1n3(8)L=-l1l2lnn1(9)其中ai=cos 0sin i-sin 0cos icos timi(10)bi=-cos 0cos isin timi(11)mi=1-P2i0(12)li=arccos Pi0-zi(13)Pi0=cos 0sin i+cos 0cos icos(Si-i+0)(14)设观测值不相关,权矩阵为单位矩阵,根据最小二乘准则函数ni=1piv2i=min(15)待估参数向量结果为X=-(ATPA)-1ATPL(16)引入 IGG方案后,这里以等价权 Pi代替经验单位权阵 P,等价权形式为Pi=pivi k0k0vipik1-vik1-k0()2k0vi k1(i=1,2,n)0vi k1(17)式中,vi=viv(i)为标准化残差;vi与 v(i)表示第 i 次量测值的残差与中误差大小;pi的初值为1;k0和 k1是抗差方案中的阈值系数,通常 k0=1.01.5,k1=2.58.0。由于权因子是残差 vi的非线性函数,等价权矩阵中就包含了平差的待估量,所以式(16)一般采用迭代法求解,表达式见式(18),其中 k 表示迭代计算次数,X(k+1)表示每次迭代的参数增量。X(k+1)=(ATPkiA)-1ATP(k)iL(18)从算法的推导过程中可以看出,引入 IGG抗差方案后将实际天文观测数据按质量划分在:保权区、降权区和拒绝区。这样做的基本原则是:充分利用保权区的有效信息;限制利用降 权 区 的 可 疑 信 息;排 除 拒 绝 区 的 有 害信息18。主体观测数据应当处于保权区范围,可有效保证天文解算效率。对于受到微小干扰影响的观测数据,比如在选星结构不良情况下的观测数据,应将其视为可疑信息划为降权区,通过改变权值以减小扰动的影响。当观测数据出现大的异常误差,比如受机械装置转动造成的明显误差,应视为有害信息划为拒绝区将其剔除。由此可见,降权区综合考虑了算法的效率和可靠性,而抗差能力由降权区和拒绝区来完成。1.2.2 阈值系数迭代更新标准化残差符合正态分布规律8,在实际算法程序中,IGG权函数的阈值系数 k0和 k1根据前一次迭代的解算结果中各时段中误差大小进行调整。通常依据经验值,将前一时段解算的中误差i从小到大排列构成序列i,找到序列i的 均 值 0,并 根 据 式(4)得 到 i86.6%和i98.7%。i86.6%=(imax-imin)86.6%+imini98.7%=(imax-imin)98.7%+imin(19)利 用 序 列 三 个 特 征 值 0、i86.6%和i98.7%根据式(19)构造新的阈值系数。一般来说,k0常取值 1.5,k1常取值 2.5。774北京测绘第37 卷 第4 期k0=i86.8%0k1=i98.7%0(20)1.2.3 算法流程引入 IGG方案,改进后的算法需要迭代求解,其步骤为:(1)使用经典最小二乘平差进行多星近似等高法的初次解算,得参数 Xi、残差 vi和中误差 v(i)。(2)根据前一次解算的残差 vi和中误差v(i),构建标准化残差vi。(3)根据前一次迭代结果中误差的比例关系调整 IGG阈值系数,利用 k0和 k1改变不同观测质量的数据权值,动态调整权矩阵结构。(4)利用调整后的权矩阵 P(k)i进行再次平差解算。(5)重复 Step2、Step3 和 Step4 三个步骤,直到 X(k+1)-X(k),为迭代收敛精度,常取经验值。具体算法流程如图 1 所示。#GN?#1N*L&*1&CA1A1z5ABBL3B.LG1E5YP=Ik1k0(P)(k)iX(k+1)X(k)A1Xiviv(i)图 1 算法流程图2 数据结果处理与分析为了验证将 IGG抗差方案引入多星近似等高法的合理性和有效性,实验小组利用某自动天文测量系统进行实地观测试验,该系统由卫星天文计时器、Leica TS60 高精度图像全站仪(标称测角精度为 0.5)、便携式计算机以及相应的附属设备和软件组成。2.1 算法应用2021 年 3 月 22 日,观测条件良好,在郑州某天文点使用天文测量系统,设置观测天顶距为50等高圈,观测 24 个时段,每时段选择 12 颗星,每颗星连续观测 8 次,并按照方位角均匀分布的原则选星。分别利用传统多星近似等高法(即单位