A Weight-Directed Edge-Preserving Filtering Method for Seismic Attributes

Abstract

The present invention discloses a seismic attribute edge-preserving filtering method based on weight guidance. In the present invention, the edge intensity values of each point in the attribute are initially calculated through a structure operator, which can preliminarily determine the approximate edge positions of each geological body in the attribute. Then, the optimal threshold method is used to perform binary segmentation on the edge intensity matrix, and the segmentation result is used as a leading matrix for subsequent filtering processing, which can distinguish the strong and weak edges of each geological body in the attribute and help with subsequent filtering processing. Secondly, using the leading matrix as a guide, non-linear filtering is performed on the area where the values jump at the edge of the geological body for edge-preserving processing to retain the jump trend, and linear filtering is performed on the area where the values transition smoothly inside the geological body to smooth out the differences between the values. This processing method can obtain an initial filtering result that retains the edge information of the geological body on the premise of suppressing noise. Finally, error analysis is performed on the initial filtering result, and this post-processing method can perform faithful restoration on the original data.

Description

Technical Field

[0001] The present invention belongs to the technical field of seismic exploration, and specifically relates to a seismic attribute edge-preserving filtering method based on weight guidance. Background Art

[0002] With the continuous improvement of seismic exploration level, the development of remaining oil within thin-bed reservoirs has entered a new stage. Extracting the continuity edges of various geological bodies from within thin-bed reservoirs is an important issue that urgently needs to be solved at present. And predicting geological body boundaries based on seismic plane attributes has gradually come into the view of researchers. However, noise inevitably exists in seismic plane attributes, which seriously affects the accuracy of boundary information prediction. Suppressing the noise in plane attribute data while maintaining geological body boundary information is of great significance for geological body boundary prediction.

[0003] However, in practical applications, the attribute values at the edges of geological bodies in seismic plane attributes tend to have jumps. In the filtering process, the value weight mainly affects the filtering effect here; within the geological body, the transition trend is relatively gentle, and in the filtering process, the spatial distance weight mainly affects the filtering effect here. Different change trends of attribute values require different types of filtering methods for filtering. Relevant research shows that in areas where the transition trend of attribute values is relatively gentle, linear filtering methods mainly considering spatial distance weights can help reduce the difference between values and reduce the influence of continuous noise; in areas where attribute values jump, the influence of attribute values needs to be considered first, increasing the weight of attribute values in the filtering process, and relevant non-linear filtering methods can retain the jump trend of values and reduce the influence of non-continuous noise. However, the existing plane attribute filtering technology methods do not have this characteristic, which limits the accuracy of geological body boundary prediction. Summary of the Invention

[0004] The purpose of the present invention is to provide a seismic attribute edge-preserving filtering method based on weight guidance to solve the above-mentioned problems.

[0005] The technical solution adopted by the present invention is as follows: A seismic attribute edge-preserving filtering method based on weight guidance includes the following steps:

[0006] Step S1: Input the seismic plane attribute At, and calculate the edge strength S of each geological body in the attribute data;

[0007] Step S2: Based on the above edge strength data S, use the optimal threshold segmentation method to obtain the threshold T, and perform binary segmentation on each point in the attribute data according to the threshold T to obtain the leading matrix Gd;

[0008] Step S3: According to the leading matrix Gd, perform classification filtering on the plane attribute At to obtain the initial filtered data At1;

[0009] Step S4: Based on the initial filtered data At1, perform error analysis according to the error coefficient θ to obtain the final filtered data At2.

[0010] In a preferred embodiment, the step S1 includes the following steps:

[0011] Step S11: Input the seismic plane attribute At and normalize it into grayscale data G in the range of 0 to 255. The calculation expression for each point is:

[0012]

[0013] where (x, y) are the horizontal and vertical coordinates, and max(At(x, y)) represents the maximum value in At;

[0014] Step S12: Calculate the horizontal edge intensity S of each geological body in the grayscale data G x , and its expression is:

[0015]

[0016] where * represents the convolution operation;

[0017] Step S13: Calculate the vertical edge intensity S of each geological body in the grayscale data G y , and its expression is:

[0018]

[0019] Step S14: Based on the horizontal edge intensity S x and the vertical edge intensity S y calculate the edge intensity S of each geological body in the attribute data. The calculation expression for each point is:

[0020]

[0021] In a preferred embodiment, the step S2 includes the following steps:

[0022] Step S21: Based on the edge intensity S of each geological body in the attribute data, calculate the global optimal binarization threshold T by the optimal threshold segmentation method;

[0023] Step S22: Segment the edge intensity data S according to the global optimal binarization threshold T to obtain the leading matrix Gd. The calculation expression for each point is:

[0024]

[0025] In a preferred embodiment, the step S21 includes the following steps:

[0026] Step S211: Based on the normalized image histogram I of S, calculate the zero-order cumulative moment A of S. The calculation expression for each gray value k is as follows:

[0027]

[0028] In the formula, I(i) represents the number of times the value i appears in S, where i and k are gray values, and i, k ∈ [0, 255];

[0029] Step S212: Based on the normalized image histogram I of S, calculate the first-order cumulative moment B of S. The calculation expression for each gray value k is as follows:

[0030]

[0031] Step S213: Based on the normalized image histogram I of S, calculate the cumulative moment M(255) of S when k = 255. The expression is as follows:

[0032]

[0033] Step S214: Calculate the average variance σ1 of each gray value in S based on A, B, and M(255). The expression is as follows:

[0034]

[0035] Step S215: Select the maximum value σ1(m) in σ1. At this time, m is the preferred global optimal threshold T. The expression is as follows:

[0036]

[0037] In a preferred embodiment, the step S3 includes the following steps:

[0038] Step S31: Construct a linear filtering template L n ;

[0039] Step S32: Construct a non-linear filtering template NL n ;

[0040] Step S33: Based on the leading matrix Gd, the linear filtering template L n and the non-linear filtering template NL n perform classification filtering on the seismic plane attribute At. Specifically, perform non-linear filtering at Gd(x, y) = 1 and linear filtering at Gd(x, y) = 0 to obtain the initial filtered data At1. The calculation expression for each point in At is as follows:

[0041]

[0042] Where At(i′,j) is a sliding window matrix with size n and At(x,y) as the target in the calculation.

[0043] In a preferred embodiment, step S4 comprises the following steps:

[0044] Step S41: Set the error coefficient θ, perform error analysis on the seismic plane attribute At and the initial filtering data At1, and perform data correction processing to obtain the final filtering result At2. The calculation expression for each point in At2 is:

[0045]

[0046] In summary, due to the adoption of the above technical solution, the beneficial effects of the present invention are:

[0047] In the present invention, the edge strength value of each point in the attribute is first calculated by the structural operator, and the approximate edge position of each geological body in the attribute can be preliminarily determined. Then, the edge strength matrix is binary segmented by the optimal threshold method, and the segmentation result is used as the leading matrix for subsequent filtering processing, which can distinguish the strong and weak edges of each geological body in the attribute, and help the subsequent filtering processing; secondly, the leading matrix is used as a guide, and nonlinear filtering is performed for the region where the value jumps at the edge of the geological body to perform edge protection processing, retain the jump trend, and linear filtering is performed in the region where the value transition is gentle inside the geological body to smooth the difference between values. This processing method can obtain the initial filtering processing result with the geological body edge information retained under the premise of suppressing noise; finally, the initial filtering result is subjected to error analysis, and this lag processing method can restore the original data to a certain extent. At the same time, it can also detect the edge of the geological body and use it as a guide to distinguish and process the internal and external information of the geological body in the attribute. This processing method can retain the edge and suppress the noise to a great extent, thereby providing a basic guarantee for the subsequent extraction and calibration of geological body edge information. DETAILED DESCRIPTION

[0048] In order to make the purpose, technical solution and advantages of the embodiments of the present invention clearer, the technical solution in the embodiments of the present invention will be clearly and completely described below in combination with the embodiments of the present invention. Obviously, the described embodiments are part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0049] Example:

[0050] A weight-guided seismic attribute edge-preserving filtering method comprises the following steps:

[0051] Step S1: Input the seismic plane attribute At, and calculate the edge strength S of each geological body in the attribute data;

[0052] Step S2: Based on the above edge strength data S, use the optimal threshold segmentation method to obtain the threshold T, and perform binary segmentation on each point in the attribute data according to the threshold T to obtain the leading matrix Gd;

[0053] Step S3: According to the leading matrix Gd, perform classification filtering on the seismic plane attribute At to obtain the initial filtered data At1;

[0054] Step S4: Based on the initial filtered data At1, perform error analysis according to the error constant σ to obtain the final filtered data At2.

[0055] The said step S1 includes the following steps:

[0056] Step S11: Input the seismic plane attribute At and normalize it to the grayscale data G of 0-255. The calculation expression for each point is:

[0057]

[0058] where (x, y) are the horizontal and vertical coordinates, and max(At(x, y)) represents the maximum value in At;

[0059] Step S12: Calculate the edge strength S of each geological body in the horizontal direction of the grayscale data G x , and its expression is:

[0060]

[0061] where * represents the convolution operation;

[0062] Step S13: Calculate the edge strength S of each geological body in the vertical direction of the grayscale data G y , and its expression is:

[0063]

[0064] Step S14: Based on the edge strength S in the horizontal direction x and the edge strength S in the vertical direction y calculate the edge strength S of each geological body in the attribute data. The calculation expression for each point is:

[0065]

[0066] The said step S2 includes the following steps:

[0067] Step S21: Based on the edge strength S of each geological body in the attribute data, calculate the global optimal binarization threshold T by the optimal threshold segmentation method;

[0068] Step S22: Segment the edge intensity data S according to the globally optimal binarization threshold T to obtain the pilot matrix Gd, and the calculation expression for each point thereof is:

[0069]

[0070] The said step S21 includes the following steps:

[0071] Step S211: Based on the normalized image histogram I of S, calculate the zero-order cumulative moment A of S, and the calculation expression for each gray value k thereof is:

[0072]

[0073] In the formula, I(i) represents the number of times the value i appears in S, i and k are gray values, and i, k ∈ [0, 255];

[0074] Step S212: Based on the normalized image histogram I of S, calculate the first-order cumulative moment B of S, and the calculation expression for each gray value k thereof is:

[0075]

[0076] Step S213: Based on the normalized image histogram I of S, calculate the cumulative moment M(255) of S when k = 255, and its expression is:

[0077]

[0078] Step S214: Calculate the average variance σ1 of each gray value in S based on A, B and M(255), and its expression is:

[0079]

[0080] Step S215: Select the maximum value σ1(m) in σ1, and at this time m is the preferred globally optimal threshold T, and its expression is:

[0081]

[0082] The said step S3 includes the following steps:

[0083] Step S31: Construct a linear filtering template L n ;

[0084] Step S32: Construct a non-linear filtering template NL n ;

[0085] Step S33: Based on the pilot matrix Gd, the linear filtering template L n and the non-linear filtering template NLn The seismic plane attribute At is subjected to classification filtering, specifically, nonlinear filtering is performed at Gd(x, y) = 1, and linear filtering is performed at Gd(x, y) = 0, to obtain the initial filtering data At1. The calculation expression for each point in At is:

[0086]

[0087] Where At(i′,j) is a sliding window matrix with size n and At(x,y) as the target in the calculation.

[0088] The step S4 comprises the following steps:

[0089] Step S41: Set the error coefficient θ, perform error analysis on the seismic plane attribute At and the initial filtering data At1, and perform data correction processing to obtain the final filtering result At2. The calculation expression for each point in At2 is:

[0090]

[0091] In the present invention, the edge strength value of each point in the attribute is first calculated by the structural operator, and the approximate edge position of each geological body in the attribute can be preliminarily determined. Then, the edge strength matrix is binary segmented by the optimal threshold method, and the segmentation result is used as the leading matrix for subsequent filtering processing, which can distinguish the strong and weak edges of each geological body in the attribute, and help the subsequent filtering processing; secondly, the leading matrix is used as a guide, and nonlinear filtering is performed for the region where the value jumps at the edge of the geological body to perform edge protection processing, retain the jump trend, and linear filtering is performed in the region where the value transition is gentle inside the geological body to smooth the difference between values. This processing method can obtain the initial filtering processing result with the geological body edge information retained under the premise of suppressing noise; finally, the initial filtering result is subjected to error analysis, and this lag processing method can restore the original data to a certain extent. At the same time, it can also detect the edge of the geological body and use it as a guide to distinguish and process the internal and external information of the geological body in the attribute. This processing method can retain the edge and suppress the noise to a great extent, thereby providing a basic guarantee for the subsequent extraction and calibration of geological body edge information.

[0092] It should be noted that in this text, relational terms such as first and second are only used to distinguish one entity or operation from another entity or operation, and do not necessarily require or imply any actual relationship or order between these entities or operations. Moreover, the term "comprising", "including" or any other variant thereof is intended to cover non-exclusive inclusion, so that a process, method, article or device comprising a series of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such process, method, article or device. Without further limitation, an element defined by the statement "comprising an..." does not exclude the presence of additional identical elements in the process, method, article or device comprising said element.

[0093] The above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit it; although the present invention has been described in detail with reference to the foregoing embodiments, those of ordinary skill in the art should understand that they can still modify the technical solutions recorded in the foregoing embodiments, or perform equivalent replacements for some of the technical features; and these modifications or replacements do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A seismic attribute edge-preserving filtering method based on weight guidance, characterized in that: It includes the following steps: Step S1: Input the seismic plane attribute At, and calculate the edge intensity S of each geological body in the attribute data; Step S2: Based on the above edge intensity S, use the optimal threshold segmentation method to obtain the threshold T, and perform binary segmentation on each point in the attribute data according to the threshold T to obtain the pilot matrix Gd; Step S3: According to the pilot matrix Gd, perform classification filtering on the seismic plane attribute At to obtain the initial filtered data At1; it includes the following steps: Step S31: Construct a linear filtering template L n ; Step S32: Construct a non-linear filtering template NL n ; Step S33: Based on the pilot matrix Gd, the linear filtering template L n and the non-linear filtering template NL n Perform classification filtering on the seismic plane attribute At. Specifically, perform non-linear filtering when Gd(x,y)=1 and linear filtering when Gd(x,y)=0 to obtain the initial filtered data At1. The calculation expression for each point in At is as follows: In the formula, At(i′,j) is a sliding window matrix with a size of n taking At(x,y) as the target point in the calculation; Step S4: Based on the initial filtered data At1, perform error analysis according to the error coefficient θ to obtain the final filtered data At2; it includes the following steps: Step S41: Set the error coefficient θ, perform error analysis on the seismic plane attribute At and the initial filtered data At1, and perform data correction processing to obtain the final filtering result At2. The calculation expression for each point in At2 is:

2. The edge-preserving filtering method for seismic attributes based on weight guidance according to claim 1, wherein: The said step S1 includes the following steps: Step S11: Input the seismic plane attribute At and normalize it into grayscale data G in the range of 0 - 255. The calculation expression for each point is: In the formula, (x,y) are the horizontal and vertical coordinates, and max(At(x,y)) represents the maximum value in At; Step S12: Calculate the horizontal edge intensity S of each geological body in the grayscale data G x , and its expression is: In the formula, * represents the convolution operation; Step S13: Calculate the vertical edge intensity S of each geological body in the grayscale data G y , and its expression is: Step S14: Based on the edge intensity S in the horizontal direction x and the edge intensity S in the vertical direction y calculate the edge intensity S of each geological body in the attribute data, and the calculation expression for each point is:

3. A seismic attribute edge-preserving filtering method based on weight guidance according to claim 1, characterized in that: The said step S2 includes the following steps: Step S21: Based on the edge intensity S of each geological body in the attribute data, calculate the global optimal binary threshold T by the optimal threshold segmentation method; Step S22: Segment the edge intensity S according to the global optimal binary threshold T to obtain the pilot matrix Gd. The calculation expression for each point is:

4. The edge-preserving filtering method for seismic attributes based on weight guidance according to claim 3, characterized in that: The said step S21 includes the following steps: Step S211: Based on the normalized image histogram I of S, calculate the zero - order cumulative moment A of S. The calculation expression for each gray value k is: In the formula, I(i) represents the number of times the value i appears in S, and i,k are gray values, i,k ∈ [0,255]; Step S212: Based on the normalized image histogram I of S, calculate the first - order cumulative moment B of S. The calculation expression for each gray value k is: Step S213: Based on the normalized image histogram I of S, calculate the cumulative moment M(255) of S when k = 255. The expression is: Step S214: Calculate the average variance σ1 of each gray value in S based on A, B, and M(255). The expression is: Step S215: Select the maximum value σ1(m) in σ1. At this time, m is the global optimal threshold T. The expression is:

Images