鲁达 《测绘学报》
构建与学术的桥梁 拉近与权威的距离
CGCS2000高精度坐标转换格网模型的建立
吕志平1 , 魏子卿2, 李军3, 郭充4
1. 信息工程大学 地理空间信息学院，河南 郑州 450052；
2. 西安测绘研究所，陕西 西安 710054；
3. 西安测绘信息技术总站，陕西 西安 710054；
4. 65051部队，辽宁 大连116031
收稿日期：2013-01-15；修回日期：2013-03-12
基金项目：国家自然科学基金(41274015)；国家863计划 (2013AA122501)
第一作者简介：吕志平(1960—)，男，博士，教授，博士生导师，主要从事大地测量和测量数据处理等领域的科研与教学工作。E-mail: ssscenter@126.com
摘要：CGCS2000启用后，用于地形图图廓线和坐标网变换的高精度坐标转换格网模型的建立成为实用中的重要问题。本文阐述我国CGCS2000高精度坐标转换格网模型的建立方法。一是讨论坐标转换最小曲率方程，该方程同时具备整体光滑和局部光滑特性，可以更好地符合地图线状要素坐标转换时的连续光滑性；二是提出格网节点迭代计算方法以及拟合推估两步极小法构建区域性格网模型的方法，构建了格网模型数据结构；三是给出全国范围CGCS2000高精度坐标转换格网模型的最优格网间距和转换精度，按48433个2000国家大地控制网点坐标构建了BJS54和XAS80到CGCS2000坐标转换格网模型，转换后的点位外符精度分别为0.26m和0.03m。
The Establishment of High Precise Coordinate Transformation Grid Model of CGCS2000
LV Zhiping1, WEI Ziqing2, LI Jun3, GUO Chong4
Abstract: After implementing CGCS2000, the establishment of the high precise coordinate transformation grid model which is mainly used to transform sheet lines and coordinate grids of topographic map becomes important to be solved in the mapping application. The establishment method of the high precise coordinate transformation grid model of CGCS2000 is proposed. First, the minimum curvature equation of coordinate transformation is discussed. And it possesses both the whole and local smooth characteristics, so it can be better in keeping with consecutive smooth characteristics for the coordinate transformation of map’s linear feature. Second, the iterative calculation method of grid nodes and the establishment method of regional grid model using the resolution of collocation by two minimization steps are advanced.The data structure of grid model is constructed. Third, the countrywide optimized grid interval and transformation precision of high precise coordinate transformation grid model of CGCS2000 are coordinate transformation grid models from BJS54 and XAS80 to the CGCS2000 have been constructed using 48433 points of 2000 national geodetic control network of China, and the positioning external accuracy of both models are 0.26m and 0.03m respectively.
Key words: CGCS2000 coordinate transformation minimum curvature method grid model data structure collocatione
1 引 言
在我国大地测量产品序列中，坐标转换格网模型是2000中国大地坐标系(CGCS2000)[1, 2, 3, 4]启用后，继重力异常格网模型、高程异常格网模型和垂线偏差格网模型之后又一新的大地测量数字格网模型产品形式。
坐标转换格网模型是一定范围内两坐标系坐标分量差值的离散化数字表达，是该范围内所有等间距格网点坐标分量差计算值的集合，在数据库中以格网数据结构的形式存储。格网数据结构是把建库的地理范围按经纬线划分成规则的梯形格网，以格网的纵横交点作为节点，存储对应于该交点的坐标分量差值。
坐标转换格网模型在国际上已有20多年的应用历史。美国1927北美基准(NAD27)至1983北美基准(NAD83)的转换最先采用了这一方法[5]，该方法的转换精度在美国主要大陆地区为±0.15m[6]，满足了测量和制图等工作的需要。后来日本[7]、新西兰[8]、澳大利亚[9]、加拿大[10]等国也相继构建了格网坐标转换模型，不同的是在构建格网模型时，根据本国的实际采用了不同的数学方法，如最小曲率法、最小二乘配置法、多元回归法等。除了经典大地坐标系与现代大地坐标系的转换外，美国还建立了NAD83向高精度区域网(HARN)转换的格网模型，由于HARN是美国各州建立的，所以几乎每个州都有一套格网模型。目前常用的GIS软件(如ArcGIS)中，内置有格网坐标转换模型，许多国家的坐标转换格网文件存储在安装目录下，用户可直接调用。在我国，对于坐标转换的数学模型已有了较深入的研究[11, 12, 13, 14, 15, 16]，但对用于地图坐标变换的坐标转换格网模型的研究和应用还刚起步。
针对我国CGCS2000启用后坐标转换应用的多类型，可按精度指标分3类构建我国坐标转换模型产品，一是用于大众化导航应用的低精度转换，坐标分量精度为±5～±10m；二是用于三维坐标相似转换的中精度转换，坐标分量精度为±0.5～±5m；三是用于地形图图廓线和坐标网转换的高精度转换，坐标分量精度不超过±0.5m。本文针对高精度转换应用类型，阐述了我国CGCS2000高精度坐标转换方法，讨论了坐标转换最小曲率方程，提出了格网节点迭代计算方法和区域性格网模型的建立方法，构建了格网模型数据结构，给出了全国范围CGCS2000高精度坐标转换格网模型的最优格网间距和转换精度。CGCS2000高精度坐标转换格网模型已在全国1∶25万和1∶5万地形图更新生产等工程中应用。
2 坐标转换格网模型的高精度
建立某区域的坐标转换格网模型就是把该区域按经纬线划分成等间距的规则格网单元，根据区域内公共点上已知的坐标分量差值，采用一定的数学模型(如最小曲率法)求得各个格网节点的坐标分量差。坐标转换格网模型能实现高精度的坐标转换，理由如下：
分析产生两坐标框架间坐标差异的原因，可分为以下3类：一是两坐标系的定义差，二是建立该两坐标框架(即大地控制网)时的局部性系统误差，三是建立该两坐标框架时的偶然误差。对于偶然误差，如控制点的坐标值一经确定，则观测偶然误差也将以系统性偏差的特性对控制测图产生影响，而这类误差是非常复杂的。文献[17]针对卫星导航数据的误差特性也作过类似的分析。因此，综合这3类偏差，要用一个数学公式构建高精度坐标转换模型是不可能的。坐标转换格网模型将数学模型和数字模型有机结合，可以高精度地拟合两坐标框架间的各种偏差，从而实现高精度的坐标转换。
3 最小曲率法坐标转换格网模型3.1 坐标转换最小曲率方程
建立格网模型的数学公式较多，有最小二乘配置法、多元回归法、相似变换、最小曲率法等[18]，单从转换精度看，这些模型相差不大。如果某一模型的精度不够，可通过改进计算策略、缩小搜索半径来提高精度。例如文献[13, 14]采用相似变换模型通过缩小分区范围也可获得很高的转换精度。但采用最小曲率法还可得到很好的光滑性。
在第2节中分析的坐标差异产生的3类原因中，定义差和大区域性的系统偏差导致的差异在地理空间是平滑改变的，曲率也是连续的，体现为整体性的光滑；小区域性的系统偏差和测量中的偶然误差导致的差异则为局部性，体现为局部性的光滑。
通过最小曲率法建立最小曲率格网曲面来拟合坐标之间的差异，能够同时满足整体光滑和局部光滑的要求。
把格网化的坐标位移值所属地理区域比作一块大的薄片，各点都对这个薄片施加垂直方向上的拉力，则薄片会发生弯曲和变形，每个点上变形的幅度取决于在这个点上施加的力，距离点较远处的变形则取决于临近的拉力，显然，点之间的弯曲是平滑的。点上的坐标位移值等同于点上施加的向上的拉力，位移表面的平滑弯曲代表着各个位置的坐标位移。
描述平面变形的双谐偏微分方程为
式中，fn为施加的力；L、B为点位的经纬度坐标；u为位移量(分别取两坐标系的经差ΔL、纬差ΔB、大地高差ΔH)。式(1)的解是一个二维空间变量的三阶多项式，要求观测点上的位移值等于观测值，即
式中，u(Lm,Bm)是观测点m处的位移量；wm是点m处的观测值，分别取m点两坐标系的经差ΔLm、纬差ΔBm、大地高差ΔHm。
对于一个平面来说，认为边界外不存在外力，且曲面在边界法线方向上弯曲率为0，即
为边界条件。
格网点(Li,Bi)上的曲率是
，取所有格网节点的曲率平方和为极小，可推导出计算格网节点的估计值ui,j的差分方程[19]
可见要求当前节点的基准位移值ui,j，需知道周围其他12个节点的基准位移值(见图 1)。
图 1 差分方程中的节点布局Fig. 1 The nodes distribution in difference equation
图选项
式(4)适用于离开边界且周围没有观测点的情况。当格网节点位于边界上或附近有观测点，则格网节点的曲率计算公式需要作相应的调整[19, 20]。
由上可知，用最小曲率法生成的插值面类似于一个通过各个数据值的具有最小弯曲量的薄弹性片。最小曲率法试图在尽可能严格地尊重数据(见式(2))的同时，生成尽可能光滑的曲面(见式(4))。式(4)的效果与分段连续样条函数一样，可以同时模拟整体和局部光滑。而采用其他数学模型(如多元回归[18]、分区七参数模型[13, 14]等)建立格网，节点值是采用该节点搜索半径范围内的观测值来计算的，转换精度虽然很好，但并不顾及节点间的光滑性。
3.2 格网节点的最小曲率法迭代计算
由式(4)可知，计算每个节点都要用到邻近的其他节点，所以为了解出差分方程，需要给所有节点赋初值，且进行迭代计算。
迭代时，从代表长波趋势的粗格网开始，在两维方向上逐次加密格网，每次加密格网都需要迭代，这样可以兼顾整体光滑和局部光滑的要求。格网间隔的选择极为重要，间隔过于稀疏会混淆数据，间距过密将不能真实地描绘区域趋势。但究竟该取多大的格网间距才算最优，只能根据多次试验结果，依据我国大地控制网的实际情况确定。
将粗格网细化的方法按如下步骤实现：
(1) 确定最优格网间距以及划分的格网数。为了减少各向异性误差的影响，将二维方向的格网间距取成一致。
(2) 确定粗格网的格网间距。为了能使用差分方程，第一次按粗格网间距划分的格网数必须大于等于4，而且能整除最优格网的格网数。
(3) 依次细化粗格网直到达到最优格网。第一次所用的粗格网节点的初值可以用离格网节点最近的观测值，或者由周围的观测值加权得到。迭代过程中加密格网产生的新格网的初始值可通过上一次细化的粗格网线性内插来获得。初始值越接近渐进可展解的值，所需迭代的总次数就越少。
3.3 格网模型的数据结构
通过格网模型构建了描述两坐标框架差异空间分布的有序数值阵列，它是表示区域D内的3个三维向量有限序列组
式中，Li、Bi是第i个格网点的平面坐标；ΔLi、ΔBi、ΔHi是两坐标系的坐标差值。当该序列中各平面位置格网呈规则排列时，其平面坐标可省略，因每一个格网点与其他相邻格网点之间的拓扑关系都隐含在该阵列的行列号当中。此时式(5)就简化为一维向量序列组
在地图和影像等产品的生产应用中，只需要平面基准的转换，而不需要大地高的转换，因而式(6)只取前两式。
从实用的角度看，采用两维方向间距相同的正方形格网在数据处理方面有很多优点，一是结构简单，除了每个格网点处的坐标分量偏差外，只需要记录一个起算点的位置坐标和格网间距，数据存储量很少；二是计算方便，可以根据该区域的原点坐标和格网间距，对相应矩阵元素的行列号经过简单的运算而获得任意格网点的平面位置；三是适合于实时服务以及大规模使用和管理，可以借助于简单的格网数据处理方法进行数据压缩处理，如行程编码法、四叉树法、多级格网法和霍夫曼码法等。
格网坐标转换模型数据文件包括3个逻辑部分：① 元数据，指描述格网坐标转换模型一般特征的数据，如名称、边界、构建单位、精度等参数；② 数据头，定义格网坐标转换模型的起点坐标、坐标类型、格网间隔、行列数等；③数据体，沿行列分布的坐标分量偏差值数字阵列。
格网坐标转换模型数据库按分块层次结构建立空间索引实现格网数据的快速查找、存取和计算操作等。按国家标准图幅大小进行块的划分，一个块就是一幅图的范围，末级子块由若干行和若干列的格网单元组成。通过行列结构索引，便可唯一确定格网坐标转换模型数据库范围内任意空间位置的坐标偏差值。由这种分块和行列层次划分形成的空间索引可以保证格网坐标转换模型数据的快速查找和无缝存取。
格网模型可采用多级格网结构，不同级格网间距不同，格网越小，所表达的连续坐标分量差值的分辨率越高。例如，将建库地区根据已知数据的分布分为多个区域，不同区域采用不同间距的格网。
4 全国范围CGCS2000坐标转换格网模型的建立4.1 最优格网间距的选取
以一幅典型的1∶50万地形图幅为例说明最优格网间距的选取试验。该图幅范围内(34°B～36°B，114°L～117°L)有410个全国天文大地网与空间大地网联合平差点，这些公共点具有1954北京坐标系(BJS54)、1980西安坐标系(XAS80)和CGCS2000坐标，数据点的平均间距约为7.2′(约13km)。按照每10个点中取1点的方式均匀提取41个点作为检核点，剩下的369个点建立模型。
按照最小曲率法建立格网模型的方法，采取不同格网间距，分析格网间距的差异对最终转换精度的影响。
表 1 不同格网间距对应的内符精度比较Tab. 1 The corresponding internal accuracy comparison of different grid interval
| 坐标转换 | 格网 间距 |
|
| BJS54 → CGCS2000 |
|
| XAS80 → CGCS2000 |
|
表选项
表 2 不同格网间距对应的外符精度比较Tab. 2 The corresponding external accuracy comparison of different grid interval
| 坐标转换 | 格网 间距 |
|
| BJS54 → CGCS2000 |
|
| XAS80 → CGCS2000 |
|
表选项
由最小曲率方程可知，理论上表 1的内符偏差值应为0，实际上由于求内插点值时采用了线性内插，故使之不为0。因而，表 1的内符偏差反映了格网点和内插点之间的坐标差值存在非线性变化的程度。
由表 2可知，转换精度与格网密度相关，在观测数据分辨率范围内，格网间距越小，转换精度越好，但格网密度超过观测数据分辨率时，受高频振荡影响，转换精度反而会变差。
表 2中，BJS54到CGCS2000的转换，格网间距为2.5′×2.5′时，转换精度已经很好，继续加密取1.5′×1.5′时已经没有改进，外符精度反而变差了。说明在2.5′×2.5′范围内，两坐标系差值的变化已比较平缓，按线性内插已基本符合实际情况。因此2.5′×2.5′的格网间距应该作为BJS54到CGCS2000转换的最佳选择。
表 2中，XAS80到CGCS2000的转换，格网间距为3.0′×3.0′时，转换精度已经很好，取2.5′×2.5′时已经没有改进，外符精度反而会变差。说明我国天文大地网经整体平差后局部精度也已提高。因此3.0′×3.0′的格网间距应该作为XAS80到CGCS2000转换的最佳选择。
由此可见，我国大地网的最优格网间距大约为公共点平均间距′)的0.3倍到0.4倍(4～5km)，在这一尺度上，任一点的大地坐标改正值可按双线性内插方法获得。
4.2 格网模型的构建
以全国天文大地网与空间大地网联合平差建立的2000国家大地控制网[21]的48433个点为观测量，建立最小曲率法坐标转换格网模型。格网间隔统一取为2.5′×2.5′。对于XAS80，因3.0′×3.0′与2.5′×2.5′的结果差别不大，为实用方便，也取为2.5′×2.5′。计算结果见表 3。
表 3 全国范围格网模型内符精度Tab. 3 The internal accuracy of the grid model in the whole country
| 坐标转换 | 格网间距 | 公共点数 | 中误差/m | |
| B分量 | L分量 | |||
| BJS54 →CGCS2000 | 2.5′×2.5′ | 48 433 | 0.062 | 0.071 |
| XAS80 →CGCS2000 | 2.5′×2.5′ | 48 433 | 0.006 | 0.012 |
表选项
4.3 格网模型的外部检核
将2000国家大地控制网48433个点按经纬度排序，每隔100点选取1点，在全国范围共选出484个点作为外部检核。剩余47949个点为观测量，建立2.5′×2.5′间隔的格网模型。将未参与构建格网模型的484个点的BJS54(或XAS80)坐标按所建格网模型转换至CGCS2000，并与已知的CGCS2000坐标比较，计算坐标转换标准偏差，结果见表 4。计算表明，我国BJS54到CGCS2000坐标转换格网模型的点位精度约为0.25m，XAS80到CGCS2000坐标转换格网模型的点位精度约为0.03m。相比可见XAS80的局部性误差比BJS54已大大减弱。
表 4 全国范围格网模型外符精度Tab. 4 The external accuracy of the grid model in the whole country
| 坐标 转换 | 格网 间距 | 公共 点数 | 外检 点数 | 中误差/m | 点位误差 /m | |
| B分量 | L分量 | |||||
| BJS54 → CGCS2000 | 2.5′×2.5′ | 48 433 | 484 | 0.124 | 0.259 | 0.259 |
| XAS80 → CGCS2000 | 2.5′×2.5′ | 48 433 | 484 | 0.017 | 0.028 | 0.028 |
表选项
5 区域性CGCS2000坐标转换格网模型的建立
全国性坐标转换格网模型不一定能满足城市区域性测绘工作的应用需求，例如：分米级精度的格网模型不能满足1∶1000地形图的坐标转换精度要求。我国还有十多万个直接用于控制测图的国家低等级控制网点和城市独立坐标系控制网点，这些点尚无CGCS2000成果，地方测绘部门通过区域网联合平差或GPS测量可以获得它们的CGCS2000坐标。公共点数量和密度的增加可以构建格网间距更密、精度更高的坐标转换格网模型。
从谱分析的观点来看，当对一个函数进行分格网取平均，则对该函数在一定程度上作了平滑，即对观测值进行了高频滤波。较大格网面积的数据意味着坐标差数据场的分辨率较低，这必然造成所确定的坐标差高频部分的丢失。区域性格网模型的构建将提高全国性格网模型在该区域的分辨率和精度。
根据区域内新增公共点构建区域格网的方法，采用文献[22]提出的拟合推估两步极小解法，第1步由全国格网点内插新增公共点上的拟合值，第2步由各公共点的剩余误差拟合协方差函数，再用拟合推估法构建区域格网。于是，区域格网节点的坐标改正量由两部分组成，一是全国性格网内插得到的量，二是残差改正量。地图数据库中的要素坐标使用区域格网模型转换得到的结果将与低等级控制点平差值的精度相匹配，相当于将地图坐标也进行了重新平差。
计算试验采用某区域(范围8°×12°)805个具有BJS54和WGS84(视为CGCS2000)坐标的公共点数据进行。均匀选取269个点视为国家等级点，模拟计算全国性格网。剩下的536个点视为城市控制点。536点中再均匀提取27个点作为外部检核点，剩下的509点进行区域格网模拟计算，具体方案如下：
试验1(全国格网模拟计算)：用269个点建立0.3°的父格网，内插27个外部检核点，统计转换精度。试验2(区域格网模拟计算)：在试验1的基础上，以509个点按拟合推估两步极小解法建立0.1°的区域格网，由区域格网内插27个外部检核点，统计转换精度。
计算中，协方差函数采用Reilly函数
式中，A为待定常数；σ2(0)为观测值方差，认为噪声Δ=0，则σ2(0)等于各点残差的平方和均值；σ(d)为两点间的协方差；d为两点间的距离。
假如区域格网节点g周围有n个数据点，则该点的协方差推估方程为
根据周围数据点的残差值vi(i=1,2，…,n)便可推估出区域格网节点g的残差值vg。
试算结果见表 5。由表可见，区域格网不仅提高了全国格网的分辨率，也提高了转换精度。由于本算例采用的是同一控制网的数据，数据间的内洽性较好，故对转换精度的提高不太明显。实际应用中，由于城市控制网与国家控制网是分别布设和计算的，故存在区域性系统偏差，其将体现在区域格网模型中，从而较好地提高转换精度。
表 5 区域格网试算外部检核精度Tab. 5 The external accuracy in regional grid calculationm
| 平均值 | 标准差 | |v| 最大值 |
| B分量 | L分量 | B分量 | L分量 | B分量 | L分量 |
| 试验1 | -0.122 | -0.093 | 0.059 | 0.029 | 0.283 | 0.146 |
| 试验2 | -0.045 | -0.029 | 0.049 | 0.027 | 0.176 | 0.068 |
表选项
目前我国已建立了2.5′×2.5′的全国范围格网模型[23, 24]，对于有更高转换精度要求的城市地区可在此基础上按上述方法加密区域格网。区域格网数据集作为全国格网的子集存储在各地分布式数据库中，如果应用范围位于区域格网内，则可根据精度要求调用该区域格网模型。
6 结 论
(1) 我国BJS54、XAS80坐标系成果向CGCS2000的转换，按不同用户需求分为高、中、低精度3类，其中，高精度转换模型的转换精度为分米级，采用数学模型与数字模型相结合的格网模型构建。格网模型可以拟合大地控制网建立时的各种误差，满足制图等高精度转换的应用需求。
(2) 以最小曲率法构建格网模型，可以同时具备整体光滑和局部光滑特性，更好地符合地图线状要素坐标转换时的连续光滑性。
(3) 格网节点最小曲率法迭代计算，先要通过试验确定最优格网间距，然后依次细化粗格网直到最优格网。格网模型的数据结构可简化为一维向量序列组，按分块层次结构建立空间索引。格网模型可采用多级格网结构，不同级格网间距不同，格网越小，所表达的连续坐标分量差值的分辨率越高。
(4) 全国范围CGCS2000坐标转换最优格网间距为2.5′×2.5′。按48433个2000国家大地控制网点坐标构建了BJS54和XAS80到CGCS2000坐标转换格网模型，转换后的点位精度分别为0.26m和0.03m，已应用于全国1∶25万和1∶5万地形图更新生产。
(5) 在全国格网模型的基础上，根据各城市公共点的密度情况构建区域格网，以提高格网模型的分辨率，获得优于分米级的转换精度，满足1∶1000以上比例尺地形图的坐标转换精度要求。区域格网以全国格网为基础采用拟合推估两步极小法构建。
【引文格式】吕志平，魏子卿，李军，等。GCS2000高精度坐标转换格网模型的建立[J].测绘学报，2013，42(6)：791-797.
精
彩
回
顾
论文推荐| 张良培：基于分裂合并的多模型拟合方法在点云分割中的应用
李德仁当选国际宇航科学院院士，集5院院士于一身！
论文推荐| 邸凯昌：视觉SLAM技术的进展与应用
院士论坛| 李德仁：遥感双院士的中国梦
论文推荐| 王密：高分辨率光学卫星影像高精度在轨实时云检测的流式计算
《测绘学报》编委金双根当选欧洲科学院院士
学术前沿| 唐新明：雷达卫星自动成图的精密干涉测量关键技术
《测绘学报》 “数字摄影测量与机器视觉专辑”在CPGIS2018北京论坛发布
听李德仁、杨元喜、龚健雅三位院士讲述“我的科研故事”
机器视觉| 晏磊：航空遥感平台通用物理模型及可变基高比系统精度评价
权威 | 专业 | 学术 | 前沿
微信投稿邮箱 | song_qi_fan@163.com
微信公众号中搜索「测绘学报」，关注我们，长按上图二维码，关注学术前沿动态。
进群请备注：姓名+单位+稿件编号