基于张量投票耦合霍夫变换的地质线性体提取方法

摘要

本发明提供了一种基于张量投票耦合霍夫变换的地质线性体提取方法，包括以下步骤：对遥感影像进行预处理，得到预处理后的遥感影像；在预处理后的遥感影像的N个多光谱波段中选择3个最优波段，得到由最优波段组合的遥感影像，其中N为大于等于3的自然数；利用高斯高通滤波对由最优波段组合的遥感影像进行锐化处理，增强线性化边缘信息；对增强线性化边缘信息的遥感影像中进行边缘检测，得到遥感影像中的所有边缘点；将遥感影像中的所有边缘点从图像坐标系转换到参数坐标系中，从参数坐标系中提取地质线性体。本发明对解译专家的知识和经验依赖更少，极大程度上缩短了处理的时间，节省了大量的人力，具有较强的实用价值和推广意义。

说明书

技术领域

本发明属于遥感地质及影像处理技术领域，具体涉及一种基于张量投票耦合霍夫变换的地质线性体提取方法。

背景技术

断裂带、断层等地质构造属于地质薄弱区，由于侵蚀等作用形成线状地貌。在地质构造的作用下，这些线状地貌在遥感影像上往往表现出明显的线状分布，称为地质线性体，其控制着地下流体(矿液、地下水和油气等)的运移和矿产资源的空间赋存，其方位与密度对分析区域构造运动趋势与活动程度有着深远的科学意义和实用价值。

从上世纪八十年代以来，国内外不少遥感地质工作者利用专家知识经验与图像处理的方法，提取了遥感影像与数字高程模型(Digital elevation model,DEM)的地质线性体，分析了区域地质构造趋势及影响程度。如Jansson等(2005)使用Landsat 7ETM+和数字高程模型提取地质线性体，并绘制了Wales东北部冰川地貌；吴婧等(2011)使用Canny边缘检测和霍夫变换提取了断裂构造信息，并使用ENVI+IDL等编程语言实现；袁小祥等(2011)利用假彩色合成、主成分变换、缨帽变换、波段比值、地貌渲染等方法，突出了地质线性体在多源遥感数据上的反差，提取了活动断层；Alaa,A.M等(2011)根据DEM数据山体阴影增强可以线性特征的特点，使用线索追踪算法提取山体阴影渲染后的影像上的地质线性体，并用B-spline曲线整合表达，最终评价了区域地质构造环境状况；Yusof等(2011)分析了高速公路周围滑坡灾害分布与地质线性体密度之间的关系；刘智荣等(2012)利用反差增强、彩色合成、方向滤波与图像融合等图像增强处理算法，提取与分析了银川活动断层信息及分布；Bahiru等(2016)使用Landsat ETM+和SRTMDEM数据提取并绘制了乌干达地区地质线性构造分布，研究了该区域金矿的分布，对矿产预测有非常重要的研究意义。

上述工作都不同程度地取得了一些成就，主要存在以下不足：

1)目视解译结果的正确性依赖于解译专家的经验和知识背景，耗时耗力且效率低。

2)计算机解译的精度与处理速度、数据源的分辨率等因素有关，影像越大，处理速度就慢；分辨率过高，则容易受到道路、土地利用界线等线性特征影响，而产生大量的错误线性边缘，也就产生了过多的噪声。过度依赖参数设置，导致方普适性差。

发明内容

针对现有技术存在的不足，本发明的目的在于，提供一种基于张量投票耦合霍夫变换的地质线性体提取方法，解决现有技术中依赖于解译专家的经验和知识背景，耗时耗力且效率低、处理速度慢、噪声大及普适性差的问题。

为了解决上述技术问题，本发明采用如下技术方案予以实现：

一种基于张量投票耦合霍夫变换的地质线性体提取方法，包括以下步骤：

步骤1，对遥感影像进行预处理，得到预处理后的遥感影像；

步骤2，在预处理后的遥感影像的N个多光谱波段中选择3个最优波段，得到由最优波段组合的遥感影像，N为大于等于3的自然数；

步骤3，利用高斯高通滤波对由最优波段组合的遥感影像进行锐化处理，增强线性化边缘信息；

步骤4，对步骤3增强线性化边缘信息的遥感影像中进行边缘检测，得到遥感影像中的所有边缘点；

步骤5，将遥感影像中的所有边缘点从图像坐标系转换到参数坐标系中，从参数坐标系中提取地质线性体；

所述图像坐标系为以遥感影像的任一角作为原点，遥感影像的水平方向为x轴，遥感影像的竖直方向为y轴；

所述参数坐标系为ρ＝x cosθ+y sinθ，其中θ和ρ为边缘点对应的参数坐标系下的极坐标。

进一步地，所述步骤1中对遥感影像进行预处理，包括：

从遥感影像库中任选一遥感影像作为当前遥感影像，该当前遥感影像上含云量低于5％，且具备有理函数模型及统计影像灰度信息；所述影像灰度信息包括每个波段的灰度方差和标准差；

对该当前遥感影像进行辐射定标、大气校正、影像裁剪处理。

进一步地，所述步骤2中在预处理后的遥感影像的N个多光谱波段中选择3个最优波段，包括：

选择预处理后的遥感影像的多光谱波段中波段指数OIF最大时所对应的3个波段作为最优波段；

其中，通过式(1)得到波段指数O_IF：

式(1)中，S_i为第i个波段的标准差；R_ij为第i个波段和第j个波段的相关系数；i≠j，i＝1,2,...,N,j＝1,2,...,N；

D(i)、D(j)分别为第i个波段和第j个波段的方差；Cov(i,j)为第i个波段和第j个波段之间的协方差。

进一步地，所述步骤3中利用高斯高通滤波对由最优波段组合的遥感影像进行锐化处理，包括：

通过式(2)对由最优波段组合的遥感影像进行高斯高通滤波，得到滤波后的遥感图像H(u,v)：

式(2)中，D(u,v)是由最优波段组合的遥感影像的频率域中点(u,v)与频率矩形中心的距离，D₀为常数。

进一步地，所述步骤4中对步骤3增强线性化边缘信息的遥感影像中进行边缘检测，得到遥感影像中的所有边缘点，包括：

步骤41，任选增强线性化边缘信息的遥感影像中的一个像素点作为当前像素点，对当前像素点作拉普拉斯算子卷积操作，得到张量矩阵T；

其中，I为增强线性化边缘信息的遥感影像，和别为遥感影像I沿x和y方向的二阶导数；其中，x、y是以遥感影像I的任一角作为原点，遥感影像I的水平方向为x轴，遥感影像I的竖直方向为y轴建立的图像坐标系中的x轴和y轴；

步骤42，将张量矩阵T进行矩阵谱分解，得到棒状分量和球状分量；

其中，为棒状分量，为球状分量；

步骤43，当(λ1-λ2)＞λ2时，所述当前像素点为边缘点；否则，所述当前像素点为非边缘点；

步骤44，重复步骤41至步骤43，直至增强线性化边缘信息的遥感影像的所有像素点都作为当前像素点为止，得到遥感影像中的所有边缘点。

进一步地，所述步骤5中将遥感影像中的所有边缘点从图像坐标系转换到参数坐标系中，从参数坐标系中提取地质线性体，包括：

步骤51，遍历参数坐标系，寻找局部最大值的点，将该局部最大值的点作为峰值点，设峰值点的坐标为(ρ，θ)，其中(ρ，θ)为遥感影像中地质线性体的斜率和截距；

步骤52，将参数坐标系下峰值点所对应的坐标转换到图像坐标系中，按照边缘点的方向及端点距离连接成线，得到地质线性体影像，即完成地质线性体的提取。本发明与现有技术相比，具有如下技术效果：

1.本发明综合了波段选择、图像增强、边界检测、线性提取等算法与规则，提供了基于张量投票耦合霍夫变换的地质线性体提取方法，相对于单纯的目视解译等方法，对解译专家的知识和经验依赖更少，极大程度上缩短了处理的时间，节省了大量的人力，具有较强的实用价值和推广意义；

2.较Canny边缘检测算法，基于张量投票的边缘检测算法既可以在边缘检测的基础上进行边界检测，也可以在灰度图像上直接用二维圆形投票域进行张量投票过程，然后进行投票解释，并根据边界特征的显著性提取边界，具有鲁棒性；

3.本发明可以处理多源遥感数据，方法的参数设置更为均衡，具有较好的普适性，对区域构造演化及板块运动方向有较强的指示作用。

附图说明

图1为本发明的基于张量投票耦合霍夫变换的地质线性体提取方法的总体流程图；

图2为本发明实例的研究区7-5-4波段合成遥感影像；

图3为本发明实例的研究区高斯高通滤波后增强影像；

图4为本发明张量投票的示意图；

图5为本发明实例的研究区张量投票后边缘检测影像；

图6为本发明霍夫变换的示意图；

图7为本发明地质线性体提取影像。

以下结合附图对本发明的具体内容作进一步详细解释说明。

具体实施方式

本发明中的遥感影像为Landsat 8OLI冬季多光谱遥感影像。

以下给出本发明的具体实施例，需要说明的是本发明并不局限于以下具体实施例，凡在本申请技术方案基础上做的等同变换均落入本发明的保护范围。

如图1所示，本实施例提供了基于张量投票耦合霍夫变换的地质线性体提取方法，包括以下步骤：

步骤1，对遥感影像进行预处理，得到预处理后的遥感影像；

具体地，本实施例从遥感影像库中任选一遥感影像作为当前遥感影像，该当前遥感影像上含云量低于5％，且具备有理函数模型及统计影像灰度信息；所述影像灰度信息包括每个波段的灰度方差和标准差；

对该当前遥感影像进行辐射定标、大气校正、影像裁剪处理。

其中，本实施例中的遥感影像取Landsat 8OLI冬季多光谱遥感影像，因为冬季影像植被覆盖较少，线性边缘在影像上突出明显。影像上云含量应低于5％，具备有理函数模型(Rational Polynomial Coefficients，RPC)及统计影像灰度信息，以便与地质资料比对，验证预处理的精度；

影像灰度统计信息主要包括，每个波段的灰度方差D_i、标准差S_i，统计信息如下表所示，用于后续分析。

表1 Landsat 8OLI多光谱遥感影像

步骤2，在预处理后的遥感影像的N个多光谱波段中选择3个最优波段，得到由最优波段组合的遥感影像，N为大于等于3的自然数；

具体地，本实施例选择预处理后的遥感影像的多光谱波段中波段指数OIF最大时所对应的3个波段作为最优波段；

其中，通过式(1)得到波段指数OIF：

式(1)中，S_i为第i个波段的标准差；R_ij为第i个波段和第j个波段的相关系数；i≠j，i＝1,2,...,N,j＝1,2,...,N；

D(i)、D(j)分别为第i个波段和第j个波段的方差；Cov(i,j)为第i个波段和第j个波段之间的协方差。

本实施例中利用Landsat 8OLI的遥感影像的6个多光谱波段，计算波段相关系数矩阵，如表2所示。

表2研究区Landsat 8OLI波段间相关系数矩阵

其中，相关系数越大说明，波段间冗余信息越多，对影像的波段之间的相关系数进行分析，不利于突出不同地质体的信息。表2中Band 5与其他波段的相关系数普遍较小，因此影像合成可优先考虑Band 5。Band 7为短波红外波段，对岩石、特定矿物反应敏感，用于区分主要岩石类型、岩石水热蚀变，探测相关粘土矿物等，其次优先考虑Band 7。

遥感影像的波段组合如表3所示，从表3中可以看出，7-5-4组合方式与7-5-6的组合方式是最大的，但Band分别与Band 4、Band 6的相关系数为0.447和0.417，因此应选择7-5-4波段合成方式。最终确定7、5、4波段分别赋予红、绿、蓝合成得到的影像，如图2，增强各种地物间的对比度。即图2为由最优波段组合的遥感影像。

表3 Landsat 8OLI波段组合和相应的OIF指数

此外，去除多余波段，目的在于抵制其余波段对合成影像的干扰，如使用全部波段，任意两个波段的相关系数一旦过高，必会造成数据冗余。

步骤3，利用高斯高通滤波对由最优波段组合的遥感影像进行锐化处理，增强线性化边缘信息；

具体地，通过式(2)对由最优波段组合的遥感影像进行高斯高通滤波，得到滤波后的遥感图像H(u,v)：

式(2)中，D(u,v)是由最优波段组合的遥感影像的频率域中点(u,v)与频率矩形中心的距离，D₀为常数。本实施例中D0为截止频率，具体取20。

本实施例利用高斯高通滤波对遥感影像线性边缘的锐化作用，处理步骤2中得到的由最优波段组合的遥感影像，突出线性边缘，获取地质线性体提取的增强线性化边缘信息的遥感影像。

其优点在于：高通滤波是一种典型的图像锐化增强算子，采用高通滤波的方法可以让高频分量顺利通过、低频分量受到抑制，也就是增强边缘特征，抑制非边缘等噪声。而高斯高通滤波就是高通滤波的一种，增强了图像的对比度，使其具有更为明显的区分信息。

如图3所示为经过高斯高通滤波后的遥感影像，山谷线和山脊线等地质地貌线在影像上被突出，呈现明显的线性地貌。

步骤4，对步骤3增强线性化边缘信息的遥感影像中进行边缘检测，得到遥感影像中的所有边缘点；

进一步的，步骤4的边界点矢量和叠加的显著性特征，定义如下：

其中，DF(L，k)显著性函数，L是曲线长度，k是曲率，σ是投票邻域范围，c是控制曲率衰减度的系数。

具体地，步骤41，任选增强线性化边缘信息的遥感影像中的一个像素点作为当前像素点，对当前像素点作拉普拉斯算子卷积操作，得到张量矩阵T；

其中，I为增强线性化边缘信息的遥感影像，和分别为遥感影像I沿x和y方向的二阶导数；其中，x、y是以遥感影像I的任一角作为原点，遥感影像I的水平方向为x轴，遥感影像I的竖直方向为y轴建立的图像坐标系中的x轴和y轴；

本实施例利用拉普拉斯算子计算张量矩阵T中的二阶导数，算子如下：

进行奇异值分解，可得：

步骤42，将张量矩阵T进行矩阵谱分解，得到棒状分量和球状分量；

其中，为棒状分量，为球状分量；

步骤43，当(λ1-λ2)＞λ2时，所述当前像素点为边缘点；否则，所述当前像素点为非边缘点；

具体地，当λ1≈λ2时，该当前像素点位于一个区域内部点或者交叉处；当λ1和λ2取值都非常小，该当前像素点判断为非边缘点。

步骤44，重复步骤41至步骤43，直至增强线性化边缘信息的遥感影像的所有像素点都作为当前像素点为止，得到遥感影像中的所有边缘点。

如图4为张量投票的示意图，其中，线性边缘点由于其灰度突变性，经过8邻域范围点对其投票，矢量均有显著变化，具体表现为显著性函数大于0，因此将其视为边缘点且被保存；非线性边缘点由于邻域范围点为同一物质，经过邻域投票，矢量和无显著变化，具体表现为显著性函数约等于0；因此将其视为非边缘点且被删除。

本实施例的步骤4利用张量表示的待提取特征经过稀疏和多尺度密集投票，叠加矢量和，非地质线性体边界点得到不同方向的投票，而矢量相互作用而抵消，边界点由于只得到某一侧的投票而得到增强，最后通过投票解释过程获取区域的边界；利用显著性函数计算边界点叠加后矢量和的显著特征，实现张量投票方法的边缘检测，得到二值黑白影像，如图5所示。

步骤5，将遥感影像中的所有边缘点从图像坐标系转换到参数坐标系中，从参数坐标系中提取地质线性体；

本实施例使用霍夫变换将边缘点从图像坐标系转换到参数坐标系中，转而求参数坐标系的峰值，记录峰值对应的位置，提取地质线性体。

图像坐标系为以遥感影像的任一角作为原点，遥感影像的水平方向为x轴，遥感影像的竖直方向为y轴；

参数坐标系为ρ＝x cosθ+y sinθ，其中θ和ρ为边缘点对应的参数坐标系下的极坐标。

具体地，步骤51，遍历参数坐标系，寻找局部最大值的点，将该局部最大值的点作为峰值点，设峰值点的坐标为(ρ，θ)，其中(ρ，θ)为遥感影像中地质线性体的斜率和截距；

本实施例使用霍夫变换将边缘点从图像坐标系转换到参数坐标系，原始的边缘点坐标转换到了以(ρ，θ)为基准的参数坐标系，不断累加。遍历ρ-θ空间，寻找出局部最大值(极值)的点，称为峰值点。

图6从左至右分别为图像坐标系下的像素点通过霍夫变换转换到参数坐标下的示意图，再通过不断累加得到图6最右边中的一个峰值点。

步骤52，将参数坐标系下峰值点所对应的坐标转换到图像坐标系中，按照边缘点的方向及端点距离连接成线，得到地质线性体影像，即完成地质线性体的提取。如图7所示为叠加霍夫变换后的该区域地质线性体分布图。