在日常的测绘工作中,不同坐标系下的大地坐标值经常需要相互转换,然而现实中由于控制点桩的建立年代

不同、精度等级不一以及标石被移动破坏等因素[1],导致参与九参数解算的公共点坐标中经常会伴有粗差,

直接影响参数的解算精度,导致坐标转换失败。

现实中用于解决大地坐标转换时公共点坐标粗差的算法较多,比如Baarda粗差数据探测法、抗差最小二乘法

及拟准检定法等。Baarda算法通过剔除粗差坐标参与九参数解算的方式来计算未知转换参数;抗差最小二乘

法采用选权迭代的方式减少粗差坐标在参数解算时的影响;拟准检定法通过特定算法发现并改正公共点粗差

坐标值,从而实现参数的精确解算。为了获得以上3种算法用于不同大地坐标转换时的抗差效果,文章通过具

体的算例得出3种上述不同抗差算法的精度对比结论,从而为工程实际提供一定的参考依据。

1 不同大地坐标系间相互转换的误差方程
不同大地坐标系间坐标的相互转换,通过公共点坐标可建立误差方程式[2]:

 
式中,X=[△X△Y△Zεxεyεzm△a△f]T为待求的坐标转换参数,(△X,△Y,△Z)为3个方向的坐标平移参

数,(εx,εy,εz)为x、y、z 3个方向的坐标轴旋转参数,m为缩放比例参数,△a为椭球的长半轴差,△f为扁

率差,未知参数共计9个。

公共点的误差方程常数项矩阵Ai以及系数矩阵li分别为[3]:

 
 
式(2)、(3)中Li为大地经度,Bi为大地纬度,Hi为大地高,a为地球长半轴,e为地球参考椭球的扁率,Mi为地球

经线圈曲率半径,Ni为地球卯酉圈曲率半径,为第一偏心率,为地球扁率,△Bi为纬度差,△Li为经度差,△Hi

为大地高差,ρ*≈206 265*。

2 坐标转换抗差算法
2.1 Baarda粗差数据探测法
Baarda数据探测法用于不同大地坐标转换是将公共点坐标粗差纳入函数模型,采用检验统计量判定并定位坐

标粗差位置,之后使含有粗差的公共点坐标不参与未知参数解算,最后采用最小二乘原理并结合余下不含粗

差的公共点解算转换九参数。[4]

解算步骤如下:

(1)将原有空间直角坐标系下的n个公共点坐标设定为固定值,令大地坐标系的坐标为等精度的独立观测值,

运用最小二乘原理计算出观测值的改正数Vi与协因数矩阵QViVi,构成u检验统计量[5]:

 
式中,σ0为母体标准差,α(比如α=0.1)为给定显著水平。

(2）查看标准正态分布表得，若，此时剔除该统计量所对应的公共点，使其不参与九参数的解算；

(3）将剩余公共点坐标按最小二乘原理进行平差计算，并按步骤（1）重新构成u统计量进行检验，直到所

有公共点均满足，即公共点坐标无粗差为止。

注：σ20可用平差后求得的方差估值代替，构成τ检验统计量，τ统计量为服从自由度f=2n-4的τ分布。

(4)将剩下未超限的(至少3个)公共点坐标代入式(1),通过最小二乘法求解出转换九参数,将转换参数代入式

(1)求解各待定点大地坐标值,即可完成大地坐标转换。

2.2 抗差最小二乘法
为了达到抵抗粗差的目的,抗差最小二乘法用于不同大地坐标转换是利用抗差权函数(如IGGIII权函数)对公

共点的各方向坐标重新进行定权[6],以降低粗差坐标值在九参数解算时的权重,之后采用迭代计算的方式计

算未知九参数。

解算步骤[7]:

(1)根据空间直角坐标各观测向量的协方差矩阵D1=δ12、D2=δ22,计算初始权矩阵PL=(D1+D2)-1;

(2)依据式(1)转换模型,采用最小二乘原理,求解初始转换参数Xi=(BTPB)-1BTPL与初始残差值

(3)结合初始残差值采用抗差权函数编写matlab程序代码选权迭代计算,直到前后两次迭代得到的转换参数

差值绝对值均小于给定的K值(按精度要求设定)时,停止迭代计算,得出最终的等价权矩阵,求出抗差最小二

乘准则下的转换九参数值;

(4)将最终满足迭代条件的九参数代入式(1)求解各控制点大地坐标值(△Bi,△Li,△Hi),完成转换。

2.3 拟准检定法
拟准检定法利用真误差与观测值间的确定关系R△=-RL,结合拟稳平差的思想,提出“拟准观测”概念。拟准

观测指认为无粗差且又未加以判定的观测值。[8]

解算步骤:

 
(2)依据单位权中误差σ赞0是否超限判断九参数估值是否合理。若合理,则停止拟准观测(公共点坐标无粗

差存在),采用直接最小二乘法求解参数;若不合理,需根据步骤(3)、(4)探测粗差坐标位置,并利用式估算粗

差值大小。

式中,R=E-J=E-A(ATPA)-1ATP为正交补投影,E为n阶单位矩阵。Cb构成方式为:假设确定的粗差个数有b个,可

得到b个n维单位向量ej=(0,…,1,0,…,0)T,即对应的第j个观测有粗差,设定第个分量为1,其余为0,则Cb=

(e1,e2,…,eb)。[9]

(3)初选拟准观测

 
式中,G参见式(5):

 
初选拟准观测时,可将初始残差中明显小于其它残差值的观测值设定为初选拟准观测。

(4)复选拟准观测

复选拟准观测是采用前次的r个拟准观测,且对复选拟准观测进行平差,计算真误差估值由大到小排序。生成

真误差拟准解序列[10]:

 
复选拟准观测在实际的操作中可以逐个或多个增加。

计算。当W3>k(如设定k=3.1)时,停止复选拟准观测,其后的观测皆为粗差观测值;反之,继续复选拟准观测,

重复步骤(4),直至满足条件W3>3.1为止。

(5)坐标转换

(1)利用式计算粗差估值;

(2)根据粗差探测位置及粗差估值修正公共点当中的粗差坐标;

(3)将改正后的公共点坐标代入式(1),采用最小二乘原理计算转换九参数,完成坐标转换。

3 算例解析
从某工程控制网当中,选出4个分布较均匀且同时具备WGS84与BJ54坐标的公共点参与九参数求解及坐标转换

的验证。4个公共点的坐标见如表1(中央子午线经度为123°00′00″)所示。

表1 已知点公共点坐标及大地高     下载原表

 表1 已知点公共点坐标及大地高
(1)公共点坐标无粗差情况下。直接利用最小二乘原理解算未知九参数分别为:

 
由于公共点坐标没有粗差,转换后经度、纬度及大地高最大差值分别为0.000 2″、0.000 007″及2 mm,满

足工程中大地坐标转换的精度要求,如表2所示。

表2 直接最小二乘法坐标转换精度(公共点坐标无粗差)     下载原表

 表2 直接最小二乘法坐标转换精度(公共点坐标无粗差)
注:表2中△B、△L、△H分别为转换坐标与已知坐标见的差值

(2)公共点坐标有粗差情况下。为便于研究,在D1点的纬度坐标上人为加入2秒粗差,分别采用直接最小二乘

法、Baarda数据探测法、抗差最小二乘法和拟准检定法进行九参数的解算及坐标转换。

(1)直接最小二乘法。利用最小二乘原理求取其转换九参数值分别为:

 
(2)Baarda数据探测法。将表1中4个控制点均作为公共点参与大地坐标转换,当D1点纬度值存在2秒粗差时,

用Baarda数据探测算法进行粗差探测与定位,可以看出D1点剔除D1点三维坐标参与九参数的计算,之后将D2

、D3、D4点利用直接最小二乘法求参,求取的转换九参数值如下:

 
(3)抗差最小二乘法。当D1点纬度坐标有2秒粗差的情况下,采用抗差估计算法求得的转换九参数值分别为:

 
(4)拟准检定法。用拟准检定法对公共点粗差坐标进行处理,改正公共点粗差坐标,最终采用改正后的公共点

坐标进行求参,求得的转换九参数值分别为:

 
通过上述九参数解算及精度对比表2看出,由于公共点坐标有2秒粗差存在,直接最小二乘法计算出的平移参

数、旋转参数及扁率变化率误差较大。如表3所示,转换后经纬度、大地高最大差值分别为0.001″、0.2″

及4 mm,转换精度差;Baarda数据探测法、抗差最小二乘法及拟准检定法解算出的转换九参数与无公共点坐

标粗差情况下的转换参数值基本一致,且检核点坐标与已知坐标间的差值很小。需要指出的是,由于式(1)转

换模型需要求解9个未知转换参数,至少需要3个公共点参与九参数的求解,该算例中若有2个公共点坐标存在

粗差,则Baarda数据探测法将剔除该2个公共点,此时只剩下2个无粗差的公共点,无法利用最小二乘法进行九

参数的解算。

表3 不同算法转换精度对比(公共点坐标存在粗差)     下载原表

 表3 不同算法转换精度对比(公共点坐标存在粗差)
注:表3中△B、△L、△H分别为转换坐标与已知坐标见的差值

4 结束语
Baarda数据探测法、抗差最小二乘及拟准检定法用于不同大地坐标间的转换均能解决公共点坐标存在粗差

所导致的转换精度差的问题。由于Baarda数据探测法存在剔除粗差公共点参与九参数求解的缺点,因此,实

际工程中若公共点个数有限,建议采用抗差最小二乘和拟准检定法抵抗粗差公共点坐标的影响,实现不同大

地坐标间的转换。