四参数仿射变换在大数据整合中的实践与探讨

吴厚清　蔡以雷

（浙江省测绘资料档案馆　浙江杭州　310000）

四参数仿射变换在大数据整合中的实践与探讨

吴厚清蔡以雷

（浙江省测绘资料档案馆浙江杭州310000）

摘要：根据布尔莎七参数模型转换原理，灵活运用四参数模型对未知坐标系的空间数据进行了转换纠正，采用迭代方式保证其转换成果的精度，能满足测绘与地理信息部门在大数据时代对获取的数据进行整合利用的需求，并为未知坐标系的数据整合提供新的思路。

关键词：数据整合；布尔莎七参数；坐标转换

1　引言

在大数据时代，测绘与地理信息部门即是数据的生产者，也是数据的整合者。在工作实践中，往往要探讨对不同部门和行业获取的数据进行整合，数据整合的前提是对不同的数据进行统一的坐标定位，也即必须进行坐标转换。

2　常用坐标系统

来自不同部门和行业的数据，通常具有不同的坐标系，基本的有大地坐标系和平面直角坐标系二大类，而大地坐标系又分地心坐标系和参心坐标系，采用不同的椭球几何参数。常见的有1984世界大地坐标系 （即WGS-84坐标系）、1954年北京坐标系、1980西安坐标系、2000国家大地坐标系等。

目前世界各国采用最广泛的高斯-克吕格投影和墨卡托投影（utm）均是正形投影（等角投影），即该投影在小区域范围内使平面图形与椭球面上的图形保持相似。为了限制长度变形，根据国际测量协会规定，将全球按一定经差分成若干带。我国采用6度分带或3度分带。

高斯平面直角坐标系一般以中央经线投影为纵轴x，赤道投影为横轴y，两轴交点即为各带的坐标原点。为了避免横坐标出现负值，在投影中规定将坐标纵轴西移500公里当作起始轴。为了区别某一坐标系统属于哪一带，通常在横轴坐标前加上带号，如（4231898m，21655933m），其中21即为带号。城建坐标多采用三度带的高斯-克吕格投影。同一坐标系下的大地坐标（即经纬度坐标b，l）与其对应的高斯平面直角坐标（x，y）有严格的转换关系。

3　坐标系统转换

关于坐标转换，首先要搞清楚转换的严密性问题，即在同一个椭球里的坐标转换都是相对严密的，而在不同的椭球之间的转换是相对近似的。例如，由1954北京坐标系的大地坐标转换到1954北京坐标系的高斯平面直角坐标是在同一参考椭球体范畴内的坐标转换，其转换过程是严密的。由WGS-84的大地坐标转换到1954北京坐标系，就属于不同椭球体间的转换。

不同椭球体间的坐标转换在局部地区的采用的常用办法是相似变换法，即利用部分分布相对合理高等级公共点求出相应的转换参数。一般而言，比较严密的是用布尔莎七参数的相似变换法，即x平移，y平移，z平移，x旋转，y旋转，z旋转，尺度变化k。要求得七参数就需要在一个地区需要3个以上的已知点，如果区域范围不大，最远点间的距离不大于30km（经验值），这可以用三参数，即x平移，y平移，z平移，而将x旋转，y旋转，z旋转，尺度变化k视为0，所以三参数只是七参数的一种特例。很多研究表明，使用七参数模型换算时，大地高对于其平面坐标换算的结果即平面坐标（x，y）影响不大；如果不考虑高程的影响，对于不同椭球体下的高斯平面直角坐标可采用四参数的相似变换法，即四参数（x平移，y平移，尺度变化m，旋转角度α）。如果用户要求的精度低于20米，在一定范围（2＇*2＇）内，就直接可以用二参数法（δb，δl）或（δx，δy）修正。但在实际操作中，这也取决于选取的公共点是否合理，并保证其足够的精度。

4四　参数仿射变换实例

对未知坐标系数据的转换，在无需考虑高程的影响下，可采用四参数仿射变换，即四参数（x平移，y平移，尺度变化k，旋转角度α）来转换未知坐标系统的数据，同时可适当增加选取公共点，来保证其足够的精度。在浙江省地理空间数据交换和共享平台数据整合实践中，获取了浙江省统计局的第三次经济普查工作中产生的浙江省范围内的建筑物地址数据（不属于2000坐标系数据），这些地址数据都存在相应的普查区工作边界，而普查区工作边界与浙江省县及乡镇级以上的行政区域境界存在高度的相关性，利用边界之间的相关性来获取全省不同区域多个同名控制点，在实际操作中，为提高转换精度，将浙江省区域范围按第三次经济普查工作划分的普查工作区为基本单元，运用四参数仿射变换并采用迭代方式，分别求得每个基本单元的转换参数，就可以得到我们需要的2000坐标系成果。

4.1技术路线

由于矢量数据都是由点组成，只要利用点数据就可以根据四参数的仿射变换的计算公式进行转换。其中最关键的是找出同名控制点在两个不同坐标系中的不同坐标。四参数的仿射变换的计算公式如下：x1=a*x0+b*y0+c，y1=m*x0+n* y0+d，其中（a,b,c,m,n,d）就是需要求的参数，其实中间是有内在关系可以省略掉两个参数的，最后需要求一个旋转角度、两个偏移参数，再加一个尺度变化。

4.2程序实现

4.2.1获取同名控制点

首先利用普查区工作边界跟浙江省县及乡镇级以上的境界的高度相关性的特点，各自取得界线相交点作为同名控制点在不同坐标系中的坐标数据组，控制点部分截图示例如图1

4.2.2获取旋转角度

根据对角两个控制点数据获取旋转角度参数。

＇〈summary〉

＇获取旋转角度

＇〈/summary〉

＇〈paramname="fromCoordPoint1"〉源点 1〈/ param〉

＇〈paramname="toCoordPoint1"〉目标点 1〈/ param〉

＇〈paramname="fromCoordPoint2"〉源点 2〈/ param〉

＇〈paramname="toCoordPoint2"〉目标点 2〈/ param〉

＇〈returns〉返回旋转角度〈/returns〉

PrivateFunctionGetRotation（fromPoint1As IPoint,toPoint1AsIPoint,fromPoint2AsIPoint, toPoint2AsIPoint）AsDouble

Dima,bAsDouble

a=（toPoint2.y-toPoint1.y）*（fromPoint2.xfromPoint1.x）-（toPoint2.x-toPoint1.x）*（fromPoint2.y-fromPoint1.y）

b=（toPoint2.x-toPoint1.x）*（fromPoint2.xfromPoint1.x）+（toPoint2.y-toPoint1.y）*（fromPoint2.y-fromPoint1.y）

If（Math.Abs（b）〉0）Then

GetRotation=Math.Tan（a/b）

Else

GetRotation=Math.Tan（0）

EndIf

EndFunction

4.2.3获取缩放比例

＇〈summary〉

＇根据对角两个控制点数据和旋转角度获取尺度变化参数。

＇〈/summary〉

＇〈paramname="fromCoordPoint1"〉源点 1〈/ param〉

＇〈paramname="toCoordPoint1"〉目标点 1〈/ param〉

＇〈paramname="fromCoordPoint2"〉源点 2〈/ param〉

＇〈paramname="toCoordPoint2"〉目标点 2〈/ param〉

＇〈paramname="rotation"〉旋转角度〈/param〉

＇〈returns〉返回缩放比例〈/returns〉

PrivateFunctionGetScale（fromPoint1AsIPoint, toPoint1AsIPoint,fromPoint2AsIPoint,toPoint2As IPoint,rotationAsDouble）AsDouble

Dima,bAsDouble

a=toPoint2.x-toPoint1.x

b=（fromPoint2.x-fromPoint1.x）*Math.Cos（rotation）-（fromPoint2.y-fromPoint1.y）*Math.Sin（rotation）

If（Math.Abs（b）〉0）Then

GetScale=a/b

Else

GetScale=0

EndIf

EndFunction

4.2.4获取X、Y偏移量

＇〈summary〉

＇根据一个控制点数据、旋转角度和尺度变化获取X、Y两个方向的偏移量参数。

＇〈/summary〉

＇〈paramname="fromCoordPoint1"〉源点 1〈/ param〉

＇〈paramname="toCoordPoint1"〉目标点 1〈/ param〉

＇〈paramname="rotation"〉旋转角度〈/param〉

＇〈paramname="scale"〉尺度变化〈/param〉

＇〈returns〉返回X方向偏移量〈/returns〉

PrivateFunctionGetXTranslation（fromPoint1As IPoint,toPoint1AsIPoint,rotationAsDouble,scale1 AsDouble）AsDouble

GetXTranslation=（toPoint1.x-scale1*（fromPoint1.x*Math.Cos（rotation）-fromPoint1.y* Math.Sin（rotation）））

EndFunction

＇〈summary〉

＇得到Y方向偏移量

＇〈/summary〉

＇〈paramname="fromCoordPoint1"〉源点 1〈/ param〉

＇〈paramname="toCoordPoint1"〉目标点 1〈/ param〉

＇〈paramname="rotation"〉旋转角度〈/param〉

＇〈paramname="scale"〉尺度变化〈/param〉

＇〈returns〉返回Y方向偏移量〈/returns〉

PrivateFunctionGetYTranslation（fromPoint1As IPoint,toPoint1AsIPoint,rotationAsDouble,scale1 AsDouble）AsDouble

GetYTranslation=（toPoint1.y-scale1*（fromPoint1.x*Math.Sin（rotation）+fromPoint1.y* Math.Cos（rotation）））

EndFunction

4.2.5调用四参数变换模型

＇说明：采用相似变换模型（四参数变换模型）

＇X=ax-by+c

＇Y=bx+ay+d

PrivateFunctionUnaffineTransformation（pPoint AsIPoint）AsIPoint

DimaAsDouble,bAsDouble

DimtPointAsIPoint

GetFourParampPoint＇得到四参数

a=m_Scale*Math.Cos（m_RotationAngle）b=m_Scale*Math.Sin（m_RotationAngle）SettPoint=NewPoint

tPoint.x=a*pPoint.x-b*pPoint.y+m_DX tPoint.y=b*pPoint.x+a*pPoint.y+m_DY SetUnaffineTransformation=tPoint EndFunction

4.3坐标转换的成果

通过对浙江省统计局全省普查边界数据的转换，进而获取普查工作区范围内的2000坐标系的建筑物地址数据。经检测，经坐标转换后的数据平均误差在3米以内，满足同比例尺地形图的精度要求。如下：图2为转换前局部数据截图，其平均误差在500多米，图3为转换后局部数据截图，平均误差在3米以内。

5　结论

经过浙江省统计局的第三次经济普查工作产生的数据转换项目，我们可以更深的理解不同坐标系之间的转换其实是两种不同的椭球参数之间的转换。其中七参数的布尔莎模型，是比较严密的转换模型，不管是四参数、三参数、二参数都是其特定情况的简化转换公式。在转换未知坐标系数据过程中，灵活套用布尔莎模型的简化转换公式，是一种比较经济实用的选择。

参考文献：

［1］陈俊勇.《测绘通报》2003年02期《世界大地坐标系1984的最新精化》

［2］程鹏飞，文汉江，成英燕，王华.《测绘学报》2009年03期《2000国家大地坐标系椭球参数与GRS80和WGS84的比较》

［3］徐登云，郝丽娟.《西部资源》2012年02期 《2000国家大地坐标系与GRS80及WGS84的比较》