CN104915640B - 一种钻孔图像结构面自动识别与参数提取方法 - Google Patents
一种钻孔图像结构面自动识别与参数提取方法 Download PDFInfo
- Publication number
- CN104915640B CN104915640B CN201510266649.2A CN201510266649A CN104915640B CN 104915640 B CN104915640 B CN 104915640B CN 201510266649 A CN201510266649 A CN 201510266649A CN 104915640 B CN104915640 B CN 104915640B
- Authority
- CN
- China
- Prior art keywords
- structural plane
- value
- signal
- tendency
- inclination angle
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Fee Related
Links
- 238000000034 method Methods 0.000 title claims abstract description 56
- 239000002131 composite material Substances 0.000 claims abstract description 30
- 238000000605 extraction Methods 0.000 claims abstract description 21
- 238000006243 chemical reaction Methods 0.000 claims abstract description 9
- 238000005192 partition Methods 0.000 claims abstract description 6
- 238000012216 screening Methods 0.000 claims abstract description 6
- 230000003044 adaptive effect Effects 0.000 claims abstract description 3
- 230000008569 process Effects 0.000 claims description 19
- 101150106671 COMT gene Proteins 0.000 claims description 16
- 230000008878 coupling Effects 0.000 claims description 16
- 238000010168 coupling process Methods 0.000 claims description 16
- 238000005859 coupling reaction Methods 0.000 claims description 16
- 238000009826 distribution Methods 0.000 claims description 16
- 230000008859 change Effects 0.000 claims description 12
- 230000000877 morphologic effect Effects 0.000 claims description 6
- 238000005553 drilling Methods 0.000 claims description 4
- 239000000203 mixture Substances 0.000 claims description 4
- 238000000205 computational method Methods 0.000 claims 1
- 238000012795 verification Methods 0.000 claims 1
- 230000006870 function Effects 0.000 description 15
- 238000004458 analytical method Methods 0.000 description 9
- 239000011435 rock Substances 0.000 description 7
- 239000000284 extract Substances 0.000 description 4
- 238000012545 processing Methods 0.000 description 4
- 230000006978 adaptation Effects 0.000 description 3
- 238000004364 calculation method Methods 0.000 description 3
- 238000005516 engineering process Methods 0.000 description 3
- 230000013011 mating Effects 0.000 description 3
- 230000008901 benefit Effects 0.000 description 2
- 238000013461 design Methods 0.000 description 2
- 238000011161 development Methods 0.000 description 2
- 230000018109 developmental process Effects 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 230000002708 enhancing effect Effects 0.000 description 2
- 238000001914 filtration Methods 0.000 description 2
- 230000001788 irregular Effects 0.000 description 2
- 230000009131 signaling function Effects 0.000 description 2
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 description 1
- 241000208340 Araliaceae Species 0.000 description 1
- 240000007594 Oryza sativa Species 0.000 description 1
- 235000007164 Oryza sativa Nutrition 0.000 description 1
- 235000005035 Panax pseudoginseng ssp. pseudoginseng Nutrition 0.000 description 1
- 235000003140 Panax quinquefolius Nutrition 0.000 description 1
- 230000004888 barrier function Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000015572 biosynthetic process Effects 0.000 description 1
- 238000004422 calculation algorithm Methods 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 238000007796 conventional method Methods 0.000 description 1
- 238000005315 distribution function Methods 0.000 description 1
- 235000008434 ginseng Nutrition 0.000 description 1
- 238000010191 image analysis Methods 0.000 description 1
- 238000003709 image segmentation Methods 0.000 description 1
- 238000012804 iterative process Methods 0.000 description 1
- 238000007726 management method Methods 0.000 description 1
- 239000000463 material Substances 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 238000004321 preservation Methods 0.000 description 1
- 235000009566 rice Nutrition 0.000 description 1
- 238000000926 separation method Methods 0.000 description 1
- 238000007619 statistical method Methods 0.000 description 1
- 238000003786 synthesis reaction Methods 0.000 description 1
- 230000009466 transformation Effects 0.000 description 1
- 230000001131 transforming effect Effects 0.000 description 1
- 230000000007 visual effect Effects 0.000 description 1
- 230000003245 working effect Effects 0.000 description 1
- 238000013316 zoning Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
- G06F2218/02—Preprocessing
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
- G06F2218/08—Feature extraction
Landscapes
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Artificial Intelligence (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Signal Processing (AREA)
- Physics & Mathematics (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Geophysics And Detection Of Objects (AREA)
- Image Analysis (AREA)
- Image Processing (AREA)
Abstract
本发明公开了一种钻孔图像结构面自动识别与参数提取方法,步骤是:A、描述特征信号:在钻孔图像中,提取每行像素点中的特征值及其对应的点坐标;B、划分结构面区域:根据合成信号的整体均值与局部均值构成自适应区域变化的阈值,并二值化合成信号以进行区域划分;C、匹配模板正弦函数:根据已知模板正弦函数,依次迭代匹配模板正弦函数;D、筛选正弦曲线:根据每一个位置构成的正弦曲线匹配值信号的最大值点或极大值点来筛选正弦曲线;E、再匹配与精确化:根据已获的位置、倾向、倾角,在特定波动范围内进行再匹配筛选;F、参数提取与转化。方法简单易行,稳定可靠,操作简便,工作效率高,一次性完成所有结构面的自动识别与参数提取。
Description
技术领域
本发明属于岩体工程与图像识别技术领域,更具体涉及一种针对钻孔图像结构面的自动识别与参数提取方法,突破了以往需要人工干预和操作的技术难题,提高了钻孔图像分析处理的工作效率,推动了钻孔图像识别技术的发展,适用于钻孔摄像中结构面的全自动化处理与智能化分析。
背景技术
在地质工程、水电工程、石油开发和地质灾害防治工程方面,常常需要了解岩石结构面的形态特征。在岩石结构面中,这些形态特征的自动化识别与倾向倾角等参数的智能提取对于勘探、工程设计与施工等具有重要的实际意义。
目前,结构面形态特征的判读与提取基本上停留在人工操作识别的基础上。对于结构面的识别,通常由人工给出结构面曲线的控制点,进而由计算机拟合出标准的正弦曲线,确定拟合曲线与图像曲线互相符合后,再根据拟合参数重新计算出结构面的几何参数。这个识别操作的过程存在较多的人为参与环节,不仅仅工作效果低,而且不同的操作人员可能会给出不同的识别结果,造成识别结果的多样性,可靠性较低。另外,通过图像处理提取结构面的方法,采用图像空间变换、图像滤波与增强、图像分割、Hough变换与参数提取等步骤可以实现结构面参数的唯一判读。然而,该方法仍然需要人为输入参数的引导和必要的人工干预,无法实现钻孔图像结构面参数的全自动化连续判读与智能识别。特别是在深至几千米的钻孔图像中,钻孔图像数据量大,结构复杂,人工操作费时费力,难以短期内完成。因此,结构面形态特征的全自动识别与参数提取具有重要的实际意义,迫切需要解决。另外,如何在钻孔图像中寻找新的特征量来辅助结构面的识别判读,也是解决这个问题的关键。
于是,本发明提出了一种全新的结构面自动识别与参数提取方法,用以解决钻孔图像中结构面全自动化判读与倾向倾角等参数的智能提取。针对利用钻孔摄像技术获取的钻孔图像,本发明所述方法充分利用钻孔图像本身的灰度特征及其分布特性,分别了提取每行的最大灰度值,最小灰度值和最大梯度变化值,并组成全新的特征量信号,用于辅助结构面的识别。再根据结构面是正弦曲线的固有特性,采用模板正弦函数依次进行匹配运算,筛选出最优最符合结构面的正弦曲线及其对应的参数,实现了钻孔图像结构面的自动识别与参数提取。
发明内容
针对背景技术存在的问题,本发明的目的是在于提供了一种钻孔图像结构面自动识别与参数提取的方法。该方法简单易行,稳定可靠,可以实现全自动化处理,无需人为因素干预,极大地提高了工作效率,一次性准确地完成率所有结构面的识别与参数的提取。
为解决上述技术问题,本发明采用如下技术方案:
一种钻孔图像结构面自动识别与参数提取方法,其步骤是:
A、描述特征信号:在钻孔图像中,提取每行像素点中的最小灰度值点,最大灰度值点,灰度梯度最大值点,分别组成一维的最小值信号MinV(i),最大值信号MaxV(i),最大梯度值信号MaxG(i);再利用MinV(i)、MaxV(i)和MaxG(i)信号组成最终特征量信号即合成信号ComS(i);
钻孔图像中像素点的灰度值最小的叫最小灰度值,灰度值最小的点叫最小灰度值点,每行中的最小灰度值点组成最小值信号;具体计算方法见具体实施方式中的公式(1)及其相关描述。
钻孔图像中像素点的灰度值最大的叫最大灰度值,灰度值最大的点叫最大灰度值点,每行中的最大灰度值点组成最大值信号;具体计算方法见具体实施方式中的公式(2)及其相关描述。
钻孔图像中该像素点的灰度值与相邻像素点的灰度值的差值记为梯度值,梯度值最大的叫最大梯度值,梯度值最大的点叫最大梯度值点,每行中的最大梯度值点组成最大梯度值信号;具体计算方法见具体实施方式中的公式(6)及其相关描述。
B、划分结构面区域:根据合成信号的整体均值与该区域的局部均值的平均值来组成自适应区域变化的阈值,用以二值化合成信号ComS(i),得到分块信号ComT(i);再根据分块信号ComT(i)中每块的宽度W和钻孔孔径D得到结构面倾角β的范围为0到2W/D度,根据特征量像素点的X坐标分布情况得到结构面倾向α的范围为α1到α2度。(具体计算方法见具体实施方式中的公式(8)到(10)及其相关描述。)
C、匹配模板正弦函数:针对步骤B所得的每个分块区域和倾向倾角范围,根据已知的模板正弦函数y(x),每一个位置(每行),在其限定的倾向倾角范围内,倾向倾角以StepL1(StepL1≥1°)度为步长,依次迭代匹配模板正弦函数;以该分块区域的灰度平均值Tgray(i)作为该块的划分阈值,根据每个像素点的灰度值大小给予不同的权值Qpoint(i,j);然后,分别统计模板正弦函数所经过区域及其附近区域像素点个数的权值总和TempNpoint(i),并选择总权值最大的模板正弦函数作为该位置结构面的正弦曲线,每个位置的最大权值构成结构面正弦曲线的匹配值信号MaxN(i);(具体计算方法见具体实施方式中的公式(11)到(13)及其相关描述。)
D、筛选正弦曲线:根据在每一个位置正弦曲线信号MaxN(i)的分布情况,选择MaxN(i)信号的最大值点所对应的正弦曲线作为该结构面的正弦曲线,并保存该正弦曲线的位置、倾向和倾角。(MaxN(i)信号中值最大的点,详情见图6,请参考具体实施方式。)
E、再匹配与精确化:根据步骤D已获得的结构面位置、倾向、倾角,在当前位置、倾向、倾角的基础上分别减一个FluctV(FluctV≥1)值到再到分别加一个Fluct值的波动范围内(即从位置y-FluctV厘米到位置y+FluctV厘米,从倾向α-FluctV度到倾向α+FluctV度,从倾角β-FluctV度到倾角β+FluctV度),依次改变其位置、倾向、倾角,并以StepL2(StepL1>Step2>0)度为步长,按照步骤C和步骤D,再一次迭代匹配正弦曲线,统计加权像素点个数,得到更优的结构面正弦曲线及其位置、倾向和倾角。
F、参数提取与转化:把步骤E得到的结构面正弦曲线位置、倾向、倾角,结构面正弦曲线位置、倾向、倾角作为该结构面的位置、倾向、倾角,再根据基准位置和图像尺寸转化为过程实践中所需要的结构面深度、倾向、倾角。
本发明与现有技术相比,具有以下优点和有益效果:
1、本发明提供了一种全新的结构面特征量描述方法。在原始的钻孔图像中,直接提取钻孔图像中每行的特征量信号(最大灰度值、最小灰度值、最大灰度梯度变化值及其所在位置的坐标,也就是特征像素点坐标)来描述结构面,规避了图像滤波、增强、锐化等效益不大的传统图像预处理方法,简化了信号处理过程,准确描述了结构面的形态,为结构面的自动识别提供了判断依据和有效条件。
2、本发明中结构面特征量的提取与利用,缩小了匹配范围和迭代过程,提高了运行速度。用特征量信号来准确划分结构面分布区域,确定结构面正弦曲线倾向倾角的分布范围与相对中心位置,极大地优化了该方法的迭代匹配过程,简化了计算量,提高了运算速度。
3、本发明直接采用正弦函数为模板,省略了传统方法的曲线拟合过程。为了充分利用先验知识(标准结构面在展开后的钻孔图像中是标准正弦曲线)来回避图像处理中盲目的目标识别。
4、本发明所得的结构面参数的准确性高,并且稳定可靠、鲁棒性好。利用模板正弦函数在所有可能的倾向倾角范围内依次进行迭代匹配操作,遍历了所有可能性,筛选出了最优的模板正弦函数,并且优中选优,步骤明确清楚,方法可靠、稳定、鲁棒性好。
5、本发明所述方法自适应性较好。在迭代匹配的过程中,阈值的选取始终兼顾局部与整体,自适应岩体结构面图像的像素点灰度变化情况,适用于不同岩体的结构面,有效地保证了结构面正弦曲线的自适应判读与目标识别。
6、本发明所述方法提取的参数实用性较好。提取的结构面参数精度可以根据需要,进行再匹配与优化,并自动转化为实际过程中所需要的数据格式。如需要倾向倾角的精度达到0.01度,那么可以设定匹配过程的StepL2步长可以设定为0.01度,然后依次递增继续匹配。
7、本发明所述方法自动化程度高,便于实现全自动化处理,无需人员干预。在人工确认选定分析区域之后,该方法可以全自动完成整个钻孔图像中所有结构面的自动识别与参数提取,并完成所有数据的保存与管理,简化了操作过程,极大地提高了工作效率。
8、本发明所述方法运行速度快,工作效率高。比如:在100米深的钻孔图像中,该方法只需要连续运行5个小时左右的时间就可以一次性完成所有结构面的自动识别与参数提取,而通常的人工判读方式则需要几天时间,并且很难保质保量地完成所有结构面的识别与参数的精确提取。
另外,本发明所述步骤均可以通过专业软件(比如Visual C++,Matlab等)的编程实现,例如通过C++语言编程,可以编写实现以上步骤的程序代码,计算机在得到钻孔图像之后依据以上步骤的程序代码即可实现图像中结构面的自动识别与参数提取。
附图说明
图1为一种钻孔图像结构面自动识别与参数提取方法的示意图。
其中:1-获取钻孔图像、2-描述和提取钻孔图像中特征信号、3-根据特征信号划分结构面区域、4-在每个划分区域内匹配模板正弦函数、5-根据匹配点数筛选出点数最多的正弦曲线、6-根据已得参数和需要在小范围内重新再匹配正弦曲线、7-提取正弦曲线的有关参数并转化实际需要;
图2为一种结构面在展开后的钻孔图像中的形态。
其中:a为地质学上的标准结构面在展开后的钻孔图像中所存现的正弦曲线;b为地质学上的不规则结构面在展开后的钻孔图像中所存现的不规则曲线。
图3为一种数据ZK1中3m到4.5m中的一段钻孔图像。
其中:有4个结构面(正弦曲线带,灰度值较暗的地方),其中1个结构面非常明显突出,表现为黑色曲线,在图像最下面;其中1个结构面比较狭窄,在上面的第二个位置。
图4为一种钻孔图像中提取出的特征信号与实际结构面的对比图。
其中:上半部分曲线上凸的地方与下半部分结构面相对应,4个凸点对应于4个结构面。
图5为一种钻孔图像中的分块信号与分块结果图。
其中:上下两条水平线之间限定结构面的范围。
图6为一种匹配模板正弦函数中倾向倾角的参考信号。
其中:方波曲线为分块信号,最右边波动较大的信号为X坐标分布信号。
图7为一种筛选信号MaxN(i)的分布情况图。
其中:图中总共显示了6个结构面区域的信息,前面4个为图2的4个结构面信息,后面两个在图中没有显示。
图8为一种钻孔图像结构面的自动识别结果(表3中编号1~4的效果图)。
具体实施方式
实施例1:
一种钻孔图像结构面自动识别与参数提取方法,在已经获取钻孔图像的前提下的具体步骤如下,并以实测钻孔图像数据为例进行详细说明。
一种钻孔图像结构面自动识别与参数提取方法,其步骤是:
(1)描述特征信号:
在钻孔图像中,理想状态下结构面在钻孔图像中是标准的正弦曲线带,非理想状态下结构面在钻孔图像中也表现为类似正弦曲线状的曲线带,如图2所示。这些正弦曲线带的灰度值分布比较明显突出,表现为暗淡或者明亮,记为目标区域,如图3所示。在非目标区域,由于孔壁的不平整和岩石材质的差异,图像中不可避免地存在许多明暗不均、大小不一的孤立光斑点。这些光斑点在很大程度上影响了目标区域的定位与识别。但是,在目标区域中结构面曲线带的灰度值明显相对连续地最黑或者相对连续地最亮,并且分布区域相对集中连续。另外,结构面曲线带在垂直方向的投影相对集中。如果结构面之间的间隔较大,则结构面曲线带在垂直方向上的投影是断续的一小段,如图4所示。
为了更好地进一步凸显这些特征,避免被其他无关点覆盖,本发明采用每行像素点灰度值的最小值,最大值来描述结构形态特征。设图像中左上角为像素坐标的起点,像素点坐标从上往下从左往右依次增加。在钻孔图像中,选择每一行像素点的最小灰度值作为该行的最小特征值,表达式为公式(1);选择每一行像素点的最大灰度值作为该行的最大特征值,表达式为公式(2);其中,MinV(i)为第i行的最小像素值,MaxV(i)为第i行的最大像素值,N为钻孔图像的宽度,是一个常数值,在本发明中N为常数1024。
MinV(i)=min[f(i,0),f(i,1),…,f(i,j),…,f(i,N-1)](i∈N) (1)
MaxV(i)=max[f(i,0),f(i,2),…,f(i,j),…,f(i,N-1)](i∈N) (2)
为了描述结构面曲线带的灰度与周边灰度的变化情况,本发明采用灰度梯度来描述。记图像为f(x,y),f(i,j)表示像素点在(i,j)处的灰度值。在3×3的邻域内,在点(i,j)处,水平方向(x轴)的梯度如公式(3)所示,垂直方向(y轴)的梯度如公式(4)所示。在大量的计算处理当中,计算机常常采用x方向梯度y方向梯度的绝对值的和作为梯度的值,如公式(5)所示。其中,Gx(i,j)为像素点(i,j)在x方向的梯度值,Gy(i,j)为像素点(i,j)在y方向的梯度值,G(i,j)为像素点(i,j)处的梯度值。
于是,每一行像素点的最大梯度值可以表示为公式(6),其中,MaxG(i)为第i行的最大梯度值。
MaxG(i)=max[G(i,0),G(i,2),…,G(i,j),…,G(i,N-1)] (6)
在钻孔图像中,每行的最小灰度值和最大灰度值代表了图像特征的两个极点,而结构面正弦曲线带就是由这些极点(钻孔图像中的黑点和亮点)组成的。这些极点一般分布在结构面正弦曲线带中,并且在垂直方向的投影相对集中(正弦曲线是水平分布的)。为了避免这些极点被非目标区域的光斑点淹没,使他们进一步突显出来,本发明采用每行的最大值MaxV(i)与最小值MinV(i)的差值作为该行的基准值,再加上该行的最大梯度值MaxG(i)或者其倍数作为该行的合成信号特征值ComS(i),其表达式如公式(7)所示。
ComS(i)=MaxV(i)-MinV(i)+λ·MaxG(i)(λ≥1) (7)
其中,ComS(i)为第i行的特征值,λ在本发明中可取1.5。在钻孔图像的整个分析区域中,最大值信号MaxV(i)、最小值信号MinV(i)和最大梯度值信号MaxG(i)构成了合成信号ComS(i)。合成信号ComS(i)用来表示该行像素点的灰度特征。该特征信号有效地表达和凸显了图像中的结构面分布情况,如图4所示。
(2)划分结构面区域:
由图4可知,合成信号ComS(i)有效地表达了结构面正弦曲线的分布区域和相对中心位置。为了更好进行计算机处理和分析,有必要对合成信号进行阈值分割和二值化处理,精确表达结构面位置。在实际钻孔图像中当中,由于钻孔深至几千米,岩石结构会发生变化从而导致钻孔图像本身的背景条件发生很大变化,故图像中的阈值选取往往需要兼顾局部与整体,用以自适应钻孔图像的整体变化。合成信号的整体阈值取整个分析区域的平均值,如公式(8)所示。另外,再取合成信号在当前分析区域中10~100行的均值作为该区域的局部阈值,本发明本次运算可以取50,如公式(9)所示。
其中,M为整个合成信号的个数,也就是整个分析区域的行数,为合成信号的整体阈值;w为该行上下偏移的一个宽度,一般取20行即可,为合成信号的局部阈值。于是合成信号的二值化方法如公式(10)所示,式中ComT(i)信号为分块信号。该信号是合成信号的二值化结果。
计算机根据分块信号ComT(i)可以很方便地实现结构面正弦曲线的定位与分析。当分块信ComT(i)非0时,记为结构面存在的区域。实际当中可能会存在一些狭小的区域,为了避免丢失一些可能存在信息,这些区域可以合并到邻近相对较大的区域当中。钻孔图像和分块信号和区域划分的结果如图5所示。其中,红色区域线为分块结果的起始和终止标志线。
该方法根据合成信号的分布特征得出了便于计算机处理的分块信号。根据分块信号,该方法有效地划分了每个结构面所在区域,也定位了每个结构面的相对中心位置。结构面所在区域的大小意味着结构面的最大倾角。结构面所在的中心位置和特征量极值点的坐标也就意味着正弦曲线的峰峰值和谷谷值的相对位置,从而得出结构面倾向的大致范围。这些信息的进一步的提取与利用,为目标区域内结构面特征曲线(正弦曲线)的识别与分析提供了判断依据,也简化了匹配运算过程。
(3)匹配模板正弦函数:
由于得到的钻孔图像是展开的,标准结构面在钻孔图像中是一条标准的正弦曲线。故本发明直接采取正弦函数为结构面正弦曲线的标准模板,对每个区域内的每一行,依次改变正弦函数的倾向倾角,计算每次正弦函数所经过的区域,并计算该区域以及附近的像素点的灰度值,统计满足该模板正弦函数阈值的像素点个数,然后选取最大个数的模板正弦函数作为该行的可能正弦曲线。即在每一个位置(每一行),在其限定的倾向倾角范围内,倾向倾角以StepL1(StepL1≥1)度为步长,依次迭代匹配模板正弦函数,然后再所有的行中筛选出最优的正弦曲线。
本实例中StepL1取值为1度。模板正弦函数如公式(11)所示。其中,D为钻孔孔径的直径大小,α为结构面的倾向,β为结构面的倾角。再根据分块信号ComT(i)中每块的宽度W和钻孔孔径D得到结构面倾角β的范围为0到2W/D度,根据特征量像素点的坐标分布情况X分布得到结构面倾向α的范围为α1到α2度,如图6所示。
根据分块信号,结构面所在区域的大小可以限定结构面的最大倾角,也就是模板正弦函数中的β角。再根据结构面所在的中心位置和特征量极值点的坐标可以初步划分出正弦曲线的峰峰值和谷谷值的相对位置,从而得出结构面倾向的大致范围,也就是模板正弦函数中的α角。这些信息的进一步的提取与利用,从而缩小模板正弦函数的倾向倾角的范围,大大减少匹配过程中的计算量,提高了算法的执行效率与有效性。
在使用标准模板正弦函数来匹配结构面时,为了判断模板函数所经过区域以及附近的像素点是否属于该结构面,本发明采用了一个阈值来划分。该阈值利用本区域像素点灰度值的平均值来表达该区域的特征,用来判断模板正弦函数所经过的像素点是否为结构面上的点,计算公式如式(12)所示。其中m为该分块区域的行数,n为该分块区域的列数,i0为该分块区域的起始行,j0位该分块区域的起始列,Tgary为判断阈值。
在实际当中,为了充分利用结构面中像素点的灰度特征。该方法给灰度值划分了几个等级,并给予了不同的权值。在统计匹配模板正弦函数像素点个数的时候,不同等级灰度值的像素点数加权的权值也不一样。像素点数的加权值函数如公式(13)所示。其中,Qpoint(i,j)为像素点(i,j)处的加权值。
在使用模板正弦函数y(x)进行像素点灰度值匹配时,需要统计满足阈值条件的总加权像素点个数。统计方法如公式(14)所示。其中,TempNpoint(i)为图像中第i行中某个倾向倾角状态下的模板正弦函数匹配到的总加权像素点个数,N在本发明中始终是1024.
由公式(14)可以得到每个模板正弦函数在每一行中的总加权像素点个数的匹配情况,选择总加权像素点个数最大的一个模板正弦函数所对应的正弦曲线作为该行的正弦曲线,并记录该正弦曲线的倾向、倾角、位置以及总加权像素点个数等参数。该正弦曲线下的总加权像素点个数随着行数的变化关系记为函数MaxN(i)。该函数MaxN(i)记录了分块区域内每行模板正弦曲线的总加权像素点数(记为匹配点数)的最大值。其中一段钻孔图像的分析结果如下表1所示。
表1结构面正弦曲线有关参数(其中一段数据)
序号 | 位置 | 倾向 | 倾角 | 匹配点数 | 序号 | 位置 | 倾向 | 倾角 | 匹配点数 |
8 | 130 | 66 | 8 | 273 | 28 | 150 | 90 | 6 | 775 |
9 | 131 | 153 | 8 | 293 | 29 | 151 | 87 | 5 | 768 |
10 | 132 | 145 | 8 | 308 | 30 | 152 | 81 | 5 | 753 |
11 | 133 | 142 | 8 | 335 | 31 | 153 | 74 | 5 | 738 |
12 | 134 | 139 | 8 | 354 | 32 | 154 | 70 | 5 | 647 |
13 | 135 | 79 | 8 | 405 | 33 | 155 | 73 | 4 | 615 |
14 | 136 | 83 | 8 | 419 | 34 | 156 | 66 | 4 | 613 |
15 | 137 | 85 | 8 | 435 | 35 | 157 | 60 | 4 | 600 |
16 | 138 | 94 | 8 | 475 | 36 | 158 | 53 | 4 | 545 |
17 | 139 | 95 | 8 | 523 | 37 | 159 | 49 | 4 | 538 |
18 | 140 | 100 | 8 | 613 | 38 | 160 | 44 | 4 | 533 |
19 | 141 | 101 | 8 | 738 | 39 | 161 | 39 | 4 | 517 |
20 | 142 | 103 | 8 | 835 | 40 | 162 | 34 | 4 | 506 |
21 | 143 | 103 | 7 | 786 | 41 | 163 | 28 | 5 | 487 |
22 | 144 | 108 | 7 | 870 | 42 | 164 | 24 | 4 | 477 |
23 | 145 | 108 | 6 | 1026 | 43 | 165 | 20 | 5 | 461 |
24 | 146 | 109 | 6 | 961 | 44 | 166 | 14 | 4 | 449 |
25 | 147 | 103 | 5 | 933 | 45 | 167 | 11 | 5 | 437 |
26 | 148 | 117 | 6 | 859 | 46 | 168 | 7 | 5 | 438 |
27 | 149 | 95 | 6 | 835 | 47 | 198 | 173 | 6 | 85 |
(4)筛选正弦曲线:
从表1可知,在分块区域内,每一行基本上都找到了一个匹配个数最多的结构面正弦曲线。在这些正弦曲线中,最符合结构面的正弦曲线就在其中。很明显,匹配点数最大的地方很可能是结构面的中心。在最大值附近,匹配点数变化最快的地方也有可能是结构面的边界区域。据此可以筛选出最适合的正弦曲线来代表当前的结构面。如图7所示,可以看出该段MaxN(i)信号存在6个区域块,每块选出一个匹配点数最大的位置作为该区域的结构面位置。该段分析区域的筛选结果表2所示。表2中编号为1的一行是表1的筛选结果,编号2-6来源与表1中未能列出的数据。
表2结构面自动识别有关参数
编号 | 像素位置 | 倾向 | 倾角 | 匹配点数 |
1 | 145 | 108 | 6 | 1026 |
2 | 213 | 111 | 6 | 437 |
3 | 298 | 119 | 6 | 636 |
4 | 387 | 118 | 8 | 2083 |
5 | 466 | 120 | 5 | 324 |
6 | 565 | 105 | 4 | 1315 |
(5)再匹配与精确化:
在上述过程中,为了减少计算量加快匹配过程,模板正弦函数的倾向倾角都可以每StepL1(StepL1=1)度为单位,依次从0度增加到360度和90度,或在倾向α(α1,α2)和倾角β[0,2W/D)的范围内依次递增循环匹配。故上述得到的结构面参数的倾向倾角的精度是1度。
为了使结构面参数更加精确,有必要在上述已得条件的基础上进行再一次匹配。根据已获得的结构面位置、倾向、倾角,在当前位置、倾向、倾角的基础上分别减一个FluctV(FluctV≥1)值到再到分别加一个Fluct值的波动范围内,即从位置y-FluctV厘米到位置y+FluctV厘米,从倾向α-FluctV度到倾向α+FluctV度,从倾角β-FluctV度到倾角β+FluctV度。然后,再依次改变其位置、倾向、倾角,并以StepL2(StepL1>Step2>0)度为步长,按照以上步骤,再一次迭代匹配正弦曲线,统计加权像素点个数,得到更优的结构面正弦曲线及其位置、倾向和倾角。
再匹配的步长可以以实际精度为基础,比如实际精度要精确到0.1,那么本实例中StepL2就以0.1度为单位,取StepL2=0.1,FluctV=2,在位置从y-FluctV厘米到y+FluctV厘米,倾向从α-FluctV度到α+FluctV度,倾角从β-FluctV度到β+FluctV度的波动范围内,依次以StepL2=0.1为步长改变其位置、倾向和倾角,按照以上步骤(步骤3和步骤4),再一次迭代匹配正弦曲线,统计加权像素点个数,筛选结构面正弦曲线,使得到的结果参数更加精确化,以满足实际要求。本发明对上述过程进行了一次再匹配与精确化操作,最后得到的结果面参数如表3所示。根据匹配参数最后画出来的结构面正弦曲线如图8所示,其中每条灰色的正弦曲线正好匹配了各自的结构面。
表3再匹配与筛选结果数据
编号 | 像素位置 | 倾向 | 倾角 | 匹配点数 |
1 | 145 | 108.7 | 5.9 | 1064 |
2 | 213 | 111.2 | 6.1 | 471 |
3 | 298 | 119.4 | 6.4 | 673 |
4 | 387 | 117.8 | 8.6 | 2194 |
5 | 466 | 120.5 | 5.3 | 352 |
6 | 565 | 104.6 | 4.7 | 1347 |
(6)参数提取与转化:
经过以上过程中,最后得到的结构面正弦曲线的位置是当前分析区域的像素坐标。该像素坐标位置经过图像基准位置(起点位置)的转化之后可以得到结构面的实际位置即深度位置。
经过以上过程中,最后得到的结构面正弦曲线的倾向等效于正弦曲线的相位。该相位的角度与钻孔图像中正北方向的位置的夹角即为地质学上结构面的倾向,并可以适当转化为人们更加便于理解接受的角度。结构面的倾向可以等于钻孔图像中正弦曲线的相位加上270度,结构面倾向的取值范围为0至360度;
经过以上过程中,最后得到的结构面正弦曲线的倾角也就是正弦曲线的幅值。该幅值与结构面倾角的转化关系如公式(15)所示。其中,D为钻孔孔径的直径大小,Aimplitude为正弦曲线的幅值,β为结构面的倾角。结构面的倾角等于钻孔图像中正弦曲线的幅值除以钻孔半径的反正切函数值,
本发明根据图像的基准位置和图像的尺寸大小与实际位置的对应关系把有关参数转化为实践过程中所需要的结构面深度、倾向、倾角,并把数据同时保存在数据库和有关数据表格中以方便查阅和理解。该转化过程满足了实际工程的需要。本发明对上述过程进行了一次参数提取与转化,最后得到的结构面有关参数如表4所示。
表4结构面自动识别结果有关参数
编号 | 深度 | 像素位置 | 匹配点数 | 倾向 | 倾角 |
1 | -3218 | 145 | 1064 | 161.3 | 47.8 |
2 | -3423 | 213 | 471 | 157.5 | 49.3 |
3 | -3678 | 298 | 673 | 150.2 | 51.7 |
4 | -3979 | 387 | 2194 | 152.6 | 56.4 |
5 | -4184 | 466 | 352 | 148.3 | 43.2 |
6 | -4476 | 565 | 1347 | 164.7 | 36.8 |
公式定义与有关符号说明:
f(i,j)表示钻孔图像中像素点在(i,j)处的灰度值,也表示像素点坐标函数。
MinV(i)为第i行的最小像素值,也为钻孔图像的最小灰度值信号分布函数。
MaxV(i)为第i行的最大像素值,也为钻孔图像的最大灰度值信号分布函数。
MaxG(i)为第i行的最大梯度值,也为钻孔图像的最大梯度值信号分布函数。
Gx(i,j)为像素点(i,j)在x方向的梯度值,也为图像中x方向的梯度函数。
Gy(i,j)为像素点(i,j)在y方向的梯度值,也为图像中y方向的梯度值函数。
G(i,j)为像素点(i,j)处的梯度值,也为图像中像素坐标点的梯度函数。
ComS(i)为第i行的合成信号值,也为钻孔图像的合成信号函数。
为钻孔图像中合成信号的整体均值。
为钻孔图像中合成信号的局部均值。
ComT(i)为第i行的分块信号值0或1,也为钻孔图像的分块信号函数。
W为当前块的宽度。
D为钻孔孔径。
α为结构面的倾向,也等同于正弦曲线的倾向。
β为结构面的倾角,也等同于正弦曲线的倾角。
StepL1为匹配模板正弦函数的步长,本实例中StepL1取值为1。
StepL2为正弦曲线再匹配过程中的步长,本实例中StepL2取值为0.1。
FluctV为位置、倾向、倾角的预设定波动范围,本实例中FluctV取值为2。
y(x)为模板正弦函数,包含了倾向倾角的限定范围。
Tgary(i)为第i行的判断阈值,也为分块区域内的判断阈值函数。
Qpoint(i,j)为像素点(i,j)处的加权值,也为整个区域内的加权函数。
TempNpoint(i)为图像中第i行中某个倾向倾角状态下的模板正弦函数匹配到的总加权像素点个数.
N在本发明中始终为常数1024。
实施例2:
下面以一般情况(没有特定钻孔数据)为例,对本发明作进一步说明。
一种基于钻孔图像的结构面全自动识别方法,包括以下步骤:
步骤1、利用钻孔摄像技术获取某地的钻孔图像。
步骤2、选定钻孔图像中的一段图像数据,或者特定的部分进行分析。
步骤3、提取该部分中每行像素点中的最小灰度值点,最大灰度值点,灰度梯度最大值点,根据公式(7)组成结构面特征量的合成信号。
步骤4、根据局部阈值与整体阈值来二值化合成信号,组成分块信号用来划分结构面区域,并根据分块信号宽度来限定结构面倾角的变化范围,根据分块信号在特征点的坐标分布情况来获得结构面倾向的变化范围。
步骤5、根据已知的模板正弦函数,在其限定的倾向和倾角范围内,以1度为步长,依次迭代匹配区域内每一行的模板正弦函数。以该分块区域的灰度平均值作为该块的划分阈值,根据该阈值统计模板正弦函数所经过位置及其附近区域所有满足阈值的像素点个数,并选择总个数最大的模板正弦函数作为该位置的结构面正弦曲线信号。
步骤6、由于区域内每一行都有一个个数最大的正弦曲线,而一个区域内一般只存在一个结构面。故在这些正弦曲线中再次选择总个数最大的正弦曲线作为该结构面的正弦曲线。
步骤7、根据前面步骤已获得的结构面位置、倾向、倾角等信息,在±2°的波动范围内,改变其位置、倾向和倾角,以0.1度为步长,按照步骤3和步骤4,再一次迭代匹配正弦曲线,统计加权像素点个数,再次选择一个更优的结构面正弦曲线。
步骤8、把以上步骤得到的正弦曲线的位置、倾向、倾角作为该结构面的位置、倾向、倾角,再根据基准位置和图像尺寸转化为过程实践中所需要的结构面深度、倾向、倾角等信息。
步骤9、把数据另存到自己的位置,或者备份到不被覆盖的存储空间去,重新选择分析区域,再进行结构面的自动识别与参数提取分析。
其中本实施例中的步骤3到步骤8与实施例1中的步骤1到步骤6基本相同。
相关概念定义与注释:本发明中所述结构面是地质学上的概念。结构面在钻孔图像上呈现正弦曲线状,称为结构面正弦曲线。结构面参数即位置、倾向、倾角分别与钻孔图像中结构面正弦曲线的位置、相位、幅值存在相互对应的关系,即:(1)结构面的位置与钻孔图像中结构面正弦曲线的位置对应,钻孔图像中正弦曲线的位置是像素坐标位置,经过图像基准位置(起点位置)的转化之后可以得到结构面的实际位置即深度位置;(2)结构面的倾向与钻孔图像中结构面正弦曲线的相位对应,结构面的倾向等于钻孔图像中正弦曲线的相位加上270度,结构面倾向的取值范围为0至360度;(3)结构面的倾角与钻孔图像中结构面正弦曲线的幅值对应,结构面的倾角等于钻孔图像中正弦曲线的幅值除以钻孔直径的反正切函数值,如附图2所示。总之,本发明所述的结构面参数(位置、倾向、倾角)和钻孔图像中正弦曲线的参数(位置、相位、幅值)与模板正弦函数的参数(位置、相位、幅值)的意义是一致的。为了描述方便和概念的模糊不清,钻孔图像中结构面正弦曲线与正弦函数的位置、相位、幅值分别统一简述为结构面的位置、倾向、倾角。另外,本发明中所述的钻孔图像等同于钻孔图像的展开图,本发明所述方法也是针对该类型的钻孔图像而言的。
Claims (2)
1.一种钻孔图像结构面自动识别与参数提取方法,其步骤是:
A、描述特征信号:在钻孔图像中,提取每一行中所有像素点灰度值的最小值和最大值以及像素点灰度梯度的最大值,并保存他们对应点的坐标,分别组成一维的最小值信号MinV(i),最大值信号MaxV(i)和最大梯度值信号MaxG(i);再利用MinV(i)、MaxV(i)和MaxG(i)信号组成最终结构面特征量信号即合成信号ComS(i);
合成信号ComS(i)采用的公式为:
ComS(i)=MaxV(i)-MinV(i)+λ·MaxG(i) (λ≥1) (1)
其中,λ为一个大于等于1的加权值参数,经过实际验证取1.5,MinV(i)为每行的最小值信号,MaxV(i)为每行的最大值信号,MaxG(i)为每行的最大梯度值信号;
B、划分结构面区域:根据合成信号ComS(i)的整体均值与该区域的局部均值的平均值来组成自适应区域变化的阈值,用以二值化合成信号ComS(i),并得到分块信号ComT(i);再根据分块信号ComT(i)中每块的宽度W和钻孔孔径D得到结构面倾角β的范围为0到2W/D度,根据合成信号特征量对应点的X坐标分布情况得到结构面倾向α的范围为α1到α2度;
其中,分块信号ComT(i)采用的公式(2)、公式(3)和公式(4)为:
其中,M为整个合成信号的个数,也是整个分析区域的行数,为合成信号的整体阈值;w为该行上下偏移的一个宽度,可以取20行,为合成信号的局部阈值;
C、匹配模板正弦函数:针对步骤B所得的每个分块区域和倾向、倾角范围,根据已知的模板正弦函数y(x),每一个位置,在其限定的倾向、倾角范围内,倾向、倾角以StepL1(StepL1≥1°)度为步长,依次迭代匹配模板正弦函数;以该分块区域的灰度平均值Tgray(i)作为该块的划分阈值,根据每个像素点的灰度值大小给予不同的权值Qpoint(i,j);然后,分别统计模板正弦函数所经过区域及其附近区域像素点个数的权值总和TempNpoint(i),并选择总权值最大的模板正弦函数作为该位置结构面的正弦曲线,每个位置的最大权值构成结构面正弦曲线的匹配值信号MaxN(i);
其中,模板正弦函数y(x)采用公式(5),分块区域灰度平均值Tgray(i)的计算方法采用公式(6),像素点灰度权值的分配关系Qpoint(i,j)采用公式(7),统计每行模板正弦函数所经过区域及其附近区域的像素点个数权值的总和TempNpoint(i)采用公式(8);
其中,y0为当前行的位置,记为正弦曲线中心的初始位置,α为结构面的倾向,β为结构面的倾角,D为钻孔的孔径,N为图像的列数,W为该块的宽度,m为该分块区域的行数,n为该分块区域的列数,i0为该分块区域的起始行,j0为该分块区域的起始列,Tgary为判断阈值;
D、筛选正弦曲线:根据在每一个位置正弦曲线匹配值信号MaxN(i)的分布情况,选择MaxN(i)信号的最大值点或极大值点所对应的正弦曲线作为该结构面的正弦曲线,并保存该正弦曲线的位置、倾向和倾角;
E、再匹配与精确化:根据步骤D已获得的结构面位置、倾向、倾角,在当前位置、倾向、倾角的基础上分别减一个FluctV(FluctV≥1)值到再到分别加一个Fluct值的波动范围内,依次改变其位置、倾向、倾角,并以StepL2(StepL1>Step2>0)度为步长,按照步骤C和步骤D,再一次迭代匹配模板正弦函数,统计加权像素点个数,得到精度更高的结构面正弦曲线及其位置、倾向和倾角;
F、参数提取与转化:把步骤E得到的结构面正弦曲线的位置、倾向、倾角作为该结构面的位置、倾向、倾角,再根据基准位置和图像尺度转化为实践过程中所需要的结构面深度、倾向、倾角。
2.根据权利要求1所述的一种钻孔图像结构面自动识别与参数提取方法,其特征在于:所述步骤D中每行筛选出一个像素点个数权值最大的模板正弦函数,并组成正弦曲线信号MaxN(i)信号;根据块内MaxN(i)的最大值或极大值再次筛选块内的结构面形态参数。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201510266649.2A CN104915640B (zh) | 2015-05-22 | 2015-05-22 | 一种钻孔图像结构面自动识别与参数提取方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201510266649.2A CN104915640B (zh) | 2015-05-22 | 2015-05-22 | 一种钻孔图像结构面自动识别与参数提取方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN104915640A CN104915640A (zh) | 2015-09-16 |
CN104915640B true CN104915640B (zh) | 2017-03-08 |
Family
ID=54084692
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201510266649.2A Expired - Fee Related CN104915640B (zh) | 2015-05-22 | 2015-05-22 | 一种钻孔图像结构面自动识别与参数提取方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN104915640B (zh) |
Families Citing this family (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN105241379A (zh) * | 2015-10-08 | 2016-01-13 | 莆田市荣兴机械有限公司 | 换气阀孔尺寸测量方法 |
CN107945272B (zh) * | 2017-11-09 | 2021-10-08 | 长江三峡勘测研究院有限公司(武汉) | 一种基于高清钻孔彩电的岩体结构面的搜索方法 |
CN108109157B (zh) * | 2017-12-18 | 2021-07-06 | 武汉大学 | 一种基于数字式全景钻孔图像的岩体评估分析方法 |
CN116579989B (zh) * | 2023-04-21 | 2024-05-14 | 中铁九局集团电务工程有限公司 | 一种基于深度相机的隧道打孔倾角修正方法 |
Family Cites Families (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
JP2622802B2 (ja) * | 1992-11-05 | 1997-06-25 | 株式会社奥村組 | 地質の調査方法 |
CN104502990A (zh) * | 2015-01-06 | 2015-04-08 | 中国地质大学(武汉) | 一种基于数码图像的隧道掌子面地质调查方法 |
-
2015
- 2015-05-22 CN CN201510266649.2A patent/CN104915640B/zh not_active Expired - Fee Related
Also Published As
Publication number | Publication date |
---|---|
CN104915640A (zh) | 2015-09-16 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN108416686B (zh) | 一种基于煤炭资源开发的生态地质环境类型划分方法 | |
CN104915640B (zh) | 一种钻孔图像结构面自动识别与参数提取方法 | |
CN107273608B (zh) | 一种油藏地质剖面图矢量化方法 | |
CN108038448A (zh) | 基于加权熵的半监督随机森林高光谱遥感影像分类方法 | |
CN106875481B (zh) | 一种三维可视化遥感影像地表分类模型的制作方法 | |
CN100595782C (zh) | 一种融合光谱信息和多点模拟空间信息的分类方法 | |
CN104620136A (zh) | 为创建储层属性模型组合多点静态和基于对象的方法的混合方法 | |
CN111475596A (zh) | 一种基于多层级轨迹编码树的子段相似性匹配方法 | |
CN104620135A (zh) | 用于编辑多点相模拟的方法 | |
CN104899965A (zh) | 一种基于清分机的多国纸币序列号识别方法 | |
CN111738055A (zh) | 多类别文本检测系统和基于该系统的票据表单检测方法 | |
CN101840582B (zh) | 一种地籍图地块的边界数字化方法 | |
CN108829711A (zh) | 一种基于多特征融合的图像检索方法 | |
CN110147778A (zh) | 稀土矿开采识别方法、装置、设备及存储介质 | |
Mohammadpour et al. | Automatic lineament extraction method in mineral exploration using CANNY algorithm and hough transform | |
CN104036294A (zh) | 基于光谱标记的多光谱遥感图像自适应分类方法 | |
Van de Voorde et al. | Extraction of land use/land cover related information from very high resolution data in urban and suburban areas | |
Lee et al. | Iterative static modeling of channelized reservoirs using history-matched facies probability data and rejection of training image | |
Chattaraj et al. | Semi-automated object-based landform classification modelling in a part of the Deccan Plateau of central India | |
CN112800590B (zh) | 一种机器学习辅助的两相流油藏随机建模的网格粗化方法 | |
CN103455816A (zh) | 一种笔画宽度提取方法、装置及一种文字识别方法、系统 | |
Čomić et al. | Morse-smale decompositions for modeling terrain knowledge | |
CN108716393B (zh) | 一种油砂sagd动用储量优选方法 | |
CN105373802A (zh) | 基于区间Type-2模糊支持向量机的场景图像分类方法 | |
CN114155440A (zh) | 一种耕地非农化的自动检测方法及系统 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
C10 | Entry into substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
C14 | Grant of patent or utility model | ||
GR01 | Patent grant | ||
CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20170308 |
|
CF01 | Termination of patent right due to non-payment of annual fee |