KR100565773B1 - Seismic imaging method with reverse-time migration using poynting vector - Google Patents

Abstract

The present invention relates to a method of imaging seismic data to improve the quality of the image in consideration of the propagation direction of the wave in the method of imaging the seismic sensing data using the inverse time structure correction. The technical task is to remove the noise generated under the conventional imaging conditions and to realize a clear seismic image, and calculate the wave propagation direction as a pointing vector to observe the wave propagated from the sound wave and back propagated. The cross-correlation is calculated only when the angle between the received waves is within a certain range, and the cross-correlation values due to the wave having the undesired reflection angle are not used as the imaging data, thereby eliminating distorted noise. It features.

Seismic wave, image, propagation direction, sound wave, observation wave, pointing vector, cross-correlation, noise reduction

Description

Seismic Imaging Method with Reverse-time Migration Using Poynting Vector}

Figure 1a is a photograph of the amplitude of the wave field S (x, y, z, t) propagated from the sound source at all points in the basement.

1b is a photograph of the amplitude of the observed wave R (x, y, z, t) propagated back at the location of the receiver;

Figure 2 is a flow chart for performing a reverse time structure correction in conventional depth structure correction before polymerization.

FIG. 3 is an image photograph of noise (white arrow) caused by cross-correlation with waves propagating along the horizontal strata by a conventional inverse time-correction method.

4 is a cross-correlation of a sound source wave or a reverse wave propagation wave reflected from the interface of a medium having a large velocity change and propagated to the earth by a reverse time structure correction method with a reverse wave propagated wave or sound wave propagating downward from the surface. Image of noise (white arrow) caused by

FIG. 5 shows an example of wave amplitude and propagation direction calculation of a wave using the same in a one-dimensional medium. FIG.

FIG. 6 is an image photograph showing values of a sound wave S (x, y, t) and a pointing vector indicating a propagation direction of a sound wave at 0.5, 1.0, and 1.5 seconds in a two-dimensional medium. FIG.

FIG. 7 is an image showing observation wave R (x, y, t) propagated in 0.5 second, 1.0 second and 1.5 second in a two-dimensional medium and a pointing vector indicating the propagation direction of the observation wave. Picture.

8A is a photographic image of conventional inverse time structure correction for 1 and 11 sound source data.

8b is a photographic image of inverse time structure correction according to the present invention for 1 and 11 sound source data.

Figure 9a is a photographic image of performing a conventional reverse time structure correction for the artificial exploration data.

9b is an image photograph of the reverse time structure correction according to the present invention for the artificial exploration data.

10 is a flowchart for imaging seismic data by performing inverse time structure correction of the present invention in consideration of the propagation direction of waves.

The present invention relates to a method for imaging seismic sensing data, and more particularly, in a method of imaging seismic sensing data using reverse time structure correction, seismic data for improving the image quality in consideration of the propagation direction of waves. To an imaging method.

In general, the imaging of seismic sensing data enables the analysis of the structure of the strata easily and easily by visualizing the data obtained during the exploration, and the seismic data is generated by the wave generated from the sound source through the strata. The underground structure is imaged based on the time and amplitude to reach.

In general, the imaging of seismic waves is based on the inaccurate initial velocity model, and the methods of realizing the best images through structural migration have been proposed. Such structural correction methods include Kirchhoff type Migration, One-way wave equation type migration or Reverse-time Migration (two-way wave equation type Migration) is proposed.

The imaging of the seismic wave is based on the inaccurate initial velocity model, and the depth model is corrected by analyzing the first image implemented from the depth pre-polymerization. By repeating this process, depth polymerization before polymerization is performed until a satisfactory underground image is obtained.

In realizing the above-described depth polymerization structure correction image, the conventional inverse temporal structure correction includes the sound wave S (x, y, z, t) propagated at the position of the sound source and the observation wave R (x propagated at the position of the receiver. Using an imaging condition of zero-lag cross correlation with , y, z, t) with the same weight "1.0" at all underground points, this imaging condition is It is expressed as (1).

In the above equation, I (x, y, z) is the structure corrected image value, Tmax is the last recording time recorded in the receiver, S (x, y, z, t) is the sound wave propagated at the location of the sound source, R (x, y, z, t) is the back-propagation in the position of sujingi observation wave, the above equation (1) ti = 0 from ti = Tmax S (x up, y, z, t) and R (x, y, z, t) , multiply the two wave fields ( S and R ) at the same time zone ti , and add them together to add the structural correction image I (x, y, z) (the form of the value is amplitude (+) , Having any size of (-)).

1 is an image photograph of correlated wave fields of inverse time structure correction in depth structure correction before polymerization. Figure 1a is a photograph of the amplitude of the sound source wave S (x, y, z, t) propagated from the sound source at all points in the basement, Figure 1b is the observed wave R (x, y, z propagated back at the position of the receiver is a photograph of the amplitude of , t) at all points in the basement. The reverse time-depth structure correction before the polymerization is performed by recording S (x, y, z, t) and R (x, y, z, t) at regular time intervals and then correlating each other with zero delay. Imaging. In the conventional structural correction, cross correlation is calculated by multiplying S (x, y, z, t) and R (x, y, z, t) by a weight of 1.0 at all x, y, z points.

On the other hand, Figure 2 is a flow chart for performing inverse time structure correction in the conventional depth-of-polymerization depth structure correction.

Referring to the drawings, a step (S10) of collecting sound source collection data that has undergone basic processing (deconvolution, polymerization rate analysis, static correction, filtering, etc.) and determining a speed model;

Calculating and storing a sound source wave S (x, y, z, t) at a predetermined time interval Δt (S20);

Determining whether time t is greater than or equal to Tmax and repeating step S20 if small (S30);

Calculating the reverse propagated R (x, y, z, t) when it is determined in step S30 that time t is greater than or equal to Tmax (S40);

Calculating cross-correlation S (x, y, z, t) * R (x, y, z, t) at all x, y, z points (S50);

Determining whether the time t is greater than or equal to Tmax and repeating steps S40 and S50 if small (S60); And

If it is determined in step S60 that the time t is greater than or equal to Tmax , the result of adding S (x, y, z, t) * R (x, y, z, t) to the value ( I (x, y, z) )) as a step (S70) for generating a corrected image of field structure section; as made.

In the above step S70, the value of the result value I (x, y, z) at a specific point x, y, z is represented as an integer or real value of plus or minus (-), and this value is implemented as a seismic data processing image. The program is imaged using a program such as SU (Seismic Unix, Colorado School of Mine) package tool, SEP (Stanford Exploration Project), Geoframe or Promax. If the value I (x, y, z) is positive, the image is dark black. If the value I (x, y, z) has a negative value, the image is expressed in a color close to white. This can be done using other common imaging methods.

However, in order to image the seismic data as described above, the conventional imaging using the inverse time-depth structure correction before polymerization cross-correlates the sound wave and the reverse propagated observation wave at all points in the basement.

1) propagation along the boundary of the strata, due to noise due to cross correlation with head waves propagating to the ground or propagating to the surface in areas where the velocity increases gradually with the depth of the strata. The quality is deteriorated. In other words, the image is distorted by waves that are not necessary. An example of such noise is shown in FIG. 3, which is a conventional method of propagating along the boundary of a horizontal strata and generating noise (white arrow) due to cross-correlation with waves propagating to the earth's surface. It is a video photograph.

2) The quality of the structure correction image is degraded by the cross-correlation noise between the reflected wave propagated to the surface and reflected from the interface of the medium with large velocity changes. An example of such a noise is illustrated in FIG. 4, which is a view of a sound wave or an observation wave propagated downward from the surface by a sound wave or an observation wave reflected from the boundary of a medium having a large velocity change by the conventional method and propagating to the ground. Image of noise (white arrow) caused by cross-correlation of.

As seen from the above, in the conventional imaging of the seismic data, the above noise is generated and the quality of the image is degraded. Therefore, there is a problem in that the stratum analysis cannot be performed smoothly. This noise is characterized by the propagation direction of the two correlated waves being nearly opposite.

The present invention is the noise generated in the imaging of the conventional acoustic wave data, that is, the head waves propagating to the surface in the region where the speed propagates along the boundary surface of the strata, or the speed is gradually increased with the depth of the strata. A clear seismic image can be realized by removing the noise due to cross correlation between the reflected wave propagated to the surface and reflected from the interface of the medium with large noise or speed change. It is a technical task.

The technical problem of the present invention can be achieved by using the seismic wave data reflected from the basement to the image through the structure correction, the rest of the wave data to make the cross-correlation value "0" so that it is not used for imaging.

That is, in the imaging of the seismic data, the noises appearing as shown in FIG. 3 and FIG. 4 are characterized in that the propagation directions of the two correlated wave fields are close to the opposite directions, but the seismic data reflected from the basement is structurally corrected. The reflected angles of these reflected waves are generally in the range of 0 to 120˚. Therefore, noise can be removed by setting the cross-correlation value to "0" (zero) in the case of having a reflection angle deviating from a certain angle (for example, 0 to 120 °) in consideration of the propagation directions of two correlated wave fields. .

As described above, the cross-correlation is calculated and used only when the angle between the sound source wave and the reverse wave propagated observation wave (received wave from the receiver) is within a certain range, and image the cross-correlation value due to the wave having the reflection angle that is not necessary. By not using it as data, distorted noise can be removed. To this end, the wave propagation direction must be considered, and the wave propagation direction can be calculated by using a pointing vector defined as Equation (2) below.

The pointing vector is defined as a vector representing the flow of energy in electromagnetics, but in the case of a seismic wave, the pointing vector can be applied and can be defined as follows. In the present invention, the pointing vector is used to obtain the propagation direction of the wave. .

는

와 같이 계산할 수 있으며,

는 파동 진폭(음향 매질의 경우에는 압력) P(x,y,z,t) 의 공간에 대한 변화율(gradient)을 나타내고,

는 P(x,y,z,t) 의 시간에 대한 변화율을 나타내는 것이며, i, j, k 는 방향성분이다.From the stomach

Is

Can be calculated as

Represents the rate of change of the wave amplitude (pressure in the case of an acoustic medium) P (x, y, z, t) over space,

Represents the rate of change of P (x, y, z, t) with respect to time, and i, j, k are aromatic components.

이고, B 지점에서는

이다. 즉, 위 곱셈에 의하여 마이너스(-) 곱하기 플러스(+) 는 마이너스(-) 값을 가지게 되는 것이나, 수학식(2)의 전단부의 마이너스(-)로 인하여 전체적으로 플러스(+) 값을 가지게 되므로 플러스(+) 방향(+ X 방향)으로 파가 진행하고 있음을 파악할 수가 있게 된다.Attached Figure 5 is a view showing an example of the amplitude and direction of propagation of the wave of the wave in one dimension using the same medium, and when the wave is propagating from left to right (+ X direction), in point of the A

At point B

to be. That is, by the above multiplication, minus (-) multiplication plus (+) has a negative value (-), but due to the minus (-) of the front end of Equation (2) as a whole has a positive value (+) It is possible to grasp that the wave is traveling in the (+) direction (+ X direction).

Therefore, the wave propagation direction has the vector value obtained as the pointing vector in Equation (2), and the reflection angle θ can be obtained by obtaining the cosθ value using the vector product. That is, the vector A and vector B = enemy A and the vector B | A | * | B | cosθ Since, the number of the obtained value cosθ and will be determined in the angle of reflection θ to the arc cosine (arc cosine).

를 나타낸 영상사진이며, 도 7은 2차원 매질에서 0.5초, 1.0초 및 1.5초에서의 관측파 R(x,y,t) 와 이 관측파의 포인팅벡터의

를 나타낸 영상사진이다.FIG. 6 and FIG. 7 are image photographs of time-phase wave and wave propagation directions in a two-dimensional medium, and FIG. 6 is a sound wave S (x, y, t) at 0.5, 1.0, and 1.5 seconds in a two-dimensional medium . And a pointing vector representing the propagation direction of this sound wave.

7 shows the observed wave R (x, y, t) and the pointing vector of the observed wave at 0.5, 1.0 and 1.5 seconds in a two-dimensional medium.

It is an image photograph.

After obtaining the vector value of the sound source wave ( S ) and the reverse wave propagated observation wave ( R ), the cosθ value is obtained from the vector product S · R = | S | * | R | cosθ of the two waves. The reflection angle θ can be obtained as a cosine value, and weighted according to the reflection angle can be used to implement imaging according to the reflection angle. For example, when the angle θ between two waves is 0 ≦ θ ≦ 120, “1.0” is 120 < θ If the weight is "0.0", the reflection angle is 120 < θ In this case, since the value becomes "0", the cross-sectional view of the structure correction image when the reflection angle is θ ≤ 120 can be obtained.

8 and 9 are photographs comparing a conventional photograph of performing reverse temporal structure correction with a photograph of performing reverse temporal structure correction in consideration of the wave direction of the present invention. The inverse time structure correction of the present invention taking into account the propagation direction of the wave calculates the cross-correlation only when the reflection angle is θ ≤ 120 and images the cross-correlation value as “0” when the reflection angle is 120 < θ .

FIG. 8A is a photographic image of a conventional reverse time structure correction for one (upper photo) and 11 (lower photo) sound sources, and FIG. 8B is one (upper photo) and 11 (lower photo). The result of performing inverse time structure correction according to the present invention on the sound source data of 9 is a structural correction image cross-sectional view of the artificial acoustic wave exploration data, 9a is a photograph by a conventional method, 9b is a photograph according to the present invention.

As can be seen in the figures of FIGS. 8 and 9, in the case of the image photograph which has performed the conventional inverse time structure correction, it can be seen that noise exists in the upper part of the reflection surface. It can be seen that the reflection surface is cleared by removing the noise.

FIG. 10 shows a flow of imaging the seismic data by performing the inverse time structure correction of the present invention in consideration of the wave propagation direction.

The present invention collects sound source data that has undergone basic processing (Deconvolution, polymerization rate analysis, static correction, filtering, etc.) and determines a speed model (S100);

Calculating and storing a sound source wave S (x, y, z, t) (S200);

Calculating and storing a Poynting Vector of the sound source wave S (x, y, z, t) (S300);

Determining whether the time t is greater than or equal to Tmax and repeating steps S200 and S300 in step S350;

Calculating the backward propagated observation wave R (x, y, z, t) when it is determined in step S350 that time t is greater than or equal to Tmax (S400);

Calculating a Poynting Vector of the observation wave R (x, y, z, t) (S500);

Calculating cross-correlation S (x, y, z, t) * R (x, y, z, t) according to the reflection angle at all x, y, z points (S600);

Determining whether the time t is greater than or equal to Tmax and repeating steps S400 to S600 if the time t is small (S650); And

If it is determined in step S650 that time t is greater than or equal to Tmax , the result of adding S (x, y, z, t) * R (x, y, z, t) ( I (x, y, z) Generating a depth structure corrected image cross-sectional view as a step (S700).

In the above step S600, the reflection angle is obtained from the vector values of the sound source wave ( S ) and the observed wave ( R ), and then the cosθ value is obtained from the vector product S · R = | S | * | R | cosθ of the two waves. can be obtained the reflection angle θ as the arc cosine of a value, the value of θ is the 0≤ θ ≤120 it indicates give a "1.0" to the weight, 120 <θ If is a weight of "0.0" to calculate the cross-correlation,

In the above step S700, the value of the result value I (x, y, z) at a specific point x, y, z is represented as an integer value of plus (+) or minus (-). Imaging is done using Seismic Unix, Colorado School of Mine) package tool, SEP (Stanford Exploration Project), Geoframe or Promax. For example , the tone of a specific color is determined based on the value I (x, y, z) is "0", and the value is positive and the absolute value is increased so that it is darker, and the value I (x, y, z) z) If this is negative and the absolute value is larger, it can be expressed thinly or imaged using other common imaging methods.

On the other hand, the imaging method according to the present invention can be usefully used for AVO (Amplitude Versus Offset) analysis.

AVO analyzes the physical characteristics of the reflection medium by observing the amplitude change of the reflected wave with the distance between the sound source and the receiver.

In general, the structural correction of seismic survey data is based on the Kirchhoff integration method and the Wave Equation method.The structural correction by the wave equation is obtained by obtaining the numerical analysis of the wave equation. By doing so, it is known that the resolution and accuracy are superior to the Kirchhoff structure correction method. In the case of Kirchhoff structural correction involving ray tracing, it is easy to obtain a structural correction image within a certain reflection angle (for example, an angle of 0 to 10 degrees, 10 to 20 degrees, or 20 to 30 degrees). Depth-structure correction in the reflection angle region) However, there is no imaging method according to the reflection angle in the structure correction using the wave equation unless the dashed line tracking is performed separately.

However, even in the structure correction using the wave equation, if the reverse propagation correction is performed by calculating the propagation direction of the wave according to the present invention, the high quality image according to the reflection angle can be obtained even by the structure correction using the wave equation. More accurate AVO analysis is possible because the more accurate amplitude can be obtained than the method.

In the above, the present invention has been described above, but the present invention is not limited to the above-described technical configuration, but may be modified and modified within the scope not departing from the spirit of the technical idea of the present invention. The spirit should also be seen as belonging to the following claims. In particular, it should be noted that the technical idea of obtaining a wave traveling direction by applying a pointing vector is a unique configuration of the present invention.

 The image implemented by the method of imaging the acoustic wave data of the present invention removes noise due to cross-correlation with head waves, and the sound source wave or observation wave reflected from the boundary surface of the medium having a large speed change The noise generated due to cross-correlation with the observation wave or sound wave propagating toward it is eliminated, thereby eliminating distortion in the elastic wave image, which leads to a remarkable effect of increasing the quality of the image. AVO (Amplitude Versus Offset) analysis has the effect of more accurate AVO analysis than Kirchhoff method with dashed line tracking.

Claims (17)

In imaging seismic data through depth compensation before polymerization using inverse time structure correction, Calculate the cross-correlation value using Imaging Condition, which is a zero-lag cross correlation of the sound wave reflected from all the underground points and the observed wave propagated from the receiver (wave received from the receiver). In consideration of the propagation directions of the two correlated waves, the weight of the wave having the reflection angle corresponding to a certain angle is made as "1", and the weight of the wave corresponding to the remaining angle is made as "0" (zero), A method of imaging seismic data by inverse time structure correction using a pointing vector, characterized in that an image is implemented as a result of multiplying this weight by the above correlation value. The method of claim 1, The propagation direction of the wave is a pointing vector

Can be calculated using

Is

Can be calculated as

Represents the rate of change of the wave amplitude (pressure in the case of an acoustic medium) P (x, y, z, t) over space,

Represents a rate of change of P (x, y, z, t) with respect to time, and i, j, k are direction components, and the method of imaging acoustic wave data by inverse time structure correction using a pointing vector. The method according to claim 1 or 2, The reflection angle of the wave making the weight "1" in the propagation direction is 0 to 120 °, and for the other angles, the weight is "0" (zero) for inverse time structure correction using a pointing vector. Method of imaging seismic data by The method according to claim 1 or 2, The reflection angle is obtained by pointing vector values of the sound wave and the observation wave, and then the reflection angle is obtained as a vector product of two waves. The method of claim 3, wherein The reflection angle is obtained by pointing vector values of the sound wave and the observation wave, and then the reflection angle is obtained as a vector product of two waves. In order to image seismic data through depth structure correction before inversion using structure correction, sound source collection data that have undergone basic processing (deconvolution, polymerization rate analysis, static correction, filtering, etc.) Collecting and determining a speed model (S10); Calculating and storing a sound source wave S (x, y, z, t) at a predetermined time interval Δt (S20); Determining whether time t is greater than or equal to Tmax and repeating step S20 if small (S30); Calculating the reverse propagated R (x, y, z, t) when it is determined in step S30 that time t is greater than or equal to Tmax (S40); Calculating cross-correlation S (x, y, z, t) * R (x, y, z, t) at all x, y, z points (S50); Determining whether the time t is greater than or equal to Tmax and repeating steps S40 and S50 if small (S60); And S (x, y, z, t) * R (x, y, z, t) when the time t is determined to be equal to or greater than Tmax in the step S60 , resulting in the value ( I (x, y, z) generating a depth structure corrected image cross-sectional view as a step (S70); Calculate the propagation direction of the wave with respect to the sound wave S (x, y, z, t) and the observed wave R (x, y, z, t ) by using the pointing vector. ( S and R ) The cross-correlation value is obtained using zero-lag cross-correlation imaging conditions, but considering the propagation direction of two correlated waves, the weight of a wave having a reflection angle corresponding to a certain angle Is set to "1", and the weights of the waves corresponding to the remaining angles are set to "0" (zero), and the pointing vector using the pointing vector is characterized by realizing the image as a result of multiplying the weight by the above correlation value. Imaging Method of Seismic Data by Inverse Time Structure Correction. The method of claim 6, The propagation direction of the wave is a pointing vector

Can be calculated using

Is

Can be calculated as

Represents the rate of change of the wave amplitude (pressure in the case of an acoustic medium) P (x, y, z, t) over space,

Represents a rate of change of P (x, y, z, t) with respect to time, and i, j, k are direction components, and the method of imaging acoustic wave data by inverse time structure correction using a pointing vector. The method according to claim 6 or 7, The reflection angle of the wave making the weight "1" in the propagation direction is 0 to 120 °, and for the other angles, the weight is "0" (zero) for inverse time structure correction using a pointing vector. Method of imaging seismic data by The method according to claim 6 or 7, The reflection angle is obtained by obtaining the pointing vector values of the sound source wave S and the observation wave R , and then calculating the reflection angle θ as a vector product of the two waves. Imaging Method of Material. The method of claim 8, The reflection angle is obtained by obtaining the pointing vector values of the sound source wave S and the observation wave R , and then calculating the reflection angle θ as a vector product of the two waves. Imaging Method of Material. In imaging seismic data through depth compensation before polymerization using inverse time structure correction, Collecting sound source data that has undergone basic processing (deconvolution, polymerization rate analysis, static correction, filtering, etc.) and determining a speed model (S100); Calculating and storing a sound source wave S (x, y, z, t) (S200); Calculating and storing a Poynting Vector of the sound source wave S (x, y, z, t) (S300); Determining whether the time t is greater than or equal to Tmax and repeating steps S200 and S300 in step S350; Calculating the backward propagated observation wave R (x, y, z, t) when it is determined in step S350 that time t is greater than or equal to Tmax (S400); Calculating a Poynting Vector of the observation wave R (x, y, z, t) (S500); Calculating cross-correlation S (x, y, z, t) * R (x, y, z, t) according to the reflection angle at all x, y, z points (S600); Determining whether the time t is greater than or equal to Tmax and repeating steps S400 to S600 if the time t is small (S650); And If it is determined in step S650 that time t is greater than or equal to Tmax , the result of adding S (x, y, z, t) * R (x, y, z, t) ( I (x, y, z) Generating a depth structure corrected image cross-sectional view (S700); and a method for imaging seismic data by inverse time structure correction using a pointing vector. The method of claim 11, In the steps S300 and S500 the pointing vector is

Can be calculated using

Is

Can be calculated as

Represents the rate of change of the wave amplitude (pressure in the case of an acoustic medium) P (x, y, z, t) over space,

Represents a rate of change of P (x, y, z, t) with respect to time, and i, j, k are direction components, and the method of imaging acoustic wave data by inverse time structure correction using a pointing vector. The method of claim 11, Calculate the propagation direction of the wave with respect to the observation wave R (x, y, z, t) back propagated with the sound source wave S (x, y, z, t) using the pointing vector. ( S (x, y, z, t) and R (x, y, z, t) ) The cross-correlation value is obtained using zero-lag cross-correlation imaging conditions, but considering the propagation direction of two correlated waves, the weight of a wave having a reflection angle corresponding to a certain angle Is set to "1", and the weights of the waves corresponding to the remaining angles are set to "0" (zero), and the pointing vector using the pointing vector is characterized by realizing the image as a result of multiplying the weight by the above correlation value. Imaging Method of Seismic Data by Inverse Time Structure Correction. The method of claim 12, Calculate the propagation direction of the wave with respect to the observation wave R (x, y, z, t) back propagated with the sound source wave S (x, y, z, t) using the pointing vector. ( S (x, y, z, t) and R (x, y, z, t) ) The cross-correlation value is obtained using zero-lag cross-correlation imaging conditions, but considering the propagation direction of two correlated waves, the weight of a wave having a reflection angle corresponding to a certain angle Is a value of "1", and the weights of the waves corresponding to the remaining angles are "0" (zero). Method of imaging seismic data by temporal structure correction. The method according to claim 13 or 14, The reflection angle of the wave making the weight "1" in the propagation direction is 0 to 120 °, and for the other angles, the weight is "0" (zero) for inverse time structure correction using a pointing vector. Method of imaging seismic data by The method according to any one of claims 11 to 14, The reflection angle is obtained by obtaining the pointing vector values of the sound source wave S and the observation wave R , and then calculating the reflection angle θ as a vector product of the two waves. Imaging Method of Material. The method of claim 15, The reflection angle is obtained by obtaining the pointing vector values of the sound source wave S and the observation wave R , and then calculating the reflection angle θ as a vector product of the two waves. Imaging Method of Material.

Images