海量遥感影像分块逆映射快速坐标转换

范 铀，伍超云，王 琨，范 冲

(1. 广东南方数码科技股份有限公司，北京 100055； 2. 中南大学地球科学与信息物理学院，湖南 长沙 410083)

海量遥感影像分块逆映射快速坐标转换

范 铀1，伍超云2，王 琨2，范 冲2

(1. 广东南方数码科技股份有限公司，北京 100055； 2. 中南大学地球科学与信息物理学院，湖南 长沙 410083)

提出了一种海量遥感影像分块逆映射快速坐标转换算法，并针对分块窗口的形状和大小对处理过程所占内存与处理效率造成的影响作了试验与分析并得出结论。该算法使得处理过程所占内存不受原影像与处理后影像的尺寸大小所限制，并且通过与主流遥感影像处理软件对比试验的结果可以看出，该算法可以完成超大幅遥感影像与分幅遥感影像的坐标转换处理，处理分幅影像时有效地避免了影像间的接缝，同时处理效率达到主流的专业遥感影像处理软件水平。

逆映射；分块；大幅影像；分幅影像；坐标转换

2008年7月1日起，我国正式启用CGCS2000。随着CGCS2000逐渐取代现有的国家参心坐标系[1-2]，不同地区可能拥有不同坐标系下的多套成果，如1954北京、1980西安坐标系下的大量遥感影像、DEM等栅格文件等。为了将数据统一到CGCS2000[3]下，坐标转换[4]的工作尤为重要。

许多商业遥感软件(如ENVI、ArcGIS等)功能强大，能够进行坐标转换的工作[5]。但在批量处理海量遥感影像数据时，需人工干预处较多，显得过于笨重；且封装性高，二次开发比较困难，与项目结合程度不高[6]；同时在分幅影像坐标转换的特殊情况下，由于坐标转换后影像的坐标发生了改变，但分幅规则并没有变，因此需将转换后的影像按照其分幅规则进行重新分幅。传统遥感影像处理软件并没有设置影像重分幅功能，一般是先对所有影像进行坐标转换，然后将转换后的影像拼接起来，再按分幅规则进行裁剪[7-8]，因此转换过程效率普遍偏低；而且当转换过程中影像出现形变而导致转换后影像留有黑边的情况时，ENVI会存在相应的问题[9]。

本文提出一种基于逆映射分块的大幅遥感影像坐标转换方法，即先对目标影像进行分块，然后将每一块映射到原影像中进行重采样赋值。该方法对分块的形状与大小作了讨论，可在不超过坐标转换平台软硬件环境限制下提高大幅影像坐标转换效率;并且在进行分幅影像的坐标转换过程中，将重分幅过程与坐标转换同时进行，可提高坐标转换效率，有效地避免分幅影像因坐标转换产生的形变而导致的接缝问题。

1 逆映射分块算法

1.1 逆映射算法

目前的商业化遥感影像处理软件对分幅影像进行坐标转换时，一般是先对所有影像进行坐标转换，然后将转换后的影像拼接起来，再按分幅规则进行裁剪，工作量繁琐，耗时较长。本文提出的逆映射影像坐标转换算法流程如图1所示，为了提高转换效率，对目标影像进行了分块处理；在对分幅影像坐标转换的同时对影像进行了重分幅处理，省去了坐标转换后影像拼接、分幅裁剪等过程，大大减少了工作量；在对目标影像进行逆映射分块赋值过程中，在像素重采样时可将因坐标转换产生形变而导致的无值像素剔除，从而避免转换后影像出现接缝的情况。

图1 逆映射算法处理分幅影像坐标转换流程

1.2 分块窗口形状与大小的选择

1.2.1 分块窗口形状的选择

在分块处理影像时，影响处理效率与占用内存的因素主要有两个：调用影像读写接口的次数与总共读写的像素缓存大小。

数字影像的排列规则为栅格排列，因此缓存读写窗口必须为矩形。ENVI/IDL处理影像时常用的方法为利用envi_get_slice函数逐行读入数据进行处理，可以理解为长为1、宽为影像列数的块[10]。这种分块方法对影像作波段运算、平移、缩放、拼接、分割等处理时，只要保证没有重叠区域就不会出现数据冗余。但是当处理过程包含旋转、变形等过程时，赋值块映射到原图后，赋值块A与其最小外接矩形A′之间就会出现一定冗余的数据，程序读写额外的冗余数据会导致读写接口调用次数增多、数据读写总量增大及冗余像素处理耗时的情况，使影像处理效率变低，因此应尽量减少冗余数据。

在坐标换带、坐标基准转换等坐标转换过程中，会出现平移、旋转、缩放等仿射变换，以及弯曲、偏扭、膨胀效应等变形，在变形过程中会产生冗余数据[11]，并且产生冗余数据的比例与赋值块的形状有一定关系。如当矩形旋转θ角后，矩形与其外接矩形的关系如图2所示，旋转矩形与其外接矩形的面积比P的计算公式为

(1)

影像处理过程中，旋转角θ已定，因此可以将其视为常数，化简后公式为

(2)

式中，C和D均为常数。此时若设矩形的长宽比为e=a/b，变形后公式为

(3)

由其关于e的一阶偏导数可得，只有在e=1时，面积比P有最小值。由此得出，当分块窗口为正方形时，其旋转后产生的冗余数据最少。

图2 矩形旋转θ角后与其外接矩形的关系

为讨论在处理前后影像变形较明显时赋值块矩形长宽比对图像处理效率的影响，本文设计了不同长宽比的赋值块下对同一幅三波段航拍影像(2610×1908像素，14.7 MB)的坐标换带试验(本文全部试验采用的计算机系统配置为CPU：Intel Core i5 3.20 GHz；内存：4 GB；系统环境：Win7 64位操作系统)。试验设计各赋值块窗口面积均相等，但长宽比不同，以排除赋值块大小对试验的影响。试验结果见表1。随着长宽比远离1∶1，过程占用内存会缓缓增加，处理用时也会增长，由于形变量并不是很大，变化并不明显，但当长宽比为1∶1时，影像处理用时最短，占用内存最少，因此处理效率最高。

表1 不同窗口形状下内存占用与耗时情况对比

1.2.2 分块窗口大小的选取

程序处理影像时，主要的时间分配在读写影像与像素处理两部分。由于像素处理部分的效率主要由处理内容而定，因此读写影像过程是影响影像处理效率的重要因素。在影像读写过程中，主要影响效率的是调用接口读写影像的次数及总共读写的像素量大小。在不考虑数据冗余的情况下，一张影像总共读写的像素量是确定的，因此分块窗口越大，每次读写的像素量越多，调用接口读写影像的次数就越少，影像处理的效率也就越高。

但是，分块窗口的大小不仅影响影像处理的效率，同时也影响系统占用内存的大小。对矩形与其外接矩形的面积比公式作θ的一阶偏导可知，当赋值块窗口为正方形时，面积比P的最大值为2。同时由于常见的需处理的影像都由R、G、B 3个波段组成，因此处理影像过程中内存缓存最少为分块窗口面积的9倍。而计算机分配内存时由于内存分配方法的不同，以及内存碎片、内存记录等原因，实际占用的内存都大于此理论值[12]，因此分块窗口越大，运行时占用内存就越高。由试验可得：该方法在不同尺寸窗口下处理影像的速度与内存占用情况见表2。

表2 不同窗口大小下内存占用与耗时情况对比

2 逆映射分块影像快速坐标转换算法

本文提出的基于逆映射分块算法的遥感影像快速坐标转换算法的具体步骤如下：

(1) 获取原影像四至坐标(如图3(a)所示O1～O4)，根据处理过程确定目标影像的范围(如图3(b)所示目标影像)，按规定的分块窗口对其进行分块(如图3(c)所示)，并以从上到下、从左到右的顺序依次对每个块进行赋值。若右边或下边的块不足分块窗口尺寸，则按剩余像素行列重新分配块尺寸(如图3(c)所示块B、C)。

(2) 在对某块A赋值时(如图3(c)所示)，首先计算其四至坐标，计算公式为

(4)

式中，BlockTX、BlockBX为该块的上下边X坐标；BlockLY、BlockRY为左右边Y坐标，将其分别组合起来就是该块的四至坐标；ImgLTX、ImgLTY为影像左上角坐标；BlockXSize、BlockYSize为该块的像素尺寸；iBlock、jBlock为该块在所有分块中的行列序号；Kx、Ky表示像素在X、Y方向上的长度，由于空间坐标系与屏幕坐标系的关系，Kx一般为负。根据该块像素尺寸建立缓存用于赋值，称为赋值块。

(3) 获取该块的四至坐标后，将其按与处理过程相逆的过程映射到原图中，并计算出其最小外接矩形A′(如图3(d)所示)的四至坐标。

(4) 遍历所有原图，找出与该矩形有覆盖的原图，计算出覆盖区域的四至坐标，并将此区域的像素值读入缓存(如图3(d)所示)。由于读入缓存的像素块用于采样取值，因此将其称为取值块。

(5) 对赋值块A的每个像素进行重采样，对于某个像素，首先计算出该像素坐标(pX,pY)，计算公式为

(5)

式中，i、j该像素在赋值块中的行列序号。计算出像素坐标后将该坐标映射到取值块上，根据其在取值块中的位置进行重采样。最邻近法[13]效率最高，但插值出来的图像质量会有损失；三次卷积法[14]能最大限度地避免图像质量损失，但是计算量巨大；综合来看，本文采用双线性插值[15]的方法，在效率和质量上都较理想。

(6) 全部采样完毕后将该赋值块的值赋到处理后影像中，并清理所有缓存，然后对下一块重复步骤(2)—步骤(5)，对其进行赋值。所有块赋值完毕后，影像处理过程结束。

3 试验与讨论

为了验证本文介绍方法的高效性，设计了一组试验(数据见表3)，将本文方法与ENVI、ERDAS、ArcGIS软件同时进行坐标转换，采样方法全部采用双线性内插法，测试其影像坐标转换效率。试验分别使用ENVI、ERDAS、ArcGIS的坐标转换功能与本文方法依次对3幅影像进行1980西安到CGCS2000的坐标转换，其中ENVI和ERDAS并未定义1980西安和CGCS2000坐标系，需自定义坐标系。转换效率评价如表4、图4—图6所示。

图3 逆映射影像分块具体步骤示例

图4 中幅影像转换前后效果对比

图5 大幅影像转换前后效果对比

图6 巨幅影像转换前后效果对比

名称尺寸大小波段类型中幅影像6575×481963MB3TIF大幅影像30544×338153.6GB3IMG巨幅影像91768×17785746.5GB3IMG

试验结果表明，本文方法在遥感影像的坐标转换过程中，与ENVI、ERDAS的处理效率相当，内存占用较多，但CPU占用较少，并且在转换不同图幅的影像时表现得更加稳定，并没有因为图幅不同而产生内存占用与CPU占用的波动，转换后影像的清晰度也与原图一致，没有出现失真的情况。

表4 影像坐标转换效率对比情况

ENVI在进行遥感影像的坐标转换过程时效率较高，在处理不同图幅的影像时CPU、内存的占用也较稳定，并且内存占用较少；只是在处理45 GB的巨幅影像时只能保存为ENVI格式，转存为IMG格式时提示文件过大无法保存。

ERDAS进行坐标转换时效率较本文方法稍低，并且在处理不同图幅影像时CPU占用与内存占用都有不同幅度的波动，稳定性较差。

ArcGIS转换遥感影像时的效率相比前3种方法慢很多，并且CPU占用较高。在转换3.6 GB的大幅影像时，用ArcGIS转换后的影像并没有像素值。

4 结 语

本文方法能完成大幅(45 GB)、分幅影像的坐标转换过程，并且可以完成大幅、分幅影像的拼接、分割、重分幅等处理过程，处理后影像清晰度与原图一致，没有出现失真现象。坐标转换过程与ERDAS、ENVI的效率相近，占用内存较大，占用CPU较小，与ArcGIS相比效率要快很多，并且内存占用与CPU占用相对稳定，不会出现太大波动。同时在对分幅影像进行坐标转换处理时，省去了坐标转换后影像拼接、分幅裁剪等过程，大大提高了分幅影像坐标转换效率，并避免了转换后影像出现接缝的情况。而且本文方法逻辑简明，适应性强，适用于遥感影像批处理。

[1] 张训虎,刘晋虎,何川,等. 2000 国家大地坐标系转换常见问题分析[J]. 测绘通报, 2016 (9): 52-55.

[2] 田桂娥，宋利杰,尹利文,等.地方坐标系与CGCS2000坐标系转换方法的研究[J]. 测绘工程, 2014,23(8): 66-69.

[3] 吕志平, 魏子卿, 李军, 等. 我国 CGCS2000 高精度坐标转换格网模型的建立[J]. 测绘学报, 2013, 42(6): 791-797.

[4] 何林, 柳林涛, 许超钤, 等. 常见平面坐标系之间相互转换的方法研究——以 1954 北京坐标系, 1980 西安坐标系, 2000 国家大地坐标系之间的平面坐标相互转换为例[J]. 测绘通报, 2014 (9): 6-11.

[5] 查东平,林辉,孙华,等. 基于GDAL的遥感影像数据快速读取与显示方法的研究[J]. 中南林业科技大学学报, 2013, 33(1):58-62.

[6] 陈谷良. 基于WebGIS的城市电子地图框架设计研究[J]. 中国西部科技, 2011,9(35):32-34.

[7] 徐健梅, ArcGIS在数字正射影像图坐标转换中的应用[J]. 长江工程职业技术学院学报, 2010, 27(3):22-23.

[8] 林旭芳. 基于 ArcGIS 平台的遥感影像快速分幅方法[J]. 测绘通报, 2014 (S2): 179-181.

[9] 崔丽华,刘善军,闻彩焕,等. 基于IDL与ENVI的MODIS遥感影像镶嵌[J]. 河北理工大学学报(自然科学版), 2009, 31(2):127-130.

[10] 张晓东, 吴正鹏, 陈楚, 等. 影像坐标转换的一体化处理研究[J]. 城市勘测, 2013(2): 116-117.

[11] 徐永明.遥感二次开发语言IDL[M].北京:科学出版社，2014.

[12] 曹元大, 张建芳. Windows 下专家系统开发工具的内存管理[J]. 电脑开发与应用, 1996,9(1):12-14.

[13] ANBARJAFARI G, DEMIREL H. Image Super Resolution Based on Interpolation of Wavelet Domain High Frequency Subbands and the Spatial Domain Input Image[J]. ETRI Journal, 2010, 32(3): 390-394.

[14] KEYS R. Cubic Convolution Interpolation for Digital Image Processing[J]. IEEE Transactions on Acoustics, Speech, and Signal Processing, 1981, 29(6): 1153-1160.

[15] MALVAR H S, HE L, CUTLER R. High-quality Linear Interpolation for Demosaicing of Bayer-patterned Color Images[C]∥Acoustics, Speech, and Signal Processing, 2004. Proceedings. IEEE International Conference on.[S.l.]：IEEE, 2004.

Massive Remote Sensing Image Fast Coordinate Transformation Algorithm of the Inverse Mapping Blocking Method

FAN You1，WU Chaoyun2，WANG Kun2，FAN Chong2

(1. South Digital Technology Company, Beijing 100055， China；2. School of Geosciences and Info-physics, Central South University, Changsha 410083， China)

This paper proposes massive remote sensing image coordinate transformation algorithm of the inverse mapping blocking method, analyzes the influence from the geometric shape and size of the blocks to the memory usage and efficiency of coordinate transforming processes with experiments then draws conclusions. This algorithm makes the memory usage of coordinate transforming processes independent from the size of the original image or the transformed image. By experiments contrast with main current remote sensing image processing software, it is observed that this algorithm is able to achieve the coordinate transforming processes for huge remote sensing images and subdivided images, in the process of subdivided images coordinate transform this algorithm is able to avoid the seams between images, and the processing efficiency of this algorithm achieves the level of main current remote sensing image processing software.

inverse mapping; blocking; huge images; subdivided images; coordinate transformation

范铀，伍超云，王琨，等.海量遥感影像分块逆映射快速坐标转换[J].测绘通报，2017(4)：53-57.

10.13474/j.cnki.11-2246.2017.0119.

2016-12-28；

2017-03-13

国家973计划(2012CB719904)

范 铀(1977—)，男，工程师，主要从事遥感数据处理方面的研究。E-mail：you.fan@southgis.com

P237

A

0494-0911(2017)04-0053-05