技术领域
[0001] 本
发明涉及地质遥感判释技术，特别是涉及一种利用图像处理和先验知识提取断裂构造信息方法。
背景技术
[0002]
铁路、公路、
水利等长大线路的线位、大型站房、附属设施等工程布置要尽量绕避或以大
角度通过断裂构造，尤其是区域性活动断裂构造，以提高工程的安全性。断裂构造已成为长大线路工程地质勘察的重点内容。
[0003] 近些年，出现了基于遥感信息的人工判释方法，用于进行断裂构造信息提取。但是，这种方法主要靠人工进行目视解译，尚存在一些问题：一方面，受判释人员的经验和主观因素影响大；另一方面，没有充分发挥出遥感技术的优势。仅靠目视判读提取断裂构造信息难度大效率低，很难满足工程建设和安全生产的要求。
发明内容
[0004] 针对现有地质遥感判释技术进行断裂构造信息提取存在的问题，本发明的目的是提供一种利用图像处理和先验知识提取断裂构造信息方法，以实现断裂构造信息的快速准确提取，革新传统的受主观因素影响大的人工判释方法。
[0005] 为实现上述目的，本发明采用的技术方案是提供一种利用图像处理和先验知识提取断裂构造信息方法，该方法用于对区域断裂构造信息的提取，其中，图像处理在遥感图像处理
软件中进行，数据综合分析在
地理信息系统软件中进行，先验知识来自区域地质资料和建立的判释标志库，用于对自动提取的线性信息进行连通和剔除，得到区域断裂构造信息，该方法包括有以下步骤：
[0006] 获取工程区高
分辨率遥感影像及
基础地质资料S1：获取工程区成像时间在两年内、成像分辨率大于5m的遥感影像数据，高分辨率遥感影像应带有有理函数
传感器模型（RPC）参数，影像获取通过代理公司或卫星地面站购买，对工程区进行现场踏勘，并收集工程区地质图和地形图资料；
[0007] 对影像进行边缘增强处理S2：对步骤S1获取的高分辨率遥感影像在遥感图像处理软件中进行预处理，利用遥感图像处理软件中的高通滤波或锐化处理，增强地物的边缘信息，将边缘信息增强结果与步骤S1获取的影像按4：6的比例
叠加，并输出成TIFF格式文件；
[0008] 利用
边缘检测算法提取地物边缘信息S3：对步骤S2得到的地物边缘增强影像利用Canny边缘检测算法进行边缘检测处理，提取地物的边缘信息；
[0009] 利用形态匹配算法提取线性信息S4：对步骤S3得到的地物边缘信息，采用先膨胀后侵蚀的方法进行数学形态学滤波处理，利用Hough变换提取线性信息，此过程利用在MatLab软件中编程实现；
[0010] 利用先验知识对线性信息进行连通和剔除S5：将步骤S2得到的地物边缘增强影像和步骤S4得到的线性信息，发送到地理信息
系统软件，进行空间叠加分析，利用先验知识对线性信息进行连通和剔除，得到断裂构造线；
[0011] 通过地理信息系统综合分析得到断裂构造信息S6：将步骤S5得到的断裂构造线发送到地理信息系统软件，进行长度的量测，划分断裂构造等级，得到区域断裂构造信息。
[0012] 本发明的效果是该提取断裂构造信息方法是利用遥感图像处理技术辅助提取断裂构造信息，与所述人工判释方法相比，能够大幅减弱判释人员经验和主观因素的影响，提高断裂构造信息提取的准确性至80%以上，为更好更快地提取断裂构造信息提供了新的技术手段。
附图说明
[0013] 图1为本发明的一种利用图像处理和先验知识提取断裂构造信息方法的
流程图；
[0014] 图2为本发明实例的研究区原始数据示意图；
[0015] 图3为本发明实例的地物边缘信息增强影像示意图；
[0016] 图4为本发明实例的边缘检测算法提取边界示意图；
[0017] 图5为本发明实例的形态滤波算法提取的线性信息示意图；
[0018] 图6（1）为本发明实例的断裂构造判释标志的河流的直角拐弯示意图；
[0019] 图6（2）为本发明实例的断裂构造判释标志的
断层沟谷与河流的突然拐弯示意图；
[0020] 图7为本发明实例的区域断裂构造线连通和剔除结果示意图；
[0021] 图8为本发明实例的区域断裂构造分布结构图。
[0022] 图中：
[0023] S1、获取工程区高分辨率遥感影像及基础地质资料
[0024] S2、对影像进行边缘增强处理
[0025] S3、利用边缘检测算法提取地物边缘信息
[0026] S4、利用形态匹配算法提取线性信息
[0027] S5、利用先验知识对线性信息进行连通和剔除
[0028] S6、通过地理信息系统综合分析得到断裂构造信息
具体实施方式
[0029] 结合附图对本发明的一种利用图像处理和先验知识提取断裂构造信息方法作进一步描述。
[0030] 如图1所示，本发明的一种利用图像处理和先验知识提取断裂构造信息方法，技术步骤包括：获取工程区高分辨率遥感影像及基础地质资料S1、对影像进行边缘增强处理S2、利用边缘检测算法提取地物边缘信息S3、利用形态匹配算法提取线性信息S4、利用先验知识对线性信息进行连通和剔除S5、通过地理信息系统综合分析得到断裂构造信息S6。
[0031] S1、获取工程区高分辨率遥感影像及基础地质资料
[0032] 获取工程区成像时间在两年内、成像分辨率大于5m的遥感影像数据，如分辨率为2.5m的SPOT5、ALOS等卫星影像、分辨率为2.1m的资源三号卫星影像等。高分辨率遥感影像应带有有理函数传感器模型（RPC）参数，影像获取可通过代理公司或卫星地面站购买。对工程区进行现场踏勘，并收集工程区地质图和地形图资料。
[0033] S2、对影像进行边缘增强处理
[0034] 对步骤S1获取的高分辨率遥感影像在遥感图像处理软件中进行预处理，利用遥感图像处理软件中的高通滤波或锐化处理，增强地物的边缘信息。将边缘信息增强结果与步骤S1获取的影像按4：6的比例叠加，并输出成TIFF格式文件。
[0035] ①遥感影像预处理
[0036] 遥感影像预处理包括影像的几何校正和镶嵌处理，可在遥感图像处理软件中实现。
[0037] 几何校正和镶嵌处理均已有成熟技术，可在ENVI、ERDAS等遥感图像处理软件中实现。
[0038] ②地物边缘信息增强
[0039] 地物边缘信息增强可采用锐化方法或高通滤波方法。
[0040] 利用锐化方法进行地物边缘信息增强，其
数据处理包括遥感图像的载入、锐化算子的选择和进行锐化处理，可在遥感图像处理软件中实现。
[0041] 利用高通滤波方法进行地物边缘信息增强，其数据处理包括遥感图像的载入、高通
滤波器的选择和进行高通滤波处理，可在遥感图像处理软件中实现。
[0042] ③地物边缘信息增强结果输出
[0043] 将得到的地物边缘信息增强结果与步骤S1获取的影像按4：6的比例叠加，并输出成TIFF格式文件，得到地物边缘信息增强影像。
[0044] S3、利用边缘检测算法提取地物边缘信息
[0045] 对步骤S2得到的地物边缘信息增强影像利用Canny边缘检测算法进行边缘检测处理，提取地物的边缘信息。主要步骤如下：
[0046] ①对影像进行平滑处理，以减少噪声影响
[0047] 利用高斯滤波器对步骤S2得到的影像进行平滑处理。首先将步骤S2得到的影像进行归一化处理，以浮点阵列的形式存储，而后对该阵列进行一维和二维高斯滤波，以减少噪声影响。
[0048] ②计算梯度的幅值和方向
[0049] 对平滑处理后的影像，利用一阶偏导的有限差分来计算梯度的幅值和方向。
[0050] ③对梯度幅值进行非极大值抑制
[0051] 对各
像素点的梯度值进行非极大值抑制，即在当前像素点的3×3邻域内，如果当前点的梯度幅值大于沿梯度方向上相邻两个像素点的梯度幅值，则认为该点是可能的边缘点，将其对应的标志
位置1。否则，则认为该点为非边缘点，将其对应的标志位清0。
[0052] ④对经过非极大值抑制处理后的图像进行双
阈值化处理，消除虚假边缘并连接断续边缘，提取地物边缘信息
[0053] 利用文献“冯珂等,一种改进的Canny边缘检测AGT算法.计算机应用与软件,2012,29(3):265-266.”的方法确定高低阈值。用高阈值与经过非极大值抑制后的图像进行比较，记录边缘点。对所有的边缘点，在8邻域内
迭代寻找大于低阈值的点，并标记为边缘点。连接断续的边缘点，得到地物边缘信息。
[0054] S4、利用形态匹配算法提取线性信息
[0055] 对步骤S3得到的地物边缘信息，采用先膨胀后侵蚀的方法进行数学形态学滤波处理，利用Hough变换提取线性信息，此过程可在MatLab中编程实现。
[0056] ①对步骤S3得到的地物边缘信息，进行数学形态滤波处理
[0057] 数学形态滤波处理主要为图像膨胀和图像侵蚀。首先，根据形态学膨胀算子对边缘点所在区域进行区域填充，在确保边界完整的同时，使图像的边缘轮廓变粗，得到图像膨胀处理结果；而后，根据形态学侵蚀算子对图像膨胀处理结果进行细化处理，使得边界更加符合实际，同时为后续边缘提取过程减小难度和错误。
[0058] ②利用hough变换提取线性信息
[0059] 利用在MatLab中编程处理的方式，对数学形态滤波结果进行Hough变换处理，提取影像线性信息。
[0060] S5、利用先验知识对线性信息进行连通和剔除
[0061] 将步骤S2得到的地物边缘增强影像和步骤S4得到的线性信息，发送到地理信息系统软件进行空间叠加分析，利用先验知识对线性信息进行连通和剔除，得到断裂构造线。
[0062] ①数据的叠加分析
[0063] 将步骤S4得到的线性信息和步骤S2得到的影像发送到地理信息系统软件，进行二者的空间叠加分析，保持步骤S4得到的线性信息叠放在步骤S2得到的影像前。
[0064] ②建立先验知识
[0065] 首先，判释人员根据野外踏勘和对地质图、地形图等区域地质资料的分析，初步认清区域的断裂构造分布特征。
[0066] 其次，根据对区域断裂构造的认识，提取区域断裂构造的判释标志。判释标志包括断层
破碎带的直接出露，断层沟谷、断层垭口、断层三角面的分布，以及河流的突然拐弯，岩性或相邻
地貌单元的突变特征。根据提取区域断裂构造的判释标志，建立判释标志库。
[0067] 最后，对区域断裂构造的初步认识和建立的判释标志库，共同构成先验知识。
[0068] ③利用先验知识对线性信息进行连通和剔除
[0069] 基于步骤S2得到的影像，依据先验知识，对步骤S4得到的线性信息进行连通和剔除。
[0070] 首先，利用先验知识，判断步骤S4得到的线性信息是否为断裂构造上的离散线段。具体过程为将步骤S4得到的线性信息和步骤S2得到的影像进行叠加分析，利用先验知识，判断步骤S4得到的线性信息是否与步骤S2得到的影像上区域断裂构造分布特征和判释标志相符。若相符则为断裂构造上的离散线段，若不符则不是断裂构造上的离散线段。
[0071] 其次，对线性信息进行连通和剔除。根据所述利用先验知识判断步骤S4得到的线性信息是否为断裂构造上的离散线段，若不是断裂构造上的离散线段，则进行直接剔除，以删除非断裂构造的线性信息；若是断裂构造上的离散线段，则根据步骤S2得到的影像，利用先验知识对所述影像上同一直线上的线性信息进行连通，形成一条完整的断裂构造线。
[0072] 最后，循环所述对线性信息进行连通和剔除过程，对整个工程区内线性信息进行连通和剔除，得到区域断裂构造线。
[0073] S6、通过地理信息系统综合分析得到断裂构造信息
[0074] 将步骤S5得到的断裂构造线发送到地理信息系统软件，进行长度的量测，划分断裂构造等级，得到区域断裂构造信息。
[0075] 首先，数据的载入与长度测量。将步骤S5得到的断裂构造线发送到地理信息系统软件。利用地理信息系统软件中的距离量算功
能量测每一条断裂构造线的长度（L），并在属性中进行标记。
[0076] 其次，划分断裂构造线等级。断裂构造等级划分标准为：L≧100km为一级断裂；50km≦L﹤100km为二级断裂；25km≦L﹤50km为三级断裂；L﹤25km为四级断裂。依据断裂构造划分标准，按照所述量测的断裂构造线长度对各断裂构造线进行等级划分，得到区域断裂构造信息。
[0077] 下面通过实例对上述方法进行验证：
[0078] S1、获取工程区高分辨率遥感影像及基础地质资料
[0079] 购买了新建浦建龙梅铁路线局部地区分辨率为2.5m的ALOS卫星影像，影像带有有理函数传感器模型（RPC）参数，研究区原始数据如图2所示。图中影像上线性信息纵横交错，地物间
色调反差较小，地物边界和纹理信息不明显。对工程区进行现场踏勘，并收集到工程区地质图和地形图资料。
[0080] S2、对影像进行边缘增强处理
[0081] ①遥感影像预处理
[0082] 基于影像带有的有理函数传感器模型（RPC）参数，利用ERDAS遥感图像处理软件对获取的影像进行了几何校正。
[0083] ②地物边缘信息增强
[0084] 在ENVI遥感图像处理软件中，利用高通滤波方法对所述影像进行地物边缘信息增强处理。
高通滤波器使用为3×3的变换核，中心值为“8”，周围像元值为“-1”。
[0085] ③地物边缘信息增强结果输出
[0086] 将得到的地物边缘信息增强结果与步骤S1获取的影像按4：6的比例叠加，并输出成TIFF格式文件，得到地物边缘信息增强影像如图3所示。由图可见，边缘增强处理后的影像地物界线清晰，边缘信息反映明显。
[0087] S3、利用边缘检测算法提取地物边缘信息
[0088] 对步骤S2得到的地物边缘信息增强影像利用Canny边缘检测算法进行边缘检测处理，提取地物的边缘信息。图4为提取的地物边缘信息结果，由图可见，地物的边界和纹理信息显示为白色线条，其他信息屏蔽为黑色。
[0089] S4、利用形态匹配算法提取线性信息
[0090] 对步骤S3得到的地物边缘信息，采用先膨胀后侵蚀的方法进行数学形态学滤波处理，在MatLab中编程实现Hough变换提取线性信息的过程（部分代码如下）。提取的线性信息结果见图5。由图可见，线性信息被很好地提取出来，其中包括断裂信息、河流信息、地物边界信息等。
[0091] Hough变换提取线性信息的部分代码如下：
[0092]
[0093]
[0094] S5、利用先验知识对线性信息进行连通和剔除
[0095] 将步骤S2得到的地物边缘增强影像和步骤S4得到的线性信息，发送到地理信息系统软件进行空间叠加分析，利用先验知识对线性信息进行连通和剔除，得到断裂构造线。
[0096] ①数据的叠加分析
[0097] 将步骤S4得到的线性信息和步骤S2得到的影像发送到地理信息系统软件，进行二者的空间叠加分析，保持步骤S4得到的线性信息叠放在步骤S2得到的影像前。
[0098] ②建立先验知识
[0099] 首先，判释人员根据野外踏勘和对地质图、地形图等区域地质资料的分析，初步认清区域的断裂构造分布特征。
[0100] 其次，根据对区域断裂构造的认识，提取区域断裂构造的判释标志。图6为提取的该区部分断裂构造判释标志，图6（1）为河流的直角拐弯，可见河流受断裂影响发生两次直角拐弯，并在沿断裂分布段呈直线展布，很好地反映出断裂构造的分布特征；图6（2）为断层沟谷与河流的突然拐弯综合形成的判释标志，断层沟谷笔直而深切，具有明显的方向性，沟谷延伸方向即为断裂构造的走向，期间河流受断裂构造影响而发生直角拐弯。
[0101] 最后，对区域断裂构造的初步认识和建立的判释标志库，共同构成先验知识。
[0102] ③利用先验知识对线性信息进行连通和剔除
[0103] 基于步骤S2得到的影像，依据先验知识，对步骤S4得到的线性信息进行连通和剔除。
[0104] 首先，利用先验知识，判断步骤S4得到的线性信息是否为断裂构造上的离散线段。具体过程为将步骤S4得到的线性信息和步骤S2得到的影像进行叠加分析，利用先验知识，判断步骤S4得到的线性信息是否与步骤S2得到的影像上区域断裂构造分布特征和判释标志相符。若相符则为断裂构造上的离散线段，若不符则不是断裂构造上的离散线段。
[0105] 其次，对线性信息进行连通和剔除。根据所述利用先验知识判断步骤S4得到的线性信息是否为断裂构造上的离散线段，若不是断裂构造上的离散线段，则进行直接剔除，以删除非断裂构造的线性信息；若是断裂构造上的离散线段，则根据步骤S2得到的影像，利用先验知识对所述影像上同一直线上的线性信息进行连通，形成一条完整的断裂构造线。
[0106] 最后，循环所述对线性信息进行连通和剔除过程，对整个工程区内线性信息进行连通和剔除，得到区域断裂构造线如图7所示。由图可见，通过联通和剔除过程处理，研究区内断裂构造分布特征被很好地呈现出来，该区断裂构造走向以北东向为主。
[0107] S6、通过地理信息系统综合分析得到断裂构造信息
[0108] 将步骤S5得到的断裂构造线发送到地理信息系统软件，进行长度的量测，划分断裂构造等级，得到区域断裂构造信息。
[0109] 首先，数据的载入与长度测量。将步骤S5得到的断裂构造线发送到地理信息系统软件。利用地理信息系统软件中的距离量算功能量测每一条断裂构造线的长度（L），并在属性中进行标记。
[0110] 其次，划分断裂构造线等级。断裂构造等级划分标准为：L≧100km为一级断裂；50km≦L﹤100km为二级断裂；25km≦L﹤50km为三级断裂；L﹤25km为四级断裂。依据断裂构造划分标准，按照所述量测的断裂构造线长度对各断裂构造线进行等级划分，得到区域断裂构造信息。根据量测，研究区内断裂构造长度（L）均小于25km，研究区内均为四级断裂，区域断裂构造分布见图8。
