周翠敏,黃映聰,3,4,5*,何月順
(1. 東華理工大學地球科學學院,330013,南昌;2.江西省數(shù)字國土重點實驗室,330013,南昌;3. 東華理工大學網(wǎng)絡(luò)與信息中心,330013,南昌;4.江西省放射性地學大數(shù)據(jù)技術(shù)工程實驗室,330013,南昌;5.核資源與環(huán)境國家重點實驗室,330013,南昌;6.東華理工大學信息工程學院,330013,南昌)
近年來,大數(shù)據(jù)與人工智能算法的引入發(fā)展了數(shù)學地球科學,開啟了大數(shù)據(jù)與人工智能的時代[1]。隨著不同傳感器、不同分辨率、不同成像方式的對地觀測技術(shù)的不斷發(fā)展,人類對地觀測能力和水平也在不斷地提高[2],但也間接地增加了影像的數(shù)據(jù)量,如NASA 地球觀測數(shù)據(jù)的數(shù)據(jù)量以每天4 TB的速度增長[3]。且遙感影像數(shù)據(jù)也呈現(xiàn)出具有體量大(Volume)、速度快(Velocity)、準確性高(Veracity)、多樣化(Variety)、價值高(Value)的“5V”特征[4]。由此遙感也呈現(xiàn)出“大數(shù)據(jù)”的特征,進入了大數(shù)據(jù)時代[2]。
急速增長的數(shù)據(jù)量為數(shù)據(jù)的存儲帶來了諸多困難。長期以來占據(jù)大量存儲空間的海量數(shù)據(jù)得不到有效的利用[5],造成了大量存儲空間的浪費。因此,需要深度挖掘遙感大數(shù)據(jù)特性[6],最大化地發(fā)揮其作用,滿足遙感大數(shù)據(jù)的應(yīng)用需求。遙感大數(shù)據(jù)利用的最終目標也在于對遙感大數(shù)據(jù)中隱藏知識的挖掘[2]。在大數(shù)據(jù)背景下,合理的解釋與應(yīng)用是最大化地發(fā)揮其作用的良好途徑[7]。目前,基于云端遙感信息挖掘應(yīng)用廣泛,如自然災害預測[8]、土地利用方式演變[9-10]、城市化研究[11]、交通行為預測[12]等。對云端遙感影像數(shù)據(jù)進行挖掘有著重要意義與價值。
遙感影像地面控制點的選取是影像處理過程中的一項重要工作。地面控制點的選取主要經(jīng)過采集地面控制點,計算誤差,選擇幾何模型,重采樣輸出、驗證校正結(jié)果等步驟。步驟較為繁瑣,耗時較長且誤差較大,校正的效率較低。隨著遙感影像數(shù)據(jù)的快速增長,廣大的遙感影像元數(shù)據(jù)為影像校正提供了數(shù)據(jù)支撐。尺度不變特征變換(Scale-invariant feature transform,SIFT)算法作為一種提取局部特征的算法,以其提取的特征點的亮度、尺度等具有不變性而得到了廣泛的應(yīng)用[13-14],為本文的研究提供了技術(shù)支撐。
在SIFT算法的基礎(chǔ)上,針對目前人工選取地面控制點方面存在的不足,以臺灣省臺中市為研究區(qū),對云端遙感影像數(shù)據(jù)進行挖掘,運用SIFT算法提取影像特征[15],基于灰度圖像[16]通過建立影像之間的相似性獲取匹配點[17-18],將其應(yīng)用于地面控制點的自動提取方面的研究。通過影像的特征進行匹配來完成影像校正,實現(xiàn)海量遙感影像數(shù)據(jù)的地面控制點的自動提取,提高地面控制點的提取效率。
SIFT(Scale-invariant feature transform,尺度不變特征轉(zhuǎn)換)算法是一種有效地描述影像局部特征的算法。該算法因其對物體的尺度、亮度、旋轉(zhuǎn)等具有不變性而廣泛應(yīng)用于圖像處理等領(lǐng)域。SIFT算法實現(xiàn)特征匹配主要有以下4個步驟。
1)尺度空間極值檢測。通過高斯微分函數(shù)來識別具有尺度、旋轉(zhuǎn)等不變的點作為潛在特征點。圖像尺度空間L(x,y,σ)可表示為尺度的高斯核函數(shù)G(x,y,σ)與圖像I(x,y)的卷積,其表達式為:
L(x,y,σ)=G(x,y,σ)×I(x,y),
其中:
式中:(x,y)代表圖像的像素位置。σ是尺度空間因子,值越小表示圖像被平滑得越少,相應(yīng)的尺度也就越小。小尺度對應(yīng)圖像的細節(jié)特征,大尺度對應(yīng)圖像的概貌特征。
2)關(guān)鍵點方向計算。根據(jù)關(guān)鍵點的分布特性,指定特征點的方向參數(shù),使算子在旋轉(zhuǎn)方向上具有不變性。
3)特征點描述。通過將特征點為中心的16×16的鄰域,分割成16個4×4的子區(qū)域,然后在每個4×4的子區(qū)域中計算生成一個8方向直方圖,通過上述計算獲得4×4×8=128維的特征描述子。
4)特征點匹配。根據(jù)圖像之間的特征點描述子來辨別特征點之間的相似性。一般可采用歐氏距離進行判斷,其距離越小,特征點就越相似。
SIFT算法作為一種成熟的方法,具有尺度、旋轉(zhuǎn)、亮度等不變性且配準精度高等優(yōu)點,對圖像的配準等具有較好的應(yīng)用性,廣泛應(yīng)用于圖像配準等方面的研究。SIFT特征是圖像的局部特征,具有尺度不變特性,在對不同分辨率的影像校正時具有優(yōu)越性。表 1 是主流圖像特征提取HOG、LBP算法與SIFT算法的對比。
表1 SIFT算法與HOG、LBP算法的比較
海量遙感數(shù)據(jù)源主要有數(shù)據(jù)量大、覆蓋范圍廣、分辨率高等特征。隨著不同傳感器、不同分辨率、不同成像方式的對地觀測技術(shù)地不斷發(fā)展,人類對地觀測能力和水平也在不斷地提高,也間接地增加了影像的數(shù)據(jù)量。本文主要基于瓦片影像數(shù)據(jù)進行實驗,瓦片影像金字塔模型是多分辨率的層次模型[19],從底層到頂層,分辨率逐漸變低,但表示的地理范圍不變[20]。隨著瓦片級別的不斷增大,分辨率增高,瓦片數(shù)據(jù)量便會增大。常見的瓦片地圖的參數(shù)如表 2所示,本文所采用的參考瓦片影像源自百度地圖服務(wù)器。
表2 谷歌地圖、百度地圖和高德地圖的瓦片參數(shù)對比
表2 (續(xù))
地理坐標系(Geographic Coordinate System),是使用三維球面定位地球表面位置,通過經(jīng)緯度對地球表面點位定位的坐標系。地理坐標系包括3個部分:角度測量單位、本初子午線和參考橢球體[21]。
本次實驗涉及到WGS-84坐標系、國測局制定的GCJ-02、百度坐標BD09 3種坐標系。由于百度地圖的坐標系是在GCJ-02基礎(chǔ)上進行了加密,為了方便數(shù)據(jù)轉(zhuǎn)換,本次實驗首先將參考影像的坐標系由BD09坐標系轉(zhuǎn)為GCJ-02坐標系,最后將GCJ-02坐標系轉(zhuǎn)為WGS-84坐標系。轉(zhuǎn)換公式如下:
BD09坐標系轉(zhuǎn)GCJ-02坐標系公式:
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:百度坐標維度,
bd_lon:百度坐標經(jīng)度。
GCJ-02坐標系轉(zhuǎn)WGS84坐標系公式:
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的緯度。
為了提高遙感影像配準的效率,本文通過SIFT算法對目標影像與參考影像分別進行特征提取與匹配。主要流程有以下幾個方面:1)對目標影像進行預處理、彩色合成、重采樣、假彩色轉(zhuǎn)灰度圖像、特征點提取等;2)在百度開放平臺下載參考影像瓦片并進行拼接投影、坐標轉(zhuǎn)換、特征提取等;3)匹配目標影像與參考影像特征點,剔除匹配錯誤點;4)組合目標圖像特征點圖像坐標與參考圖像特征點地理坐標為地面控制點文件,實現(xiàn)地面控制點坐標的自動提取。具體匹配流程如圖 1所示。
圖1 自動提取地面控制點流程圖
本實驗的硬件平臺為Intel(R) Core(TM) i7-9700KF CPU,頻率為3.60 GHz,內(nèi)存為16 G。編程平臺為PyCharm2021,操作環(huán)境為64位的Windows10操作系統(tǒng),編程語言為Python3.8,引用庫為GDAL、OpenCV 2.0。
待校正影像數(shù)據(jù)來源于地理空間數(shù)據(jù)云(http://www.gscloud.cn/search),研究區(qū)為臺灣省臺中市分辨率為30 m的Landsat8影像,格式為“TIFF”。通過對 Landsat8 影像進行輻射校正、大氣校正、影像拼接與裁剪等預處理。為減少計算量,提高計算速度,本次研究截取臺中地區(qū)的數(shù)據(jù)(表3、圖2)作為待校正目標影像,并將分辨率重采樣為64 m/pixel。
圖2 研究區(qū)目標影像示意圖
表3 目標影像參數(shù)
參考影像獲取是根據(jù)研究區(qū)的四至范圍和百度地圖瓦片編碼規(guī)則,計算參考影像的行列號,并生成瓦片影像的url地址。根據(jù)瓦片影像的url地址在百度地圖瓦片影像服務(wù)器上下載,然后,利用瓦片行列號所對應(yīng)的位置進行圖像拼接。通過瓦片行列號起止值計算影像的四至地理范圍并將其轉(zhuǎn)換為WGS84坐標系對應(yīng)值,通過計算經(jīng)向和緯向像元分辨率,構(gòu)造仿射變換矩陣。最后,利用Python中的GDAL庫(柵格空間數(shù)據(jù)轉(zhuǎn)換庫)將拼接后的影像寫入仿射變換矩陣和投影信息。參考影像及參考影像參數(shù)如(表4、圖3)所示。
表4 研究區(qū)參考影像數(shù)據(jù)特征一覽表
(西至119.7925°,東至122.1465°,南至21.6798°,北至25.4759°)
由于彩色圖像的信息量過大,而進行圖像識別時,其實只需要使用灰度圖像里的信息就已足夠,所以為了提高運算速度需要將彩色圖像轉(zhuǎn)為灰度圖像;同時,在對圖像處理時,需要計算圖像的紋理信息,所以就必須要將假彩圖像轉(zhuǎn)成灰度圖像來對影像特征點進行提取。其次,由于741波段組合具有兼容中紅外、近紅外及可見光波段信息的優(yōu)勢,圖面色彩豐富,層次感好,具有極為豐富的地質(zhì)信息和地表環(huán)境信息;而且清晰度高,干擾信息少,地質(zhì)可解譯程度高的優(yōu)點。所以基于地學遙感研究的經(jīng)驗做法,最終選取741波段對影像進行假彩色合成,最后將假彩色影像轉(zhuǎn)換為灰度圖像。
根據(jù)下載獲取到的遙感影像數(shù)據(jù),使用SIFT算子計算研究區(qū)域遙感影像的特征點,獲取關(guān)鍵點和描述符,使用SIFT算法檢測角點。首先,參數(shù)的選擇主要是先找到每個像素的中心點,按照最接近中心的特征點的原則進行均勻的選取。通過將數(shù)組的描述轉(zhuǎn)化為矩陣,分別對參考影像和目標影像的2個矩陣進行交叉運算。經(jīng)計算,結(jié)果如表 5所示,得到參考影像的特征點49 465個,耗時13.745 s,目標影像的特征點37 704個,耗時8.861 s。目標影像的特征點及坐標值疊加顯示在目標影像上的效果如圖4所示。同理,即可得到參考影像的特征點,根據(jù)參考影像的特征點,即可得到對應(yīng)的理論坐標及地理信息。將得到的參考影像特征點疊加在影像上的顯示效果如圖 5所示。所用到的特征點提取readFeturue函數(shù)式如下:
圖4 研究區(qū)目標遙感影像疊加特征點示意圖
圖5 研究區(qū)參考影像疊加特征點示意圖
表5 特征點提取結(jié)果表
sift=cv2.xfeatures2d.SIFT_create(),
img=cv2.imread(file,cv2.IMREAD_GRAYSCALE),
kp,des=sift.detectAndCompute(img,None)。
獲取到目標影像和參考影像圖像特征后,通過定義FLANN匹配器進行相似特征點的匹配及錯誤匹配特征點的剔除。將參考影像與目標影像的數(shù)組矩陣進行交叉運算,計算目標特征點與參考特征點的相似度,得到特征點提取的點對,即GCPs文件。本文共匹配749對特征點,匹配耗時139.17 s,匹配結(jié)果如圖 6所示,其主要實現(xiàn)過程如下:1)獲取目標影像特征點坐標及特征點描述符;2)獲取參考影像特征點坐標及特征點描述符;3)計算目標影像與參考影像特征點描述符匹配度;4)剔除錯誤匹配特征點;5)計算特征匹配點對中參考圖像特征點所在地理坐標;6)組合目標圖像特征點圖像坐標與參考圖像特征點地理坐標為地面控制點文件。
圖6 目標影像與參考影像特征點匹配結(jié)果示意圖
在完成目標影像和整個參考影像特征匹配后,根據(jù)目標影像特征點圖像坐標及參考影像圖像坐標及其對應(yīng)的地理坐標,生成適合ArcGIS影像配準的地面控制點文件(共獲取749個地面控制點)(表6)和地理信息文件(表7),該控制點文件可直接在Arcmap的地理配準相應(yīng)模塊導入,提取的所有地面控制點文件可顯示在ArcGIS平臺上(如圖7)。將控制點文件導入到ArcGIS中進行計算之后,得到匹配結(jié)果誤差(如表6)。校正完成之后,實際圖像就得到理論的地理坐標值,即具備地理投影信息。根據(jù)表6的殘差結(jié)果,直觀地表現(xiàn)了該方法在此方面應(yīng)用的優(yōu)越性,計算精度較高,精度可達99%以上。
圖7 控制點顯示在ArcGIS平臺上的示意圖
表6 地面控制點的坐標計算結(jié)果及誤差(部分)
表7 校正完成之后生成的投影參數(shù)
1)遙感影像的校正主要有2種方式:一種方法是從影像到地圖的校正,參考具有標準地圖投影的地圖,采集地面控制點,確定輸入像元坐標和該點對應(yīng)的地圖坐標之間的幾何關(guān)系來完成校正。另一種方法是從影像到影像的校正,在2幅影像中選取同名點,使2幅影像匹配達到一致。
本研究嘗試采用網(wǎng)絡(luò)云端提供的海量遙感數(shù)據(jù)進行配準研究,選用百度平臺提供的瓦片數(shù)據(jù)進行影像校正,以提高影像校正的效率。通過對百度瓦片進行編碼轉(zhuǎn)換,將分辨率轉(zhuǎn)換為瓦片級別,利用云端影像數(shù)據(jù)實現(xiàn)了目標影像與參考影像的智能匹配,表明使用云端海量遙感瓦片數(shù)據(jù)進行影像校正是可行的。
2)本研究主要對全圖進行特征點的提取與匹配,選擇大于目標影像范圍的影像作為參考影像,以確保目標影像的區(qū)域落在參考影像區(qū)域內(nèi)。因此,參考影像特征點要多于目標影像特征點。本次實驗共提取參考影像特征點49 465個,提取目標影像特征點37 704個,最終成功提取地面控制點749個,實現(xiàn)了地面控制點的自動化提取,且精度高達99%以上。但當提取的特征點較多時,可能會增加計算量。如果通過優(yōu)化SIFT算法,只提取固定地物的特征,可能會減少特征點的提取數(shù)量來降低計算量,提高提取效率。例如:只提取圖像固定地物(標志性建筑、道路、河道等)的特征,通過只提取長期固定不變的地物特征,減少特征點的提取,降低計算量,提高特征提取的速度與效率,這將是未來研究工作需要解決的問題。
3)SIFT特征是圖像的局部特征,具有尺度不變特性。因此,本研究選取分辨率為64 m的影像作為目標影像,選取128 m分辨率的影像作為參考影像,分別對2個影像進行特征提取并成功實現(xiàn)特征匹配,說明利用云端瓦片影像可以實現(xiàn)用低分辨率影像校正較高分辨率的影像。
在遙感影像校正過程中,通過使用云端瓦片影像進行特征提取與匹配,能夠?qū)崿F(xiàn)影像智能化校正,而且精度較高。此外,影像校正不受分辨率的影響,可以用較低分辨率的影像校正較高分辨率的影像。