影像坐标转换的一体化处理研究

张晓东，吴正鹏，陈楚，王琳

(天津市测绘院，天津 300381)

影像坐标转换的一体化处理研究

张晓东∗，吴正鹏，陈楚，王琳

(天津市测绘院，天津 300381)

实际工作中，常存在各种坐标系相互转换的需求。在传统数据生产流程中，产生了过多不必要的中间成果，同时影响数据处理效率。本文针对天津市特定情况，为提高数据处理效率，减少对硬盘空间的需求，研究一种新的影像数据坐标转换方法，实现一体化处理。

坐标转换；多项式纠正；一体化

1 引 言

天津市目前常用的坐标系有1990年天津市任意直角坐标系(以下简称90系)、1980西安坐标系(以下简称80系)、2000国家大地坐标系、WGS-84坐标系，以满足规划、国土、测绘等不同部门的需要。其中90系是为满足经济的迅猛发展和城市基本建设的需要而形成的自己一套独立坐标系统[1]。在实际生产工作中，常存在90系转80系，80系转2 000国家大地坐标系等需求。对于遥感影像而言，数据量大，在坐标转换过程中带来硬盘空间不足，效率低下的问题。针对上述生产中的问题及天津市特定情况，本文以90系(1∶2 000)转80系(1∶2 000)为例，介绍影像坐标转换的一体化快速实现方法。

2 传统坐标转换方法

天津市地理位置约为东经116°42′～118°03′，其中央经线为117°，坐标转换过程中不存在3°带或6°带跨带问题。因各坐标系参数严格保密，为便于常规生产，一般是先获取不同坐标系之间的数量相当的同名点对，根据多项式纠正实现影像坐标转换。同名点对的获取方法如下:分别在X方向和Y方向上每隔一定距离均匀生成标准点位，根据这些标准点位，从当地测绘主管部门生成不同坐标系之间的对应坐标点对。

以90系转80系为例，目前各作业部门通常做法如下:

(1)在ERDAS软件中拼接90系影像；

(2)将拼接后的90系影像采用二次多项式方法纠正到80系；

(3)对纠正后的80系影像分幅编号。

在以上步骤中，存在以下问题:若数据量过大，拼接后的90系影像和纠正得到的80系影像占用大量磁盘空间；拼接影像耗时过长，且这一过程属于不必要时间；采用分步操作，作业员需要人工干预过多，不符合一体化数据处理思想。

3 一体化坐标转换方法

因90系和80系具有标准的分幅编号，根据地理坐标可以找到相应的图幅编号，根据图幅编号同样可以找到对应像素的地理坐标，两种坐标系1∶2 000的分幅方式定义如下:

90系1∶2 000分幅方法:采用40 cm×50 cm矩形分幅，图号采用天津市地形图系列编号办法，即每幅1∶10 000图分为25幅1∶2 000图，编号在1∶10 000图幅号后由1～25顺序编排，如300-95-2[1]。

80系1∶2 000分幅方法:采用50 cm×50 cm分幅方式，四角点按整千米数进行分幅。地形图编号采用图廓西南角坐标千米数编号法，北方向坐标在前，东方向坐标在后，取至1 km(如4345-543)，如图1所示。

图1 分幅方式说明

提出一体化坐标转换思路方法后，通常想到的做法如图2所示，即对于给定的80系影像编号，遍历其每一个像素，根据其地理坐标及坐标转换二次多项式系数计算相应的90系坐标，判断该点位于哪一幅90系影像，打开该影像并通过双线性内插计算该点的像素值。

图2 坐标转换的通常方法

上述方法存在问题如下:在图2中第三步计算该坐标位于90系影像的分幅编号名称，并打开对应影像，对于每一个像素都需要进行上述操作，即频繁打开并关闭影像，严重影响数据处理效率。对于得到一副1∶2 000图80系影像，其计算时间长达半个多小时。为减少影像打开次数，提高数据处理效率，其改进后的坐标转换一体化数据处理流程如图3所示。

图3 影像坐标转换一体化数据处理流程

首先，预先建立与影像像幅大小一样的标志数组，并初始化为0；逐一遍历80系影像，判断该处标志数组是否为0，如果为0则根据其地理坐标计算相应的90系影像坐标，并计算该坐标所处的90系影像分幅编号名称，打开该90系影像；为防止频繁打开该90系影像，需要找到与该90系影像相关的所有像素，因此，再次遍历80系影像的每一个像素，如果标志数组相应位置为0，则根据其地理坐标计算相应的90系坐标，如果该坐标位于已经打开的90系影像，则通过双线性内插[2]计算该处像素值，并将标志数组相应位置置为1。上述流程涉及关键技术如下:

(1)影像坐标转换的仿射变换系数计算

因为天津市面积相对较小，坐标转换过程中不存在跨带问题，因此80系和90系之间的转换可以近似看作为平移、缩放、旋转、仿射、偏扭、弯曲等基本形变的合成，可以用以下多项式表达[3，4]。

上式多项式系数的计算根据一定数量的同名点对通过最小二乘法[5]计算得到。

(2)90系图幅号与坐标之间的转换①已知图幅号计算其左下角坐标

已知90系1∶2 000图的编号为y10000-x10000-index，则其所在图幅的左下角坐标计算如下: xGeo=x10000+[(index-1)%5]×1000；yGeo=y10000+[4-(index-1)/5]×800；②已知坐标计算其所在图幅号

已知某点地理坐标为(xGeo，yGeo)，求在90系1∶2 000图的编号y10000-x10000-index，用伪代码表示过程如下:

index=int(xGeo/1 000)；

x10 000=index-index%5；

index=int(yGeo/1 000)；

y10 000=index-index%4；

int yID=5-int((yGeo-y10 000∗1 000))/800；

int xID=int((xGeo-x10 000∗1 000))/1 000+1；

index=(yID-1)∗5+xID；

4 结 语

本文对上述流程总结并通过程序实现，目前已经应用于实际生产中。坐标转换之后的影像位置与ERDAS纠正之后的影像位置完全吻合。对于20 000幅DOM，采用常规作业方法，由于需要人工干预，且对硬盘空间需求较大，从90系转到80系需要将近1个月的时间。采用本文提出的方法，其过程是全自动化，中间不需要任何人工干预，其数据处理时间约60 h，有效减少了数据处理时间，提高了整个作业效率。

由于天津市是一个直辖市，其区域面积相对各省面积小很多，且不存在3°带或6°带跨带问题，因此采用二次多项式纠正进行坐标转换能够保证其精度。

[1] 天津市测绘院.天津市基础地理信息要素数据字典[R]. 2011.

[2] 冈萨雷斯.数字图像处理(第二版)[M].北京：电子工业出版社，2008.

[3] 张剑清，潘励，王树根.摄影测量学[M].武汉：武汉大学出版社，2004.

[4] 张祖勋，张剑清.数字摄影测量学[M].武汉：武汉大学出版社，2007.

[5] 武汉大学测绘学院.误差理论与测量平差基础[M].武汉：武汉大学出版社，2003.

The Research on Integrated Processing of Coordinate Transformation

Zhang Xiaodong，Wu Zhengpeng，Chen Chu，Wang Lin
(Tianjin Institute of Surveying and Mapping，Tianjin 100381，China)

In the practicalwork，there is a demand for the transformation of different coordinate systems..However，there are some redundant results in the traditional data processing which affects the processing efficiency.Aiming at the specific situation of Tianjin，the integrated processing of coordinate transformation is put forward in order to improve the efficiency and reduce the demand for the hard disk space.

Coordinate transform；Polynomial rectification；Integration

1672-8262(2013)02-116-02

P226+.3

A

2012—07—11

张晓东(1968—)，男，高级工程师，主要从事有关航空航天摄影测量工作。