阅读说明：本技术 一种岩层孔隙率分布计算方法 (Rock stratum porosity distribution calculation method ) 是由 王少锋 皮滋滋 尹江江 周子龙 于 2021-08-02 设计创作，主要内容包括：本发明公开了一种岩层孔隙率分布计算方法,包括输出地层运动原始三彩图像；对地层运动原始三彩图像进行预处理,生成对应的三维图像数字矩阵；生成岩层裂隙层间分布图；沿岩层方向遍历岩层裂隙层间分布图,得到岩层裂隙信息提取区域图；并绘制岩层裂隙层孔隙率对应的分布曲线图；对岩层裂隙层间分布图进行全图像范围遍历,生成岩层孔隙率分布云图。本发明提供的这种岩层孔隙率分布计算方法,利用计算机对岩层图像的处理,减少了人工运算出现的遗漏或疏忽,避免了传统孔隙率计算方法中实验造成的有关误差,且在一定程度上加快了运算速度,从数据上直观呈现了岩层孔隙率在一定走向上的分布；同时减少了其他图像干扰,计算更加准确、高效。(The invention discloses a rock stratum porosity distribution calculation method which comprises the steps of outputting original three-color images of stratum movement; preprocessing the original three-color image of the stratum movement to generate a corresponding three-dimensional image digital matrix; generating a stratum fracture interlayer distribution map; traversing the stratum fracture interlayer distribution map along the stratum direction to obtain a stratum fracture information extraction area map; drawing a distribution curve chart corresponding to the porosity of the rock stratum fracture layer; and traversing the stratum fracture interlayer distribution map in a full-image range to generate a stratum porosity distribution cloud map. According to the rock stratum porosity distribution calculation method provided by the invention, the computer is used for processing the rock stratum images, so that omission or negligence caused by manual operation is reduced, related errors caused by experiments in the traditional porosity calculation method are avoided, the operation speed is accelerated to a certain extent, and the distribution of the rock stratum porosity in a certain trend is visually presented from data; meanwhile, the interference of other images is reduced, and the calculation is more accurate and efficient.)

1. A rock stratum porosity distribution calculation method is characterized by comprising the following steps:

s1, simulating a stratum fracture enrichment area, and outputting a corresponding stratum movement original three-color image;

s2, preprocessing the original three-color image of the stratum movement to generate a corresponding three-dimensional image digital matrix;

s3, removing redundant rock stratum pixel information by analyzing the three-dimensional image digital matrix to obtain a rock stratum fracture interlayer distribution map;

s4, traversing a stratum fracture interlayer distribution map along the stratum direction to obtain a stratum fracture representative pixel point and a stratum fracture information extraction area map; drawing a distribution curve chart corresponding to the porosity of the rock stratum fracture layer;

and S5, setting a full-image-range traversal window, performing full-image-range traversal on the stratum fracture interlayer distribution diagram, and generating a stratum porosity distribution cloud chart along the vertical and horizontal directions of the full-image-range traversal window.

2. The method for calculating porosity distribution of a rock formation according to claim 1, wherein the step S1 specifically includes: simulating the stratum of the stratum fracture enrichment area by using a discrete unit method, carrying out parameterization processing on the contact between a rock matrix and the rock to generate a stratum condition parameter, and outputting a corresponding stratum motion original tricolor image1 according to the stratum condition parameter; the original color map stores pixel values in a color mode of three color [ r, g, b ] channels.

3. A method of calculating a porosity distribution in a rock formation as claimed in claim 2, wherein said formation condition parameters include: the mass density, the elastic property and the plastic property of the block body are used for determining the mechanical property of the material; the rigidity, the shearing rigidity, the tensile strength, the friction angle and the cohesion are used for determining the mechanical characteristics of the contact between the blocks; the size of the formation block.

4. A rock formation porosity distribution calculating method according to claim 3, wherein in step S2, the preprocessing includes: generating a corresponding three-dimensional image digital matrix Array1 from the stratum movement original three-color image1, wherein the digital matrix Array1 reflects three color [ r, g, b ] channels in a three-dimensional form, respectively displays pixel positions in the corresponding channels in a two-dimensional Array form [ · ] i ] j, and reflects three original color pigment values through pixel values; wherein each pixel location contains three color channels, the red, green and blue tristimulus pixel values are distributed in the pixel location in the form of r, g, b, with r, g, b e [0,255 ].

5. The method for calculating porosity distribution of a rock formation according to claim 4, wherein the step S3 specifically includes the steps of:

A1. removing redundant upper grid lines in the stratum movement original tricolor image1 output in the step S1;

A2. the unnecessary lower grid lines in the stratum movement original tricolor image1 output in step S1 are removed.

6. The method for calculating porosity distribution of a rock formation according to claim 5, wherein the step A1 specifically includes: and taking the lower grid lines along the upper pixel rows and all columns as a reference, and setting:

w₁＝I(i,j,1)/I(i,j,2)

w₂＝I(i,j,2)/I(i,j,3)

wherein, w₁Expressing the ratio of the pixel value corresponding to the red channel of the pixel point to the pixel value of the green channel; w is a₂Expressing the ratio of the pixel value corresponding to the green channel of the pixel point to the pixel value of the blue channel; i (I, j,1) represents the pixel value of the pixel point in the red channel; i (I, j,2) represents the pixel value of the pixel point in the green channel; i (I, j,3) represents the pixel value of the pixel point in the blue channel;

if the grid line edge satisfies the upper setting range, I (I, j,1) is set to I (I, j,2) is set to I (I, j,3) is set to 0.

7. The method of calculating porosity distribution of a rock formation according to claim 6, wherein the step A2 includes setting, along the lower pixel row, all columns as references, the following network lines:

sum＝I(i,j,1)+I(i,j,2)+I(i,j,3)

if the grid line is satisfied with the downward setting range, let I (I, j,1) be I (I, j,2) be I (I, j,3) be 0; i (I, j,1) represents the pixel value of the pixel point in the red channel; i (I, j,2) represents the pixel value of the pixel point in the green channel; i (I, j,3) represents the pixel value of the pixel point in the blue channel; sum represents the sum of pixel values of the pixel points in red, green and blue channels.

8. The method for calculating porosity distribution of rock formation according to claim 7, wherein in step S4, the method for generating the formation fracture information extraction area image3 includes: traversing pixels in the stratum fracture interlayer distribution image2, wherein the form of a pixel two-dimensional array is [. cndot. ] i ] j, summing pixel values [ r, g, b ] of three color channels of pixel points, setting a threshold value T, and calculating:

sum＝I(i,j,1)+I(i,j,2)+I(i,j,3)

wherein, I (I, j,1) represents the pixel value of the pixel point in the red channel; i (I, j,2) represents the pixel value of the pixel point in the green channel; i (I, j,3) represents the pixel value of the pixel point in the blue channel; sum represents the sum of pixel values of the pixel points in red, green and blue channels; when sum > T, I (I, j,1) ═ I (I, j,2) ═ I (I, j,3) ═ 255; and when sum is less than or equal to T, I (I, j,1) ═ I (I, j,2) ═ I (I, j,3) ═ 0, and through traversal and calculation, an image formed by the region where the pixel point with the pixel value of 0 is the formation fracture information extraction region image 3.

9. The method for calculating porosity distribution of rock formation according to claim 8, wherein in step S4, the method for generating a distribution graph image4 corresponding to the porosity of fracture layer of rock formation comprises: setting the size of a traversal window as s x k, the step length as s, wherein s represents the length of the s x k window, k represents the width of the s x k window, and k takes a fixed value and represents the original width of a rock stratum fracture information extraction region image3 image for a rock stratum fracture information extraction region image 3; in an s x k traversal window, the porosity size Pr of a stratum fracture layer is as follows:

Pr＝n/N

wherein n is the number of pixels of the stratum fracture of the S x k traversal window; and N is the total number of the traversal pixel points of the traversal window.

10. The method of calculating porosity distribution in rock formation according to claim 9, wherein said step S5 includes: the window scans a stratum fracture information extraction area image3, the size of a traversal window is set to be m x h, and the step length is taken to be m; traversing the rock stratum fracture information extraction area image3 along the horizontal direction, and then traversing the rock stratum fracture information extraction area image3 along the vertical direction with the step length being m; in an m x h traversal window, the porosity size Pr' of a stratum fracture layer is as follows:

Pr'＝n'/N'

wherein n' is the number of pixel points of the stratum fracture of the traversal window m x h; n' represents the total number of the traversal pixel points of the m x h traversal window; and calculating the rock stratum fracture rate Pr' in each m x h traversal window, visualizing the porosity rate of each window through progressive colors, and gradually displaying the porosity rate represented by the traversal windows at different positions to obtain a distributed cloud image5 of the rock stratum porosity rate under the rock stratum fracture information extraction area image3 standard.

Technical Field

The invention belongs to the field of geotechnical engineering, and particularly relates to a rock stratum porosity distribution calculation method.

Background

Nowadays, with the improvement of the infrastructure level, the continuous development of geological development engineering faces many safety problems in geotechnical engineering. For example, a goaf after mining mineral resources is prone to collapse due to improper maintenance; engineering excavation causes water penetration accidents. The reasons for the above safety problems are generally that the formation has a large number of fractures in engineering practice, which changes the fluid flow properties of the formation; the developmental status of the fracture is generally calculated by porosity. Therefore, the method has important application significance in researching and calculating the porosity distribution of the stratum in rock engineering projects.

The porosity of the rock stratum is the ratio of the sum of various holes and fractures in the stratum to the total volume of the rock mass, is an important physical property of the rock stratum, and is important data information for researching rock dynamics. A plurality of researches on the calculation of the porosity distribution of the rock stratum are carried out at home and abroad. Generally, the method for calculating the porosity of the rock mass is to perform a caving experiment on the rock mass under the condition of neglecting the internal porosity of the rock, wherein the caving rock does not deform, and the porosity can be expressed by a formula: n is Vf/(Vs + Vf), and n is the porosity of the formation; vf is the volume of the rock mass pore, and is expressed by the difference of the volume of the rock mass before and after the caving; vs is the total volume of rock behind the rock mass fracture. In addition, measurements are made in the laboratory using rock samples or drill cores, and measurements of porosity values at a location in the formation are made using sonic and neutron density logs. Meanwhile, the prior art rock stratum pore distribution calculation method is as follows: CN107402176B discloses a method and a device for determining fracture porosity, wherein the method comprises: obtaining an outcrop sample and a rock core sample of a target area; performing three-dimensional laser scanning on the outcrop sample, and establishing a reservoir geological model of a target area; determining the extension length of the fracture through a reservoir geological model of a target area; determining the opening degree of the crack, the filling coefficient of the crack and the penetration coefficient of the crack according to the core sample; and determining the fracture porosity of the target area according to the extension length of the fracture, the opening degree of the fracture, the filling coefficient of the fracture and the penetration coefficient of the fracture. According to the scheme, the outcrop sample with a good representation effect is used for establishing a geological model of the reservoir, and the cracks are quantitatively determined; and the fracture porosity is accurately determined by combining the core data, so that the technical problems of low accuracy and poor accuracy of the determined fracture porosity in the conventional method are solved, and the technical effect of quickly and accurately quantitatively determining the fracture porosity is achieved. The CN105866002B invention discloses an accurate oil shale-containing nuclear magnetic resonance porosity testing method, which can accurately represent the physical property characteristics of an oil shale-containing reservoir through testing the nuclear magnetic resonance T2 spectral distribution and the porosity of saturated kerosene shale with different waiting time and echo intervals, and comparing the waiting time and the echo intervals with the helium porosity to calibrate the waiting time and the echo intervals, and the nuclear magnetic resonance T2 spectral distribution during the optimal waiting time and the echo intervals. The invention attributes the factors influencing the testing precision of the nuclear magnetic resonance porosity of the oil-containing shale to the type of the sample saturated fluid and testing parameters (waiting time and echo interval), and eliminates the influence of the hydration of the shale sample by saturated kerosene; calibrating the optimal waiting time and echo interval by analyzing the T2 spectrum distribution change and helium porosity; by adopting the method, the optimal test parameters capable of accurately representing the porosity of the oil-bearing shale can be obtained by testing the nuclear magnetic resonance T2 spectrum distribution and the porosity of the saturated kerosene shale with different waiting time and echo intervals.

The crack porosity calculation method and the nuclear magnetic resonance porosity test method in the prior art can analyze and calculate the crack and porosity distribution of the formation or the shelter to a certain extent, but the technology has some limitations. For example, the fracture porosity calculation method can determine multiple property parameters such as the extension length and the penetration coefficient of a fracture by sampling a target area and performing laser scanning, but the laser scanning method lacks accurate research on the overall view of the formation, lacks comprehensive generalization and has higher requirements on experimental instruments and operation; the nuclear magnetic resonance method compares the resonance spectrum distribution and the porosity with corresponding parameters of a calibration object to obtain the oil saturation and the porosity of the shale reservoir, but a certain error exists in the nuclear magnetic resonance method instead of a direct measurement means; meanwhile, the prior art is complex in mode and low in operation speed.

Disclosure of Invention

The invention aims to provide a rock stratum porosity distribution calculation method which can reduce calculation steps and has the advantages of wide application range, high efficiency and strong reliability.

The rock stratum porosity distribution calculation method provided by the invention comprises the following steps:

s1, simulating a stratum fracture enrichment area, and outputting a corresponding stratum movement original three-color image;

s2, preprocessing the original three-color image of the stratum movement to generate a corresponding three-dimensional image digital matrix;

s3, removing redundant rock stratum pixel information by analyzing the three-dimensional image digital matrix to obtain a rock stratum fracture interlayer distribution map;

s4, traversing a stratum fracture interlayer distribution map along the stratum direction to obtain a stratum fracture representative pixel point and a stratum fracture information extraction area map; drawing a distribution curve chart corresponding to the porosity of the rock stratum fracture layer;

and S5, setting a full-image-range traversal window, performing full-image-range traversal on the stratum fracture interlayer distribution diagram, and generating a stratum porosity distribution cloud chart along the vertical and horizontal directions of the full-image-range traversal window.

The step S1 specifically includes: simulating the stratum of the stratum fracture enrichment area by using a discrete unit method, carrying out parameterization processing on the contact between a rock matrix and the rock to generate a stratum condition parameter, and outputting a corresponding stratum motion original tricolor image1 according to the stratum condition parameter; the original color map stores pixel values in a color mode of three color [ r, g, b ] channels.

The formation condition parameters include: the mass density, the elastic property and the plastic property of the block body are used for determining the mechanical property of the material; the rigidity, the shearing rigidity, the tensile strength, the friction angle and the cohesion are used for determining the mechanical characteristics of the contact between the blocks; the size of the formation block.

In step S2, the preprocessing includes: generating a corresponding three-dimensional image digital matrix Array1 from the stratum movement original three-color image1, wherein the digital matrix Array1 reflects three color [ r, g, b ] channels in a three-dimensional form, respectively displays pixel positions in the corresponding channels in a two-dimensional Array form [ · ] i ] j, and reflects three original color pigment values through pixel values; wherein each pixel location contains three color channels, the red, green and blue tristimulus pixel values are distributed in the pixel location in the form of r, g, b, with r, g, b e [0,255 ].

The step S3 specifically includes the following steps:

A1. removing redundant upper grid lines in the stratum movement original tricolor image1 output in the step S1;

A2. the unnecessary lower grid lines in the stratum movement original tricolor image1 output in step S1 are removed. The step a1 specifically includes: and taking the lower grid lines along the upper pixel rows and all columns as a reference, and setting:

w₁＝I(i,j,1)/I(i,j,2)

w₂＝I(i,j,2)/I(i,j,3)

wherein, w₁Expressing the ratio of the pixel value corresponding to the red channel of the pixel point to the pixel value of the green channel; w is a₂Expressing the ratio of the pixel value corresponding to the green channel of the pixel point to the pixel value of the blue channel; i (I, j,1) represents the pixel value of the pixel point in the red channel; i (I, j,2) represents the pixel value of the pixel point in the green channel; i (I, j,3) represents the pixel value of the pixel point in the blue channel;

if the grid line edge satisfies the upper setting range, I (I, j,1) is set to I (I, j,2) is set to I (I, j,3) is set to 0.

Step a2 includes setting the following network lines along the lower pixel rows, all columns as the reference:

sum＝I(i,j,1)+I(i,j,2)+I(i,j,3)

if the grid line is satisfied with the downward setting range, let I (I, j,1) be I (I, j,2) be I (I, j,3) be 0; i (I, j,1) represents the pixel value of the pixel point in the red channel; i (I, j,2) represents the pixel value of the pixel point in the green channel; i (I, j,3) represents the pixel value of the pixel point in the blue channel; sum represents the sum of pixel values of the pixel points in red, green and blue channels.

In step S4, the method for generating the formation fracture information extraction area image3 includes: traversing pixels in the stratum fracture interlayer distribution image2, wherein the form of a pixel two-dimensional array is [. cndot. ] i ] j, summing pixel values [ r, g, b ] of three color channels of pixel points, setting a threshold value T, and calculating:

sum＝I(i,j,1)+I(i,j,2)+I(i,j,3)

wherein, I (I, j,1) represents the pixel value of the pixel point in the red channel; i (I, j,2) represents the pixel value of the pixel point in the green channel; i (I, j,3) represents the pixel value of the pixel point in the blue channel; sum represents the sum of pixel values of the pixel points in red, green and blue channels; when sum > T, I (I, j,1) ═ I (I, j,2) ═ I (I, j,3) ═ 255; and when sum is less than or equal to T, I (I, j,1) ═ I (I, j,2) ═ I (I, j,3) ═ 0, and through traversal and calculation, an image formed by the region where the pixel point with the pixel value of 0 is the formation fracture information extraction region image 3.

In the step S4, the method for generating the distribution graph image4 corresponding to the porosity of the formation fracture layer includes: setting the size of a traversal window as s x k, the step length as s, wherein s represents the length of the s x k window, k represents the width of the s x k window, and k takes a fixed value and represents the original width of a rock stratum fracture information extraction region image3 image for a rock stratum fracture information extraction region image 3; in an s x k traversal window, the porosity size Pr of a stratum fracture layer is as follows:

Pr＝n/N

wherein n is the number of pixels of the stratum fracture of the S x k traversal window; and N is the total number of the traversal pixel points of the traversal window.

The step S5 includes: the window scans a stratum fracture information extraction area image3, the size of a traversal window is set to be m x h, and the step length is taken to be m; traversing the rock stratum fracture information extraction area image3 along the horizontal direction, and then traversing the rock stratum fracture information extraction area image3 along the vertical direction with the step length being m; in an m x h traversal window, the porosity size Pr' of a stratum fracture layer is as follows:

Pr'＝n'/N'

wherein n' is the number of pixel points of the stratum fracture of the traversal window m x h; n' represents the total number of the traversal pixel points of the m x h traversal window; and calculating the rock stratum fracture rate Pr' in each m x h traversal window, visualizing the porosity rate of each window through progressive colors, and gradually displaying the porosity rate represented by the traversal windows at different positions to obtain a distributed cloud image5 of the rock stratum porosity rate under the rock stratum fracture information extraction area image3 standard.

According to the rock stratum porosity distribution calculation method provided by the invention, the computer is utilized to process the rock stratum images, so that omission or negligence which possibly occurs in the manual calculation of the porosity is reduced, related errors caused by related experiments in the traditional porosity calculation method are avoided, the calculation speed of the rock stratum porosity distribution is greatly accelerated to a certain extent, and the distribution of the rock stratum porosity in a certain trend is visually presented from the data; meanwhile, the interference of other images is reduced, and the calculation is more accurate and efficient.

Drawings

FIG. 1 is a schematic flow diagram of the process of the present invention.

Fig. 2 is a schematic diagram of a stratum movement original image according to an embodiment of the invention.

FIG. 3 is a schematic illustration of a formation fracture zonal profile of an embodiment of the present invention.

FIG. 4 is a schematic diagram of a formation fracture information extraction zone in accordance with an embodiment of the present invention.

FIG. 5 is a graph illustrating a formation fracture distribution according to an embodiment of the present invention.

FIG. 6 is a cloud chart illustrating the porosity of a formation according to an embodiment of the present invention.

Detailed Description

FIG. 1 is a schematic flow chart of the method of the present invention: the rock stratum porosity distribution calculation method provided by the invention comprises the following steps:

s1, simulating a stratum fracture enrichment area, and outputting a corresponding stratum movement original three-color image;

step S1 is specifically that a stratum of a stratum fracture enrichment area is simulated by using a discrete unit-based method, and a corresponding stratum motion original tricolor image1 is output according to stratum condition parameters; the original three-color image1 of the stratum movement highlights qualitative performance of rock stratum fracture expansion, such as fracture development and distribution positions, fracture density parameters and the like; the original color map stores pixel values in a color mode of three color [ r, g, b ] channels.

When the formation movement is simulated by using UDEC numerical software based on a discrete unit method, parameterization processing is carried out on rock matrixes (blocks) and rock contacts (joints and rock layers in the rock layers) to generate rock layer condition parameters, wherein the rock layer condition parameters comprise: mass density, elastic property, plastic property and the like of the block body so as to determine the mechanical property of the material; rigidity, shear rigidity, tensile strength, friction angle, cohesion and the like to determine the mechanical characteristics of the contact between blocks; formation bulk size, etc.

S2, preprocessing the original three-color image of the stratum movement;

the pretreatment process specifically comprises the following steps: the stratum movement original three-color image1 is used for generating a corresponding three-dimensional image digital matrix Array1, the digital matrix Array1 reflects three color [ r, g, b ] channels in a three-dimensional form, pixel positions are displayed in the corresponding channels in a two-dimensional Array form [ · ] i ] j respectively, and three original color pigment values are reflected through pixel values.

Wherein each pixel location contains three color channels, the red, green and blue tristimulus pixel values are distributed in the pixel location in the form of r, g, b, with r, g, b e [0,255 ].

S3, detecting an original three-color image1 of stratum movement by analyzing a three-dimensional image digital matrix Array1, extracting rock stratum fracture boundary position information, and removing redundant rock stratum pixel information to obtain a rock stratum fracture interlayer distribution map image 2;

the concrete method for extracting the rock stratum fracture boundary position information comprises the following steps of removing redundant upper and lower grid lines in the original tricolor image1 of the stratum motion output in the step S1, wherein the redundant upper and lower grid lines represent redundant interlayer information, so that the fracture interlayer information is more prominent, the computer calculation burden is reduced, and the calculation speed is accelerated.

A1. And taking the lower grid lines along the upper pixel rows and all columns as a reference, and setting:

w₁＝I(i,j,1)/I(i,j,2)

w₂＝I(i,j,2)/I(i,j,3)

wherein, w₁Representing pixel points I [ I, j]The ratio of the pixel value corresponding to the red channel to the pixel value of the green channel; w is a₂Representing pixel points I [ I, j]The ratio of the pixel value corresponding to the green channel to the pixel value corresponding to the blue channel. I (I, j,1) represents the pixel value of the pixel point in the red channel; i (I, j,2) represents the pixel value of the pixel point in the green channel; and I (I, j,3) represents the pixel value of the pixel point in the blue channel.

If the grid line meets the upper set range, making I (I, j,1) I (I, j,2) I (I, j,3) 0; in this embodiment, the grid line edge is set to have a range of 0.9 < w₁,w₂Less than 1.1, i.e. 0.9 and 1.1 are the upper and lower bounds of the corresponding ratio of the pixel channel value in this embodiment, respectively;

A2. and (3) setting the following data lines along the lower pixel row and all columns as reference:

sum＝I(i,j,1)+I(i,j,2)+I(i,j,3)

if the grid line is satisfied with the downward setting range, let I (I, j,1) be I (I, j,2) be I (I, j,3) be 0; the grid line is set to a downward range of 130 < sum < 170 in this embodiment; sum represents the sum of pixel values of the pixel points in red, green and blue channels.

The extraction of the interlayer fracture characteristic information is realized through the steps A1 and A2, and a stratum fracture interlayer distribution map image2 is obtained.

S4, traversing a stratum fracture interlayer distribution map image2 along the stratum direction, setting a threshold T, selecting a T value according to an actual value of a fracture pixel point of a stratum fracture interlayer distribution map image2, and selecting 8 in the embodiment; calculating, identifying and extracting the pixel points of which the sum of the RGB values Isum is smaller than a set threshold value T to obtain a rock stratum fracture representative pixel point gray [ i, j]And a formation fracture information extraction area image 3; representing the rock stratum fracture by pixel gray [ i, j]Summing to obtain the number n of the traversal window rock stratum fracture pixel points₁Calculating the porosity of a rock stratum fracture layer; drawing a distribution curve graph image4 corresponding to the porosity of the fracture layer of the rock stratum;

in step S4, the specific method of generating the formation fracture information extraction area image3 is: traversing pixels in the stratum fracture interlayer distribution image2, wherein the two-dimensional array of the pixels is in the form of [ ·]_i]_jThree color channel pixel values [ r, g, b ] for a pixel point]Summing, setting a threshold value T, and calculating:

sum＝I(i,j,1)+I(i,j,2)+I(i,j,3)

when sum > T, I (I, j,1) ═ I (I, j,2) ═ I (I, j,3) ═ 255; and when sum is less than or equal to T, I (I, j,1) ═ I (I, j,2) ═ I (I, j,3) ═ 0, and through traversal and calculation, an image formed by the region where the pixel point with the pixel value of 0 is the formation fracture information extraction region image 3.

The specific method for drawing the distribution curve image4 corresponding to the porosity of the fracture layer of the rock stratum comprises the following steps: setting the size of a traversal window as s x k, the step length as s, wherein s represents the length of the s x k window, k represents the width of the s x k window, and k takes a fixed value and represents the original width of a rock stratum fracture information extraction region image3 image for a rock stratum fracture information extraction region image 3; in an s x k traversal window, the porosity size Pr of a stratum fracture layer is as follows:

Pr＝n/N

wherein n is the number of pixels of the stratum fracture of the S x k traversal window; and N is the total number of the traversal pixel points of the traversal window. And visualized by a computer program, thereby obtaining a layer crack distribution graph image4 of a unit length in the horizontal direction.

S5, setting a full-image traversal window, and performing full-image range traversal on the stratum fracture interlayer distribution image2, wherein the full-image range traversal window moves the window pixel by pixel along the vertical and horizontal directions according to a set step length m; after the full-image-range traversal of the stratum fracture interlayer distribution image2 is completed, the stratum fracture interlayer distribution image can be regarded as a stratum porosity distribution cloud image 5.

Step S5 specifically includes window scanning the formation fracture information extraction area image3, and setting the traversal window size m × h; the step length is m, where m is 1 in this embodiment; traversing the rock stratum fracture information extraction area image3 along the horizontal direction, and then traversing the rock stratum fracture information extraction area image3 along the vertical direction with the step length being m; in an m x h traversal window, the porosity size Pr' of a stratum fracture layer is as follows:

Pr'＝n'/N'

wherein n' is the number of pixel points of the stratum fracture of the traversal window m x h; n' represents the total number of the traversal pixel points of the m x h traversal window; and calculating the rock stratum fracture rate Pr' in each m x h traversal window, visualizing the porosity rate of each window through progressive colors, and progressively displaying the porosity rate represented by the traversal windows at different positions, thereby obtaining a distributed cloud image5 of the rock stratum porosity rate under the rock stratum fracture information extraction area image3 standard.

In this embodiment, a coal formation is used as an example of the application of the method of the present invention. And mining a certain coal seam in certain coal mines in Shanxi by adopting an underground longwall coal mining method, wherein the coal mine is concentrated in a No. 4 coal seam closest to the earth surface. The thickness of the No. 4 coal seam is 6.5m, the cutting height of a coal mining machine is 3.0m, and the thickness of top coal caving is 3.5 m. The width of longwall mining is 200 meters. The geological conditions of the No. 4 coal bed and the overlying strata are relatively simple, no fault or hydrological influence exists, and the structures of the coal bed and the overlying strata are basically horizontal structures.

The method comprises the following application steps:

step 1, as shown in fig. 2, is a schematic diagram of a stratum movement original image according to an embodiment of the present invention. And (3) simulating the stratum movement of the stratum fracture enrichment area by using a UDEC numerical software based on a discrete unit method, wherein the sizes of blocks of a2 nd key layer (a main key layer) and a1 st key layer are respectively designed to be 42.6m and 23.5m in the UDEC model. And simultaneously, the actually measured lithology and geometric properties of the coal-series stratum above the No. 4 coal seam and the mechanical properties of the complete rock matrix in the stratum are obtained. Outputting a corresponding stratum movement original image as shown in figure 2, wherein the image shows qualitative performance of stratum fracture development, such as fracture distribution positions and fracture density parameters; and the pigment values of the RGB three channels.

Step 2, inputting the stratum movement original image1 in the step 1 into a computer program, reading and storing the original image, generating a corresponding three-dimensional image digital matrix Array1, wherein the digital matrix Array1 reflects three color channels in three dimensions: r, g and b, respectively displaying pixel positions in a two-dimensional array form [. cndot. ] i ] j form in corresponding channels, and reflecting three original color pixel values of red, green and blue through pixel values;

step 3, as shown in fig. 3, is a schematic diagram of a formation fracture interval distribution diagram according to an embodiment of the present invention. The concrete method for extracting the stratum movement original image and the position information of the rock stratum fracture boundary thereof by analyzing the stratum movement original image matrix Array in the step 2 comprises the following steps: and (3) removing redundant upper and lower grid lines (representing redundant interlayer information) in the stratum motion original image output in the step (1), aiming at making fracture interlayer information more prominent, reducing the calculation load of a computer and accelerating the calculation speed. And taking the lower grid lines along the upper pixel rows and all columns as a reference, and setting:

w₁＝I(i,j,1)/I(i,j,2)

w₂＝I(i,j,2)/I(i,j,3)

wherein, w₁Expressing the ratio of the pixel value corresponding to the red channel of the pixel point to the pixel value of the green channel; w is a₂Expressing the ratio of the pixel value corresponding to the green channel of the pixel point to the pixel value of the blue channel; i (I, j,1) represents the pixel value of the pixel point in the red channel; i (I, j,2) represents the pixel value of the pixel point in the green channel; i (I, j,3) represents the pixel value of the pixel point in the blue channel;

satisfy 0.9 < w₁,w₂If < 1.1, let I (I, j,1) be I (I, j,2) be I (I, j,3) be 0;

and (3) setting the following data lines along the lower pixel row and all columns as reference:

sum＝I(i,j,1)+I(i,j,2)+I(i,j,3)

if 130 < sum < 170 is satisfied, I (I, j,1) ═ I (I, j,2) ═ I (I, j,3) ═ 0 is assumed.

The interlayer fracture characteristic information is extracted through the above formula, and a stratum fracture interlayer distribution diagram is obtained.

Step 4, ergodically scanning image pixels distributed among stratum fracture layers in the step 3, wherein the form of a pixel two-dimensional array is [ [. cndot. ] i ] j, summing pixel values [ r, g and b ] of three color channels of the pixels, setting a threshold value T, selecting the T value according to actual values of fracture pixel points of the stratum fracture distribution diagram, and taking the image as 8:

sum＝I(i,j,1)+I(i,j,2)+I(i,j,3)

wherein, when sum > T, I (I, j,1) ═ I (I, j,2) ═ I (I, j,3) ═ 255 is set; when sum is less than or equal to T, I (I, j,1) ═ I (I, j,2) ═ I (I, j,3) ═ 0 is set. Fig. 4 is a schematic diagram of a formation fracture information extraction area according to an embodiment of the present invention. Through traversal and calculation, an image formed by the region where the pixel point with the pixel value of 0 is located is a stratum fracture information extraction region map.

In step 4, fig. 5 is a graph illustrating a formation fracture distribution according to an embodiment of the present invention. The method for drawing the rock stratum fracture distribution curve comprises the steps of setting the size of a traversal window to be 2 x 10 and taking the step length to be 2 for the rock stratum fracture information extraction regional graph. In a traversal window (2 x 10), the porosity size Pr of a stratum fracture layer is Pr-N/N;

wherein N is the number of the traversing window (2 x 10) rock stratum fracture pixel points, and N is the total number of the traversing window (2 x 10) traversing pixel points. The data results are visualized by a computer program, from which a plot of the formation fracture distribution per unit length in the horizontal direction as shown in fig. 5 can be obtained.

And 5, scanning a stratum fracture information extraction area graph shown in the graph 4 through a window, extracting stratum fracture area characteristics according to the step 4, setting the size of a traversal window to be 5 x 10, taking the step length to be 1, and traversing the stratum fracture information extraction area shown in the graph 4 along the horizontal direction and the vertical direction by taking the step length to be 1. In a traversal window (5 x 10), the porosity dimension Pr' of the formation fracture layer is:

Pr'＝n'/N'

wherein N 'is the number of the pixel points of the formation fracture in the traversal window (5 × 10), and N' is the total number of the pixel points of the traversal window (5 × 10). Fig. 6 is a cloud chart illustrating the porosity of a rock formation according to an embodiment of the present invention. And calculating the porosity of the rock stratum fracture layer in each traversal window, visualizing the porosity of each window through progressive colors by a computer program, and gradually displaying the porosity represented by the traversal windows at different positions, thereby obtaining a distribution cloud picture of the rock stratum porosity as shown in fig. 6 under the standard of a rock stratum fracture information extraction area picture.

And analyzing and calculating the porosity distribution of the rock stratum by applying the method to the No. 4 coal seam of the coal mine rock stratum. The calculation results are shown in FIGS. 2 to 6.

As can be seen from fig. 2 to 6, a rock stratum porosity distribution calculation method based on digital image recognition is applied to the coal seam No. 4 of the coal mine rock stratum, and UEDC stratum simulation, digital matrix extraction, rock stratum fracture boundary information extraction, rock stratum fracture interlayer distribution diagram acquisition, rock stratum fracture information extraction area diagram acquisition, porosity distribution curve diagram and cloud diagram drawing are successfully performed on the coal rock stratum. The method can efficiently and quickly calculate the porosity of the rock layer, and the whole operation implementation process accords with the actual condition of the working rock body.