基于ArcGIS与C#语言的省级历史遥感影像数据坐标转换应用研究

韩军 于亚杰 郭玮

摘 要：充分整合历史遥感影像数据，挖掘其服务价值对自然资源管理、经济社会发展具有十分重要的意义。本文以河北省自然资源厅2022—2023年度“全省历史遥感影像开放”项目为例，对河北省近10年遥感影像成果进行收集分析，利用布尔莎七参数模型，采用C#语言编制七参数计算工具，以ArcGIS10.8作为数据处理平台，完成全省原有参心坐标系影像数据向2000国家大地坐标系转换，评定其转换精度，为全省历史遥感影像数据融合及应用系统建设统一了数据基准。

关键词：遥感影像；ArcGIS；C#；CGCS2000；坐标转换；布尔莎七参数

在科技革新与社会效应的驱动下，地理信息数据日益丰富，应用前景广泛［1］。为进一步提升测绘地理信息数据的服务价值，促进资源共享，河北省自然资源厅启动了“全省历史遥感影像开放”项目。由于历史原因，历年遥感影像数据坐标系不统一，无法保障整个项目质量。本文利用C#语言开发坐标转换软件，通过主流的测绘地理信息数据处理平台，实现大量多源遥感影像数据高效、快速、便捷地完成坐标转换。

1 坐标转换模型

本项目作业范围覆盖整个河北省，面积约为188800km2，最终成果为CGCS2000下大地坐标（BLH）。数据源主要有两种情况：一是源坐标系均为2000国家大地坐标系，但中央子午线不同，可采用高斯反算实现平面坐标的向大地坐标的转换；二是数据源为不同椭球基准，则采用布尔沙七参数模型实现。该模型为三维模型，在空间直角坐标系中，两坐标系之间存在严密的转换模型，不存在模型误差和投影变形误差，适用于椭球面3°及以上的省级及全国范围控制点坐标转换［24］。

布尔沙七参数模型包括3个平移参数（X0，Y0，Z0）、3个旋转参数（εX、εY、εz，也被称为3个欧勒角）和1个尺度参数（m）。本次采用布尔莎模型七参数法实现原参心坐标系到CGCS2000的转换。如式（1）、式（2）所示［56］：

布尔莎七参数公式为：

X2Y2Z2=（1+m）X1Y1Z1+0εZ-εY-εZ0εXεY-εX0X1Y1Z1+X0Y0Z0（1）

一般展开式形式是：

X2=X0+（1+m）X1+εZY1-εYZ1Y2=Y0+（1+m）Y1-εZX1+εXZ1Z2=Z0+（1+m）Z1+εYX1-εXY1（2）

上述公式中，X0、Y0、Z0、εX、εY、εZ、m即为七参数，（X1，Y1，Z1）代表原有参心坐标系成果，（X2，Y2，Z2）代表CGCS2000坐標系下成果。转换参数求取采用最小二乘法，至少需要3个公共点。当公共点个数较多时，可根据测量平差原理列出观测值的误差方程式，组成并解算法方程，求得转换参数并评估精度。

2 已有数据资料

2.1 遥感影像数据

本项目主要数据来源为历年全省地理国情监测、1∶10000基础测绘更新、全省航空摄影遥感资料获取与处理等项目数据成果，具体情况如下表所示。

2.2 控制点数据

本项目选取了河北省自然资源档案馆提供的均匀覆盖全省的国家A级、B级及省C级GPS点共计238个，作为坐标转换参数计算依据。通过附合性检查，精度良好。

3 技术路线与方案

3.1 坐标转换流程

河北省横跨38、39、40共3个分带，因此涉及不同类型的坐标转换。由于分区的原因，不同子午线邻带接边处易出现转换盲区，以分带38（114°）向39（117°）转换为例，首先对38（114°）的影像进行坐标转换，同时在两分带接边处，选取39（117°）的影像换带至38（114°），与中央子午线为117°的影像一并进行坐标转换，从而消除接边漏洞。由于数据源存在多种分辨率，接边处存在不满幅的问题，需要采用数据融合方法解决。以分辨率0.8m、1.0m为例，首先将两套影像数据按本文所述流程分别进行坐标转换，然后将1.0m影像融合至0.8m影像，完成坐标转换，对转换后的影像进行坐标精度、接边精度等方面检查，合格后方可开展下一步工序。

针对不同椭球转换，以1980西安坐标系向CGCS2000坐标系转换为例（1954年北京坐标系转换相同），主要流程如下：

（1）利用选取的均匀覆盖全省的国家、省级高等级GPS点共计238个参与转换参数计算与检核，其中选定200个点为求参点，38个点为检核点。

（2）分别对控制点1980西安坐标系、CGCS2000坐标系数据进行高斯反算，得到基于相应椭球下的大地坐标BLH，通过数学模型将转换至空间直角坐标XYZ，然后采用布尔沙七参数模型，基于最小二乘法求取转换参数，并计算重合点坐标残差。

（3）利用计算得到的七参数，将待转换数据全部转换至CGCS2000坐标系下。

3.2 采用C#开发坐标转换软件

利用BursaWolf七参数转换模型，通过QR分解解决了矩阵求逆过程中出现的数值不稳定问题，按照最小二乘法基于Visual Studio 2015平台实现坐标系之间的转换。软件界面如下图所示。

主要代码如下：

private void button3_Click（object sender，EventArgs e）

{

double［，］S=new double［7，1］；

if（newX ！=null）

{

Matrix Matrix=new Matrix（）；

double［，］B=new double［3*length，7］；for（int i=0；i<3*length；i++）for（int j=0；j<7；j++）B［i，j］=0；

double［，］l=new double［3*length，1］；

for（int i=0；i

{

B［3*i，0］=1；B［3*i+1，1］=1；B［3*i+2，2］=1；

B［3*i，3］=oldX［i，0］；B［3*i+1，3］=oldX［i，1］；B［3*i+2，3］=oldX［i，2］；

B［3*i，5］=-oldX［i，2］；B［3*i+1，6］=-oldX［i，0］；B［3*i+2，4］=-oldX［i，1］；

B［3*i，6］=oldX［i，1］；B［3*i+1，4］=oldX［i，2］；B［3*i+2，5］=oldX［i，0］；

l［3*i，0］=newX［i，0］；l［3*i+1，0］=newX［i，1］；l［3*i+2，0］=newX［i，2］；

}

para=Matrix.MultiplyMatrix（Matrix.Athwart（Matrix.MultiplyMatrix（Matrix.Transpose（B），B）），Matrix.MultiplyMatrix（Matrix.Transpose（B），l））；S=para；

xtextBox.Text=para［0，0］.ToString（）；ytextBox.Text=para［1，0］.ToString（）；ztextBox.Text=para［2，0］.ToString（）；

mtextBox.Text=Convert.ToString（（para［3，0］-1）*1000000）；

extextBox.Text=Convert.ToString（para［4，0］/para［3，0］*206265）；

eytextBox.Text=Convert.ToString（para［5，0］/para［3，0］*206265）；

eztextBox.Text=Convert.ToString（para［6，0］/para［3，0］*206265）；

}

本軟件支持空间直角坐标系与大地坐标系互换，通过源坐标、目标坐标的空间直角坐标系计算布尔沙七参数，可实现坐标文件的批量转换。

3.3 基于ArcGIS平台实现数据坐标转换

坐标投影转换坐标投影转换实际上就是将源坐标投影信息删除，然后重新给其定义与目标投影体系相同的坐标系统［7］。首先，启动ArcGIS Desktop平台，添加待转换遥感影像数据，由于影像数据量较大，因此需要选择创建金字塔。其次，在Arc Toolbox中使用数据管理工具的投影和变换工具创建自定义地理转换（Create Custom Geographic Transformation），设置地理变换名称（Geographic Transformation Name），依次输入源地理坐标系（Import Geographic Coordinate System），输出目标地理坐标系（Output Geographic Coordinate System）和自定义地理变换方法（Custom Geographic Transformation），设置3.2中计算出的七参数。最后，选择新建的地理（坐标）转换方式，进行投影转换。

4 精度要求与评价

4.1 精度要求

本项目最终目的是将遥感影像数据制作电子地图，通过瓦片形式实现公众服务。依据规范要求，遥感影像平面位置中误差规定为：平地、丘陵地区不大于0.5mm（图上），山地、高山地不大于0.75mm（图上），以明显地物点平面位置中误差的2倍为其最大误差。影像拼接数字影像图应与相邻影像图接边，接边误差不大于2个像元。

4.2 精度评价

坐标转换精度采用内符合和外符合精度评价，依据计算转换参数的重合点残差中的误差评估坐标转换精度，残差小于3倍点位中误差的点位精度满足要求。经计算，本项目内符合转换各点残差最小为：±0.002m，最大为±0.091m，内符合中误差为±0.032m；外符合转换各点残差最小为：±0.012m，最大为±0.131m，外符合中误差为±0.058m，完全满足本项目精度要求。

将新得到的数据与目标数据、全省行政界线、基础线划图等参考图层进行叠加，发现两者数据吻合情况良好，没有发现变形与偏移；与源数据相比，属性信息没有发生任何丢失，满足了质量要求。

结语

本文通过结合项目实例，采用ArcGIS与C#相结合，实现了省级历史遥感影像数据坐标转换的方法和流程，建立了严密、统一的坐标转换关系模型，在保证成果质量的情况下，提高了工作效率。经检验，本次工作所采取的转换方法与流程切实可行，对大范围地理信息数据坐标转换问题提供了良好的解决方案。

参考文献：

［1］贾继鹏，厉芳婷，侯爱羚.基于多源数据融合的基础地理信息数据动态更新［J］.地理空间信息，2021，44（01）：4143+4.

［2］徐仕琪，张晓帆，周可法，等.关于利用七参数法进行WGS84和BJ54坐标转换问题的探讨［J］.测绘与空间地理信息，2007，30（5）：3338.

［3］谢艳玲，夏正清.基于ArcGIS的1980西安坐标系到2000国家坐标系的转换研究［J］.测绘与空间地理信息，2014，37（12）：220224.

［4］成英燕，程鹏飞，秘金钟，等.大尺度空间域下1980西安坐标系与WGS84坐标系转换方法研究［J］.测绘通报，2007（12）：58.

［5］孔祥元，郭际明，刘宗泉.大地测量学基础［M］.武汉：武汉大学出版社，2010.

［6］郭英起，唐彬，张秋江，等.基于空间直角坐标系的高精度坐标转换方法研究［J］.大地测量与地球动力学，2012，32（3）：125128.

［7］谭玲，秦龙，毛莉莉.基于ArcGIS的DOM坐标系统转换［J］.北京测绘，2017（3）：111113.

基金项目：本论文为河北省自然资源厅2023年度财政项目，立项编号：13000023P006CA410209K

作者简介：韩军（1976— ），男，汉族，河北石家庄人，本科，高级工程师，航空摄影室主任，主要从事测绘地理信息数据采集处理及系统研发工作。