周翠敏,黄映聪,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个影像进行特征提取并成功实现特征匹配,说明利用云端瓦片影像可以实现用低分辨率影像校正较高分辨率的影像。
在遥感影像校正过程中,通过使用云端瓦片影像进行特征提取与匹配,能够实现影像智能化校正,而且精度较高。此外,影像校正不受分辨率的影响,可以用较低分辨率的影像校正较高分辨率的影像。