WO2021238340A1 - Improved markov chain monte carlo two-dimensional rock section reconstruction method and system - Google Patents

Abstract

An improved Markov chain Monte Carlo two-dimensional rock section reconstruction method and system. The porosity of a rock section electron microscope scanning image is calculated; neighborhood systems to be used and conditional probabilities corresponding to different combinations of the neighborhood systems are determined; iterative evaluation from right to left is carried out on even-numbered rows, iterative evaluation from left to right is carried out on odd-numbered rows, the first pixel point and the last pixel point in each row are obtained via pixel points in the previous row and in the same column and by means of two vertical neighborhoods, and a reconstruction path is obtained; according to an anisotropic characteristic exhibited by a rock section image, considering that the conditional probabilities are also anisotropic, the conditional probability of a binarized image in each direction is scanned along the reconstruction path; and a reconstructed image is obtained by means of Monte Carlo calculation and assignment. In the method, a conditional probability scanning direction is matched with a reconstruction path, and a random number is selected to reconstruct an image, such that the influence and effect of all directions are equalized, thereby improving the degree of fitting of a variation function of the reconstructed image and that of an original image while anisotropy is exhibited.

Description

Improved Markov chain Monte Carlo two-dimensional rock slice reconstruction method and system

This application claims the priority of a Chinese patent application filed with the Chinese Patent Office on May 28, 2020, the application number is 202010468814.3, and the invention title is "Improved Markov Chain Monte Carlo 2D Rock Slice Reconstruction Method and System" , Its entire content is incorporated into this application by reference.

Technical field

The invention belongs to the technical field of computer image processing, and relates to an improved Markov chain Monte Carlo two-dimensional rock slice reconstruction method and system.

Background technique

Image reconstruction is widely used in various fields such as medicine and geology due to its high convenience and low cost, and its purpose is to perform feature prediction and structural analysis of unknown areas. In order to make the characteristics of the reconstructed image more suitable for reality, choosing a suitable imaging method is a direct factor that affects the reconstruction effect.

In the prior art, Kejian Wu’s article "An Efficient Markov Chain Model for the Simulation of Heterogeneous Soil Structure" systematically explained the MCMC method, using the idea of Markov chain to select an appropriate neighborhood system from a fixed The direction traverse scans the rock electron microscope slice image, the distribution law is statistically analyzed to obtain the conditional probability, and then the Monte Carlo algorithm is used to calculate the random number to solve the calculation problem, and the quality of the reconstruction result is judged by the fitting of the variogram function. However, the path will continue to change directions during the reconstruction process, and in order to make the result more closely fit the original image, weights are set for the conditional probability, which directly changes the statistical information, and the fitting effect of the variogram cannot be quantitatively judged. Therefore, when the anisotropy of the rock image is widespread, the improved method mentioned in the present invention is more reasonable.

Summary of the invention

Aiming at the defects and deficiencies in the prior art, the present invention provides an improved Markov chain Monte Carlo two-dimensional rock slice reconstruction method and system, aiming at the widespread feature of rock image anisotropy , And the current existing technology does not consider the conditional probability caused by this feature and there is also an anisotropy problem.

In order to achieve the above objectives, the present invention adopts the following technical solutions:

An improved Markov chain Monte Carlo 2D rock slice reconstruction method, including the following steps:

Step 1: Binarize the original image to calculate the porosity:

For the original image with a size of m*n, the threshold of the image is obtained by the maximum variance between classes, and the original image is binarized according to the threshold, and the total number of pixels with different pixel values is counted. The porosity is calculated as follows:

In the binarized image, a skeleton with a pixel value of 0 represents a skeleton, and a pore with a pixel value of 1 is represented; in the above formula, num ₁ is the total number of pixels with a pixel value of 1, and num ₀ is a pixel with a pixel value of 0 The total number of;

Step 2: Determine the neighborhood system and conditional probability calculation formula:

Determine the neighborhood system to be used on the basis of the binarized image obtained in step 1, list all possible combinations of the neighborhood system, and find the probability that the current pixel is 0 or 1 under each combination condition ; The conditional probability calculation formula is:

Probability that the pixel is 0 under the neighborhood condition:

The probability that the pixel is 1 under the neighborhood condition:

In the above formula, N is a certain combination in the domain system, x _r is the pixel in the combination, x _(i,j) is the pixel currently sought, and num _r→0 is the combination of x _{(i ,j)} is the number of times that is 0, num _r→1 is _{the number of times that x (i,j)} is 1;

Step 3: Determine the reconstruction path and the conditional probability scan direction:

Set an empty image with the same size as the original image. In the empty image, iteratively evaluate from right to left on even-numbered rows, and iteratively evaluate from left to right on odd-numbered rows. The first and last pixels of each row are determined by it The pixel points in the same column in the previous row are obtained through the vertical two neighborhoods to obtain the reconstruction path; according to the reconstruction path, scan the conditional probability of each direction of the binarized image obtained in step 1;

Step 4: Use Monte Carlo calculation and assignment:

Assign the first pixel in the binarized image to the first pixel of the empty image, and then use the horizontal two-neighbors to calculate and assign the first row from left to right, and then the last pixel in the second row through the vertical Two-neighborhood, the last pixel in the first row is used as a condition to calculate the pixel value, and then starting from the n-1th pixel in the second row, the left-to-four-neighborhood is used for matching calculation, and the third row is iterated from left to right , Repeat the previous operation for the following line breaks, until the last pixel is determined;

Step 5: The porosity satisfies the condition to stop the loop, and the reconstructed image is obtained;

_Regarding the porosity as the condition for the end of the cycle, when the porosity P _{KXD of the} reconstructed image and the original image porosity p kxd satisfy the following formula, the cycle is stopped, and the reconstructed image is obtained;

p _kxd -0.05≤P _KXD ≤p _{kxd +0.05} (4)

Otherwise, reset the empty image and repeat steps 4 and 5 until formula (4) is satisfied, and finally a reconstructed image is obtained.

The invention also includes the following technical features:

Specifically, the method for evaluating the reconstructed image obtained in step 5 is as follows: respectively calculating the variation value of the binarized image and the reconstructed image;

The calculation formula of the variation value is as follows:

Among them, a is the step length, that is, the distance between two points, which is selected from 1 to 50 here, N(a) is the number of two points with a distance of a, x is the pixel value, and i is the pixel point coordinate;

Calculate the average distance between the variogram of the binarized image and the reconstructed image. The calculation formula is as follows:

In the formula, R(a) is the variation value of the reconstructed image, r(a) is the variation value of the binarized image; the data represents the fit of the variation curve of the binarized image and the reconstructed image, Take the difference between the reconstructed image and the binarized image at the same step length, and accumulate the difference under the non-synchronization length to get the average value, the smaller the value, the better the fit.

The present invention also provides a rock slice reconstruction system based on Markov chain Monte Carlo, including:

The porosity calculation module is used to obtain the threshold of the original image based on the maximum variance between clusters, and binarize the original image according to the threshold. The binarized image represents the skeleton with a pixel value of 0 and the pore with a pixel value of 1 , Count the total number of pixels with pixel values of 0 and 1, and calculate the porosity;

The domain system and conditional probability determination module is used to determine the neighborhood system based on the binary image, list all possible combinations of the neighborhood system, and find the probability that the pixel value is 0 or 1 under each combination condition;

Reconstruction path and conditional probability determination module in each direction. It is used to evaluate iteratively from right to left in even rows in an empty image, and iteratively evaluate from left to right in odd rows, the first and last of each row The pixel point is obtained from the pixel point in the same column of its previous row through the vertical two neighborhoods to determine the reconstruction path; according to the reconstruction path, scan the conditional probability of each direction of the binarized image obtained in step 1;

The calculation and assignment module is used to iteratively calculate the value of each pixel in the empty image using Monte Carlo;

Obtain a reconstructed image module, which is used to use porosity as a condition for the end of the cycle to obtain a reconstructed image;

The result fit evaluation module is used to characterize the fit of the variation curve of the binarized image and the reconstructed image through data, and take the difference between the reconstructed image and the binarized image at the same step size, Accumulate the difference values under non-synchronization length to get the average value, the smaller the value, the better the fit.

Compared with the prior art, the present invention has the following beneficial technical effects:

The MCMC reconstruction method proposed by the present invention, based on the anisotropy of the rock image, considers that the conditional probability also has the characteristics of anisotropy. The conditional probability of this direction is scanned from different directions. In the iterative evaluation process, the conditional probability is scanned The direction is matched with the reconstruction path, and random numbers are selected to reconstruct the image, and the influence of each direction is equalized. While characterizing the anisotropy, it improves the fit of the variogram between the reconstructed image and the original image.

Attached drawings

Figure 1 is a schematic diagram B of original image A and binarization scanned by a rock electron microscope;

Figure 2 is a schematic diagram of the horizontal and vertical two neighborhoods; (gray is the neighborhood system pixels)

Figure 3 is a schematic diagram of the left and right four neighborhoods; (gray is the pixels of the neighborhood system)

Figure 4 is a schematic diagram of a serpentine reconstruction path;

Figure 5 is a schematic diagram of the prior art method and the imaging result of the present invention;

Figure 6 is a schematic diagram of a variogram fitting curve;

Figure 7 is a flow chart of the present invention.

The specific embodiments of the present invention will be described in detail below. It should be understood that the specific embodiments described here are only used to illustrate and explain the present invention, and are not used to limit the present invention.

Detailed ways

The present invention aims at the characteristics of anisotropy in the scanning images of rock slices by electron microscopy, uses Markov chain ideas and Monte Carlo sampling methods, and takes into account the anisotropy shown in conditional probabilities in different directions due to the characteristic information of the rock image. And use iterative calculation method for imaging display.

Example 1:

As shown in Figures 1 to 7, this embodiment provides an improved Markov chain Monte Carlo two-dimensional rock slice reconstruction method. Figure 7 is a flowchart of the present invention. This method performs conditional probability scanning. Considering the anisotropy, match the conditional probability scanning direction with the reconstruction path, and select random numbers to reconstruct the image, equalize the influence of each direction, and characterize the anisotropy while improving the change between the reconstructed image and the original image. The degree of fit of the difference function; specifically includes the following steps:

Step 1: Binarize the original image to calculate the porosity:

The left A of Figure 1 is the original image. For the original image of size m*n, the threshold of the image is obtained by the maximum inter-class variance, and the original image is binarized according to the threshold. The right B of Figure 1 is the binarized rock image. Statistics The total number of pixels with different pixel values, and the porosity is calculated as follows:

In the binarized image, a skeleton with a pixel value of 0 represents a skeleton, and a pore with a pixel value of 1 is represented; in the above formula, num ₁ is the total number of pixels with a pixel value of 1, and num ₀ is a pixel with a pixel value of 0 Specifically, in this embodiment, the porosity of the binarized rock image calculated from the right B of Figure 1 is 17.68%;

Step 2: Determine the neighborhood system and conditional probability calculation formula:

Determine the neighborhood system to be used on the basis of the binarized image obtained in step 1, list all possible combinations of the neighborhood system, and find the probability that the current pixel is 0 or 1 under each combination condition ; The conditional probability calculation formula is:

The probability that the pixel is zero under the neighborhood condition:

The probability that the pixel is 1 under the neighborhood condition:

In the above formula, N is a certain combination in the domain system, x _r is the pixel in the combination, x _(i,j) is the pixel currently sought, and num _r→0 is the combination of x _{(i ,j)} is the number of times that is 0, num _r→1 is _{the number of times that x (i,j)} is 1;

In this embodiment, Figure 2 shows the horizontal and vertical two neighborhoods respectively, and its conditional probability calculation formula is:

Among them, N ₂ is a certain combination of two neighborhoods, x _r is the pixel in the combination, x _(i,j) is the pixel currently sought, and num _r→0 is x _{(i, j)} is the number of times 0, num _r→1 is _{the number of times x (i, j)} is 1; the expression on the left calculates the probability that the pixel is 0 under the condition of the neighborhood, and the expression on the right calculates the probability of the pixel under the condition of the neighborhood The probability that the pixel is 1;

Figure 3 shows the left and right four neighborhoods, and its conditional probability calculation formula is:

Among them, N ₄ is a certain combination of two neighborhoods, x _r is the pixel in the combination, x _(i,j) is the pixel currently sought, and num _r→0 is x _{(i, j)} is the number of times 0, num _r→1 is _{the number of times x (i, j)} is 1; the expression on the left calculates the probability that the pixel is 0 under the condition of the neighborhood, and the expression on the right calculates the probability of the pixel under the condition of the neighborhood The probability that the pixel is 1; the related explanation is similar to the second neighborhood;

Step 3: Determine the reconstruction path and the conditional probability scan direction:

Figure 4 is the principle diagram of the serpentine path. An empty image with the same size as the original image is set. In the empty image, the even rows are evaluated iteratively from right to left, and the odd rows are evaluated iteratively from left to right. The first and last pixel of the row are obtained from the pixels in the same column of the previous row through the vertical two neighborhoods to obtain the reconstruction path; according to the reconstruction path, scan the binarized image obtained in step 1 for each direction Conditional Probability;

Step 4: Use Monte Carlo calculation and assignment:

Assign the first pixel in the binarized image to the first pixel of the empty image, and then use the horizontal two-neighbors to calculate and assign the first row from left to right, and then the last pixel in the second row through the vertical Two-neighborhood, the last pixel in the first row is used as a condition to calculate the pixel value, and then starting from the n-1th pixel in the second row, the left-to-four-neighborhood is used for matching calculation, and the third row is iterated from left to right , Repeat the previous operation for the following line breaks, until the last pixel is determined;

Step 5: The porosity satisfies the condition to stop the loop, and the reconstructed image is obtained;

_Regarding the porosity as the condition for the end of the cycle, when the porosity P _{KXD of the} reconstructed image and the original image porosity p kxd meet the following formula, stop the cycle, otherwise, reset the empty image and repeat steps 4 and 5;

p _kxd -0.05≤P _KXD ≤p _{kxd +0.05} (4)

Specifically, FIG. 5A is a reconstructed image using the conditional probability obtained by scanning the binarized image once in the prior art; FIG. 5B is a reconstructed image using the same conditional probability as the reconstruction direction in this embodiment; both The porosity of the image is 14.21% and 19.36%, both of which have reached the condition for the end of the cycle.

Step 6: Calculate the variation values of the binarized image and the reconstructed image respectively;

The calculation formula of the variation value is as follows:

Among them, a is the step length, that is, the distance between two points, which is selected from 1 to 50 here, N(a) is the number of two points with a distance of a, x is the pixel value, and i is the pixel point coordinate;

Calculate the average distance between the variogram of the binarized image and the reconstructed image. The calculation formula is as follows:

In the formula, R(a) is the variation value of the reconstructed image, r(a) is the variation value of the binarized image; the data represents the fit of the variation curve of the binarized image and the reconstructed image, Take the difference between the reconstructed image and the binarized image at the same step length, and accumulate the difference under the non-synchronization length to get the average value, the smaller the value, the better the fit.

Figure 6 is the variogram of the three images, the punctuation is the original image variability, the red line is the prior art Figure 5A, the blue line is the embodiment of Figure 5B, you can see that the blue line is more fitted to the punctuation, they are different from the punctuation The average distance between the two is 0.0208 (FIG. 5A) and 0.0073 (FIG. 5B), respectively. In the reconstruction process, the conditional probability in the same direction as the reconstruction path is used, that is, the imaging effect of the method of the present invention is better.

Example 2:

This embodiment provides a rock slice reconstruction system based on Markov chain Monte Carlo, including:

The porosity calculation module is used to obtain the threshold of the original image based on the maximum variance between clusters, and binarize the original image according to the threshold. The binarized image represents the skeleton with a pixel value of 0 and the pore with a pixel value of 1 , Count the total number of pixels with pixel values of 0 and 1, and calculate the porosity;

The domain system and conditional probability determination module is used to determine the neighborhood system based on the binarized image, list all possible combinations of the neighborhood system, and calculate the probability that the pixel value is 0 or 1 under each combination condition;

Reconstruction path and conditional probability determination module for each direction. It is used to evaluate iteratively from right to left in even rows in an empty image, and iteratively evaluate from left to right in odd rows, the first and last of each row The pixel point is obtained from the pixel point in the same column of its previous row through the vertical two neighborhoods to determine the reconstruction path; according to the reconstruction path, scan the conditional probability of each direction of the binarized image obtained in step 1;

The calculation and assignment module is used to use Monte Carlo to iteratively calculate and determine the pixel value of each point of the empty image;

Obtain a reconstructed image module, which is used to use porosity as a condition for the end of the cycle to obtain a reconstructed image;

The result fit evaluation module is used to characterize the fit of the variation curve of the binarized image and the reconstructed image through data, and take the difference between the reconstructed image and the binarized image at the same step size, Accumulate the difference values under non-synchronization length to get the average value, the smaller the value, the better the fit.

Claims (3)

1. An improved Markov chain Monte Carlo two-dimensional rock slice reconstruction method, which is characterized in that it specifically includes the following steps:

   Step 1: Binarize the original image to calculate the porosity:

   For the original image with a size of m*n, the threshold of the image is obtained by the maximum variance between classes, and the original image is binarized according to the threshold, and the total number of pixels with different pixel values is counted. The porosity is calculated as follows:

   

   In the binarized image, a skeleton with a pixel value of 0 represents a skeleton, and a pore with a pixel value of 1 is represented; in the above formula, num ₁ is the total number of pixels with a pixel value of 1, and num ₀ is a pixel with a pixel value of 0 The total number of;

   Step 2: Determine the neighborhood system and conditional probability calculation formula:

   Determine the neighborhood system to be used on the basis of the binarized image obtained in step 1, list all possible combinations of the neighborhood system, and find the probability that the current pixel is 0 or 1 under each combination condition ; The conditional probability calculation formula is:

   The probability that the pixel is zero under the neighborhood condition:

   

   The probability that the pixel is 1 under the neighborhood condition:

   

   In the above formula, N is a certain combination in the domain system, x _r is the pixel in the combination, x _(i,j) is the pixel currently sought, and num _r→0 is the combination of x _{(i ,j)} is the number of times that is 0, num _r→1 is _{the number of times that x (i,j)} is 1;

   Step 3: Determine the reconstruction path and the conditional probability scan direction:

   Set an empty image with the same size as the original image. In the empty image, iteratively evaluate from right to left on even-numbered rows, and iteratively evaluate from left to right on odd-numbered rows. The first and last pixels of each row are determined by it The pixel points in the same column in the previous row are obtained through the vertical two neighborhoods to obtain the reconstruction path; according to the reconstruction path, scan the conditional probability of each direction of the binarized image obtained in step 1;

   Step 4: Use Monte Carlo calculation and assignment:

   Assign the first pixel in the binarized image to the first pixel of the empty image, and then use the horizontal two-neighbors to calculate and assign the first row from left to right, and then the last pixel in the second row through the vertical Two-neighborhood, the last pixel in the first row is used as a condition to calculate the pixel value, and then starting from the n-1th pixel in the second row, the left-to-four-neighborhood is used for matching calculation, and the third row is iterated from left to right , Repeat the previous operation for the following line breaks, until the last pixel is determined;

   Step 5: The porosity satisfies the condition to stop the loop, and the reconstructed image is obtained;

   _Regarding the porosity as the condition for the end of the cycle, when the porosity P _{KXD of the} reconstructed image and the original image porosity p kxd satisfy the following formula, the cycle is stopped, and the reconstructed image is obtained;

   p _kxd -0.05≤P _KXD ≤p _{kxd +0.05} (4)

   Otherwise, reset the empty image and repeat steps 4 and 5 until formula (4) is satisfied, and finally a reconstructed image is obtained.
2. The improved Markov chain Monte Carlo two-dimensional rock slice reconstruction method according to claim 1, characterized in that the reconstructed image obtained in step 5 is evaluated, and the method is as follows: respectively calculating the binarized image And the variation value of the reconstructed image;

   The calculation formula of the variation value is as follows:

   

   Among them, a is the step length, that is, the distance between two points, which is selected from 1 to 50 here, N(a) is the number of two points with a distance of a, x is the pixel value, and i is the pixel coordinate;

   Calculate the average distance between the variogram of the binarized image and the reconstructed image. The calculation formula is as follows:

   

   In the formula, R(a) is the variation value of the reconstructed image, r(a) is the variation value of the binarized image; the data represents the fit of the variation curve of the binarized image and the reconstructed image, Take the difference between the reconstructed image and the binarized image at the same step length, and accumulate the difference under the non-synchronized length to get the average value. The smaller the value, the better the fit.
3. An improved Markov chain Monte Carlo rock slice reconstruction system, which is characterized in that it includes:

   The porosity calculation module is used to obtain the threshold of the original image based on the maximum variance between clusters, and binarize the original image according to the threshold. The binarized image represents the skeleton with a pixel value of 0 and the pore with a pixel value of 1 , Count the total number of pixels with pixel values of 0 and 1, and calculate the porosity;

   The domain system and conditional probability determination module is used to determine the neighborhood system based on the binarized image, list all possible combinations of the neighborhood system, and calculate the probability that the pixel value is 0 or 1 under each combination condition;

   Reconstruction path and conditional probability determination module for each direction. It is used to evaluate iteratively from right to left in even rows in an empty image, and iteratively evaluate from left to right in odd rows, the first and last of each row The pixel point is obtained from the pixel point in the same column of its previous row through the vertical two neighborhoods to determine the reconstruction path; according to the reconstruction path, scan the conditional probability of each direction of the binarized image obtained in step 1;

   The calculation and assignment module is used to iteratively calculate the value of each pixel in the empty image using Monte Carlo;

   Obtain a reconstructed image module, which is used to use porosity as a condition for the end of the cycle to obtain a reconstructed image;

   The result fit evaluation module is used to characterize the fit of the variation curve of the binarized image and the reconstructed image through data, and take the difference between the reconstructed image and the binarized image at the same step size, Accumulate the difference values under non-synchronization length to get the average value, the smaller the value, the better the fit.

Images