徐振恒,潘蔚,张元涛

(核工业北京地质研究院 遥感信息与图像分析技术国家级重点实验室,北京 100029)

我国幅员辽阔,不同地区的地形地貌差异明显,地质调查与矿产勘查任务艰巨。遥感技术因其具备覆盖范围广、获取数据快、受地面条件限制少等优势,在国民经济和国防建设等领域取得了广泛的应用[1--5]。但是,以往的卫星遥感应用研究大多是基于国外的遥感数据,部分数据价格昂贵,抬高了遥感在国内各领域应用的门槛,限制了基层单位、特别是民营企业应用和发展遥感技术[6]。为了改变这一局面,近年来我国大力发展了国产遥感卫星观测系统,截至目前为止,我国已经发射了200余颗地球遥感卫星,现如今已经形成了高分、资源、环境/实践、风云等遥感系列卫星观测体系,并有天绘、高景、北京、珠海、吉林、珞珈、三极等小遥感卫星系列[7]。

但是国产遥感卫星在一些性能方面与国外卫星仍存在一定差别,比如光谱分辨率。而高分系列卫星却有着较高的空间分辨率,如GF2(高分二号)--PMS最高达0.8 m,GF1(高分一号)--PMS最高也可达2 m,较高空间分辨率意味着获取的遥感影像具有更加丰富的纹理信息[7]。所以,如何利用这些丰富的纹理信息,并结合一些遥感算法来充分发挥国产高分数据的优势,弥补国产卫星光谱分辨率的不足,最终达到实际应用的目的具有重要的理论和实际价值。本文以高分卫星系列中的高分一号PMS遥感数据为例,针对我国热液铀成矿区的典型岩石种类,使用“地形结构--岩性组分模型”进行了岩性识别研究[8--10]。该模型是由潘蔚等于2009年提出,并利用ETM 遥感影像对不同地区的目视解译难以区分的二长花岗岩与沉积变质岩进行了有效的识别区分[10]。

1 实验数据与研究对象

1.1 研究区

图1 重点研究区遥感图像与取样位置示意图Fig.1 GF1--PMS image and the sample location of study area

研究区位于扬子地台与华南褶皱系的交接部位,涉及湖南、江西和广东3省,其中诸广山岩体是重点研究区(图1)。诸广山岩体是我国重要的产铀岩体之一,位于华南褶皱系的华夏地块西侧,区域上处于闽赣粤加里东隆起与湘桂粤北海西--印支坳陷的结合部位,北东向的吴川--四会断裂及北西向的汕头--安仁深断裂穿过本区[11--13]。

1.2 数据源

研究使用到的数据包括GF1--PMS传感器的全色波段遥感影像、区域地质图以及DEM 数 据。其 中,GF1--PMS 全 色 波 段(0.45～0.90μm)影像的空间分辨率为2 m[7],数据来自于“国家航天局对地观测与数据中心高分辨率对地观测网格平台”。该平台可下载高分系列遥感卫星与资源系列遥感卫星的高质量遥感影像,具有遥感数据类型丰富、数据量大、更新速度快、免费下载等优点。区域地质图为1∶50万,从国家基础地理信息平台获取,经MapGIS格式转换后,经属性编辑形成。DEM 数据下载自EarthData平台,空间分辨率30 m。

1.3 岩性样本图像选取

为了识别与铀成矿关系密切的岩性,我们从重点研究区岩性地层中选择下寒武统沉积岩--1、中下寒武统浅变质岩--1和加里东期花岗岩--1作为识别对象,使用的1 000×1 000像元作为裁剪窗口选取样本图像(图2)。

为了检验算法的普适性,研究中同种岩性单元要分别在两个不同的地区进行裁剪,尽可能选择相距遥远的地区,以避免距离过近的局部相似性的影响。为了与重点研究区进行对比,在重点研究区外的湖南和广东两省,对应选取了地质图上岩性相同的3个图像样本:下寒武统沉积岩--2、中下寒武统浅变质岩--2和加里东期花岗岩--2,6 个样本区的GF1--PMS 图像作为验证对象(图3,表1)。

图2 主要成矿区地层及岩性分布示意图Fig.2 Schematic map of strata and lithological units in the main metallogenic area

图3 GF1--PMS图像样本与对应岩性Fig.3 GF1--PMS images of different lithological type

表1 样本图像的岩性特征Table 1 Lithological feature of the image sample

2 研究原理与方法

研究中使用“地形结构--岩性组分模型”提取遥感影像岩石纹理信息中的岩性特征组分,并使用α--f(α)多重分形谱对岩性特征组分进行描述(图4)。

2.1 地形结构--岩性组分模型

研究中使用潘蔚等在2009年提出的“地形结构--岩性组分模型”来对各样本遥感影像的岩性特征进行提取[10]。该模型基于光学遥感图像的形态--组分模型,使用小波变换并结合DEM 的分形维数来确定最佳小波分解级别,从而有效分离出遥感影像的地形特征组分与岩性特征组分[14--15]。该模型的原理为:因为地形变化是连续的,可以用小波分解的低频成分表示;岩性组分是跳跃不连续的,可以用小波分解的高频成分表示。又由于地形具有良好的自相似性,所以可以使用分形维数对其进行有效的描述,那么最能反映遥感影像地形特征组分的小波分解低频成分的分形维数,理论上应该接近于对应的DEM 的分形维数,故可据此确定小波分解的最佳级别,从而达到有效分离遥感影像的地形特征组分与岩性特征组分的目的[16--19]。小波分解后,高频成分图像具有水平、垂直、对角线3个方向,而岩石纹理的发育与岩性、产状和经受的应力作用有关,三者的空间不均匀性也会造成岩石纹理方向分布不均匀性,因此利用高频成分进行岩性特征组分表示时,需要考虑3个方向的影响。

图4 遥感影像岩性特征组分提取流程Fig.4 Extraction process for lithology feature of remote sensing image

本研究基于该模型,依据流程进行了遥感影像岩性特征组分的提取。处理流程主要通过MATLAB 编译实现,其中分形维数使用的是MATLAB的Fraclab2.2工具箱进行计算。首先,输入高分影像灰度图并进行浮点化、归一化操作;其次,对灰度图进行二维多级小波变换,并对各级小波变换低频分量进行重构,得到各级低频成分;然后,调用Fraclab 工具箱中的“数盒子法(Box--Counting Method)”分别计算上一步中得到的重构后的各级小波变换低频成分以及对应区域DEM 的分形维数;接下来,通过对比二者分形维数之间差值的绝对值,取低频成分分形维数与DEM 分形维数差值最小的分解级别作为最佳分解级别;最后,保留最佳分解级别对应低频成分作为地形特征组分,重构该级别3方向(水平、垂直、对角线)高频分量,并取重构后3方向高频成分均值为岩性特征组分[16--18]。

2.2 多重分形谱与岩性分类

多重分形谱的描述方法较为常用的为α--f(α)谱与q--Dq谱[20],本研究使用前者进行计算,其常用于对图像纹理的分析计算。α--f(α)多重分形谱最经典的计算方法是基于小波变换模极大值的奇异谱计算法[21],奇异指数α的定义为:对于函数f(x)上一点x0,设k＞0(k为常数)且α(x0)＞n,如果存在多项式Pn(x－x0),使得在x0的某个邻域内的任意一点x,满足|f(x)－Pn(x－x0)|≤k|x－x0|α,则称α(x0)是函数在点x0处的奇异指数。

α--f(α)多重分形谱可以用来描述构成复合分形体的具有不同奇异程度的所有单个分形体的维数分布和变化特征,其计算步骤为:先计算每个点的奇异指数α;再计算由所有相同奇异值α的点构成的集合的分形维数;最终求得每个奇异指数所代表的点集的分形维数,就可以得到与不同奇异指数对应的一系列分形维数,即多重分形谱。实验中使用Fraclab 工具箱中的“Multifractal Spectra”工具来具体实现α--f(α)多重分形谱的计算。

利用上述方法,计算了6个样本图像的岩性组分图像的α--f(α)谱。通过对多重分形谱特征的分析,来研究岩性分类识别的方法,探索多重分形谱开展岩性分类识别的数学描述模型和可能的参数。

3 结果分析与讨论

图5 各样本图像岩性特征组分的α--f(α)多重分形谱(－0.15＜α＜0.8)Fig.5 Multifractalα--f(α)spectrum(－0.15＜α＜0.8)of lithological component of sample images

图6 各样本图像岩性特征组分的α--f(α)多重分形谱(－0.15＜α＜0)Fig.6 Multifractalα--f(α)spectrum(－0.15＜α＜0)of lithological component of sample images

通过对各岩性样本平均高频成分(水平、垂直、对角线3方向平均加权)多重分形谱的观察可以得到:在α＜0区间内,尤其是－0.15＜α＜0,3种岩性之间表现出较好的区分性,可以相互区别。对于同一岩性,其多重分形谱也是相似的,具体表现为其在数值上(f(α))与变化趋势(曲线斜率)上相近似(图5,图6)。对于相同的f(α),α的值从大到小依次为:下寒武统＞中下寒武统＞加里东期,其中加里东期与其余两者的区别十分明显,下寒武统与中下寒武统数值上相对较为接近,但也有明显差别,中下寒武统整体上要大于下寒武统,且其两者曲线走势不同。多重分形谱识别结果与目视解译结果相同,加里东期遥感影像细节纹理特征多呈发散树枝状,明显区别于下寒武统与中下寒武统的遥感影像纹理特征。

在α＞0区间,加里东期表现最好,不仅自身曲线相似而且明显区别于其余两者,中下寒武统α较小时表现较好,但随着α的增大,曲线走势逐渐分离,而下寒武统表现最差,不仅自身曲线相分离,且与中下寒武统相混淆。

综上所述,多重分形谱在α＜0区间具有岩性识别能力,可有效识别不同岩性单元。

4 结论

本文使用了国产高分一号遥感卫星PMS传感器影像数据在我国南方热液铀成矿区进行了岩性识别应用研究。结果为:通过“地形结构--岩性组分模型”提取的高分一号遥感影像岩性特征组分的多重分形谱,可以对不同类型岩性单元进行有效区分,其区分度主要表现在多重分形谱α＜0区间内。不同岩性单元遥感影像纹理特征具有不同的优势方向,一般情况下需要综合考虑各个方向的纹理特征才能够表征整体的岩性特征。可见,国产高分一号遥感卫星数据在结合纹理算法模型的基础上,可以有效弥补光谱分辨率的不足,为我国地质资源勘查应用提供一种新的方法与补充。