周翠敏，黃映聰,3,4,5*，何月順

(1. 東華理工大學(xué)地球科學(xué)學(xué)院，330013，南昌；2.江西省數(shù)字國(guó)土重點(diǎn)實(shí)驗(yàn)室，330013，南昌；3. 東華理工大學(xué)網(wǎng)絡(luò)與信息中心，330013，南昌；4.江西省放射性地學(xué)大數(shù)據(jù)技術(shù)工程實(shí)驗(yàn)室，330013，南昌；5.核資源與環(huán)境國(guó)家重點(diǎn)實(shí)驗(yàn)室，330013，南昌；6.東華理工大學(xué)信息工程學(xué)院，330013，南昌)

0 引言

近年來(lái)，大數(shù)據(jù)與人工智能算法的引入發(fā)展了數(shù)學(xué)地球科學(xué)，開(kāi)啟了大數(shù)據(jù)與人工智能的時(shí)代[1]。隨著不同傳感器、不同分辨率、不同成像方式的對(duì)地觀測(cè)技術(shù)的不斷發(fā)展，人類(lèi)對(duì)地觀測(cè)能力和水平也在不斷地提高[2]，但也間接地增加了影像的數(shù)據(jù)量，如NASA 地球觀測(cè)數(shù)據(jù)的數(shù)據(jù)量以每天4 TB的速度增長(zhǎng)[3]。且遙感影像數(shù)據(jù)也呈現(xiàn)出具有體量大(Volume)、速度快(Velocity)、準(zhǔn)確性高(Veracity)、多樣化(Variety)、價(jià)值高(Value)的“5V”特征[4]。由此遙感也呈現(xiàn)出“大數(shù)據(jù)”的特征，進(jìn)入了大數(shù)據(jù)時(shí)代[2]。

急速增長(zhǎng)的數(shù)據(jù)量為數(shù)據(jù)的存儲(chǔ)帶來(lái)了諸多困難。長(zhǎng)期以來(lái)占據(jù)大量存儲(chǔ)空間的海量數(shù)據(jù)得不到有效的利用[5]，造成了大量存儲(chǔ)空間的浪費(fèi)。因此，需要深度挖掘遙感大數(shù)據(jù)特性[6]，最大化地發(fā)揮其作用，滿足遙感大數(shù)據(jù)的應(yīng)用需求。遙感大數(shù)據(jù)利用的最終目標(biāo)也在于對(duì)遙感大數(shù)據(jù)中隱藏知識(shí)的挖掘[2]。在大數(shù)據(jù)背景下，合理的解釋與應(yīng)用是最大化地發(fā)揮其作用的良好途徑[7]。目前，基于云端遙感信息挖掘應(yīng)用廣泛，如自然災(zāi)害預(yù)測(cè)[8]、土地利用方式演變[9-10]、城市化研究[11]、交通行為預(yù)測(cè)[12]等。對(duì)云端遙感影像數(shù)據(jù)進(jìn)行挖掘有著重要意義與價(jià)值。

遙感影像地面控制點(diǎn)的選取是影像處理過(guò)程中的一項(xiàng)重要工作。地面控制點(diǎn)的選取主要經(jīng)過(guò)采集地面控制點(diǎn)，計(jì)算誤差，選擇幾何模型，重采樣輸出、驗(yàn)證校正結(jié)果等步驟。步驟較為繁瑣，耗時(shí)較長(zhǎng)且誤差較大，校正的效率較低。隨著遙感影像數(shù)據(jù)的快速增長(zhǎng)，廣大的遙感影像元數(shù)據(jù)為影像校正提供了數(shù)據(jù)支撐。尺度不變特征變換(Scale-invariant feature transform，SIFT)算法作為一種提取局部特征的算法，以其提取的特征點(diǎn)的亮度、尺度等具有不變性而得到了廣泛的應(yīng)用[13-14]，為本文的研究提供了技術(shù)支撐。

在SIFT算法的基礎(chǔ)上，針對(duì)目前人工選取地面控制點(diǎn)方面存在的不足，以臺(tái)灣省臺(tái)中市為研究區(qū)，對(duì)云端遙感影像數(shù)據(jù)進(jìn)行挖掘，運(yùn)用SIFT算法提取影像特征[15]，基于灰度圖像[16]通過(guò)建立影像之間的相似性獲取匹配點(diǎn)[17-18]，將其應(yīng)用于地面控制點(diǎn)的自動(dòng)提取方面的研究。通過(guò)影像的特征進(jìn)行匹配來(lái)完成影像校正，實(shí)現(xiàn)海量遙感影像數(shù)據(jù)的地面控制點(diǎn)的自動(dòng)提取，提高地面控制點(diǎn)的提取效率。

1 基于SIFT算法的圖像配準(zhǔn)

1.1 SIFT算法原理

SIFT(Scale-invariant feature transform，尺度不變特征轉(zhuǎn)換)算法是一種有效地描述影像局部特征的算法。該算法因其對(duì)物體的尺度、亮度、旋轉(zhuǎn)等具有不變性而廣泛應(yīng)用于圖像處理等領(lǐng)域。SIFT算法實(shí)現(xiàn)特征匹配主要有以下4個(gè)步驟。

1)尺度空間極值檢測(cè)。通過(guò)高斯微分函數(shù)來(lái)識(shí)別具有尺度、旋轉(zhuǎn)等不變的點(diǎn)作為潛在特征點(diǎn)。圖像尺度空間L(x,y,σ)可表示為尺度的高斯核函數(shù)G(x,y,σ)與圖像I(x,y)的卷積，其表達(dá)式為：

L(x,y,σ)=G(x,y,σ)×I(x,y)，

其中：

式中：(x,y)代表圖像的像素位置。σ是尺度空間因子，值越小表示圖像被平滑得越少，相應(yīng)的尺度也就越小。小尺度對(duì)應(yīng)圖像的細(xì)節(jié)特征，大尺度對(duì)應(yīng)圖像的概貌特征。

2)關(guān)鍵點(diǎn)方向計(jì)算。根據(jù)關(guān)鍵點(diǎn)的分布特性，指定特征點(diǎn)的方向參數(shù)，使算子在旋轉(zhuǎn)方向上具有不變性。

3)特征點(diǎn)描述。通過(guò)將特征點(diǎn)為中心的16×16的鄰域，分割成16個(gè)4×4的子區(qū)域，然后在每個(gè)4×4的子區(qū)域中計(jì)算生成一個(gè)8方向直方圖，通過(guò)上述計(jì)算獲得4×4×8=128維的特征描述子。

4)特征點(diǎn)匹配。根據(jù)圖像之間的特征點(diǎn)描述子來(lái)辨別特征點(diǎn)之間的相似性。一般可采用歐氏距離進(jìn)行判斷，其距離越小，特征點(diǎn)就越相似。

SIFT算法作為一種成熟的方法，具有尺度、旋轉(zhuǎn)、亮度等不變性且配準(zhǔn)精度高等優(yōu)點(diǎn)，對(duì)圖像的配準(zhǔn)等具有較好的應(yīng)用性，廣泛應(yīng)用于圖像配準(zhǔn)等方面的研究。SIFT特征是圖像的局部特征，具有尺度不變特性，在對(duì)不同分辨率的影像校正時(shí)具有優(yōu)越性。表 1 是主流圖像特征提取HOG、LBP算法與SIFT算法的對(duì)比。

表1 SIFT算法與HOG、LBP算法的比較

1.2 海量遙感數(shù)據(jù)源的特征

海量遙感數(shù)據(jù)源主要有數(shù)據(jù)量大、覆蓋范圍廣、分辨率高等特征。隨著不同傳感器、不同分辨率、不同成像方式的對(duì)地觀測(cè)技術(shù)地不斷發(fā)展，人類(lèi)對(duì)地觀測(cè)能力和水平也在不斷地提高，也間接地增加了影像的數(shù)據(jù)量。本文主要基于瓦片影像數(shù)據(jù)進(jìn)行實(shí)驗(yàn)，瓦片影像金字塔模型是多分辨率的層次模型[19]，從底層到頂層，分辨率逐漸變低，但表示的地理范圍不變[20]。隨著瓦片級(jí)別的不斷增大，分辨率增高，瓦片數(shù)據(jù)量便會(huì)增大。常見(jiàn)的瓦片地圖的參數(shù)如表 2所示，本文所采用的參考瓦片影像源自百度地圖服務(wù)器。

表2 谷歌地圖、百度地圖和高德地圖的瓦片參數(shù)對(duì)比

表2 (續(xù))

1.3 影像坐標(biāo)關(guān)系的轉(zhuǎn)換

地理坐標(biāo)系(Geographic Coordinate System)，是使用三維球面定位地球表面位置，通過(guò)經(jīng)緯度對(duì)地球表面點(diǎn)位定位的坐標(biāo)系。地理坐標(biāo)系包括3個(gè)部分：角度測(cè)量單位、本初子午線和參考橢球體[21]。

本次實(shí)驗(yàn)涉及到WGS-84坐標(biāo)系、國(guó)測(cè)局制定的GCJ-02、百度坐標(biāo)BD09 3種坐標(biāo)系。由于百度地圖的坐標(biāo)系是在GCJ-02基礎(chǔ)上進(jìn)行了加密，為了方便數(shù)據(jù)轉(zhuǎn)換，本次實(shí)驗(yàn)首先將參考影像的坐標(biāo)系由BD09坐標(biāo)系轉(zhuǎn)為GCJ-02坐標(biāo)系，最后將GCJ-02坐標(biāo)系轉(zhuǎn)為WGS-84坐標(biāo)系。轉(zhuǎn)換公式如下：

BD09坐標(biāo)系轉(zhuǎn)GCJ-02坐標(biāo)系公式：

x=bd_lon-0.0065，

y=bd_lat-0.006，

z=math.squt(x*x+y*y)-0.00002*math.sin(y*math.sin(y*selt.x_pi)，

theat=math.atan2(y,x)-0.000003*math.cos(x*selt.x_pi)，

gg_lng=z*math.cos(theta)，

gg_lat=z*math.sin(theta)，

bd_lat:百度坐標(biāo)維度，

bd_lon:百度坐標(biāo)經(jīng)度。

GCJ-02坐標(biāo)系轉(zhuǎn)WGS84坐標(biāo)系公式：

dlat=(dlat*180)/((selt.a*(1-self.ee))/(magic*squtmagic)*self.pi,

dlng=(dlng*180)/(self.a/sqrtmagic*math.cos(radlat)*self.pi,

mglat=lat+dlat,

mglng=lng+dlng,

lng:GCJ02的經(jīng)度，

lat:GCJ02的緯度。

1.4 控制點(diǎn)坐標(biāo)提取流程

為了提高遙感影像配準(zhǔn)的效率，本文通過(guò)SIFT算法對(duì)目標(biāo)影像與參考影像分別進(jìn)行特征提取與匹配。主要流程有以下幾個(gè)方面：1)對(duì)目標(biāo)影像進(jìn)行預(yù)處理、彩色合成、重采樣、假彩色轉(zhuǎn)灰度圖像、特征點(diǎn)提取等；2)在百度開(kāi)放平臺(tái)下載參考影像瓦片并進(jìn)行拼接投影、坐標(biāo)轉(zhuǎn)換、特征提取等；3)匹配目標(biāo)影像與參考影像特征點(diǎn)，剔除匹配錯(cuò)誤點(diǎn)；4)組合目標(biāo)圖像特征點(diǎn)圖像坐標(biāo)與參考圖像特征點(diǎn)地理坐標(biāo)為地面控制點(diǎn)文件，實(shí)現(xiàn)地面控制點(diǎn)坐標(biāo)的自動(dòng)提取。具體匹配流程如圖 1所示。

圖1 自動(dòng)提取地面控制點(diǎn)流程圖

2 實(shí)驗(yàn)與分析

本實(shí)驗(yàn)的硬件平臺(tái)為Intel(R) Core(TM) i7-9700KF CPU，頻率為3.60 GHz，內(nèi)存為16 G。編程平臺(tái)為PyCharm2021，操作環(huán)境為64位的Windows10操作系統(tǒng)，編程語(yǔ)言為Python3.8，引用庫(kù)為GDAL、OpenCV 2.0。

2.1 實(shí)驗(yàn)數(shù)據(jù)

待校正影像數(shù)據(jù)來(lái)源于地理空間數(shù)據(jù)云(http://www.gscloud.cn/search)，研究區(qū)為臺(tái)灣省臺(tái)中市分辨率為30 m的Landsat8影像，格式為“TIFF”。通過(guò)對(duì) Landsat8 影像進(jìn)行輻射校正、大氣校正、影像拼接與裁剪等預(yù)處理。為減少計(jì)算量，提高計(jì)算速度，本次研究截取臺(tái)中地區(qū)的數(shù)據(jù)(表3、圖2)作為待校正目標(biāo)影像，并將分辨率重采樣為64 m/pixel。

圖2 研究區(qū)目標(biāo)影像示意圖

表3 目標(biāo)影像參數(shù)

參考影像獲取是根據(jù)研究區(qū)的四至范圍和百度地圖瓦片編碼規(guī)則，計(jì)算參考影像的行列號(hào)，并生成瓦片影像的url地址。根據(jù)瓦片影像的url地址在百度地圖瓦片影像服務(wù)器上下載，然后，利用瓦片行列號(hào)所對(duì)應(yīng)的位置進(jìn)行圖像拼接。通過(guò)瓦片行列號(hào)起止值計(jì)算影像的四至地理范圍并將其轉(zhuǎn)換為WGS84坐標(biāo)系對(duì)應(yīng)值，通過(guò)計(jì)算經(jīng)向和緯向像元分辨率，構(gòu)造仿射變換矩陣。最后，利用Python中的GDAL庫(kù)(柵格空間數(shù)據(jù)轉(zhuǎn)換庫(kù))將拼接后的影像寫(xiě)入仿射變換矩陣和投影信息。參考影像及參考影像參數(shù)如(表4、圖3)所示。

表4 研究區(qū)參考影像數(shù)據(jù)特征一覽表

(西至119.7925°，東至122.1465°，南至21.6798°，北至25.4759°)

由于彩色圖像的信息量過(guò)大，而進(jìn)行圖像識(shí)別時(shí)，其實(shí)只需要使用灰度圖像里的信息就已足夠，所以為了提高運(yùn)算速度需要將彩色圖像轉(zhuǎn)為灰度圖像；同時(shí)，在對(duì)圖像處理時(shí)，需要計(jì)算圖像的紋理信息，所以就必須要將假彩圖像轉(zhuǎn)成灰度圖像來(lái)對(duì)影像特征點(diǎn)進(jìn)行提取。其次，由于741波段組合具有兼容中紅外、近紅外及可見(jiàn)光波段信息的優(yōu)勢(shì)，圖面色彩豐富，層次感好，具有極為豐富的地質(zhì)信息和地表環(huán)境信息；而且清晰度高，干擾信息少，地質(zhì)可解譯程度高的優(yōu)點(diǎn)。所以基于地學(xué)遙感研究的經(jīng)驗(yàn)做法，最終選取741波段對(duì)影像進(jìn)行假彩色合成，最后將假彩色影像轉(zhuǎn)換為灰度圖像。

2.2 影像特征提取

根據(jù)下載獲取到的遙感影像數(shù)據(jù)，使用SIFT算子計(jì)算研究區(qū)域遙感影像的特征點(diǎn)，獲取關(guān)鍵點(diǎn)和描述符，使用SIFT算法檢測(cè)角點(diǎn)。首先，參數(shù)的選擇主要是先找到每個(gè)像素的中心點(diǎn)，按照最接近中心的特征點(diǎn)的原則進(jìn)行均勻的選取。通過(guò)將數(shù)組的描述轉(zhuǎn)化為矩陣，分別對(duì)參考影像和目標(biāo)影像的2個(gè)矩陣進(jìn)行交叉運(yùn)算。經(jīng)計(jì)算，結(jié)果如表 5所示，得到參考影像的特征點(diǎn)49 465個(gè)，耗時(shí)13.745 s，目標(biāo)影像的特征點(diǎn)37 704個(gè)，耗時(shí)8.861 s。目標(biāo)影像的特征點(diǎn)及坐標(biāo)值疊加顯示在目標(biāo)影像上的效果如圖4所示。同理，即可得到參考影像的特征點(diǎn)，根據(jù)參考影像的特征點(diǎn)，即可得到對(duì)應(yīng)的理論坐標(biāo)及地理信息。將得到的參考影像特征點(diǎn)疊加在影像上的顯示效果如圖 5所示。所用到的特征點(diǎn)提取readFeturue函數(shù)式如下：

圖4 研究區(qū)目標(biāo)遙感影像疊加特征點(diǎn)示意圖

圖5 研究區(qū)參考影像疊加特征點(diǎn)示意圖

表5 特征點(diǎn)提取結(jié)果表

sift=cv2.xfeatures2d.SIFT_create()，

img=cv2.imread(file,cv2.IMREAD_GRAYSCALE)，

kp,des=sift.detectAndCompute(img,None)。

2.3 影像特征點(diǎn)匹配

獲取到目標(biāo)影像和參考影像圖像特征后，通過(guò)定義FLANN匹配器進(jìn)行相似特征點(diǎn)的匹配及錯(cuò)誤匹配特征點(diǎn)的剔除。將參考影像與目標(biāo)影像的數(shù)組矩陣進(jìn)行交叉運(yùn)算，計(jì)算目標(biāo)特征點(diǎn)與參考特征點(diǎn)的相似度，得到特征點(diǎn)提取的點(diǎn)對(duì)，即GCPs文件。本文共匹配749對(duì)特征點(diǎn)，匹配耗時(shí)139.17 s，匹配結(jié)果如圖 6所示，其主要實(shí)現(xiàn)過(guò)程如下：1)獲取目標(biāo)影像特征點(diǎn)坐標(biāo)及特征點(diǎn)描述符；2)獲取參考影像特征點(diǎn)坐標(biāo)及特征點(diǎn)描述符；3)計(jì)算目標(biāo)影像與參考影像特征點(diǎn)描述符匹配度；4)剔除錯(cuò)誤匹配特征點(diǎn)；5)計(jì)算特征匹配點(diǎn)對(duì)中參考圖像特征點(diǎn)所在地理坐標(biāo)；6)組合目標(biāo)圖像特征點(diǎn)圖像坐標(biāo)與參考圖像特征點(diǎn)地理坐標(biāo)為地面控制點(diǎn)文件。

圖6 目標(biāo)影像與參考影像特征點(diǎn)匹配結(jié)果示意圖

2.4 實(shí)驗(yàn)結(jié)果分析

在完成目標(biāo)影像和整個(gè)參考影像特征匹配后，根據(jù)目標(biāo)影像特征點(diǎn)圖像坐標(biāo)及參考影像圖像坐標(biāo)及其對(duì)應(yīng)的地理坐標(biāo)，生成適合ArcGIS影像配準(zhǔn)的地面控制點(diǎn)文件(共獲取749個(gè)地面控制點(diǎn))(表6)和地理信息文件(表7)，該控制點(diǎn)文件可直接在Arcmap的地理配準(zhǔn)相應(yīng)模塊導(dǎo)入，提取的所有地面控制點(diǎn)文件可顯示在ArcGIS平臺(tái)上(如圖7)。將控制點(diǎn)文件導(dǎo)入到ArcGIS中進(jìn)行計(jì)算之后，得到匹配結(jié)果誤差(如表6)。校正完成之后，實(shí)際圖像就得到理論的地理坐標(biāo)值，即具備地理投影信息。根據(jù)表6的殘差結(jié)果，直觀地表現(xiàn)了該方法在此方面應(yīng)用的優(yōu)越性，計(jì)算精度較高，精度可達(dá)99%以上。

圖7 控制點(diǎn)顯示在ArcGIS平臺(tái)上的示意圖

表6 地面控制點(diǎn)的坐標(biāo)計(jì)算結(jié)果及誤差(部分)

表7 校正完成之后生成的投影參數(shù)

3 討論

1)遙感影像的校正主要有2種方式：一種方法是從影像到地圖的校正，參考具有標(biāo)準(zhǔn)地圖投影的地圖，采集地面控制點(diǎn)，確定輸入像元坐標(biāo)和該點(diǎn)對(duì)應(yīng)的地圖坐標(biāo)之間的幾何關(guān)系來(lái)完成校正。另一種方法是從影像到影像的校正，在2幅影像中選取同名點(diǎn)，使2幅影像匹配達(dá)到一致。

本研究嘗試采用網(wǎng)絡(luò)云端提供的海量遙感數(shù)據(jù)進(jìn)行配準(zhǔn)研究，選用百度平臺(tái)提供的瓦片數(shù)據(jù)進(jìn)行影像校正，以提高影像校正的效率。通過(guò)對(duì)百度瓦片進(jìn)行編碼轉(zhuǎn)換，將分辨率轉(zhuǎn)換為瓦片級(jí)別，利用云端影像數(shù)據(jù)實(shí)現(xiàn)了目標(biāo)影像與參考影像的智能匹配，表明使用云端海量遙感瓦片數(shù)據(jù)進(jìn)行影像校正是可行的。

2)本研究主要對(duì)全圖進(jìn)行特征點(diǎn)的提取與匹配，選擇大于目標(biāo)影像范圍的影像作為參考影像，以確保目標(biāo)影像的區(qū)域落在參考影像區(qū)域內(nèi)。因此，參考影像特征點(diǎn)要多于目標(biāo)影像特征點(diǎn)。本次實(shí)驗(yàn)共提取參考影像特征點(diǎn)49 465個(gè)，提取目標(biāo)影像特征點(diǎn)37 704個(gè)，最終成功提取地面控制點(diǎn)749個(gè)，實(shí)現(xiàn)了地面控制點(diǎn)的自動(dòng)化提取，且精度高達(dá)99%以上。但當(dāng)提取的特征點(diǎn)較多時(shí)，可能會(huì)增加計(jì)算量。如果通過(guò)優(yōu)化SIFT算法，只提取固定地物的特征，可能會(huì)減少特征點(diǎn)的提取數(shù)量來(lái)降低計(jì)算量，提高提取效率。例如：只提取圖像固定地物(標(biāo)志性建筑、道路、河道等)的特征，通過(guò)只提取長(zhǎng)期固定不變的地物特征，減少特征點(diǎn)的提取，降低計(jì)算量，提高特征提取的速度與效率，這將是未來(lái)研究工作需要解決的問(wèn)題。

3)SIFT特征是圖像的局部特征，具有尺度不變特性。因此，本研究選取分辨率為64 m的影像作為目標(biāo)影像，選取128 m分辨率的影像作為參考影像，分別對(duì)2個(gè)影像進(jìn)行特征提取并成功實(shí)現(xiàn)特征匹配，說(shuō)明利用云端瓦片影像可以實(shí)現(xiàn)用低分辨率影像校正較高分辨率的影像。

4 結(jié)論

在遙感影像校正過(guò)程中，通過(guò)使用云端瓦片影像進(jìn)行特征提取與匹配，能夠?qū)崿F(xiàn)影像智能化校正，而且精度較高。此外，影像校正不受分辨率的影響，可以用較低分辨率的影像校正較高分辨率的影像。