格网内插法坐标转换①

张成校，顾和和

（中国矿业大学 环境与测绘学院，江苏 徐州221008）

0 引 言

2008年7月，2000国家大地坐标系（CGCS2000）正式使用，基于1954北京坐标系（BJS54）和1980西安坐标系（XAS80）的各种坐标成果需要转换到CGCS2000。我国使用的1954年北京坐标系（BJS54）和1980年西安坐标系（XAS80），由于受当时技术手段的限制，存在较大的系统误差和累积误差，控制网间存在着扭曲变形；而基于GPS、VLBI、SLR等空间技术标定的CGCS2000精度均匀，不存在明显的误差积累且控制网间具有较好的一致性。为使转换后的坐标与新坐标系具有较好的一致性，同时满足各种比例尺地形图特别是大比例尺地形图的坐标转换精度需求，需要研究高精度坐标转换方法。日本、美国、澳大利亚［1－5]等国普遍采用格网法作为其大地坐标的主要的转换模型，完成了新旧坐标系的转换，取得了很好的转换精度。如日本采用克里金内插法，生成30s×45s的格网数据，转换精度达到±0.02 m。美国采用最小曲率内插法生成格网，全国分为多个区域，精度在±0.05～±0.5m.澳大利亚采用最小二乘配置内插格网节点改正数，全国分为多个区域，AGD84与GAD94之间转换精度达到±0.04～±0.05m，AGD66与 GAD94之间转换精度达到±0.03～±0.23m.

1 格网内插法原理

格网内插法坐标转换的基本思想［4－8]是将大的转换区域划分成小的格网单元，利用两个系统间离散点的坐标差，采用一定的内插方法计算具有一定间隔的格网节点的坐标差，利用格网节点上的坐标差内插其它任意点上的坐标差，实现不同坐标系坐标的变换。格网内插法原理如图1.

图1 格网内插法原理图

1.1 建立格网

将转换区域划分为具有一定间隔的规则格网，利用该区域内的若干个具有新旧两套坐标的公共点的坐标差值，使用数学模型求得各个规则格网节点对应的转换到新坐标系的改正量。

1.2 格网内插

根据待求点所处的格网单元，提取格网文件对应的4个格网节点的坐标改正量，采用双线性内插模型进行插值计算。原理如图2所示。

图2 由格网节点内插某一点的坐标改正量

双线性内插的公式为

其中：

t1，t2，t3，t4分别为内插点（三角形表示）周围的4个格网节点的坐标改正量；（B1，L1）、（B2，L2）、（B3，L3）、（B4，L4）分别为4个格网节点所在旧坐标系的经纬度坐标；（Bp，Lp）为待求点的经纬度坐标；z为待求点的坐标改正量。

1.3 计算步骤

1）按照一定的经纬度间隔，转换区域被划分成规则的格网根据公共点的坐标差，用最小二乘回归拟合一个简单的平面模型：ax＋by＋c＝z（x，y）；

2）用平面回归模型拟合值减去已知值得到一组残差；

3）用数学模型将格网节点周围一定搜索半径内公共点的残差拟合到格网节点上，得到格网节点上的残差值；

4）格网节点上的平面回归模型拟合值加残差值，得到格网节点上对应的新旧坐标系坐标改正量；

5）经双线性内插得到格网区域内任意一点的坐标改正量，加到旧坐标上就得到了转换后的新坐标，从而完成坐标转换。

2 格网内插法的数学模型

求取格网节点改正量的数学模型可以有多种选择，这里简单地介绍了几种模型。

2.1 加权平均模型

加权平均模型［6－7]根据周围数据点对于中心点的远近分配不同的权重，离中心点越近的数据点对中心点的影响越大，占的权重就越大，反之占的权重就越小。一般形式为

其中：X是中心点处的估值；xi是中心点周围数据点；wi是数据点对应的权。权函数有多种形式，这里取wi＝，di是数据点到中心点的距离，dmin是各数据点到中心点的最小距离。

2.2 克里金内插模型

克里金内插法［9－10]以由区域化随机变量z（x）的变异函数γ（x，y）或期望μ（x）＝E［z（x）]和协方差函数c（x，y）量化的随机模型为基础，由未观测点x0附近x1、x2、…、xn处的观测值zi＝z（xi）（i＝1，…，n）内插z（x）在x0处的值z（x0），得到z（x0）的最优无偏估计量＾z（x0）。

克里金内插法的一般公式为

式中：z（xi）为采样点xi处的观测值；x0是一个未采样点；wi（x0）为赋给观测值z（xi）的权，其估算服从无偏条件

最小的条件下得到。

普通克里金法的应用条件是：区域化随机变量满足内蕴假设条件，即随机变量z（x）的数学期望为常数；有足够的观测量确定变异函数，即随机变量的变异函数γ（x，y）为已知。

2.3 最小曲率模型

最小曲率法［7]通过将具有观测值的区域格网化，在尽可能尊重观测值的前提下，以总曲率最小的原则生成最平滑的曲面。格网点（xi，yi）上的曲率为

3 算例比较

由于CGCS2000坐标数据比较少，因此，只取某地区124个具有北京54坐标和CGCS2000坐标的点进行转换比较，其中均匀的提取其中的14个作为检核点；为了便于格网的划分，将该区域取为矩形，纬度跨度为3°，经度跨度为4°.分别用Bursa模型和最小曲率模型进行坐标转换，对转换结果进行比较。

方案1：使用Bursa模型计算待求点的坐标转换值，即计算得到的CGCS2000坐标值；将其与已知的CGCS2000坐标值进行比较得到相应的残差，如表1所示。

方案2：利用最小二乘曲率模型构建格网并内插待求点得到一组残差，如表2所示。

两种方案转换的结果将按残差平均值、标准偏差、绝对值最大值这三项指标进行比较。为便于分析，经纬度残差均换算成弧长，以m为单位。

表1 Bursa模型的残差统计

表2 最小曲率模型的残差统计

从表1看出，方案1的残差平均值和标准差较大，最大残差绝对值接近3m，残差的数量级基本在0.1m以上；从内部残差统计来看，绝对值小于0.1m的不到10%，说明模型内符合不好；从外部残差统计来看，绝对值小于0.1m的不到7%，说明模型外符合不好。这反映了我国天文大地网存在较大局部系统差。如果整个区域采用Bursa模型进行坐标转换，将不能顾及天文大地网的局部扭曲和累积误差，不能充分的拟合新旧坐标系间的坐标差异。

从表2看出，方案2转换结果的各项指标都明显好于方案1，从内部残差统计来看，标准差降到了厘米级，残差绝对值在0.3m以内，系统差已基本被消除，残差绝对值小于0.1m所占的比例超过了93%，说明最小曲率模型内符合精度较好。从外部残差来看，残差平均值和标准差都比方案1大大减小，绝对值最大值控制在0.6m以内，残差绝对值小于0.1m所占的比例超过了80%，说明最小二乘曲率模型的外符合精度也较好。总之，最小曲率模型较好地拟合了新旧坐标系间的坐标差异，有效地控制了天文大地网的局部扭曲和积累误差。

4 结 论

我国CGCS2000刚刚启用，也涉及到新旧坐标的转换问题。相对于采用Bursa模型进行坐标转换，采用格网内插转换方法转换精度高，适合用来实现高精度坐标转换。采用格网内插法坐标转换需要较多的公共点，公共点数量较多，则转换效果较好；若公共点的数量较少，则转换效果不明显。

［1]施建平.美国等国家大地坐标系的更新转换和维持［J].科协论坛，2009（7）：112－113.

［2]DEWHURST W T.NADCON－The application of minimum curvature derived surfaces in the transformation of positional data from the north american datum of 1927to the north american datum of 1983［R].NOAA Technical Memorandum NOSNGS－50，1990.

［3]MIKIO T.Coordinate transformation software“TKY2JGD”from Tokyo datum to a geocentric reference system，Japanese geodetic datum 2000［J].Times of Geographical Survey Institute，2001（97）：31－57.

［4]COLLIER P.Developmemt of Australia’s national GDA94transformation grid［R].Consultant’s Report to the Intergovernmental Committee On Surveying and Mapping，2002.

［5]李克恭，牛岸英，张斌才，等.格网法实现基础测绘成果向CGCS2000的整体转换［J].矿山测量，2010（4）：68－70.

［6]郭 充，吕志平，李 岩，等.基于格网的坐标转换方法［J].信息工程大学学报，2010，11（2）：166－169.

［7]郭 充，吕志平，徐基刚.几种模型在格网坐标转换应用中的比较与分析［J].测绘通报，2009（5）：38－41.

［8]郭 充，吕志平，于兴超，等.基于Bursa加权模型的格网坐标转换［J].测绘科学技术学报，2009，26（2）：86－88.

［9]施建平，杨华忠.克里金内插法实现坐标转换的应用研究［J].海洋测绘，2010，30（4）：30－32.

［10]张景雄.空间信息的尺度、不确定性与融合［M].武汉：武汉大学出版社，2008：83－129.