海量遙感影像分塊逆映射快速坐標轉換

范 鈾，伍超云，王 琨，范 沖

(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