周翠敏,黃映聰,3,4,5*,何月順
(1. 東華理工大學地球科學學院,330013,南昌;2.江西省數字國土重點實驗室,330013,南昌;3. 東華理工大學網絡與信息中心,330013,南昌;4.江西省放射性地學大數據技術工程實驗室,330013,南昌;5.核資源與環境國家重點實驗室,330013,南昌;6.東華理工大學信息工程學院,330013,南昌)
近年來,大數據與人工智能算法的引入發展了數學地球科學,開啟了大數據與人工智能的時代[1]。隨著不同傳感器、不同分辨率、不同成像方式的對地觀測技術的不斷發展,人類對地觀測能力和水平也在不斷地提高[2],但也間接地增加了影像的數據量,如NASA 地球觀測數據的數據量以每天4 TB的速度增長[3]。且遙感影像數據也呈現出具有體量大(Volume)、速度快(Velocity)、準確性高(Veracity)、多樣化(Variety)、價值高(Value)的“5V”特征[4]。由此遙感也呈現出“大數據”的特征,進入了大數據時代[2]。
急速增長的數據量為數據的存儲帶來了諸多困難。長期以來占據大量存儲空間的海量數據得不到有效的利用[5],造成了大量存儲空間的浪費。因此,需要深度挖掘遙感大數據特性[6],最大化地發揮其作用,滿足遙感大數據的應用需求。遙感大數據利用的最終目標也在于對遙感大數據中隱藏知識的挖掘[2]。在大數據背景下,合理的解釋與應用是最大化地發揮其作用的良好途徑[7]。目前,基于云端遙感信息挖掘應用廣泛,如自然災害預測[8]、土地利用方式演變[9-10]、城市化研究[11]、交通行為預測[12]等。對云端遙感影像數據進行挖掘有著重要意義與價值。
遙感影像地面控制點的選取是影像處理過程中的一項重要工作。地面控制點的選取主要經過采集地面控制點,計算誤差,選擇幾何模型,重采樣輸出、驗證校正結果等步驟。步驟較為繁瑣,耗時較長且誤差較大,校正的效率較低。隨著遙感影像數據的快速增長,廣大的遙感影像元數據為影像校正提供了數據支撐。尺度不變特征變換(Scale-invariant feature transform,SIFT)算法作為一種提取局部特征的算法,以其提取的特征點的亮度、尺度等具有不變性而得到了廣泛的應用[13-14],為本文的研究提供了技術支撐。
在SIFT算法的基礎上,針對目前人工選取地面控制點方面存在的不足,以臺灣省臺中市為研究區,對云端遙感影像數據進行挖掘,運用SIFT算法提取影像特征[15],基于灰度圖像[16]通過建立影像之間的相似性獲取匹配點[17-18],將其應用于地面控制點的自動提取方面的研究。通過影像的特征進行匹配來完成影像校正,實現海量遙感影像數據的地面控制點的自動提取,提高地面控制點的提取效率。
SIFT(Scale-invariant feature transform,尺度不變特征轉換)算法是一種有效地描述影像局部特征的算法。該算法因其對物體的尺度、亮度、旋轉等具有不變性而廣泛應用于圖像處理等領域。SIFT算法實現特征匹配主要有以下4個步驟。
1)尺度空間極值檢測。通過高斯微分函數來識別具有尺度、旋轉等不變的點作為潛在特征點。圖像尺度空間L(x,y,σ)可表示為尺度的高斯核函數G(x,y,σ)與圖像I(x,y)的卷積,其表達式為:
L(x,y,σ)=G(x,y,σ)×I(x,y),
其中:
式中:(x,y)代表圖像的像素位置。σ是尺度空間因子,值越小表示圖像被平滑得越少,相應的尺度也就越小。小尺度對應圖像的細節特征,大尺度對應圖像的概貌特征。
2)關鍵點方向計算。根據關鍵點的分布特性,指定特征點的方向參數,使算子在旋轉方向上具有不變性。
3)特征點描述。通過將特征點為中心的16×16的鄰域,分割成16個4×4的子區域,然后在每個4×4的子區域中計算生成一個8方向直方圖,通過上述計算獲得4×4×8=128維的特征描述子。
4)特征點匹配。根據圖像之間的特征點描述子來辨別特征點之間的相似性。一般可采用歐氏距離進行判斷,其距離越小,特征點就越相似。
SIFT算法作為一種成熟的方法,具有尺度、旋轉、亮度等不變性且配準精度高等優點,對圖像的配準等具有較好的應用性,廣泛應用于圖像配準等方面的研究。SIFT特征是圖像的局部特征,具有尺度不變特性,在對不同分辨率的影像校正時具有優越性。表 1 是主流圖像特征提取HOG、LBP算法與SIFT算法的對比。

表1 SIFT算法與HOG、LBP算法的比較
海量遙感數據源主要有數據量大、覆蓋范圍廣、分辨率高等特征。隨著不同傳感器、不同分辨率、不同成像方式的對地觀測技術地不斷發展,人類對地觀測能力和水平也在不斷地提高,也間接地增加了影像的數據量。本文主要基于瓦片影像數據進行實驗,瓦片影像金字塔模型是多分辨率的層次模型[19],從底層到頂層,分辨率逐漸變低,但表示的地理范圍不變[20]。隨著瓦片級別的不斷增大,分辨率增高,瓦片數據量便會增大。常見的瓦片地圖的參數如表 2所示,本文所采用的參考瓦片影像源自百度地圖服務器。

表2 谷歌地圖、百度地圖和高德地圖的瓦片參數對比

表2 (續)
地理坐標系(Geographic Coordinate System),是使用三維球面定位地球表面位置,通過經緯度對地球表面點位定位的坐標系。地理坐標系包括3個部分:角度測量單位、本初子午線和參考橢球體[21]。
本次實驗涉及到WGS-84坐標系、國測局制定的GCJ-02、百度坐標BD09 3種坐標系。由于百度地圖的坐標系是在GCJ-02基礎上進行了加密,為了方便數據轉換,本次實驗首先將參考影像的坐標系由BD09坐標系轉為GCJ-02坐標系,最后將GCJ-02坐標系轉為WGS-84坐標系。轉換公式如下:
BD09坐標系轉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:百度坐標經度。
GCJ-02坐標系轉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的經度,
lat:GCJ02的緯度。
為了提高遙感影像配準的效率,本文通過SIFT算法對目標影像與參考影像分別進行特征提取與匹配。主要流程有以下幾個方面:1)對目標影像進行預處理、彩色合成、重采樣、假彩色轉灰度圖像、特征點提取等;2)在百度開放平臺下載參考影像瓦片并進行拼接投影、坐標轉換、特征提取等;3)匹配目標影像與參考影像特征點,剔除匹配錯誤點;4)組合目標圖像特征點圖像坐標與參考圖像特征點地理坐標為地面控制點文件,實現地面控制點坐標的自動提取。具體匹配流程如圖 1所示。

圖1 自動提取地面控制點流程圖
本實驗的硬件平臺為Intel(R) Core(TM) i7-9700KF CPU,頻率為3.60 GHz,內存為16 G。編程平臺為PyCharm2021,操作環境為64位的Windows10操作系統,編程語言為Python3.8,引用庫為GDAL、OpenCV 2.0。
待校正影像數據來源于地理空間數據云(http://www.gscloud.cn/search),研究區為臺灣省臺中市分辨率為30 m的Landsat8影像,格式為“TIFF”。通過對 Landsat8 影像進行輻射校正、大氣校正、影像拼接與裁剪等預處理。為減少計算量,提高計算速度,本次研究截取臺中地區的數據(表3、圖2)作為待校正目標影像,并將分辨率重采樣為64 m/pixel。

圖2 研究區目標影像示意圖

表3 目標影像參數
參考影像獲取是根據研究區的四至范圍和百度地圖瓦片編碼規則,計算參考影像的行列號,并生成瓦片影像的url地址。根據瓦片影像的url地址在百度地圖瓦片影像服務器上下載,然后,利用瓦片行列號所對應的位置進行圖像拼接。通過瓦片行列號起止值計算影像的四至地理范圍并將其轉換為WGS84坐標系對應值,通過計算經向和緯向像元分辨率,構造仿射變換矩陣。最后,利用Python中的GDAL庫(柵格空間數據轉換庫)將拼接后的影像寫入仿射變換矩陣和投影信息。參考影像及參考影像參數如(表4、圖3)所示。

表4 研究區參考影像數據特征一覽表

(西至119.7925°,東至122.1465°,南至21.6798°,北至25.4759°)
由于彩色圖像的信息量過大,而進行圖像識別時,其實只需要使用灰度圖像里的信息就已足夠,所以為了提高運算速度需要將彩色圖像轉為灰度圖像;同時,在對圖像處理時,需要計算圖像的紋理信息,所以就必須要將假彩圖像轉成灰度圖像來對影像特征點進行提取。其次,由于741波段組合具有兼容中紅外、近紅外及可見光波段信息的優勢,圖面色彩豐富,層次感好,具有極為豐富的地質信息和地表環境信息;而且清晰度高,干擾信息少,地質可解譯程度高的優點。所以基于地學遙感研究的經驗做法,最終選取741波段對影像進行假彩色合成,最后將假彩色影像轉換為灰度圖像。
根據下載獲取到的遙感影像數據,使用SIFT算子計算研究區域遙感影像的特征點,獲取關鍵點和描述符,使用SIFT算法檢測角點。首先,參數的選擇主要是先找到每個像素的中心點,按照最接近中心的特征點的原則進行均勻的選取。通過將數組的描述轉化為矩陣,分別對參考影像和目標影像的2個矩陣進行交叉運算。經計算,結果如表 5所示,得到參考影像的特征點49 465個,耗時13.745 s,目標影像的特征點37 704個,耗時8.861 s。目標影像的特征點及坐標值疊加顯示在目標影像上的效果如圖4所示。同理,即可得到參考影像的特征點,根據參考影像的特征點,即可得到對應的理論坐標及地理信息。將得到的參考影像特征點疊加在影像上的顯示效果如圖 5所示。所用到的特征點提取readFeturue函數式如下:

圖4 研究區目標遙感影像疊加特征點示意圖

圖5 研究區參考影像疊加特征點示意圖

表5 特征點提取結果表
sift=cv2.xfeatures2d.SIFT_create(),
img=cv2.imread(file,cv2.IMREAD_GRAYSCALE),
kp,des=sift.detectAndCompute(img,None)。
獲取到目標影像和參考影像圖像特征后,通過定義FLANN匹配器進行相似特征點的匹配及錯誤匹配特征點的剔除。將參考影像與目標影像的數組矩陣進行交叉運算,計算目標特征點與參考特征點的相似度,得到特征點提取的點對,即GCPs文件。本文共匹配749對特征點,匹配耗時139.17 s,匹配結果如圖 6所示,其主要實現過程如下:1)獲取目標影像特征點坐標及特征點描述符;2)獲取參考影像特征點坐標及特征點描述符;3)計算目標影像與參考影像特征點描述符匹配度;4)剔除錯誤匹配特征點;5)計算特征匹配點對中參考圖像特征點所在地理坐標;6)組合目標圖像特征點圖像坐標與參考圖像特征點地理坐標為地面控制點文件。

圖6 目標影像與參考影像特征點匹配結果示意圖
在完成目標影像和整個參考影像特征匹配后,根據目標影像特征點圖像坐標及參考影像圖像坐標及其對應的地理坐標,生成適合ArcGIS影像配準的地面控制點文件(共獲取749個地面控制點)(表6)和地理信息文件(表7),該控制點文件可直接在Arcmap的地理配準相應模塊導入,提取的所有地面控制點文件可顯示在ArcGIS平臺上(如圖7)。將控制點文件導入到ArcGIS中進行計算之后,得到匹配結果誤差(如表6)。校正完成之后,實際圖像就得到理論的地理坐標值,即具備地理投影信息。根據表6的殘差結果,直觀地表現了該方法在此方面應用的優越性,計算精度較高,精度可達99%以上。

圖7 控制點顯示在ArcGIS平臺上的示意圖

表6 地面控制點的坐標計算結果及誤差(部分)

表7 校正完成之后生成的投影參數
1)遙感影像的校正主要有2種方式:一種方法是從影像到地圖的校正,參考具有標準地圖投影的地圖,采集地面控制點,確定輸入像元坐標和該點對應的地圖坐標之間的幾何關系來完成校正。另一種方法是從影像到影像的校正,在2幅影像中選取同名點,使2幅影像匹配達到一致。
本研究嘗試采用網絡云端提供的海量遙感數據進行配準研究,選用百度平臺提供的瓦片數據進行影像校正,以提高影像校正的效率。通過對百度瓦片進行編碼轉換,將分辨率轉換為瓦片級別,利用云端影像數據實現了目標影像與參考影像的智能匹配,表明使用云端海量遙感瓦片數據進行影像校正是可行的。
2)本研究主要對全圖進行特征點的提取與匹配,選擇大于目標影像范圍的影像作為參考影像,以確保目標影像的區域落在參考影像區域內。因此,參考影像特征點要多于目標影像特征點。本次實驗共提取參考影像特征點49 465個,提取目標影像特征點37 704個,最終成功提取地面控制點749個,實現了地面控制點的自動化提取,且精度高達99%以上。但當提取的特征點較多時,可能會增加計算量。如果通過優化SIFT算法,只提取固定地物的特征,可能會減少特征點的提取數量來降低計算量,提高提取效率。例如:只提取圖像固定地物(標志性建筑、道路、河道等)的特征,通過只提取長期固定不變的地物特征,減少特征點的提取,降低計算量,提高特征提取的速度與效率,這將是未來研究工作需要解決的問題。
3)SIFT特征是圖像的局部特征,具有尺度不變特性。因此,本研究選取分辨率為64 m的影像作為目標影像,選取128 m分辨率的影像作為參考影像,分別對2個影像進行特征提取并成功實現特征匹配,說明利用云端瓦片影像可以實現用低分辨率影像校正較高分辨率的影像。
在遙感影像校正過程中,通過使用云端瓦片影像進行特征提取與匹配,能夠實現影像智能化校正,而且精度較高。此外,影像校正不受分辨率的影響,可以用較低分辨率的影像校正較高分辨率的影像。