抗差估计在多星近似等高法中的应用
栗志昊 詹银虎 王成龙
(战略支援部队信息工程大学, 河南 郑州 450001)
摘 要：传统天文测量中,利用多星近似等高法测定天文经纬度时默认平差解算过程中的权矩阵为单位阵,在一定程度上影响了观测定位精度。为了综合考虑已知的误差项和尚未确定的影响因素,利用抗差估计调整观测数据的权值来降低异常观测量对定位结果的影响。为此,提出将IGGⅢ抗差估计方案引入平差解算,根据前一次解算结果的中误差大小动态调整方案阈值系数,实现权矩阵结构的动态迭代更新,并推导出详细的解算公式。通过实地观测并结合仿真试验验证了理论可行性,能够在高精度全站仪观测条件下将天文定位精度提升10%以上,并且证明抗差方案可以有效削弱不良选星空间结构对于解算结果的影响。
0 引言
光学天文大地测量技术是大地测量技术的重要组成部分,在大地水准面模型验证、高精度垂线偏差测量、地震等自然灾害预报、地球自转参数监测与预报以及深空探测等领域有着广泛的应用[1-2]。目前,配备图像全站仪的视频测量机器人已经实现了高精度、集成化,并且在向全自动化无人干预、智能化阶段推进[3-4]。而在天文定位算法方面,目前广泛采用的是利用多星近似等高法同时测定天文经纬度[5],该理论将大气垂直折光模型误差作为系统参数同天文经纬度一同解算[6],有效削弱了大气垂直折射模型误差对天文定位的影响并大大提高了观测效率。但是现阶段成熟的解算系统中,均认为观测数据质量差异较小,平差解算时默认权矩阵为单位阵,在一定程度上影响了解算精度。
为了区分不同质量的天文观测数据对于解算结果的贡献度,本文将IGGⅢ抗差估计方案[7]引入平差解算程序中,利用前一次迭代的解算结果中各时段中误差大小来确定下一次解算过程中IGGⅢ权函数的阈值系数k0和k1[8-9],实现权矩阵结构的动态调整,兼顾算法效率的同时提升抗差能力。通过处理实地观测数据得到结论,与传统多星近似等高法相比,该算法一定程度上提升了解算精度,并且选星空间结构越差时算法应用效果越好。
1 基于抗差估计的多星近似等高法
1.1 多星近似等高法[2]
根据球面几何中的边余弦公式[10],可以得到天顶距z(或高度角h),待测点天文经纬度λ、φ,观测恒星的赤经α、赤纬δ和时角t之间的关系,公式为
cos z=sin φsin δ+cos φcos δcos t
(1)
(2)
式中,时角t与观测恒星赤经α、待测点天文经度λ和观测瞬间格林尼治真恒星时S有关,具体关系式为
t=S-α+λ
(3)
观测天体的赤经α、赤纬δ可通过查询星表并进行恒星视位置计算得到,观测天体的天顶距z(或高度角h)则通过仪器直接测量得到,观测瞬间格林尼治真恒星时S可以由观测瞬间的UTC时刻加改正换算解得[11]。由于大气折射干扰对于天文观测结果影响较大,在正式解算前需要对高度数据进行大气折射改正。因此,式(1)只含有待测点天文经纬度φ、λ两个未知数,理论上只需要观测2颗恒星即可解算。
在实际天文测量中,考虑到大气折射改正不彻底和人仪差的影响,天顶距z精度较差,无法用于解算,为此引入天顶距误差Δz[2]。天顶距误差主要与大气折射残留差Δρ和由仪器指标差引起的量测高度差Δ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]。但是,实际天文观测量不会严格服从正态分布规律。当观测环境恶劣,观测数据中存在粗差时,采用最小二乘法进行解算会对精度产生影响。根据最小二乘准则,
是残差的平方函数,若天文观测值中掺入异常粗差值,单个观测值的偏差就会导致定位结果产生较大偏差,而为了得到
的极值,必然会迁就异常值,所以传统最小二乘法不具有抗干扰性[16-17]。
基于此,本文将IGGⅢ抗差估计方案应用于多星近似等高法,目的是检测并剔除异常观测值以获得可靠的定位结果。方案可以为离群值分配零权重,通过降低权重承载误差值较大的观测量,并充分利用高精度观测信息,得到最优的解算结果。
1.2.1 解算公式推导
为提高天文观测定位精度,应观测多颗星并对每颗星进行多次测量,取平差结果。在式(1)的基础上方程式可以写为
(4)
式中,Δφ为测站纬度φ与测站概略纬度φ0之差,即:Δφ=φ-φ0;Δλ为测站经度λ与测站概略经度λ0之差,即:Δλ=λ-λ0;其中测站概略纬度φ0与测站概略经度λ0从卫星计时器的导航数据中提取得到(精度要求不高);赤经αi、赤纬δi、观测瞬间格林尼治真恒星时Si和天顶距zi获取方式已经在1.1节中介绍,vi表示残差。因此,方程中只需解Δφ、Δλ和Δz三个参数,具体过程如下。
设观测值数量为n,设置λ0、φ0作为经度和纬度的初值以便利用级数展开使得方差线性化,得到
(5)
式中,V表示n维残差向量;A表示n×t阶系数矩阵;
为3维参数向量;L为n维观测向量,具体构成为
(6)
(7)
(8)
(9)
其中
(10)
(11)
(12)
li=arccos Pi0-zi
(13)
Pi0=cos φ0sin δi+cos φ0cos δicos(Si-αi+λ0)
(14)
设观测值不相关,权矩阵为单位矩阵,根据最小二乘准则函数
(15)
待估参数向量结果为
(16)
引入IGGⅢ方案后,这里以等价权
代替经验单位权阵P,等价权形式为
(17)
式中,
为标准化残差;vi与σv(i)表示第i次量测值的残差与中误差大小;pi的初值为1;k0和k1是抗差方案中的阈值系数,通常k0=1.0～1.5,k1=2.5～8.0。
由于权因子是残差vi的非线性函数,等价权矩阵中就包含了平差的待估量,所以式(16)一般采用迭代法求解,表达式见式(18),其中k表示迭代计算次数,
表示每次迭代的参数增量。
(18)
从算法的推导过程中可以看出,引入IGGⅢ抗差方案后将实际天文观测数据按质量划分在:保权区、降权区和拒绝区。这样做的基本原则是:充分利用保权区的有效信息;限制利用降权区的可疑信息;排除拒绝区的有害信息[18]。
主体观测数据应当处于保权区范围,可有效保证天文解算效率。对于受到微小干扰影响的观测数据,比如在选星结构不良情况下的观测数据,应将其视为可疑信息划为降权区,通过改变权值以减小扰动的影响。当观测数据出现大的异常误差,比如受机械装置转动造成的明显误差,应视为有害信息划为拒绝区将其剔除。由此可见,降权区综合考虑了算法的效率和可靠性,而抗差能力由降权区和拒绝区来完成。
1.2.2 阈值系数迭代更新
标准化残差符合正态分布规律[8],在实际算法程序中,IGGⅢ权函数的阈值系数k0和k1根据前一次迭代的解算结果中各时段中误差大小进行调整。
通常依据经验值,将前一时段解算的中误差σi从小到大排列构成序列[σi],找到序列[σi]的均值σ0,并根据式(4)得到[σi]86.6%和[σi]98.7%。
(19)
利用序列三个特征值σ0、[σi]86.6%和[σi]98.7%根据式(19)构造新的阈值系数。一般来说,k0常取值1.5,k1常取值2.5。
(20)
1.2.3 算法流程
引入IGGⅢ方案,改进后的算法需要迭代求解,其步骤为:
(1)使用经典最小二乘平差进行多星近似等高法的初次解算,得参数
残差vi和中误差σv(i)。
(2)根据前一次解算的残差vi和中误差σv(i),构建标准化残差
(3)根据前一次迭代结果中误差的比例关系调整IGGⅢ阈值系数,利用k0和k1改变不同观测质量的数据权值,动态调整权矩阵结构。
(4)利用调整后的权矩阵
进行再次平差解算。
(5)重复Step2、Step3和Step4三个步骤,直到
为迭代收敛精度,常取经验值。
具体算法流程如图1所示。
图1 算法流程图
2 数据结果处理与分析
为了验证将IGGⅢ抗差方案引入多星近似等高法的合理性和有效性,实验小组利用某自动天文测量系统进行实地观测试验,该系统由卫星天文计时器、Leica TS60高精度图像全站仪(标称测角精度为0.5″)、便携式计算机以及相应的附属设备和软件组成。
2.1 算法应用
2021年3月22日,观测条件良好,在郑州某天文点使用天文测量系统,设置观测天顶距为50°等高圈,观测24个时段,每时段选择12颗星,每颗星连续观测8次,并按照方位角均匀分布的原则选星。分别利用传统多星近似等高法(即单位权方案)、阈值系数采用经验值的抗差方案(k0=1.5和k1=2.5)和阈值系数采用迭代更新的抗差方案三种算法分别处理同一观测数据进行解算,解算结果如图2和表1所示。
表1 三种解算方案的结果标准差
图2 三种解算方案效果评估
图2中的三幅图分别表示应用不同算法,24个时段数据解算下的经度、纬度和大气折射残差结果对比(经纬度结果只保留秒级单位数据)。为了验证算法的应用效果,本次试验利用表征观测稳定度的定位结果标准差作为评价指标来验证算法有效性,结合表1中列出的部分统计结果可得结论:
(1)相较于引入IGGⅢ方案,阈值系数采用经验值的抗差方案(k0=1.5和k1=2.5),对于精度没有明显改进;并且IGGⅢ方案对于削弱大气折射改正残余误差也没有明显效果。
(2)采用1.2.2节方案,迭代更新IGGⅢ方案的阈值系数,经度结果提升2.75%,纬度结果提升5.38%,由于当天观测条件较好且仪器精度高造成结果提升不明显,但是从图表中可以观察到时段5和时段21数据结果得到明显改善。
(3)针对改善较为明显的时段5和时段21,分析原始观测数据,发现这两个时段选星空间结构明显劣于时段6和时段11,绘制这四个时段的选星方位角结果如图3所示。结果证明,在实际天文观测定位过程中,选星分布的均匀程度对于定位解算结果影响很大,要求所选测星要尽量满足均匀分布的原则。从时段5、21与时段6、12的对比情况来看,引入抗差方案可以在一定程度上削弱较差的选星空间分布对于结果的不良影响。
图3 不同时段下的选星空间分布(每时段均观测96次)
2.2 算法效果
为了进一步验证抗差方案对于削弱不良选星空间分布对于解算影响的效果,于2.1节中进行的试验次日,在相同观测条件下连续观测32个时段,每时段选择12颗星,每颗星连续观测8次。并在原始数据基础上,对时段1、3、4、7、11、13、14、16、18、24、32的水平角原始观测数据中加入随机噪声(劣化所选测星的空间分布结构),再利用传统多星近似等高法(即单位权方案)和阈值系数采用迭代更新的抗差方案两种算法分别解算,解算结果见图4和表2。
表2 阈值系数迭代更新的IGGⅢ方案解算精度提升效果
图4 两种解算方案效果评估
结合表2列出的相关统计信息可得出以下结论:
(1)本次实验方案相较于单位权解算,纬度精度提升了16.86%,经度精度提升了13.04%,解算效果得到明显提高。
(2)从改善明显的时段数据来看,加入噪声的时段定位结果都得到有效改善,证明抗差方案可以削弱不良选星空间分布对于解算结果的影响。
3 结束语
为解决多星近似等高法中由观测条件引入的粗差影响,引入IGGⅢ抗差估计方案通过调整权阵结构削弱粗差带来的影响,并以代表观测结果稳定度的定位结果标准差作为评价指标。最后通过实验表明以下两个结论:
(1)采用阈值系数迭代更新的抗差方案,可以有效提高多星近似等高法的定位结果精度。
(2)引入抗差估计方案后,可以有效削弱不良选星空间分布对于解算结果的影响,在推动天文测量向高精度、智能化方向发展方面具有参考意义。
参考文献
[1] 张超,詹银虎,王若璞,等. 光学天文大地测量技术发展评述[J]. 测绘科学技术学报, 2021, 38(4): 331-336, 342.
[2] 张超,郑勇,夏治国. 多星近似等高法同时测定天文经纬度的研究[J]. 测绘通报, 2006 (11): 3-5, 15.
[3] 时春霖. 基于视频测量机器人的自动天文测量技术研究[D].郑州:信息工程大学, 2018.
[4] 王若璞. 新型大地天文测量系统工程化研究[D]. 郑州:信息工程大学, 2004.
[5] 张超. 基于电子经纬仪的天文测量系统及应用研究[D]. 郑州:信息工程大学, 2009.
[6] 邸凯昌,王镓,邢琰,等. 深空探测车环境感知与导航定位技术进展与展望[J]. 测绘学报, 2021, 50(11): 1457-1468.
[7] 牛国军,陈芳. 抗差估计(IGGI方案)在粗差探测中的应用[J]. 西部探矿工程, 2005 (8): 64-66.
[8] 张恩政,唐宁敏,陈刚,等. 基于改进IGG3权函数距离误差模型的工业机器人标定[J]. 中国机械工程, 2021, 32(13): 1539-1546.
[9] SARACENO G, AGOSTINELLI C, GRECO L. Robust estimation for multivariate wrapped models[J]. METRON, 2021,79(2): 225-240.
[10] 张盼盼,郭春喜,马艳鸽. GNSS高程拟合中的粗差探测[J]. 北京测绘, 2018, 32(8): 982-986.
[11] 刘丽丽,林雪原,郁丰,等. 一种SINS/CNS/GNSS组合导航滤波算法[J]. 大地测量与地球动力学, 2021, 41(7): 676-681.
[12] 陈张雷,李崇辉,郑勇,等. 天文定位中几何精度衰减因子最小值分析[J]. 测绘学报, 2019, 48(7): 879-888.
[13] 陈熙源, 柳笛, 刘晓,等. 一种适用于弹道导弹INS/CNS/GNSS组合导航系统的最优数据融合方法:CN201910907858.9[P]. 2022.
[14] 周玟龙,詹银虎,张鹤,等. 基于GDOP贡献值递推的自动天文选星算法[J]. 测绘科学技术学报, 2021, 38(3): 240-245.
[15] LI L,ZHANG G Z,LIU J Z, et al. Estimation of fracture density and orientation from azimuthal elastic impedance difference through singular value decomposition[J]. Petroleum Science, 2021, 18(6): 1675-1688.
[16] 卢庆凯,王家峰. 浅探抗差估计在水准数据处理中的应用[J]. 科学技术创新, 2021(34): 9-11.
[17] 段顺利,孙伟,吴增林. 抗差自适应EKF在INS/GNSS紧组合中的应用[J]. 电子科技大学学报, 2019, 48(2):216-220.
[18] 吴江飞,杨元喜.相关GPS基线向量网的抗差估计[J]. 测绘学报, 2001, 30(3): 247-251.
引文格式：栗志昊,詹银虎,王成龙.抗差估计在多星近似等高法中的应用[J].北京测绘,2023,37(4):475-480.
基金项目： 国家自然科学基金(42074013)
作者简介：栗志昊(1999— ),男,河北石家庄人,本科,研究方向为天文导航。E-mail：531662918@qq.com
通信作者： 詹银虎,E-mail: oscardad@163.com
转自：“测绘学术资讯”微信公众号
如有侵权，请联系本站删除！
