A method for constructing low-frequency model of wave impedance inversion based on deep fusion of seismic multi-attributes

Abstract

The present invention discloses a method for constructing a low-frequency model for wave impedance inversion based on deep fusion of multiple seismic attributes. The method includes establishing a low-frequency model for well interpolation, then performing post-stack inversion to obtain a bandpass wave impedance inversion volume, and extracting other attribute volumes: a three-instantaneous attribute volume, a time-windowed frequency attribute volume, and a frequency-divided seismic volume. Subsequently, pseudo-well curves are extracted, and the wave impedance volume is obtained and a residual volume is generated through deep learning, ultimately obtaining a highly accurate low-frequency model. Through deep learning, the above method integrates multiple geologically significant attribute volumes, effectively reducing the uncertainty introduced by well interpolation, improving the accuracy of the low-frequency model, and providing a good foundation for high-precision seismic wave impedance inversion.

Description

Technical Field

[0001] The present invention relates to the technical field of seismic exploration for oil and gas, and more specifically, to a method for constructing a low-frequency model of wave impedance inversion based on deep fusion of seismic multi-attributes. Background Art

[0002] Currently, inversion is widely used for qualitative or quantitative reservoir prediction during the exploration phase and for reserve calculation, well pattern deployment, and reservoir dynamic monitoring during the development phase. It is an essential method for characterizing subsurface reservoirs. With the continuous advancement of exploration and development, structurally simple oil and gas reservoirs are becoming increasingly rare, and the subsurface landscape is becoming increasingly complex, requiring increasingly sophisticated oil and gas field research. Seismic impedance inversion converts seismic data into an impedance form that can be directly compared with well logging data. It is an important tool for detailed research on oil and gas field structure, reservoirs, and fluids. The construction of a low-frequency model is a key step in seismic inversion, and its accuracy can significantly impact the inversion results, especially in complex structural geological environments. The low-frequency model is derived by interpolating and extrapolating well logging information across the entire data volume, based on the constraints of seismic interpretation horizons and sedimentary patterns. This model compensates for the missing low-frequency information in the seismic data to a certain extent.

[0003] Conventional low-frequency models are constructed by directly interpolating well logs within a stratigraphic framework. These models lack both low- and high-frequency components, leading to significant uncertainty in the inversion results. Conventional low-frequency models lack geological and physical significance, and the interpolation results carry significant uncertainty, leading to strong ambiguity in the inversion results. This results in low precision and hinders detailed research on oil and gas fields. Summary of the Invention

[0004] The purpose of the present invention is to overcome the defects of poor accuracy and large uncertainty of the existing wave impedance inversion low-frequency model, and to provide a wave impedance inversion low-frequency model construction method based on deep fusion of seismic multi-attributes. This method can effectively improve the accuracy of the seismic wave impedance inversion low-frequency model, thereby obtaining high-precision inversion results.

[0005] In order to solve the above technical problems, the technical solution adopted by the present invention is:

[0006] The present invention provides a method for constructing a low-frequency model of wave impedance inversion based on deep fusion of seismic multi-attributes, comprising the following steps:

[0007] Step S1: constructing a stratigraphic framework using the marker layers picked up from the seismic data, and spatially interpolating the filtered logging data based on the stratigraphic framework to establish a low-frequency model of well interpolation wave impedance;

[0008] Step S2: performing post-stack inversion using the low-frequency model established in step S1 to obtain a bandpass impedance inversion volume;

[0009] Step S3, extracting other attribute bodies: three-moment attribute body, time window frequency attribute body and frequency-divided earthquake body;

[0010] Step S4, extracting pseudo-well curves: extracting pseudo-well curves of the well interpolation low-frequency model, the bandpass wave impedance inversion volume and other attribute volumes generated in steps S1 to S3 at the well point;

[0011] Step S5: using the pseudo-well curve extracted in step S4, obtaining the wave impedance body through deep learning;

[0012] Step S6: Generate a residual volume using the wave impedance volume obtained in step S5;

[0013] Step S7, obtaining the final low-frequency model: adding the wave impedance volume obtained in step S5 and the residual volume obtained in step S6 to obtain the final high-precision low-frequency model.

[0014] The wave impedance in step S1 is the ratio of the pressure acting on a certain area to the particle flow rate perpendicular to this area per unit time (i.e., the area multiplied by the particle vibration velocity) when the seismic wave propagates in the medium. It has the meaning of resistance and is calculated using the following formula:

[0015] imp=v×ρ

[0016] Where imp is the wave impedance, v and ρ are the velocity of the seismic wave and the medium density, respectively.

[0017] Furthermore, step S1 includes the following steps:

[0018] Step S11: filtering the well logging wave impedance curve using a low-pass filter to retain effective frequency components below 10 Hz;

[0019] Step S12: Use the inverse distance weighted method as the spatial interpolation method for well logging data. The specific calculation formula is:

[0020]

[0021] Among them, (x, y) is the coordinate of the interpolation point, (x i ,y i ) is the coordinate of the discrete point, f(x, y) is the interpolation function, n is the total number of discrete points, i and j are the discrete point numbers, From (x, y) to (x i ,y i ) horizontal distance, p is a constant greater than 0, called the weighted power index (usually, p = 2), Z i 、Zj are the vertical heights of the i-th and j-th discrete points respectively.

[0022] Furthermore, step S2 includes the following steps:

[0023] Step S21, extracting target layer statistical wavelet using original seismic data;

[0024] Step S22: using the well interpolation low-frequency model established in step S1, combined with statistical wavelets, seismic data, and well logging data, to perform post-stack impedance inversion;

[0025] Step S23: performing band-pass filtering on the wave impedance inversion result to obtain a band-pass wave impedance inversion volume.

[0026] Furthermore, the frequency components retained by the bandpass filtering in step S23 are 5 Hz to 60 Hz.

[0027] Furthermore, the three instantaneous attributes in step S3 are instantaneous amplitude, instantaneous phase and instantaneous frequency, and the specific calculation formulas are as follows:

[0028]

[0029]

[0030]

[0031] where A(t), φ(t) and ω(t) are the instantaneous amplitude, instantaneous phase and instantaneous frequency, respectively; s(t) is the seismic trace; h(t) is the Hilbert transform of s(t); t is time; and d is the difference operator.

[0032] Furthermore, the time window frequency attribute in step S3 is a main frequency attribute and an average frequency attribute, which specifically includes the following steps:

[0033] Step S31: Convert the seismic volume from the time domain to the frequency domain through Fourier transform. The formula of Fourier transform is:

[0034]

[0035] Where ω is frequency, t is time, e -iωt represents a complex function;

[0036] Step S32: Set a time window, calculate the dominant frequency or average frequency within the time window, and place the value in the middle of the time window. Then, move backward one step at a time until the complete individual is traversed, and the dominant frequency attribute body and the average frequency attribute body can be obtained.

[0037] Furthermore, the frequency division seismic volume in step S3 is obtained by frequency division of the seismic volume, and is evenly divided into 6 frequency division volumes within the effective seismic frequency band. Of course, the number of frequency division volumes can also be increased or decreased accordingly according to actual conditions. The frequency division algorithm is wavelet frequency division, and the specific calculation formula is as follows:

[0038]

[0039] Among them, g(t) is the basic wavelet, s(t) is the seismic trace, a is a non-zero real number as a scale factor, and b is a real number as a translation factor. It means taking the complex conjugate of g(t).

[0040] Furthermore, step S5 includes the following steps:

[0041] Step S51: low-pass filter the wave impedance curve to filter out frequency components above 60 Hz;

[0042] Step S52: Deep learning is performed using the pseudo-well curves extracted in step S4 and the filtered wave impedance curves to obtain a deep learning network between the pseudo-well curves and the filtered wave impedance curves. This network is then applied to the data volume generated in steps S1 to S3 to obtain a wave impedance volume. The wave impedance curve is obtained by multiplying the well logging velocity curve by the well logging density curve.

[0043] Furthermore, the deep learning algorithm is a deep feedforward neural network, and its forward propagation formula is as follows:

[0044]

[0045]

[0046] in, is the output of the jth neuron in layer l, is the value of the i-th neuron in the l+1 layer before being activated by the function, is the weight between the jth neuron in layer l and the ith neuron in layer l+1, is the bias, f is the nonlinear activation function, and i and j are the numbers of neurons in each layer.

[0047] Furthermore, step S6 includes the following steps:

[0048] Step S61, extracting the pseudo-well curve of the wave impedance body obtained in step S5;

[0049] Step S62: Subtract the target curve filtered in step S5 from the wave impedance pseudo-well curve to obtain a residual curve;

[0050] Step S63: Based on the stratigraphic grid in step S1, a residual volume is obtained by using a residual curve spatial interpolation algorithm. The residual curve spatial interpolation algorithm is an inverse distance weighted method, and the specific calculation formula is as follows:

[0051]

[0052] Among them, (x, y) is the coordinate of the interpolation point, (x i ,y i ) is the coordinate of the discrete point, f(x, y) is the interpolation function, n is the total number of discrete points, i, j are the discrete point numbers, From (x, y) to (x i ,y i ) horizontal distance, p is a constant greater than 0, called the weighted power index, Z i 、Z j are the vertical heights of the i-th and j-th discrete points respectively.

[0053] Compared with the existing technology, the beneficial effect of the present invention lies in that, based on deep learning, by introducing bandpass inversion results, the frequency components of the low-frequency model are effectively supplemented, and a variety of seismic attribute bodies with geological characteristics are integrated, which is equivalent to adding constraints of geological knowledge. While reducing the multi-solution of seismic wave impedance inversion, it also effectively improves the accuracy of seismic wave impedance inversion. BRIEF DESCRIPTION OF THE DRAWINGS

[0054] Figure 1 This is a flow chart of the method for constructing a low-frequency model using wave impedance inversion according to the present invention;

[0055] Figure 2 This is the stratigraphic grid diagram of the low-frequency model of the study area in Example 1;

[0056] Figure 3 This is a low-frequency model diagram of the well interpolation in Example 1;

[0057] Figure 4 This is the final low-frequency model diagram in Example 1. DETAILED DESCRIPTION

[0058] The present invention is further described below with reference to specific embodiments. The accompanying drawings are for illustrative purposes only and are not to be construed as limiting this patent. Certain components in the accompanying drawings may be omitted, enlarged, or reduced in size to better illustrate the embodiments, and do not represent the actual dimensions of the product. It is understandable to those skilled in the art that certain well-known structures and their descriptions may be omitted from the accompanying drawings.

[0059] The same or similar reference numerals in the drawings of the embodiments of the present invention correspond to the same or similar components; in the description of the present invention, it should be understood that if there are terms such as "front", "rear", "left", "right", etc. indicating an orientation or positional relationship, they are based on the orientation or positional relationship shown in the drawings. This is only for the convenience of describing the present invention and simplifying the description, and does not indicate or imply that the device or element referred to must have a specific orientation, be constructed and operate in a specific orientation. Therefore, the terms describing the positional relationship in the drawings are only used for illustrative purposes and cannot be understood as limiting this patent. For ordinary technicians in this field, the specific meanings of the above terms can be understood according to specific circumstances. In addition, in the present invention, descriptions such as "first", "second", etc. are only for descriptive purposes and cannot be understood as indicating or implying their relative importance or implicitly indicating the number of technical features indicated. Therefore, the features defined as "first" or "second" can explicitly or implicitly include at least one such feature.

[0060] Example 1:

[0061] This example uses the construction of a low-frequency model using seismic impedance inversion in a gas field in the Yinggehai Basin as an example. The gas reservoir in the study area is a non-classical slope-break gravity flow depositional system. Reconstructed by later water flow, the sand body distribution and contact relationships are very complex. Conventional low-frequency models have significant uncertainty, resulting in low inversion accuracy and hindering detailed description of the gas reservoir. Figure 1 This embodiment provides a method for constructing a low-frequency model using wave impedance inversion based on deep fusion of multiple seismic attributes. This method effectively improves the accuracy of the low-frequency model using seismic wave impedance inversion, laying a good foundation for subsequent inversion research and detailed description of gas reservoirs. The method includes the following steps:

[0062] Step S1, refer to Figure 2 , using the marker layers picked up from the seismic data to build a stratigraphic grid (in order to comprehensively consider the time and accuracy of model establishment, the vertical sampling rate of the stratigraphic grid is generally smaller in the target layer and appropriately larger in the non-target layer), and based on the stratigraphic grid, the filtered logging data are spatially interpolated, see Figure 3 , establish a low-frequency model of well interpolation wave impedance;

[0063] Step S2: performing post-stack inversion using the low-frequency model to obtain a bandpass impedance inversion volume;

[0064] Step S3, extracting other attribute bodies: three-moment attribute body, time window frequency attribute body and frequency-divided earthquake body;

[0065] Step S4, extracting pseudo-well curves: extracting pseudo-well curves of the well interpolation low-frequency model, the bandpass wave impedance inversion volume and other attribute volumes generated in steps S1 to S3 at the well point;

[0066] Step S5: using the pseudo-well curve to obtain a wave impedance body through deep learning;

[0067] Step S6: Generate a residual volume using the wave impedance volume obtained in the previous step;

[0068] Step S7, obtain the final low-frequency model: add the wave impedance body obtained in step S5 and the residual body obtained in step S6, and finally obtain the following Figure 4 High-precision low-frequency model shown.

[0069] Compared with the conventional low-frequency model, which brings great uncertainty to the inversion results, this embodiment, based on deep learning, effectively supplements the frequency components of the low-frequency model by introducing bandpass inversion results, and also integrates a variety of seismic attribute bodies with geological characteristics, which is equivalent to adding constraints of geological knowledge. While reducing the multi-solution of seismic wave impedance inversion, it also effectively improves the accuracy of seismic wave impedance inversion.

[0070] Example 2:

[0071] Based on the first embodiment, the step S1 performs filtering on the logging wave impedance curve by low-pass filtering, retaining only the frequency components below 10 Hz; after filtering, the inverse distance weighted method is used as the spatial interpolation method of the logging data, and the specific calculation formula is:

[0072]

[0073] Among them, (x, y) is the coordinate of the interpolation point, (x i ,y i ) is the coordinate of the discrete point, f(x, y) is the interpolation function, n is the total number of discrete points, i, j are the discrete point numbers, From (x, y) to (x i ,y i ) horizontal distance, p is a constant greater than 0, called the weighted power index, Z i 、Z j are the vertical heights of the i-th and j-th discrete points respectively.

[0074] Step S2 is to perform post-stack inversion to obtain a bandpass wave impedance inversion volume: first, the statistical wavelet of the target layer is extracted using the original seismic data; then, the well interpolation low-frequency model established in step S1 is used to combine the statistical wavelet, seismic data and logging data to perform post-stack wave impedance inversion; finally, the wave impedance inversion result is bandpass filtered, and the frequency components retained by the bandpass filter are 5 Hz to 60 Hz to obtain a bandpass wave impedance inversion volume.

[0075] Example 3:

[0076] Based on the first or second embodiment, step S3 extracts other attribute bodies: three instantaneous attribute bodies, time window frequency attribute body, and frequency-divided seismic body. The three instantaneous attributes are instantaneous amplitude, instantaneous phase, and instantaneous frequency, and the corresponding calculation formulas are as follows:

[0077]

[0078]

[0079]

[0080] where A(t), φ(t) and ω(t) are the instantaneous amplitude, instantaneous phase and instantaneous frequency, respectively; s(t) is the seismic trace; h(t) is the Hilbert transform of s(t); t is time; and d is the difference operator.

[0081] The time window frequency attributes are the main frequency attribute and the average frequency attribute. The seismic volume needs to be converted from the time domain to the frequency domain through Fourier transform. The formula of Fourier transform is:

[0082]

[0083] Where ω is frequency, t is time, e -iωt Represents a complex function.

[0084] Then, a time window is set, the dominant frequency or average frequency within the time window is calculated, and the value is placed in the middle of the time window. Then, a step is moved backward in sequence until the complete individual is traversed, and the dominant frequency attribute body and the average frequency attribute body can be obtained.

[0085] The frequency division seismic volume is obtained by frequency division of the seismic volume, which is divided into 6 frequency division volumes within the effective seismic frequency band. Of course, the number of frequency division volumes can also be increased or decreased according to the actual situation. The frequency division algorithm is wavelet frequency division, and the specific calculation formula is as follows:

[0086]

[0087] Among them, g(t) is the basic wavelet, s(t) is the seismic trace, a is a non-zero real number as a scale factor, and b is a real number as a translation factor. It means taking the complex conjugate of g(t).

[0088] Example 4:

[0089] Based on any of the above embodiments, in this embodiment, step S5 obtains the wave impedance body through deep learning, which first performs low-pass filtering on the wave impedance curve to filter out frequency components above 60 Hz, and then uses the pseudo-well curve extracted in step S4 and the filtered wave impedance curve to perform deep learning to obtain a deep learning network between the pseudo-well curve and the filtered wave impedance curve, and applies the network to the data body generated by steps S1 to S3 to obtain the wave impedance body.

[0090] The deep learning algorithm is a deep feedforward neural network, and its forward propagation formula is as follows:

[0091]

[0092]

[0093] in, is the output of the jth neuron in layer l, is the value of the i-th neuron in the l+1 layer before being activated by the function, is the weight between the jth neuron in layer l and the ith neuron in layer l+1, is the bias, f is the nonlinear activation function, and i and j are the numbers of neurons in each layer.

[0094] Next, the target curve after filtering in step S5 is subtracted from the pseudo-well curve of the wave impedance volume obtained in step S5 to obtain a residual curve. Then, based on the stratigraphic grid in step S1, a residual volume is obtained by using a residual curve spatial interpolation algorithm; the residual curve spatial interpolation algorithm is an inverse distance weighted method, and the specific calculation formula is as follows:

[0095]

[0096] Among them, (x, y) is the coordinate of the interpolation point, (x i ,y i ) is the coordinate of the discrete point, f(x, y) is the interpolation function, n is the total number of discrete points, i, j are the discrete point numbers, From (x, y) to (x i ,y i ) horizontal distance, p is a constant greater than 0, called the weighted power index, Z i 、Z j are the vertical heights of the i-th and j-th discrete points respectively.

[0097] Finally, the wave impedance volume obtained in step S5 is added to the residual volume obtained in step S6 to obtain a final high-precision low-frequency model.

[0098] In the specific contents of the above-mentioned specific implementation methods, the various technical features can be combined in any non-contradictory manner. In order to make the description concise, not all possible combinations of the above-mentioned technical features are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this specification.

[0099] Obviously, the above embodiments of the present invention are merely examples for the purpose of clearly illustrating the present invention, and are not intended to limit the embodiments of the present invention. Those skilled in the art will appreciate that other variations or modifications can be made based on the above description. It is not necessary and impossible to enumerate all embodiments here. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of the present invention shall be included within the scope of protection of the claims of the present invention.

Claims

1. A method for constructing a low-frequency model of wave impedance inversion based on deep fusion of seismic multi-attributes, characterized in that: The following steps are involved: Step S1: constructing a stratigraphic framework using the marker layers picked up from the seismic data, and spatially interpolating the filtered logging data based on the stratigraphic framework to establish a low-frequency model of well interpolation wave impedance; Step S2: performing post-stack inversion using the model established in step S1 to obtain a bandpass impedance inversion volume; Step S3, extracting other attribute bodies: three-moment attribute body, time window frequency attribute body and frequency-divided earthquake body; Step S4, extracting pseudo-well curves: extracting pseudo-well curves of the well interpolation low-frequency model, the bandpass wave impedance inversion volume and other attribute volumes generated in steps S1 to S3 at the well point; Step S5: using the pseudo-well curve extracted in step S4, obtaining a wave impedance body through deep learning; Step S6: Generate a residual volume using the wave impedance volume obtained in step S5; Step S7, obtaining the final low-frequency model: adding the wave impedance volume obtained in step S5 and the residual volume obtained in step S6 to obtain the final low-frequency model; Wherein, the step S5 includes the following steps: Step S51: low-pass filter the wave impedance curve to filter out frequency components above 60 Hz; Step S52: Perform deep learning using the pseudo-well curve extracted in step S4 and the filtered wave impedance curve to obtain a deep learning network between the pseudo-well curve and the filtered wave impedance curve, and apply the network to the data volume generated in steps S1 to S3 to obtain a wave impedance volume.

2. The method for constructing a low-frequency model based on wave impedance inversion and deep fusion of multiple seismic attributes according to claim 1 is characterized in that: The step S1 comprises the following steps: Step S11: The marker layer is a seismic event axis with stable reflection in the seismic data and traceable interpretation in the entire area; Step S12: filtering the well logging wave impedance curve using a low-pass filter to retain frequency components below 10 Hz; Step S13: Use the inverse distance weighted method as the spatial interpolation method for well logging data. The specific calculation formula is: Where (x, y) is the coordinate of the interpolation point, (x i ,y i ) is the coordinate of the discrete point, f(x, y) is the interpolation function, n is the total number of discrete points, i, j are the discrete point numbers, From (x, y) to (x i ,y i ) horizontal distance, p is a constant greater than 0, called the weighted power index, Z i , Z j are the vertical heights of the i-th and j-th discrete points respectively.

3. The method for constructing a low-frequency model based on wave impedance inversion and deep fusion of seismic multi-attributes according to claim 1 is characterized in that: The step S2 comprises the following steps: Step S21, extracting target layer statistical wavelet using original earthquake data; Step S22: using the well interpolation low-frequency model established in step S1, combined with statistical wavelets, seismic data, and well logging data, to perform post-stack impedance inversion; Step S23: performing band-pass filtering on the wave impedance inversion result to obtain a band-pass wave impedance inversion volume.

4. The method for constructing a low-frequency model based on wave impedance inversion and deep fusion of seismic multi-attributes according to claim 3 is characterized in that: The frequency components retained by the bandpass filtering in step S23 are 5 Hz to 60 Hz.

5. The method for constructing a low-frequency model based on wave impedance inversion and deep fusion of seismic multi-attributes according to claim 1 is characterized in that: The three instantaneous attributes in step S3 are instantaneous amplitude, instantaneous phase and instantaneous frequency. The specific calculation formulas are as follows: where A(t), φ(t) and ω(t) are the instantaneous amplitude, instantaneous phase and instantaneous frequency, respectively; s(t) is the seismic trace; h(t) is the Hilbert transform of s(t); t is time; and d is the difference operator.

6. The method for constructing a low-frequency model based on wave impedance inversion and deep fusion of seismic multi-attributes according to claim 5 is characterized in that: The time window frequency attributes in step S3 are the main frequency attribute and the average frequency attribute, which specifically includes the following steps: Step S31: Fourier transform Convert the seismic volume from the time domain to the frequency domain; where ω is the frequency, e -iωt represents a complex function; Step S32: Set a time window, calculate the dominant frequency or average frequency within the time window, and place the value in the middle of the time window. Then, move backward one step at a time until the complete individual is traversed, and the dominant frequency attribute body and the average frequency attribute body can be obtained.

7. The method for constructing a low-frequency model based on wave impedance inversion and deep fusion of multiple seismic attributes according to claim 5 is characterized in that: The frequency division seismic volume in step S3 is obtained by frequency division of the seismic volume. The seismic volume is divided into 6 frequency division volumes within the effective seismic frequency band. The frequency division algorithm is wavelet frequency division. The specific calculation formula is as follows: Where g(t) is the basic wavelet, a is a non-zero real number as a scale factor, and b is a real number as a translation factor. It means taking the complex conjugate of g(t).

8. The method for constructing a low-frequency model based on wave impedance inversion and deep fusion of seismic multi-attributes according to claim 1 is characterized in that: The deep learning network is a deep feedforward neural network, and its forward propagation formula is as follows: Where, is the output of the jth neuron in layer l, is the value of the i-th neuron in the l+1 layer before being activated by the function, is the weight between the jth neuron in layer l and the ith neuron in layer l+1, is the bias, f is the nonlinear activation function, and i and j are the numbers of neurons in each layer.

9. The method for constructing a low-frequency model based on wave impedance inversion and deep fusion of multiple seismic attributes according to claim 1, characterized in that: The step S6 comprises the following steps: Step S61, extracting the pseudo-well curve of the wave impedance body obtained in step S5; Step S62: Subtract the wave impedance curve filtered in step S5 from the wave impedance pseudo-well curve to obtain a residual curve; Step S63: Based on the stratigraphic grid in step S1, a residual volume is obtained by using a residual curve spatial interpolation algorithm. The residual curve spatial interpolation algorithm is an inverse distance weighted method, and the specific calculation formula is as follows: Where (x, y) is the coordinate of the interpolation point, (x i ,y i ) is the coordinate of the discrete point, f(x,y) is the interpolation function, n is the total number of discrete points, i, j are the discrete point numbers, From (x,y) to (x i ,y i ) horizontal distance, p is a constant greater than 0, called the weighted power index, Z i , Z j are the vertical heights of the i-th and j-th discrete points respectively.

Images