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

(1. 廣東南方數(shù)碼科技股份有限公司，北京 100055； 2. 中南大學(xué)地球科學(xué)與信息物理學(xué)院，湖南 長沙 410083)

海量遙感影像分塊逆映射快速坐標(biāo)轉(zhuǎn)換

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

(1. 廣東南方數(shù)碼科技股份有限公司，北京 100055； 2. 中南大學(xué)地球科學(xué)與信息物理學(xué)院，湖南 長沙 410083)

提出了一種海量遙感影像分塊逆映射快速坐標(biāo)轉(zhuǎn)換算法，并針對分塊窗口的形狀和大小對處理過程所占內(nèi)存與處理效率造成的影響作了試驗與分析并得出結(jié)論。該算法使得處理過程所占內(nèi)存不受原影像與處理后影像的尺寸大小所限制，并且通過與主流遙感影像處理軟件對比試驗的結(jié)果可以看出，該算法可以完成超大幅遙感影像與分幅遙感影像的坐標(biāo)轉(zhuǎn)換處理，處理分幅影像時有效地避免了影像間的接縫，同時處理效率達(dá)到主流的專業(yè)遙感影像處理軟件水平。

逆映射；分塊；大幅影像；分幅影像；坐標(biāo)轉(zhuǎn)換

2008年7月1日起，我國正式啟用CGCS2000。隨著CGCS2000逐漸取代現(xiàn)有的國家參心坐標(biāo)系[1-2]，不同地區(qū)可能擁有不同坐標(biāo)系下的多套成果，如1954北京、1980西安坐標(biāo)系下的大量遙感影像、DEM等柵格文件等。為了將數(shù)據(jù)統(tǒng)一到CGCS2000[3]下，坐標(biāo)轉(zhuǎn)換[4]的工作尤為重要。

許多商業(yè)遙感軟件(如ENVI、ArcGIS等)功能強大，能夠進(jìn)行坐標(biāo)轉(zhuǎn)換的工作[5]。但在批量處理海量遙感影像數(shù)據(jù)時，需人工干預(yù)處較多，顯得過于笨重；且封裝性高，二次開發(fā)比較困難，與項目結(jié)合程度不高[6]；同時在分幅影像坐標(biāo)轉(zhuǎn)換的特殊情況下，由于坐標(biāo)轉(zhuǎn)換后影像的坐標(biāo)發(fā)生了改變，但分幅規(guī)則并沒有變，因此需將轉(zhuǎn)換后的影像按照其分幅規(guī)則進(jìn)行重新分幅。傳統(tǒng)遙感影像處理軟件并沒有設(shè)置影像重分幅功能，一般是先對所有影像進(jìn)行坐標(biāo)轉(zhuǎn)換，然后將轉(zhuǎn)換后的影像拼接起來，再按分幅規(guī)則進(jìn)行裁剪[7-8]，因此轉(zhuǎn)換過程效率普遍偏低；而且當(dāng)轉(zhuǎn)換過程中影像出現(xiàn)形變而導(dǎo)致轉(zhuǎn)換后影像留有黑邊的情況時，ENVI會存在相應(yīng)的問題[9]。

本文提出一種基于逆映射分塊的大幅遙感影像坐標(biāo)轉(zhuǎn)換方法，即先對目標(biāo)影像進(jìn)行分塊，然后將每一塊映射到原影像中進(jìn)行重采樣賦值。該方法對分塊的形狀與大小作了討論，可在不超過坐標(biāo)轉(zhuǎn)換平臺軟硬件環(huán)境限制下提高大幅影像坐標(biāo)轉(zhuǎn)換效率;并且在進(jìn)行分幅影像的坐標(biāo)轉(zhuǎn)換過程中，將重分幅過程與坐標(biāo)轉(zhuǎn)換同時進(jìn)行，可提高坐標(biāo)轉(zhuǎn)換效率，有效地避免分幅影像因坐標(biāo)轉(zhuǎn)換產(chǎn)生的形變而導(dǎo)致的接縫問題。

1 逆映射分塊算法

1.1 逆映射算法

目前的商業(yè)化遙感影像處理軟件對分幅影像進(jìn)行坐標(biāo)轉(zhuǎn)換時，一般是先對所有影像進(jìn)行坐標(biāo)轉(zhuǎn)換，然后將轉(zhuǎn)換后的影像拼接起來，再按分幅規(guī)則進(jìn)行裁剪，工作量繁瑣，耗時較長。本文提出的逆映射影像坐標(biāo)轉(zhuǎn)換算法流程如圖1所示，為了提高轉(zhuǎn)換效率，對目標(biāo)影像進(jìn)行了分塊處理；在對分幅影像坐標(biāo)轉(zhuǎn)換的同時對影像進(jìn)行了重分幅處理，省去了坐標(biāo)轉(zhuǎn)換后影像拼接、分幅裁剪等過程，大大減少了工作量；在對目標(biāo)影像進(jìn)行逆映射分塊賦值過程中，在像素重采樣時可將因坐標(biāo)轉(zhuǎn)換產(chǎn)生形變而導(dǎo)致的無值像素剔除，從而避免轉(zhuǎn)換后影像出現(xiàn)接縫的情況。

圖1 逆映射算法處理分幅影像坐標(biāo)轉(zhuǎn)換流程

1.2 分塊窗口形狀與大小的選擇

1.2.1 分塊窗口形狀的選擇

在分塊處理影像時，影響處理效率與占用內(nèi)存的因素主要有兩個：調(diào)用影像讀寫接口的次數(shù)與總共讀寫的像素緩存大小。

數(shù)字影像的排列規(guī)則為柵格排列，因此緩存讀寫窗口必須為矩形。ENVI/IDL處理影像時常用的方法為利用envi_get_slice函數(shù)逐行讀入數(shù)據(jù)進(jìn)行處理，可以理解為長為1、寬為影像列數(shù)的塊[10]。這種分塊方法對影像作波段運算、平移、縮放、拼接、分割等處理時，只要保證沒有重疊區(qū)域就不會出現(xiàn)數(shù)據(jù)冗余。但是當(dāng)處理過程包含旋轉(zhuǎn)、變形等過程時，賦值塊映射到原圖后，賦值塊A與其最小外接矩形A′之間就會出現(xiàn)一定冗余的數(shù)據(jù)，程序讀寫額外的冗余數(shù)據(jù)會導(dǎo)致讀寫接口調(diào)用次數(shù)增多、數(shù)據(jù)讀寫總量增大及冗余像素處理耗時的情況，使影像處理效率變低，因此應(yīng)盡量減少冗余數(shù)據(jù)。

在坐標(biāo)換帶、坐標(biāo)基準(zhǔn)轉(zhuǎn)換等坐標(biāo)轉(zhuǎn)換過程中，會出現(xiàn)平移、旋轉(zhuǎn)、縮放等仿射變換，以及彎曲、偏扭、膨脹效應(yīng)等變形，在變形過程中會產(chǎn)生冗余數(shù)據(jù)[11]，并且產(chǎn)生冗余數(shù)據(jù)的比例與賦值塊的形狀有一定關(guān)系。如當(dāng)矩形旋轉(zhuǎn)θ角后，矩形與其外接矩形的關(guān)系如圖2所示，旋轉(zhuǎn)矩形與其外接矩形的面積比P的計算公式為

(1)

影像處理過程中，旋轉(zhuǎn)角θ已定，因此可以將其視為常數(shù)，化簡后公式為

(2)

式中，C和D均為常數(shù)。此時若設(shè)矩形的長寬比為e=a/b，變形后公式為

(3)

由其關(guān)于e的一階偏導(dǎo)數(shù)可得，只有在e=1時，面積比P有最小值。由此得出，當(dāng)分塊窗口為正方形時，其旋轉(zhuǎn)后產(chǎn)生的冗余數(shù)據(jù)最少。

圖2 矩形旋轉(zhuǎn)θ角后與其外接矩形的關(guān)系

為討論在處理前后影像變形較明顯時賦值塊矩形長寬比對圖像處理效率的影響，本文設(shè)計了不同長寬比的賦值塊下對同一幅三波段航拍影像(2610×1908像素，14.7 MB)的坐標(biāo)換帶試驗(本文全部試驗采用的計算機(jī)系統(tǒng)配置為CPU：Intel Core i5 3.20 GHz；內(nèi)存：4 GB；系統(tǒng)環(huán)境：Win7 64位操作系統(tǒng))。試驗設(shè)計各賦值塊窗口面積均相等，但長寬比不同，以排除賦值塊大小對試驗的影響。試驗結(jié)果見表1。隨著長寬比遠(yuǎn)離1∶1，過程占用內(nèi)存會緩緩增加，處理用時也會增長，由于形變量并不是很大，變化并不明顯，但當(dāng)長寬比為1∶1時，影像處理用時最短，占用內(nèi)存最少，因此處理效率最高。

表1 不同窗口形狀下內(nèi)存占用與耗時情況對比

1.2.2 分塊窗口大小的選取

程序處理影像時，主要的時間分配在讀寫影像與像素處理兩部分。由于像素處理部分的效率主要由處理內(nèi)容而定，因此讀寫影像過程是影響影像處理效率的重要因素。在影像讀寫過程中，主要影響效率的是調(diào)用接口讀寫影像的次數(shù)及總共讀寫的像素量大小。在不考慮數(shù)據(jù)冗余的情況下，一張影像總共讀寫的像素量是確定的，因此分塊窗口越大，每次讀寫的像素量越多，調(diào)用接口讀寫影像的次數(shù)就越少，影像處理的效率也就越高。

但是，分塊窗口的大小不僅影響影像處理的效率，同時也影響系統(tǒng)占用內(nèi)存的大小。對矩形與其外接矩形的面積比公式作θ的一階偏導(dǎo)可知，當(dāng)賦值塊窗口為正方形時，面積比P的最大值為2。同時由于常見的需處理的影像都由R、G、B 3個波段組成，因此處理影像過程中內(nèi)存緩存最少為分塊窗口面積的9倍。而計算機(jī)分配內(nèi)存時由于內(nèi)存分配方法的不同，以及內(nèi)存碎片、內(nèi)存記錄等原因，實際占用的內(nèi)存都大于此理論值[12]，因此分塊窗口越大，運行時占用內(nèi)存就越高。由試驗可得：該方法在不同尺寸窗口下處理影像的速度與內(nèi)存占用情況見表2。

表2 不同窗口大小下內(nèi)存占用與耗時情況對比

2 逆映射分塊影像快速坐標(biāo)轉(zhuǎn)換算法

本文提出的基于逆映射分塊算法的遙感影像快速坐標(biāo)轉(zhuǎn)換算法的具體步驟如下：

(1) 獲取原影像四至坐標(biāo)(如圖3(a)所示O1～O4)，根據(jù)處理過程確定目標(biāo)影像的范圍(如圖3(b)所示目標(biāo)影像)，按規(guī)定的分塊窗口對其進(jìn)行分塊(如圖3(c)所示)，并以從上到下、從左到右的順序依次對每個塊進(jìn)行賦值。若右邊或下邊的塊不足分塊窗口尺寸，則按剩余像素行列重新分配塊尺寸(如圖3(c)所示塊B、C)。

(2) 在對某塊A賦值時(如圖3(c)所示)，首先計算其四至坐標(biāo)，計算公式為

(4)

式中，BlockTX、BlockBX為該塊的上下邊X坐標(biāo)；BlockLY、BlockRY為左右邊Y坐標(biāo)，將其分別組合起來就是該塊的四至坐標(biāo)；ImgLTX、ImgLTY為影像左上角坐標(biāo)；BlockXSize、BlockYSize為該塊的像素尺寸；iBlock、jBlock為該塊在所有分塊中的行列序號；Kx、Ky表示像素在X、Y方向上的長度，由于空間坐標(biāo)系與屏幕坐標(biāo)系的關(guān)系，Kx一般為負(fù)。根據(jù)該塊像素尺寸建立緩存用于賦值，稱為賦值塊。

(3) 獲取該塊的四至坐標(biāo)后，將其按與處理過程相逆的過程映射到原圖中，并計算出其最小外接矩形A′(如圖3(d)所示)的四至坐標(biāo)。

(4) 遍歷所有原圖，找出與該矩形有覆蓋的原圖，計算出覆蓋區(qū)域的四至坐標(biāo)，并將此區(qū)域的像素值讀入緩存(如圖3(d)所示)。由于讀入緩存的像素塊用于采樣取值，因此將其稱為取值塊。

(5) 對賦值塊A的每個像素進(jìn)行重采樣，對于某個像素，首先計算出該像素坐標(biāo)(pX,pY)，計算公式為

(5)

式中，i、j該像素在賦值塊中的行列序號。計算出像素坐標(biāo)后將該坐標(biāo)映射到取值塊上，根據(jù)其在取值塊中的位置進(jìn)行重采樣。最鄰近法[13]效率最高，但插值出來的圖像質(zhì)量會有損失；三次卷積法[14]能最大限度地避免圖像質(zhì)量損失，但是計算量巨大；綜合來看，本文采用雙線性插值[15]的方法，在效率和質(zhì)量上都較理想。

(6) 全部采樣完畢后將該賦值塊的值賦到處理后影像中，并清理所有緩存，然后對下一塊重復(fù)步驟(2)—步驟(5)，對其進(jìn)行賦值。所有塊賦值完畢后，影像處理過程結(jié)束。

3 試驗與討論

為了驗證本文介紹方法的高效性，設(shè)計了一組試驗(數(shù)據(jù)見表3)，將本文方法與ENVI、ERDAS、ArcGIS軟件同時進(jìn)行坐標(biāo)轉(zhuǎn)換，采樣方法全部采用雙線性內(nèi)插法，測試其影像坐標(biāo)轉(zhuǎn)換效率。試驗分別使用ENVI、ERDAS、ArcGIS的坐標(biāo)轉(zhuǎn)換功能與本文方法依次對3幅影像進(jìn)行1980西安到CGCS2000的坐標(biāo)轉(zhuǎn)換，其中ENVI和ERDAS并未定義1980西安和CGCS2000坐標(biāo)系，需自定義坐標(biāo)系。轉(zhuǎn)換效率評價如表4、圖4—圖6所示。

圖3 逆映射影像分塊具體步驟示例

圖4 中幅影像轉(zhuǎn)換前后效果對比

圖5 大幅影像轉(zhuǎn)換前后效果對比

圖6 巨幅影像轉(zhuǎn)換前后效果對比

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

試驗結(jié)果表明，本文方法在遙感影像的坐標(biāo)轉(zhuǎn)換過程中，與ENVI、ERDAS的處理效率相當(dāng)，內(nèi)存占用較多，但CPU占用較少，并且在轉(zhuǎn)換不同圖幅的影像時表現(xiàn)得更加穩(wěn)定，并沒有因為圖幅不同而產(chǎn)生內(nèi)存占用與CPU占用的波動，轉(zhuǎn)換后影像的清晰度也與原圖一致，沒有出現(xiàn)失真的情況。

表4 影像坐標(biāo)轉(zhuǎn)換效率對比情況

ENVI在進(jìn)行遙感影像的坐標(biāo)轉(zhuǎn)換過程時效率較高，在處理不同圖幅的影像時CPU、內(nèi)存的占用也較穩(wěn)定，并且內(nèi)存占用較少；只是在處理45 GB的巨幅影像時只能保存為ENVI格式，轉(zhuǎn)存為IMG格式時提示文件過大無法保存。

ERDAS進(jìn)行坐標(biāo)轉(zhuǎn)換時效率較本文方法稍低，并且在處理不同圖幅影像時CPU占用與內(nèi)存占用都有不同幅度的波動，穩(wěn)定性較差。

ArcGIS轉(zhuǎn)換遙感影像時的效率相比前3種方法慢很多，并且CPU占用較高。在轉(zhuǎn)換3.6 GB的大幅影像時，用ArcGIS轉(zhuǎn)換后的影像并沒有像素值。

4 結(jié) 語

本文方法能完成大幅(45 GB)、分幅影像的坐標(biāo)轉(zhuǎn)換過程，并且可以完成大幅、分幅影像的拼接、分割、重分幅等處理過程，處理后影像清晰度與原圖一致，沒有出現(xiàn)失真現(xiàn)象。坐標(biāo)轉(zhuǎn)換過程與ERDAS、ENVI的效率相近，占用內(nèi)存較大，占用CPU較小，與ArcGIS相比效率要快很多，并且內(nèi)存占用與CPU占用相對穩(wěn)定，不會出現(xiàn)太大波動。同時在對分幅影像進(jìn)行坐標(biāo)轉(zhuǎn)換處理時，省去了坐標(biāo)轉(zhuǎn)換后影像拼接、分幅裁剪等過程，大大提高了分幅影像坐標(biāo)轉(zhuǎn)換效率，并避免了轉(zhuǎn)換后影像出現(xiàn)接縫的情況。而且本文方法邏輯簡明，適應(yīng)性強，適用于遙感影像批處理。

[1] 張訓(xùn)虎,劉晉虎,何川,等. 2000 國家大地坐標(biāo)系轉(zhuǎn)換常見問題分析[J]. 測繪通報, 2016 (9): 52-55.

[2] 田桂娥，宋利杰,尹利文,等.地方坐標(biāo)系與CGCS2000坐標(biāo)系轉(zhuǎn)換方法的研究[J]. 測繪工程, 2014,23(8): 66-69.

[3] 呂志平, 魏子卿, 李軍, 等. 我國 CGCS2000 高精度坐標(biāo)轉(zhuǎn)換格網(wǎng)模型的建立[J]. 測繪學(xué)報, 2013, 42(6): 791-797.

[4] 何林, 柳林濤, 許超鈐, 等. 常見平面坐標(biāo)系之間相互轉(zhuǎn)換的方法研究——以 1954 北京坐標(biāo)系, 1980 西安坐標(biāo)系, 2000 國家大地坐標(biāo)系之間的平面坐標(biāo)相互轉(zhuǎn)換為例[J]. 測繪通報, 2014 (9): 6-11.

[5] 查東平,林輝,孫華,等. 基于GDAL的遙感影像數(shù)據(jù)快速讀取與顯示方法的研究[J]. 中南林業(yè)科技大學(xué)學(xué)報, 2013, 33(1):58-62.

[6] 陳谷良. 基于WebGIS的城市電子地圖框架設(shè)計研究[J]. 中國西部科技, 2011,9(35):32-34.

[7] 徐健梅, ArcGIS在數(shù)字正射影像圖坐標(biāo)轉(zhuǎn)換中的應(yīng)用[J]. 長江工程職業(yè)技術(shù)學(xué)院學(xué)報, 2010, 27(3):22-23.

[8] 林旭芳. 基于 ArcGIS 平臺的遙感影像快速分幅方法[J]. 測繪通報, 2014 (S2): 179-181.

[9] 崔麗華,劉善軍,聞彩煥,等. 基于IDL與ENVI的MODIS遙感影像鑲嵌[J]. 河北理工大學(xué)學(xué)報(自然科學(xué)版), 2009, 31(2):127-130.

[10] 張曉東, 吳正鵬, 陳楚, 等. 影像坐標(biāo)轉(zhuǎn)換的一體化處理研究[J]. 城市勘測, 2013(2): 116-117.

[11] 徐永明.遙感二次開發(fā)語言IDL[M].北京:科學(xué)出版社，2014.

[12] 曹元大, 張建芳. Windows 下專家系統(tǒng)開發(fā)工具的內(nèi)存管理[J]. 電腦開發(fā)與應(yīng)用, 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，F(xiàn)AN 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

范鈾，伍超云，王琨，等.海量遙感影像分塊逆映射快速坐標(biāo)轉(zhuǎn)換[J].測繪通報，2017(4)：53-57.

10.13474/j.cnki.11-2246.2017.0119.

2016-12-28；

2017-03-13

國家973計劃(2012CB719904)

范 鈾(1977—)，男，工程師，主要從事遙感數(shù)據(jù)處理方面的研究。E-mail：you.fan@southgis.com

P237

A

0494-0911(2017)04-0053-05