WO2012064839A2 - System and method for investigating sub-surface features of a rock formation - Google Patents

Abstract

A method and system for investigating rock formations outside a borehole are provided. The method includes generating a first acoustic wave at a first frequency by a first acoustic source; and generating a second acoustic wave at a second frequency by a second acoustic source. The first and the second acoustic sources are arranged within a localized area of the borehole. The first and the second acoustic waves intersect in an intersection volume outside the borehole. The method further includes receiving a third acoustic wave at a third frequency, the third shear acoustic wave returning to the borehole due to a non-linear mixing process in a non-linear mixing zone within the intersection volume at a receiver arranged in the borehole. The third frequency is equal to a difference between the first frequency and the second frequency.

Description

SYSTEM AND METHOD FOR INVESTIGATING

SUB-SURFACE FEATURES OF A ROCK FORMATION

CROSS REFERENCE TO RELATED APPLICATION

This application claims the benefit of U.S. patent application Serial No. 61/413,172, filed on November 12, 2010, the entire contents of which is incorporated herein by reference.

GOVERNMENT RIGHTS

This invention was made with Government support under Cooperative Research and Development Agreement (CRADA) Contract Number DE-AC52- 06NA25396 awarded by the United States Department of Energy. The Government may have certain rights in this invention.

FIELD

The present invention relates generally to seismic interrogation of rock formations and more particularly to creating three-dimensional images of non-linear properties and/or the compressional to shear velocity ratio in a region remote from a borehole using a combination of sources in a borehole, and receiving and analyzing a resultant third wave formed by a mixing process.

BACKGROUND

Acoustic interrogation of subsurface features tends to be limited by the size and power of practical sources, and in practice, the output of down hole acoustic transducers is limited by the power transmission capabilities of the wireline cable. High frequency signals have a relatively short penetration distance, while low frequency signals generally require large sources, clamped to the borehole wall, to maximize energy transfer to the formation and minimize unwanted signals within the well bore. Currently, acoustic borehole tools are designed with acoustic sources in the borehole to detect returning acoustic waves that are propagating along the borehole walls or scattered by inhomogeneities of linear properties of rock formations surrounding the borehole. U.S. Patent No. 7,301,852 to Leggett, III et al. discloses a Logging While Drilling tool, designed to detect rock formation boundaries. The tool uses two acoustic source arrays emitting two acoustic waves from a borehole that generate a third wave by non-linear mixing in the rock formation surrounding the borehole at the location of intersection of the acoustic waves. The third wave continues forward and interacts linearly with heterogeneities in the subsurface properties. The third wave is scattered by the heterogeneities in the subsurface properties, and the scattered signal is detected by sensors in the logging tool. U.S. Patent No. 7,301,852 does not discuss detecting the third wave directly but rather the signal that is scattered by the heterogeneities in the rock formation. U.S. Patent No.

7,301,852 merely uses the resultant scattered wave to detect rock formation boundaries.

Attempts have been made to characterize the non-linear properties of a formation in the area of oil and gas prospecting from boreholes, but each has its own limitations. For example, U.S. Patent No. 5,521,882 to D'Angelo et al. discloses an acoustic tool designed to record with pressure receivers non-linear waves generated by non-linear mixing of two waves. The non-linear waves propagate along the borehole wall with limited penetration into the surrounding rock formation and refract back into the well bore fluid. The indication of non-linearity is utilized to provide an indication of the relative consolidation of the formation surrounding the borehole.

U.S. Patent No. 5,521,882 does not discuss measuring non-linear characteristics of a rock formation away from the borehole. U.S. Patent No. 6, 175,536 by Khan discloses a method to estimate the degree of non-linearity of earth formations from spectral analysis of seismic signals transmitted into the earth formations from a first borehole and received in a second borehole. The method in U.S. Patent No. 6,175,536 determines from the spectral analysis the presence of a frequency at a receiver located at the second borehole representing a sum or a difference of two selected frequencies of the transmitted seismic signals generated by two sources located at the first borehole. U.S. Patent No. 6, 175,536 does not discuss measuring non-linear characteristics of a rock formation in a remote region of a borehole where the receiver and the sources are located in one borehole. In light of these prior attempts, there is a need for a system and method for characterizing non-linear properties in a remote region from a borehole.

SUMMARY

An aspect of the present disclosure is to provide a method for investigating rock formations outside a borehole. The method includes generating a first acoustic wave at a first frequency by a first acoustic source; and generating a second acoustic wave at a second frequency by a second acoustic source. The first and the second acoustic sources are arranged within a localized area of the borehole. The first and the second acoustic waves intersect in an intersection volume outside the borehole.

The method further includes receiving a third acoustic wave at a third frequency, the third shear acoustic wave returning to the borehole due to a non-linear mixing process in a non-linear mixing zone within the intersection volume at a receiver arranged in the borehole. The third frequency is equal to a difference between the first frequency and the second frequency.

Another aspect of the present disclosure is to provide a system for investigating rock formations outside a borehole. The system includes a first acoustic source configured to generate a first acoustic wave at a first frequency; and a second acoustic source configured to generate a second acoustic wave at a second frequency. The first and the second acoustic sources are arranged within a localized area of the borehole. The first and the second acoustic waves intersect in an intersection volume outside the borehole. The system further includes a receiver arranged in the borehole and configured to receive a third acoustic wave at a third frequency, the third shear acoustic wave returning to the borehole due to a non-linear mixing process in a non- linear mixing zone within the intersection volume. The third frequency is equal to a difference between the first frequency and the second frequency.

In some aspects of the present disclosure, the first acoustic wave and the second acoustic wave include coded acoustic signals. In other aspects of the present disclosure, the first and second acoustic sources are configured to generate conical broadcast acoustic signals. In yet further aspects of the present disclosure, there is provided a system and method for imaging non-linear properties of the rock formation, or determining a compressional velocity, determining a shear velocity, or determining a velocity ratio between the compressional velocity and the shear velocity in the rock formation, or any combination thereof.

In yet another aspect of the present disclosure, there is provided a system for investigating non-linear properties of a rock formation around a borehole. The system includes a first sub-system configured to perform data acquisition, control and recording of data; a second subsystem in communication with the first sub-system and configured to perform non-linearity and velocity preliminary imaging; a third subsystem in communication with the first subsystem and configured to emit controlled acoustic broadcasts and receive acoustic energy; a fourth subsystem in communication with the first subsystem and the third subsystem and configured to generate a source signal directed towards the rock formation; and a fifth subsystem in communication with the third subsystem and the fourth subsystem and configured to perform detection of signals representative of the non-linear properties of the rock formation.

These and other objects, features, and characteristics of the present invention, as well as the methods of operation and functions of the related elements of structure and the combination of parts and economies of manufacture, will become more apparent upon consideration of the following description and the appended claims with reference to the accompanying drawings, all of which form a part of this specification, wherein like reference numerals designate corresponding parts in the various Figures. It is to be expressly understood, however, that the drawings are for the purpose of illustration and description only and are not intended as a definition of the limits of the invention. As used in the specification and in the claims, the singular form of "a", "an", and "the" include plural referents unless the context clearly dictates otherwise.

BRIEF DESCRIPTION OF THE DRAWINGS

Figure 1 shows a configuration for creating three-dimensional images of nonlinear properties in a region remote from a borehole, in accordance with various aspects of the disclosure;

Figure 2 shows another configuration for creating three-dimensional images of non-linear properties in a region remote from a borehole, in accordance with aspects of the disclosure; Figure 3 shows yet another configuration for creating three-dimensional images of non-linear properties in a region remote from a borehole, in accordance with aspects of the disclosure;

Figure 4 shows a flow chart for creating three-dimensional images of non- linear properties in a region remote from a borehole, in accordance with various aspects of the disclosure;

Figures 5a, 5b and 5c shows a numerical simulation of the first selection rule for a beam-beam interaction listed in Table 1 when the two primary waves are beams;

Figure 6 illustrates the geometry of the generation of the difference frequency third wave by non-linear mixing of two primary acoustic waves as governed by the non-linear mixing selection rule;

Figure 7 shows an application of aspects of the present disclosure for imaging using a beam and broad beam or plane wave;

Figure 8 shows an example configuration for a borehole-based system for remote mapping of non-linear properties and/or Vp/Vs ratio of rock formations using non-collinear acoustic mixing, in accordance with an aspect of the present disclosure;

Figure 9a shows the configuration of figure 8 for the purposes of identifying vectors representing the broadcast and scattered acoustic waves;

Figure 9b shows a vector representation of the non-collinear acoustic mixing of Figure 9a;

Figure 10a shows a representative dependence of mixing coefficient on the plane wave frequency ratio for a range of mixing zone compressional velocity and shear velocity ratios Vp/Vs, in accordance with aspects of the present disclosure;

Figure 10b shows a representative dependence of convergence angles with plane wave frequency ratio that honor the selection rules for the P+P SV interaction;

Figure 10c shows a representative dependence of scattering angles with plane wave frequency ratio that honor the selection rules for the P+P SV interaction;

Figures 11a to 11c show example results of a numerical simulation of non- collinear interaction of plane waves in a non-linear medium leading to the generation of a scattered wave that returns to the borehole, in accordance with aspects of the present disclosure;

Figure 12a shows an example representation of the directions and the times of flight for the primary and scattered acoustic waves in accordance with aspects of the present disclosure;

Figures 12b-12d show an example of a simulated signal from a first acoustic source, a simulated signal from second acoustic source, and a simulated template signal, according to an embodiment of the present disclosure;

Figure 13a shows a position of a first acoustic source and a second acoustic source and a receiver array, according to an embodiment of the present disclosure;

Figures 13b to 13d show how correlation of the received series of pulses with a modeled template signal results in identification of the signal's arrival time at the receiver array, in accordance with an embodiment of the present disclosure;

Figure 14a shows a position of a first acoustic source and a second acoustic source and a receiver array, according to an embodiment of the present disclosure;

Figures 14b to 14d show the effect of transmitting a series of coded pulses and using correlation techniques to improve signal to noise ratio, in accordance with aspects of the invention;

Figures 15a and 15b show an example non-collinear mixing arrangement in a toroid around the borehole at intersection of two coaxial cones, in accordance with aspects of the present disclosure;

Figures 16a and 16b show an example non-collinear mixing arrangement between two intersecting coaxial cones, in accordance with aspects of the present disclosure;

Figure 17a shows an example single well arrangement with a cranked rathole where there is a complete intersection of lower cone with upper cone, in accordance with aspects of the present disclosure;

Figure 17b shows another example single well arrangement with a cranked rathole where a lower transmitter emits energy near perpendicular to a borehole axis, in accordance with aspects of the present disclosure; Figure 18a shows an example of a vertical well and sidetrack with receivers in the vertical part of the well, in accordance with aspects of the present disclosure;

Figure 18b shows another example of a vertical pilot hole and horizontal sidetrack with receivers in the sidetrack, in accordance with aspects of the present disclosure; and

Figure 19 is a schematic system diagram of a system for performing survey design, data acquisition, data processing and imaging, in accordance with aspects of the present disclosure.

DETAILED DESCRIPTION

Figure 1 shows one of several possible configurations for creating three- dimensional images of non-linear properties and the compressional to shear velocity ratio in a region remote from a borehole in accordance with various aspects of the disclosure. First acoustic source 105 is arranged in borehole 1 10 to generate a steerable primary beam of acoustic energy at a first frequency fi. Second acoustic source 1 15 is also arranged in borehole 110 to generate a steerable primary beam of acoustic energy at a second frequency By way of a non-limiting example, both first acoustic source 105 and second acoustic source 1 15 may be a phased array of sources and may be configured to generate either compressional or shear steerable beams. In the present disclosure, the term "acoustic" can refer to P, SV or SH acoustic mode.

As shown in Figure 1, first acoustic source 105 is arranged on first tool body 120 and second acoustic source 115 is arranged on second tool body 125. However, the disclosure is not so limiting as first tool body 120 and second tool body 125 may also be arranged together on a common tool body (not shown). Tool bodies 120 and

125 are arranged to be independently moveable within bore hole 110 in at least two degrees of freedom including translation along the longitudinal axis 150 of borehole 1 10 and rotation 155 in azimuth about the longitudinal axis of borehole 1 10. First acoustic source 105 may be arranged above or below second acoustic source 115 in borehole 1 10. Tool bodies 120 and 125 may be arranged on a conveyed logging tool

(not shown) within borehole 110. For a given azimuth orientation of first acoustic source 105 and second acoustic source 115, the beam generated by second acoustic source 1 15 and the beam generated by first acoustic source 105 are configured such that the beams converge and intercept in mixing zones 130 remote from borehole 1 10. By a combination of independently steering the beams and changing the separation between the sources 105, 1 15, the mixing zones 130 move in a plane defined by the beams and the longitudinal borehole axis 150, while controlling the angle of interception. The distance of mixing zones 130 from borehole 1 10 can range from near the edge of borehole 1 10 to about 300 meters into the surrounding subsurface rock formation. By way of a non-limiting example, the phase difference and/or start time differences between adjacent elements in the source array 105, 115 referred to in the above paragraphs may be modified to focus the acoustic energy of the primary beams at a particular mixing zone 130.

The non-linear properties of the earth at the location between the two waves result in the generation of a third elastic wave. The third elastic wave is a result of a three-wave mixing process that occurs in nonlinear materials, in this case, rock formations. In this process, two converging non-collinear waves of different frequencies, /; and , also called primary waves, mix to form additional waves at the harmonic and intermodulation frequencies fi - f2, fi + ,

and 2xf₂, etc. The strength of the third wave is a function the non-linearity of the rocks in the mixing zones. By way of a non-limiting example, when a primary compressional (P) wave with a frequency fi and a primary shear (SV) wave with a frequency cross or intersect in a non-linear medium, a third compressional (P) or shear (SV) wave is generated with a frequency fi -f2.

Under propagation selection rules, the third wave propagation vector is co- planar with the propagation vectors of the two primary waves. Certain combinations of angle of intersection, /;// ratio and compressional to shear velocity ratio result in a third elastic wave with frequency / -/ propagating in a specific angle relative to the primary beams back to the borehole 1 10.

Sensor or receiver array 135 is arranged at a specific location in borehole 1 10 to detect the third wave returning to the borehole 110. In one embodiment, as shown for example in the Figure 1, sensor array 135 comprises more than one sensor, arranged as an array of sensors on sensor tool body 140 and separate from tool bodies 120 and 125. Sensor 135 is configured to be independently moveable within bore hole 110 along the longitudinal axis 150 of borehole 1 10. For example, sensor tool body 140 can be arranged below tool bodies 120 and 125 or arranged above and below tool bodies 120 and 125. In some embodiments, sensor tool body 140 can be connected to either or both tool bodies 120 and 125.

The third wave is detected at borehole 1 10 by sensor array 135. Figure 2 shows an arrangement similar to Figure 1, wherein receiver 135 includes three component geophone 145 clamped to the borehole walls. The resultant signal is decomposed by processing into its elevation and azimuth in order to add redundancy to the system by determining the direction of the incoming third wave arrival.

In one embodiment, a first processor or controller can be provided and configured to execute machine-readable instructions (not shown) to perform various processing tasks, such as controlling source firing and compressing or filtering the data recorded by sensor array 135. In one embodiment, the first processor can be arranged within the borehole 110. In one embodiment, a second processor can be provided and configured to execute machine-readable instructions (not shown) to assist the first processor or perform different processing tasks than the first processor. For example, the second processor may perform part or all processing activities in creating the three-dimensional images. A transmitter or transceiver (not shown) may be arranged in borehole 1 10 to transmit data up-hole through a wireline cable (not shown). In one embodiment, the second processor can be, for example, arranged outside the borehole.

At a given depth along the borehole of one of the sources 105, 115, sweeping the beams in elevation at constant relative bearing to spatially scan the mixing zone in a plane passing through the borehole axis, rotating the sources azimuthally to rotationally scan the mixing region and moving the whole assembly along borehole 1 10, results in scanning a 3D volume of mixing zones around the borehole for nonlinear properties. With sources 105, 1 15 and sensor array 135 located on independent tool bodies, redundancy in the data can be obtained and the depth of investigation can be varied. In this way, a 3D volume of the rocks surrounding the borehole can be interrogated for non-linear properties and a 3D image of non-linear properties can be processed and computed from the returned signals, i.e., signals detected by sensor array 135.

Figure 3 shows another arrangement for creating three-dimensional images of non-linear properties in a region remote from borehole 110 in accordance with another embodiment of the present invention. The arrangement of Figure 3 is similar to the arrangement in Figure 2, with the primary difference being that the sources are arranged in borehole 1 10 to produce elastic waves (e.g., spherical waves) instead of steerable directional beams. With reference to Figure 3, first acoustic source 305 is arranged in borehole 1 10 on first tool body 320 to generate a first elastic wave of acoustic energy at a first frequency fi. Second acoustic source 315 is arranged in borehole 1 10 on second tool body 325 to generate a second elastic wave of acoustic energy at a second frequency First and second elastic waves produced by sources 305, 315 are arranged to intersect away from borehole 110 at various mixing zones 130. Receiver 145 is arranged within borehole 110 to receive a third wave that is produced in the mixing zones 130 by the three-wave mixing process discussed above, and further discussed below. Since the waves produced by sources 305, 315 are essentially non-directional, mixing between the waves occurs simultaneously in the entire area of mixing zones 130, that also extends out of the plane of the Figure, and receiver 145 tends to have directional characteristics. By way of a non-limiting example, a three component geophone array may be used for this purpose. The resultant signal is decomposed by processing into multiple arrival signals at a range of elevations and azimuths and travel times. Given the locations of sources 305 and 315 and the receivers 145, the travel times and directions of each decomposed directional arrival, there is sufficient information to apply selection rules described in the following paragraphs to determine a unique mixing zone where the third wave was generated. This unique mapping allows the construction of a three dimensional (3D) image from the properties of the received signal.

Figure 4 shows a flow chart for a method of creating three-dimensional images of non-linear properties and the compressional to shear velocity ratio in a region remote from a borehole using a conveyed logging tool, according to an embodiment of the present invention. The method begins at 405 where a first acoustic source is arranged in the borehole to generate a steerable beam elastic energy at a first frequency and a second acoustic source is arranged in the borehole to generate a steerable beam of elastic energy at a second frequency. The steerable beams at the first and second frequency are arranged to intersect at a location away from the borehole. As such, the second beam is generated at the same azimuth as the first beam, but at a different elevation relative to the longitudinal axis of the borehole.

The method continues at 410 where a third elastic wave is received at the borehole by a sensor array. As discussed above, the third elastic wave is created by a mixing process, with a frequency equal to a difference between the first and second frequencies and a direction of propagation towards the borehole. At 415, a mixing location away from the borehole is determined from the arrangement of the first and second acoustic sources and properties of the third wave, by recourse to the selection rules. At 420, three-dimensional images are created of the non-linear properties using data recorded by repeating the generating of step 405, the receiving of step 410 and the determining of step 415 at a plurality of azimuths, elevations and longitudinal locations within the borehole. In cases of compressional-shear interaction, the received signals are analyzed in step 425 for the compressional/shear velocity (Vp/Vs) ratio as discussed in the above paragraphs. At 430, the non-linear properties are transformed to physical reservoir properties such as fluid saturation, effective stress, fracture density and mineralogy.

In some aspects of the present disclosure, the first and second acoustic sources may be beam or cylindrical, or spherical wave sources, and the sensor array may be any combination of non-directional single component sensors and three component geophones. Alternative permutations of the component parts offer different degrees of redundancy in signal processing and imaging.

In the special case where a primary compressional (P) wave with a frequency i and a primary shear (S) wave with a frequency cross each other, in a non-linear medium, a third P or S wave is generated with the frequency fi - f₂. If the primary P and S waves are beams with wave vectors ki and k₂, respectively, and the non-linear formation property is uniform, the kinematics of wave interaction requires the resulting third wave to be a plane wave with wave vector k3 that obeys the selection rule ki - k₂ = k3. The selection rule imposes a very tight restriction on the permissible crossing angles for the primary waves and a specific propagation direction of the third wave. The general kinematic theory for non-linear mixing of two linear plane waves and the selection rules and amplitude responses have contributions from Jones and Kobett (1963), Rollins, Taylor et al. (1964) and later by Korneev, Nihei and Myer (1998), all of which are hereby incorporated by reference in their entirety, who also provide specific relationships between non-linear parameters of the mixing medium and the non-linear mixing signal strength. For example, Equation 53 and 54 of Korneev, Nihei and Myer show that the mixing strength of P and SV (vertically polarized shear) plane waves is proportional to a specific combination of non-linear parameters of the rocks.

The selection rules governing the nonlinear interaction of two elastic plane waves can be used as guidance for the interaction of two elastic beams. These plane wave selection rules dictate that the following six nonlinear interactions produce backscattered waves.

Table 1 - Selection Rules Governing Non-Linear Interaction of Two Elastic Plane Waves. In this table, and elsewhere in this document,/; is greater than f₂.

Figures 5a, 5b and 5c shows a numerical simulation of selection rule 1 of Table 1 when the two primary waves are beams of a beam-beam interaction. A 25 kHz compressional beam, shown in Figure 5a, and a 18 kHz shear beam, shown in Figure 5b, mix to form a third beam, shown in Figure 5c, with frequency 7 kHz = 25 kHz - 18 kHz. In this example, in accordance with the plane wave predictions, a third back propagating P beam with frequency (fi -fi) at an angle of 133° to the P( ) wave is generated by nonlinear mixing in the region where the P( ) and SV( ¾ beams overlap.

The kinematics of non-linear interactions of beams results in the generation of specific combinations of wave vectors and frequencies. The third wave returns at a specific travel time, and with specific frequencies j and wave vectors k₃ such

- f₂ and k₃=ki - k₂. For a combination of fi, f₂, ki and k₂, there is a well-defined propagation wave vector k₃ of the third wave in the same plane, defined by ki and k₂. There is a direct correspondence between the signal detected at a particular receiver position and the location where the non-linear mixing of the two primary waves ki and k₂ takes place. The signal strength of the receiver would be proportional to the strength of the non-linearity of the rocks in the mixing zone, among other factors, and reach a maximum for a receiver lying on vector k₃. Therefore, the signal strength at the receivers can be geometrically mapped onto the non-linearity of the rocks along the beam trajectory as shown, for example, in Figure 1.

The geometrical theory of wave propagation indicates that the beam generated in each mixing zone would arrive at the borehole at a specific receiver defined by the geometry of the three wave vectors ki, k₂ and k₃, after a specific time delay. The strength of the returning signal at a specific location in the borehole at a particular time is dependent on the degree of non-linearity of the interaction location.. Hence, a time image of the relative strength of the non-linear properties of the rocks along the beam can be constructed. The amplitude or magnitude of a returned signal at the receivers can be itself indicative of certain petrophysical properties of the mixing zone. If the beam and plane wave are scanned in azimuth and elevation while preserving the convergence angle, a localized circumferential and radial 3D image of non-linear properties of rocks surrounding the borehole can be obtained. By moving the entire assembly up and down the borehole, repeated 3D images of non-linear properties of rocks surrounding the borehole can be obtained. By making weighted stacks of these repeated images, a final image of non-linear properties of rocks surrounding the entire borehole can be constructed through subsequent computer processing. In addition, if the sources and the receivers are part of three separate tool bodies, one or two can be moved while the third one is fixed (for example, the sources are fixed while the receiver tool body is moved up and down). Alternatively, several descents into the well may be made with different spacing between the tool bodies.

For non-linear mixing between an elastic beam and a broader beam (quasi plane wave), the selection rule is relaxed. Third waves of frequency fi -f2, centered around the wave vector k₃ = ki - k2, are generated continuously along the primary beam if the beam width is about ten wavelengths of the third wave. The resulting signal strength for fii =fi ^~f2 is a function of the average non-linear properties of the mixing region, the average ratio of velocity of fi propagation and average velocity for f2 propagation (noting that beams with frequencies fi and may be compressional or shear), the volume of the mixing zone and the geometry of the mixing. This function can be computed for various mixing modes. For example, the signal strength for a particular mixing mode such as compressional wave P for fi and SV for f₂ is given by

where U is the displacement amplitude of the third wave received at the borehole, A i is the longitudinal polarization of the compressional wave and B2 is the transverse polarization of the shear wave, β is a function of the A, B and C parameters of Landau and Lifschitz representing the non-linearity of the rocks in the mixing zone, v is the volume of the mixing zone, L is the distance from mixing zone to the receiver. F is the geometric form factor of order 1 which is dependent on the geometry of the incident beams and can be numerically computed for the particular geometry. Δ is a selection rule form factor which is a numerically computable function of the wave vectors ki, k2 and k₃ and is only significant if the interaction geometry honors the selection rules. The subscript PSyP in the formula refers to compressional-shear interaction generating a compressional wave.

In accordance with certain aspects of this disclosure, an image of the compressional to shear velocity ratio may be constructed as follows. When one of the sources generates a compressional wave (P-wave) with frequency fi and the other source generates an SV-wave with frequency and both waves are steered towards a specific intersection volume, the propagation direction of the third compressional wave (P-wave) with difference frequency =fi -f2 is controlled by the average in situ Vp/V s velocity ratio of the rock in the mixing zone as governed by the selection rules as shown in Figure 6. From the measurements of the signal in the three component receiver array 145 on Figure 2 or Figure 3, the direction of this third wave can be determined and thereby, the in situ Vp/V s of the mixing zone can be computed. If the beam and plane wave are scanned in azimuth and elevation while preserving the necessary convergence angle, a localized circumferential and radial 3D image of in situ Vp/Vs ratio of rocks surrounding the borehole can be obtained. By moving the entire assembly up and down the borehole, repeated 3D images of in situ Vp/Vs of rocks surrounding the borehole may be obtained. By making weighted stack of these repeated images, a final image of in situ Vp/Vs of rocks surrounding the entire borehole can be constructed through subsequent computer processing. Alternatively, several descents into the well may be made with different fixed spacing between the tool bodies.

In some aspects of this disclosure, an alternative determination of Vp/Vs ratio is achieved through scanning the ratio of the frequencies fi to of the primary beams. Figure 6 illustrates the geometry of the interaction of two beams such as those generated in the configuration of Figure 1, that may be analyzed using the vector mathematics and trigonometry. The lengths kj and of vectors ki and k2 are defined by the ratio of their corresponding frequencies and velocities. As shown in Figure 6, the returning angle φ is a function of fi/f2, Vp/Vs ratio and the intersection angle Θ of the two primary beams. In addition, the physical selection rules only permit the generation of a third wave at specific combinations of _/}/ , Vp/Vs ratio and angle of interception Θ, such as the example illustrated on Figure 5.

Using the symbol r for the Vp/Vs ratio and the terms defined on Figure 6, the magnitude of vector k₃ is given by the vector sum of ki and -k₂, that is fc$■ |i ^— ί½ | " ^ί½^{£ϊ an}d also by the cosine rule that states k = kl + k - 2k_xk₂ cos Θ . Combining the two equations, and substituting filVp for ki and fJVs for leads to a statement of the geometric conditions imposed by the selection rules. The quadratic equation r — 2r cos # - + 2 = 0 may be solved

for r, the Vp/V s ratio of the mixing zone. This leads to a non-limiting alternative method for measuring in situ Vp/Vs ratio of a particular mixing region by the following sequence: a) record a standard sonic waveform log to determine Vp and Vs near the borehole to acquire data to estimate the phase differences between adjacent elements in a phased source array to steer the beams at the approximate convergence angle for the geometry of the planned measurement; b) steer the P and SV sources to converge at a controlled angle Θ and mix at a particular region in space surrounding the borehole; c) vary f₂ while fixing fj and measure the amplitude of the received signal at the difference frequency fi -f₂ at the sensors in the borehole; d) identify the frequency at which the signal each receiver in the array reaches a maximum amplitude strength; and e) determine angles Θ and φ from the geometry of the sources and receivers. By sweeping the beams in elevation, rotating in azimuth, and moving the entire assembly up and down the borehole and repeating the above procedure, the Vp/Vs ratio of a 3D volume around the borehole is interrogated and thereby 3D images of in situ Vp/V s ratio of rocks surrounding the borehole may be obtained.

The methods described above provide that the frequency difference fi -f₂ is very specific, allowing for spectral analysis to enhance the signal to noise ratio of the measurements. Moreover, if both frequencies fi and f₂ are simultaneously chirped proportionally, the resulting difference frequency signal fi - f₂ would also be a well defined chirped signal. The time-varying code may include one or more of a variation in amplitude, a variation in frequency, and/or a variation in phase of the first, the second, or both the first and the second beams or waves. The third difference wave can be broadband if one of the primary frequencies is swept through a range of frequencies while their frequency ratio is fixed. Thus, the resulting third beam f₂ -fi will be swept across a wide frequency range, while preserving the same direction. This allows for improvement in signal to noise by standard auto-correlation of the chirped or coded signal.

Since the wave vector k₃ = ki - k₂ is well defined, the signal to noise discrimination of the recorded third wave from receivers 135 can be enhanced further by employing three-component receivers in the borehole. For example, the signals from the three components can be tuned to specific directivity by a technique, such as, hodogram analysis.

In some aspects of the present disclosure, the signal to noise ratio can be improved by repeating the above steps and using an inverse polarity (180 degrees out of phase) source signals and adding the results together. The returning difference frequency signal will add coherently as its amplitude is proportional to the product of the amplitudes of the two primary waves and therefore will not reverse polarity when the polarity of the primary source is reversed. On the other hand, any linear noise generated by the primary sources in the system will reverse polarity and thus cancel upon addition.

Alternative methods can be devised with various non-exclusive combinations of beams and waves. By way of a non-limiting example, a method to generate images by computer processing of acoustic and seismic signals includes the following steps. First, the method performs spectral analysis of the frequency content of the recorded third wave and applicable selection rules of the difference frequency signal in order to isolate the third wave signal generated by the non-linear mixing process. In the case that the sensors include three component geophones, determine the direction of the third wave impinging on the borehole using orientation techniques. The method continues by analyzing the amplitude of the recorded third wave as a function of frequency ratios of the primary mixing waves and determining the mixing location where the third wave signals originated, from the selection rules of non-collinear mixing in non-linear media, the wavenumbers of the first and second beams and the third wave and the locations of the two beam sources and the sensor array. The method continues by constructing seismograms determined by cross-correlation of the received signals with chirped transmitter signals for each source-receiver combination. The method continues by performing three dimensional time or depth imaging to the entire data set, in order to obtain three dimensional images of the nonlinear properties of the formation surrounding a borehole in either or both of time and distance. The methods for generating images from seismograms are known, for example, Hill et al., which is hereby incorporated by reference, have provided the general methodology for the special case of imaging from beams.

Another non-limiting alternative imaging method is illustrated in Figure 7, which shows the case of interactions of a narrow 705 and a broad (wide) beam 710. Given a smooth background model of Vp and Vs of the investigated volume, application of the selection rules enables the geometric mapping of the energy detected at a receiver location 735 on to mixing zones 730 along the narrow beam. A time image of the non-linear property can thus be constructed along the narrow beam. By rotating in azimuth and moving the assembly along the borehole, a three dimensional time image can be constructed of a volume centered on the borehole. Successive repetition of the measurement at different beam elevations, and altering the fjfi frequency ratio yields a series of three dimensional time images. This redundancy in imaging permits further refinement of the smooth background model and a three dimensional spatial image.

Non-linear parameters of rocks have been found to be related to a number of important hydrocarbon reservoir parameters, such as variations with gas, oil and water saturation, effective stress, fracture density and mineralogical content. In certain aspects of this disclosure, the 3D images of non-linear properties constructed by the above method are transformed to provide quantitative information on the distribution of these properties around the borehole at the time of recording. In addition, sequential repetitions of this method are used to detect changes in reservoir properties over time for reservoir monitoring purposes.

The recordings of received waveforms are processed to generate an image of the non-linear characteristics of the formation. The directivity of the beam and the time of flight may fix the locations where scattered waves are generated, thus, distinguishing this method from normal sonic imaging techniques using conventional non-directional monopole and dipole sources.

By way of a non-limiting example, when a primary compressional (P) wave with a frequency fi and a primary compressional (P) wave with a frequency cross in a non-linear medium, and the selection rules are honored, a third shear (SV) wave can be generated with a frequency fi -f₂. This particular configuration can be used for creating 3D images of the Vp/Vs velocity ratio and non-linear properties of the rock formations around the borehole for many ranges of distance of investigation from a borehole. This particular example of compressional non-linear mixing to generate a shear wave (i.e., P + P ^ SV) will be used to describe a number of new concepts, methodologies, processes and systems for measurement and analysis purposes in the following paragraphs . These are equally applicable to P+SV SV or any permutation of non-linear mixing of two acoustic compressional or shear waves to generate a third wave.

Figure 8 shows an example configuration for a borehole-based system for remote mapping of non-linear properties and Vp/V s velocity ratio of rock formations using non-collinear acoustic mixing in accordance with an aspect of the present disclosure. Two primary acoustic beams, for example compressional (P) waves, from upper and lower arrays of transmitters 801 and 802 located in borehole 800, are directed into the rock formation surrounding the borehole. The transmitter arrays can be oriented such that acoustic energy propagates at any azimuth

809 and <fh 811 and elevation CC\ 810 and <¾ 812 relative to the borehole axis. For suitable elevation angles <¾, <¾,and azimuth angles

and φι, the two primary P beams propagating through the rock formation intersect with convergence angle Θ 804 at an mixing zone 805, remote from borehole 800. The convergence angle Θ is defined as the angle between the directions of the two converging beams, represented on Figure 8 as the lines joining the two transmitters 801 and 802 to the mixing zone 805. If the rock formation at the point of intersection has non-linear properties, a secondary (S) shear wave SV (e.g., a shear wave polarized in the plane defined by the axes of the two intersecting compressional waves) is generated due to non-linear interaction. The secondary shear wave propagates in a direction defined by the selection rules, represented by scattering angle, ψ, 806. The scattering angle ψ, is defined as the angle between the axis of the acoustic wave from the lower transmitter and the axis of the scattered wave. In the configuration shown, energy 807 returns to the borehole and is recorded at receiver or array of receivers 803.

As discussed above, the conditions suitable for generation of a secondary shear wave can be inferred from the selection rules that can be derived by the conservation of energy and conservation of momentum. The secondary wave S must obey either of the following conditions

k₃ = ki - k₂ (3) or

k₃ = ki + k₂ (5) where ki, k₂ and k₃ are wave vectors. The first frequency condition, where = fi -f₂ , is of particular interest for investigating properties of rock formation near a borehole. As shown in Figure 9, the conditions (2) and (3) can be represented by the formation of the wave-vector triangles and satisfied by the following relationships (6), (7) and (8).

| k₃ | = 2_rt| ; - , | / Vs = ( |/, -/ | /Vp ) X (Vp/Vs) (6) I ki I = 2π/,/νρ (7) | k₂ | = 27i£/Vp (8)

Using trigonometry in the vector diagram of Figure 9b, it can be shown that these conditions can be met when equations (9) and (10) as satisfied.

sin (ψ) = (f₂/fi) sin (Θ) / (1 -/ /})

(10)

Since the Vp/Vs velocity ratio is in the range 1.5 to 3.0 for many sedimentary rocks, there are combinations of convergence angle Θ and frequency ratio fl/f2 (denoted herein as d) that permit the generation of a secondary shear wave SV that propagates back to the borehole in the configuration shown in Figure 9a.

As discussed above, the behavior of acoustic energy generated by non-linear interaction of intersecting non-collinear planar waves can be calculated. The selection rules define a set of permitted interactions. The P+P→SV interaction is additional to the partial listing on Table 1. These various interactions use certain combinations of convergence angle, frequency ratio and scattering angle that depend on whether the converging waves are compressional or shear, on the Vp/Vs velocity ratio of the material at the interaction location, and differ for interactions generating sum ( }+ ¾ and difference ifi-fi) frequency energy. The geometry presented in the following examples is based on the

interaction, and analogs could equally be presented for other permitted interactions generating difference or sum frequency resonance, for example, but not limited to, interactions including P( })+SV½) SV( }- ,) and P( })+SV½) P( - ₂).

For example, considering P^+P^^ SV ifi-fi), Figure 10a shows mixing coefficient W as a function of the frequency ratio of the two sources. The mixing coefficient W, which is a measure of the amplitude of the efficiency of the conversion generating the scattered wave, is given by equation (1 1).

2(Vs/Vp)

(11)

where D depends on the Lame coefficients λ and μ and is proportional to for a given mixing zone, as defined in equation (12).

D = -

4π(λ + 2μ)

(12)

Θ is the convergence angle of the two primary beams at the mixing zone, and m is a scaling factor related to the Landau-Lifshitz non-linear constants A and B, as expressed in equation (13),

A _n

m =— + B

2

(13)

Hence, m is constant for a given mixing zone.

Figures 10a, 10b and 10c show the dependence of mixing coefficient W, convergence angle and scattering angle on Vp/Vs velocity ratio in the range from 1.5 to 2.0 at the mixing location, assuming a representative value (from Korneev, Nihei and Myer 1998) for m of -3660 GPa. Figure 10a shows the dependence of mixing coefficient on the plane wave frequency ratio for a range of mixing zone Vp/Vs ratios. Figures 10b and 10c show the corresponding convergence and scattering angles that honor the selection rules for the P+P SV interaction. As can be understood from Figure 10a, the mixing coefficient reaches a maximum at a ratio of the second frequency f2 to the first frequency fl equal to about 0.7. In addition, from Figure 10b, it can be appreciated that when the ratio of the second frequency f2 to the first frequency fl is equal to about 0.7, the convergence angle is in the range between about 30 deg. and about 40 deg. Moreover, it can also be appreciated that when the ratio of the second frequency f2 to the first frequency fl is equal to about 0.7, a scattering angle of the returning wave relative to a direction the first wave is equal to about 40 deg.

Figures 1 1a to 11c show example results of numerical simulations of non- collinear interaction of acoustic beams in a non-linear medium leading to the generation of a scattered wave that returns to the borehole. In Figures 11a to 1 1c, borehole 1100 includes lower transmitter 1 101, upper transmitter 1102 and receiver array 1 103. Acoustic energy in the form of a compressional beam generated by lower transmitter 1 101 and upper transmitter 1102 converges at a distance from borehole at convergence angle 1104. Receiver array 1103 is arranged to receive scattered wave 1 107 which is produced by the interaction of the acoustic beams as defined by the selection rules, described above and again below, at scattering angle 1 106. In figures 1 1a, l ib and 11c, distance along the borehole in meters is shown versus radial distance away from borehole in meters. In Figure 1 1a and b, the volumetric strains associated with the compressional acoustic energy of the two beams produced by the two sources are shown. In Figure 11c, the shear strain associated with the scattered shear wave 1 107 is shown.

In the following paragraphs a coded scheme that can be used to enhance or extract measured scattered acoustic waves that originate from non-linear mixing of primary acoustic signals in a mixing zone within a rock formation around the borehole is described. Measurements of scattered acoustic waves generated from non-linear acoustic phenomena in rock formations away from the borehole can be enhanced by broadcasting coded primary acoustic signals, recording the returning signal from the non-linear interaction and subsequently using a waveform recognition method and/or a band pass filtering method based on the forecast properties of the non-linear signal. One non-limiting example is to correlate the recorded signal with a template representing a forecast of its timing and frequency content derived from the broadcast parameters in accordance with the selection rules. The result of this correlation represents an acoustic pulse traveling past the borehole. For example, in the case of the P + P ^SV interaction, the returning energy resulting from the nonlinear interaction appears to travel along the borehole at a velocity equal to the shear wave velocity of the formation divided by the cosine of the angle between its propagation direction and the borehole axis. When applied to a system such as that illustrated in Fig. 8, the coding and correlation method improves the detection of the weak non-linear signal and thereby enhances the construction of 3D images of nonlinear properties and Vp/V s velocity ratio in the volume probed by the measurement

One non-limiting implementation for utilizing coding and signal correlation for the enhancement of signals generated by non-collinear acoustic mixing in a non- linear medium is described in the following paragraphs. A detailed description for a more general implementation of the coding scheme can be found in U.S. Patent Application entitled "System and Method for Generating Micro-Seismic Events and Characterizing Properties of a Medium With Non-Linear Acoustic Interactions," with the Attorney Docket No. 091230-0397046, filed concurrently with the present application, the entire contents of which is incorporated herein by reference.

Referring to Figs. 12a through 12d, source 1201 starts to broadcast a first coded or modulated time train ui (t), consisting of a sequence of N acoustic pulses, at time t =

0. The nth acoustic pulse has frequency f_n and an amplitude envelope El_n(t-T_n) of limited time duration, where n=l,2...N and T_n is the broadcast time of the nth pulse. The time separation between sequential pulses is variable. In some embodiments, the time separation between sequential pulses is much longer than the time duration of the individual pulses and the pulses do not overlap. A second coded or modulated time train, u₂(t), is broadcast from source 1202. This second coded or modulated time train consists of a sequence of N sequential acoustic pulses and starts at time t = δ, where δ is a start time difference between a start time of a broadcast of the first coded train and a start time of a broadcast of the second coded train. In one embodiment, the start time difference can be understood as a time delay between the broadcast of the first coded train and the second coded train. As it can be appreciated, the broadcast of the first coded train can be delayed relative to the broadcast of the second coded train or vice versa. In the present description it is often referred to δ as being a "time delay". However, δ should be interpreted broadly to be a "time difference" as the second pulse sequence may start before the first. The nth acoustic pulse has frequency d* f_n and an amplitude envelope E2_n(t-(T_n+ δ)) of limited time duration, where n=l,2...N. The frequency ratio between corresponding pulses in the two trains is fixed at d. T_n+ δ is the broadcast time of the nth pulse. The amplitude envelopes El and E2 of the first and second coded signal trains, respectively, can be different or the same. Examples of the time coded signal trains are shown in Figs. 12b-12c. These coded signals can be represented mathematically by following formulas (14) and (15). In the present description the symbol " * " is used for as a multiplication operator. wi( ⁼∑ E\n(t - Tn) * exp (i 2π* β * {t -Tn )} exp(7 ζ^η)

n

(14)

u₂ (t - S) =∑ E2n(t - (Tn + S)) * exp(72π * d * /» * (t - (Tn + S))) * exp(z ζ_η) n

(15) where ζ_η is a phase of each pulse n and εχρ(ίζ_η) is a phase term of each pulse n.

The respective signal envelopes El_n(t-Tn) and E2_n(t-(Tn+d)) can have any shape or form such as a Gaussian form, etc. Similarly, although the modulated signals within the envelopes are expressed in equations (14) and (15), they can be modeled by other mathematical formulas. When the time delay δ is equal to the difference in travel times tl and t2 (that is to say δ = tl - t2) from transmitters 1201 and 1202 to the mixing zone 1204, the acoustic energy of corresponding pulses n from the two broadcast trains arrives simultaneously at the mixing zone and, if the convergence angle, frequency ratio and Vp/V s ratio at the mixing location are accord with the selection rules' criteria, generates a third series of scattered acoustic pulses with dominant frequency (1 - d *f_n, equal to the difference between the frequencies in the two primary pulses, f_n and d *f_n. This third wave, denoted w₃, recorded at the receiver, inherits the coding of the two primary signals and therefore may be expressed as equation (16).

« 3(

E3n(t - {T_n + Τ))* εχρ(ϊ2π * (l - d)* f_n * (t - {T_n + Γ)))· exp(i£) n

(16) where ζ_η is a phase of each pulse n and εχρ(ίζ_η) is a phase term of each pulse n of the third wave.

E3_n(i) is the resulting amplitude envelope due the mixing of the primary pulses and T is the total travel time from the source 1 to the recording receiver via the center of the mixing zone as further explained below.

In one embodiment, signal enhancement can be accomplished using a correlation technique well known in seismic processing industry to extract a relevant part of the non-linear interaction and travel time information in measured signal m.

The technique involves the construction of a template signal u_s that has the coded signal form as expressed in equation (17). u_s{t) = X Wn(t - T_H νφ(ί2π * g (/„) * (t - T„)) * exp(i£)

n

(17) where ζ_η is a phase of each pulse n and εχρ(ίζ_η) is a phase term of each pulse n of the template signal.

W_n is some suitably chosen envelope function and g(f_n) is some suitably chosen function of frequency^ . The selection of an appropriate function g(f_n) may be based on the shape of the expected modulated signal within the measured signal ¾ to achieve the best non-linear signal extraction. For example, g(f_n) can be selected to be (1 - d) * f_n such that u_s is equal to

^Ms( ⁼∑ Wn t - T_n )exp(z^'2 r (1 - d)(t - T_n ))exp(z^'^) as shown in Fig. 12d. However, n

other functions g(f_n) are also within the scope of the present invention. To extract the signals from non-linear interactions in wj, a cross-correlation between 113 and u_s can be performed to obtain the correlated signal M where M is mathematically defined by equation (18).

(18)

It can be shown mathematically that the resulting correlated signal M is a sharp band-limited spike with a bandwidth including all the frequencies (l-d) * fn for n=l,2,...,N when the number of pulses N is large. The term "band-limited" spike is used herein to refer to a spike signal having a limited frequency bandwidth. In one embodiment, increasing the broadcast duration by increasing the number of pulses N enhances signals generated by non-linear interaction in the correlated signal M while more effectively suppressing signals generated by linear interactions and other noise. It should be noted that the correlation technique using the coded signal pattern is one of many ways to extract and enhance signals generated by non-linear interactions. Alternative signal processing techniques including pattern recognition or frequency band filtering could equally be used for signal extraction and enhancement.

The correlated measured signal M has the following properties. First, the correlated signal contains a sharp band-limited spike, corresponding to the non-linear interaction at the mixing zone, only if the time delay δ between first and second primary coded signals is equal to the difference between the travel time ti from the first acoustic source 1201 to the mixing zone 1204 and the travel time ¾ from the second acoustic source 1202 to mixing zone 1204, i.e. δ = ti - t₂. If this condition is not met, the correlated signal is highly suppressed. Second, if the condition δ = ti - t₂ is met, the band-limited spike occurs on the correlated signal M(t) at the time T which is equal to a sum of the travel time from the first primary acoustic source to the mixing zone and the travel time from the mixing zone 1204 to the receiver within the receiver array 1203, i.e. T = t; + 1₃ = δ + 1₂ + 1₃. Third, increasing the duration of the coded signal train, i.e., increasing the number of pulses N in the broadcast train, improves discrimination of signal from noise_e, because the noise does not have the form of the template signal u_s. Numerical simulation results of a case where the lower transmitter 1201 and upper transmitter 1202 emit two coded signal trains consisting of sequential acoustic pulses with Gaussian envelopes are shown in Fig. 13a- 13d. In this non-limiting example, coded signals of 12 pulses are used with frequency pairs (574Hz, 373Hz) (624Hz, 405.6Hz) (700Hz,455Hz) (757Hz,492Hz) (802Hz, 521.3Hz)

(870.5Hz,566Hz) (947Hz, 615.5Hz) (1000Hz,650Hz) (1 120Hz, 728Hz) (1218Hz, 792Hz) (1324Hz, 861Hz) (1440Hz,936Hz). The frequency ratio f₂/fi between the pairs is a constant 0.65. The start time delay δ between the two signal trains is chosen to be equal to (tl- t2). Numerical simulation of the non-linear interaction due to the broadcast of the two coded wave trains ui and ¾ is performed on a computer. The simulated signals due to non-linear interaction from a broadcast of two sequential pulses received and recorded at six of 1 10 receivers, indexed from 1 to 1 10, in the non-limiting example of a receiver array 1203 are shown in Fig. 13(b). The template u_s for the returning coded signal is shown in Fig. 13(c). The result of correlation between the template signal with the recorded signal at each receiver is shown in Fig.

13(d). The resulting correlated signal at each receiver shown in Fig. 13(d) shows a very sharp band-limited spike. This sharp spike occurs at the time T = tl + 13 where tl is the travel time from the source 1 to the center of the mixing zone and t3 is the travel time from the center of the mixing zone to the receiver. The time delays or "move- out" across the receiver array are as though the acoustic energy is traveling along the borehole with an apparent velocity equal to the shear wave velocity divided by the cosine of the angle between the direction of the returning wave and the borehole axis.

The numerical simulation shown in Figs. 13a-13d clearly illustrates the power and utility of the coding scheme when used in conjunction with the system of measurement in Figs. 12a-12d or Fig. 8. It allows for computer processing of the recorded signals at the receivers to generate correlated records that contain band- limited spike signals with strength proportional to the strength of the non-linear interaction at the mixing zone 1204. The arrival time T of the band-limited spike is equal to the total travel time from source 1201 to mixing zone 1204 and back to the borehole at receiver 1203. The amplitude of the band-limited spikes vary with the receiver position with a maximum occurring at a particular receiver, the location of which is dependent on the scattering angle ψ 1206 of the non-linear interaction at mixing zone 1204. The scattering angle ψ is dependent on the properties of the rocks, e.g. Vp/Vs velocity ratio, at the mixing zone 1204. It should be noted that this result is a characteristic of the coding scheme and of the measurement system shown in Fig. 12a-12d or Fig. 8. The use of Gaussian envelopes and coded signals in conjunction with templates are non-limiting examples used for the purpose of illustrating the coding scheme and its characteristics. Variants of ui,U2 and u_s can be considered in order to optimize the performance of the correlation process in term of resolution and signal to noise ratio in response to various considerations imposed by field applications.

In some aspects of the present disclosure, coded acoustic signals in the primary acoustic beam can also be used to enhance the amplitude and focusing of the non-linear signal returning to the borehole, and to improve signal detection sensitivity and signal to noise ratio. Figures 14a-14d show an example of an application of the coded signal scheme to a noisy time series signal generated by numerical simulation. The noisy times series signal simulates a signal returning to the borehole as a result of non-linear interaction. White Gaussian noise with an amplitude 10% larger than the amplitude of the non-linear interaction signal is added to the time-series signal produced by the numerical simulation of wave propagation in a non-linear model before the correlation with the coding template is applied. The configuration is the same as that shown in Figs. 12a-12d and Figs. 13a-13d. Fig. 14(b) shows the simulated receiving signal containing noise recorded at 6 receivers of the receiver array 1403. Fig. 14(d) shows the signal retrieved from the noisy signals (in this case the simulated noisy signals) on the same receivers when correlating with the coded template u_s (t) of 12 pulses shown in Fig. 14(c). The coding scheme is thus shown to effectively extract the signal from the non-linear interaction and minimize the noise, a useful characteristic for field applications.

The preceding text described how a combination of broadcast coding and signal processing may be used to improve the detection and to determine the amplitude of an acoustic wave signal generated by non-collinear mixing in a nonlinear medium and to determine the acoustic travel time from sources to receiver array via the mixing zone and to infer the mixing zone's non-linearity and Vp/Vs ratio. The following paragraphs describe a non-limiting implementation of the coding and signal processing method in the context of a borehole based measurement system. Given an initial compressional and shear velocity model of the formations around the borehole, techniques such as ray tracing may be used to estimate the location of a mixing zone, 805 on figure 8, corresponding to chosen inclinations and azimuths (¾, ( 2, φι and φ₂, 809 to 812) and also the convergence and scattering angles (Θ and ψ, 804 and 806). This information is used to predict the time delay δρ and frequency ratio dp required for the pulses from the two transmitters to arrive simultaneously at the mixing zone and generate a sequence of scattered pulses by non-collinear interaction that arrive at a receiver, and to predict the total travel time Tp from the first acoustic source to a receiver via the mixing zone. Systematically scanning δ and d around their predicted values, and correlating the template u_s(t) with the recorded signal U3(t) results in a suite of results M (t, δ, d) for each receiver element of the array 803. As discussed in previous paragraphs, the correlated signal M (t, δ, d) for each receiver element contains a band-limited spike corresponding to the non-linear interaction in the mixing zone 805 if selection rules permit non-linear interaction. A search is made in the space of (t, δ, d) on M(t, δ, d) to identify the location (T_NL, 6NL, d^) of the band-limited spike corresponding to non-linear interaction at mixing zone 805, which should be in the vicinity of (Tp, δρ, dp). Differences between (T_NL, 6NL, d^) and (Γ_Ρ, δρ, dp) are indications of deviation of the propagation velocity model from the true propagation characteristics of the rock formation. These differences are then used to update the propagation velocity model by velocity tomographic inversion or other velocity updating methods. Given the relationships between frequency ratio d, Vp/Vs ratio, convergence angle and scattering angle from equation 9 and figure 10, the Vp/Vs ratio of the rock formation at the mixing zone is then calculated from the observed value dNL, the convergence angle and the scattering angle at the mixing zone - the latter two quantities can be calculated by raytracing or other numerical methods from the updated propagation velocity model. The amplitude value M (T_NL, 6NL, C NL) contains information about the non-linear mixing strength of the rock formation at the mixing zone 805. It can be used for 3D-imaging of non-linear properties of the rock formation as described in the discussion below.

The ability to process non-collinear acoustic mixing signals by correlation with a template signal to identify the combination of parameters associated with a particular mixing zone, that is M (T_NL, 6NL, C NL), T_NL, 6NL and dNL, that can be subsequently used for the determination of Vp/Vs ratio, propagation model Vp and Vs and nonlinear properties is a direct consequence of the coding and correlating protocol. This unique characteristic underpins the imaging scheme that follows.

The discussion in the paragraphs below will focus on one non-limiting example of the use of coded acoustic signals for creation of the 3D images of Vp/Vs ratio and non-linear properties of rock formation surrounding the borehole. However, the use of coded acoustic signals has broader applications beyond those related to geology and petrophysical applications including the areas of non-destructive testing and medical imaging.

By way of a non-limiting example, one measurement and data processing protocol to enhance signal to noise ratio is discussed in the following paragraphs for the system described in Fig. 8.

For the sake of simplicity, the operating protocol is described in the context of a non-attenuative rock formation. However, the following measurement and data processing protocol can equally be applied to an attenuative rock formation. The effect of attenuation in the rock formation is to displace the origin of the scattered wave by a predictable amount related to the formation Q.

First, for a given ratio frequency ratio d of the two coded signals as described in the above paragraphs, the coded primary acoustic beam signals from first acoustic energy source 801 and second acoustic energy source 802 of the measurement system described in Fig. 8 are transmitted into the earth. A time delay δ is maintained between two coded signals of source 801 and source 802. The measurement geometry honors the selection rules, and the time delay is such that energy from the two sources arrives substantially simultaneously at the mixing zone. In one embodiment, each component of the three component geophone at receiver 803 at location z3 detects and records the returning acoustic waves. This measurement is denoted m(t, δ, d, z3) . Alternatively or in addition a hydrophone may also be provided to detect a pressure signal in the returning acoustic waves. Second, the coded primary acoustic signals are transmitted from the first and the second acoustic energy sources as described in the above paragraphs but with polarity reversed, that is, with the phase shifted by 180 degrees. The signal recorded at receiver 803 is denoted m. (t, δ, d, z3).

Third, the two signals m(t,5,d,z3) and m.(t, δ, d, z3) are added together to form the combined signal, which can be denoted as mm(t, δ, d, z3). Because the signals m(t, δ, d, z3) and m.(t, δ, d, z3) have opposite polarity, signals from linear interaction in the rock formation will be cancelled out by the addition of m(t, δ, d, z3) and m.(t, δ, d, z3). However, the non-linear responses of the earth will add coherently since the amplitude of the non-linear responses is proportional to the product of the amplitudes of the two primary signals and therefore will not reverse polarity when the polarity of both primary signals is reversed. Therefore, mm(t,5,d,z3) would essentially contain a signal from non-linear interaction of the rock formation.

Fourth, a time variant band-pass filter can be applied on the obtained signal mm(t,5,d,z3), so as to keep a narrow band around the expected bandwidth of the signal. The bandwidth of the obtained signal is determined from the frequency differences and bandwidths of the two primary broadcast signals.

Fifth, a cross-correlation of the filtered signal mm(t,5,d,z3) with the template coded signal is performed as described in the above paragraphs to obtain the pulsed signal which can be denoted as mmc(t, δ, d, z3).

Sixth, hodogram analysis can be applied to the three component data obtained from the three uniaxial sensors of the receiver and/or applied to the pressure signal detected by the hydrophone. These data may be used to analyze any of the possible modes P, SH and SV and can be transformed to obtain separate measurements of any SV, SH and P arrivals, denoted mmcr(t, δ, d, z3).

Seventh, the above six steps can be repeated multiple times with different broadcast coded signals and the collection of mmcr(t,5,d,z3) signals can be stacked to improve signal to noise. The resulting stacked signal is the signal record M(t,5,d,z3) for each of the SV, SH and P arrivals. For example, in a P+P to SV mode, the SV mode is detected by the receiver while in a P+SV to P mode, the detected signal at the receiver would be a P mode. Eighth, the above seven steps can be repeated for a sweep through for a sequence of many values δ and d to obtain the entire set of M(t, δ, d, z3). The described measurement and processing protocol allow for the construction of the measurement signals (t, δ, d, z3) with high signal to noise ratio at the receiver arrays.

The measurements M (t, δ, d, z3) can be repeated for many transmitter locations zl, z2 and many receiver array locations as the transmitter arrays and the receiver arrays can be moved independently. As the primary acoustic beams from source 801 and 802 can be steered independently for any azimuth angles φΐ, φ2 and elevation angles al, a2, the measurements M(t,8,d,z3) are also repeated for many angles φΐ, φ2, al and α2. These repeated and multiple measurements may contain multiple redundant signals generated by non- linear interactions in the earth for many values of zl, z2, z3, φΐ, φ2, al and al. The redundancy allows for additional signal to noise enhancements by signal processing on computers and for creation of 3D images of rock properties around the borehole.

It should be noted that the steps described above for the measurement and processing protocol can be re-ordered or eliminated in various permutations as warranted. Furthermore, there are many additional signal processing techniques familiar to those who are experienced with the art of seismic signal processing, e.g. multidimensional filtering, time moveout analysis and stacking. These additional techniques can be added to the measurement and data processing protocol described in the above paragraphs to improve the quality of the recorded data and processed images. The correlated acoustic signal M(t, δ, d, z3) from non-linear interactions has many properties rooted in the coding methodology and the selection rules.

In the following paragraphs, an imaging method and workflow that exploits these properties to construct 3D images of the non-linear properties and Vp/Vs velocity ratio of an earth volume and determines other rock properties such as Vp and Vs is discussed. A non limiting example of the workflow is discussed below.

Referring to the measurement system described in Fig. 8 with travel time notations denoted in Figs. 12a- 12d, and the discussion in the above paragraphs, the correlated signal record M (t, δ, d, z3) after measurement and processing at a receiver location z3 for one particular transmitter location and beam angle (zl, z2, φΐ, φ2, al and α2) will contain a band-limited spike at the travel time T= tl+ t3 if the following conditions are satisfied:

a) The transmitted beams intersect and interact non-linearly at mixing zone

805.

b) The time difference δ is equal to the travel time difference tl— 12.

c) The selection rules are obeyed, i.e., the frequency ratio d used

coding scheme obeys the condition of equation (9).

(19)

where Vp/V s is the compressional to shear velocity ratio at the mixing zone 805 and Θ is the convergence angle between the first and second transmitted beams. The Vp/Vs ratio and convergence angle Θ can be calculated by raytracing from the beam geometry, location and beam direction parameters (zl, z2, z3, φΐ , φ2, al and α2) and a model of compressional velocity Vp and shear velocity Vs of the rock volume being investigated.

The above conditions dictate that the non-linear interaction at mixing zone 805 contributes to the record M (t, δ, d, z3) as a band-limited spike signal at a single point in the (t, δ, d) space for each position z3 of receiver 803. The observed location of the band-limited spike signal on the record M (t, δ, d, z3) can be denoted as (TNL, ½L_> NL, 3). The amplitude of this band-limited spike is a function of the non-linear interaction strength, β, at mixing zone 805.

A combining process (e.g., a stacking process) as commonly used in the seismic industry and borehole acoustic waveform analysis can be performed on the signal records M (t, δ, d, z3). The progressive time delay or "moveout" of the band- limited spikes on the signal records M (t, δ, d, z3) can be analyzed as a function of z3. For example, in a stacking process, the signal records can be stacked to obtain a stacked record Ms(t, δ, d, z3r) corresponding to a chosen reference location z3r. The reference location is selected such that the plurality of signal records M (t, δ, d, z3) can be stacked appropriately. This stacking process enhances the signal to noise ratio and improves the detectability of the location (TNL, &NL_> ^NL, ^z3r) of the band- limited spike originating from the non-linear interaction on the record M s(t, δ, d, z3r). Although a stacking process is described herein, other combination techniques can also be used.

Given an initial Vp and Vs velocity model, which can be estimated from the borehole well logs and assumptions about lateral continuity of the rock formation properties away from the borehole, the trajectories and travel times of acoustic beams with elevation angle al and azimuth angle φΐ from transmitter 801 at position zl and elevation angle al and azimuth angle φ2 from transmitter 802 at position z2 can be calculated by ray tracing or other numerical modeling techniques.

Therefore, the location of the mixing zone 805 can be located in space by ray tracing or other acoustic numerical modeling techniques from the parameters zl, z2, z3r, φΐ, φ2, al and α2 if the beams intersect. The predicted travel time tip, t2p and t3p between transmitter locations 801 and 802, reference receiver z3r and mixing zone 805 can also be calculated based on position and velocity using ray tracing or other acoustic numerical modeling techniques.

The predicted pulse arrival time T_p = tlp+t3p and time difference δ_ρ = tip - t2p can be predicted from the compressional velocity Vp and shear velocity Vs model and then compared with the arrival time T L and difference time 5NL of the observed band-limited spike on the actual record Ms(t, δ, d, z3r). If the Vp and Vs velocity model correctly approximates the true Vp and Vs velocity of the rock formation, the predicted times (Tp, δ_ρ ) are equal to observation times (TNL, 6NL). If there are differences between the predicted times and the observed times, these differences can be used to update the propagating Vp and Vs velocity model to minimize the differences and achieve consistency between the modeled and observed data. Various iterative velocity tomographic inversion methods that are familiar to those experienced in the art of imaging in the seismic processing industry can be used to update the propagating Vp and Vs model.

Given the Vp and Vs propagation model obtained through the above tomographic velocity inversion step, the trajectories of acoustic beams with elevation angle al and azimuth angle φΐ from transmitter 801 at position zl and elevation angle a2 and azimuth angle φ2 from transmitter 802 at position z2 can be calculated by ray tracing or other numerical modeling techniques For a given azimuth angle φΐ, there will be an azimuth angle φ2 for which the source beams will intersect at a mixing zone 805 for which the location and convergence angle Θ can be calculated by ray tracing or other numerical modeling techniques from the parameters zl, z2, z3r, φΐ, φ2, al and α2. The amplitude of the pulse at the point (TNL, 6NL, dNL) on the record Ms(t, δ, d, z3r) can then be mapped to the spatial coordinates of the mixing zone 805. Since dNL must obey Equation (19), the Vp/Vs velocity ratio at the mixing zone 805 can be calculated from Equation (19) and mapped on to its spatial position. By repeating the above mapping step for a range of values of parameters zl, z2, z3, φΐ, φ2, al and α2, a 3D image of the non-linear interaction strength β and a 3D image of velocity ratio Vp/V s can be constructed. The velocity Vp/V s ratio obtained from the above mapping method using selection rules is an alternative method to using the ratio of the propagation velocity Vp and Vs obtained from tomographic inversion of travel times.

The geometric mapping step discussed in the above paragraphs is only one example of many imaging techniques that can be used for 3D imaging for non-linear property and Vp/Vs velocity ratio from the signal records M (t, δ, d, z3) for many values of parameters zl, z2, φΐ, φ2, al and α2. Other advanced imaging techniques such as Kirchhoff imaging, beam imaging, wave equation imaging used in the seismic industry can be adapted for the 3D imaging of non- linear property and Vp/Vs velocity ratio. For example, the three dimensional image of the propagation compressional velocity Vp, the three dimensional image of the propagation of shear velocity Vs, the three dimensional image of the ratio of compressional velocity and shear velocity Vp/V s, or the three dimensional image of non-linear properties of a rock formation, or any combination of two or more thereof can be performed using Kirchhoff imaging, beam imaging or wave equation imaging. In addition, in one embodiment, the value of propagation compressional velocity or value of shear velocity or both can also be determined using tomographic velocity inversion or full wave form inversion or by iterative imaging in combination with tomographic velocity inversion or full wave form inversion.

As it can be appreciated, determining a value of a parameter can be different from imaging the parameter. Indeed, an image of a parameter may only contain relative values of the parameter and does not necessarily provide the information on the absolute value of the parameter. Therefore, obtaining an image of the velocity ratio Vp/Vs may be different from determining the value of the velocity ratio Vp/Vs. Determining a value of the velocity ratio from the image of the velocity ratio may require additional information.

For investigation of non- linear properties and the ratio of compressional to shear velocities farther into the rock formation from a borehole, sources of lower frequency, of order 500Hz to lOkHz, may be needed since acoustic energy of lower frequency can penetrate farther into the rock formations before being attenuated to a non-detectable level. Lower frequency acoustic energy in the 500Hz to 10 kHz frequency range has a wavelength much bigger than the borehole diameter. In these circumstances it is difficult to control the azimuth directions of the acoustic wave broadcast from sources deployed in the borehole. The configuration shown in Fig. 8 with two primary beams of acoustic energy emanating from a borehole may not be suitable for probing and imaging of non-linear properties at large distance greater than hundreds of feet from the borehole. In view of this, a different low frequency (500 Hz to 10 kHz) system and method is needed to maximize the distance of investigation away from the borehole. The following paragraphs describe a non-limiting example of such a system and method to generate three-dimensional images of non-linear acoustic properties and Vp/Vs velocity ratio using low frequency sources. Although the following description will refer to the frequency range 500 Hz to 10 kHz, the system and method described can also be applicable and beneficial in a higher frequency ranges, for example, 10 kHz to 500 kHz.

Figure 15a shows a system of two transmitters and a receiver or a receiver array that are arranged in the borehole to detect non-collinear mixing in a rock volume around the borehole. In one embodiment, the upper transmitter 1502 includes a linear array of transmitters which can be clamped to or undamped from the borehole. The lower transmitter 1501 includes a linear array of transmitters which can also be clamped or undamped. The receiver array 1508 includes a clamped three component receiver or a clamped three component receiver array. The transmitters and receivers can be moved together or independently. In one embodiment, array elements of transmitters 1501 and

1502 can be arranged, for example using phase control, to broadcast acoustic energy into two cones (lower acoustic cone 1504 produced by transmitter 1501 and upper acoustic cone 1505 produced by transmitter 1502), with axes collinear with transmitters 1501 and 1502. In a phase controlled system, the cone angles depend on the phase difference between transmitter elements and the rock formation velocity. The intersection of the two conical broadcasts of acoustic energy is a toroidal shaped intersection volume 1506. Where the selection rules are honored, scattered energy 1507 is generated by non-linear interaction around intersection volume center 1506. The scattered energy originating from non-linear interaction is recorded at the receiver or array of receivers 1508.

Because the vertical locations and the elevation angles of the conical broadcasts of the transmitters 1501 and 1502 are controllable, the spatial location (distance from the borehole and vertical location) of their toroidal intersection volume 1506 may be controlled and scanned over a three dimensional rock volume around the borehole 1500. To the extent that the convergence angles of the two cones and the frequency ratios of the two sources can be arranged to honor the selection rules, the scattered signals due to non- linear mixing at mixing zone 1506 that are recorded at receiver 1508 contain information of the non-linear properties of the formation at the intersection between the two cones. The data recorded by an array of three component geophones 1508 may be analyzed to determine the azimuth and elevation to its origin. The signal coding methodology described in the above paragraphs and the measurement and processing protocol described above can be applied to the system depicted in figure 15(a). The correlated signal M (t, δ, d, z3) obtained from the latter system and configuration is composed of the linear superposition of all the pulses generated by the non-linear interaction at all intersection volume segments 1509-1, 2, 3, etc. to 1509-k, that cover the entire circumference of the toroid intersection volume 1506 as shown in Fig. 15(b). As discussed in the above paragraphs, the contribution of the non-linear interaction of the wave by a segment 1509-k, corresponding to the interaction or mixing zone, to the resulting correlated signal M (t, δ, d, z3) at a receiver at position z3 in receiver array 1508 is a pulse with a travel time equal to the time of flight from transmitter 1501 to segment 1509-k and from 1509-k to receiver if the following conditions are met: a) the frequency ratio d, convergence angle Θ and Vp/V s velocity ratio at the particular segment 1506-k obeys the selection rule condition of equation (19),

b) the time delay δ between the coded source transmitters signals is equal to the difference between the times of flight from transmitter 1501 to 1509-k and from transmitter 1502 to 1509-k.

In other words, the pulses generated by non-linear mixing at the intersection volume segments 1509-1 to 1509-k in the correlated signal M (t, δ, d, z3) are distributed over a range of (t, δ, d) for each receiver at position z3. The contribution of each segment 1509-k can be mapped to a point in (t, δ, d) space in the correlated signal M(t, δ, d, z3). This property, together with the trajectory information of the signal obtained from three component receivers allows for geometric mapping of the signal amplitude at the point (t, δ, d) in correlated signal M (t, δ, d, z3) to the spatial locations of the intersection volume segments 1509-1 to 1509-k within the toroid intersection volume using the imaging methods discussed in the above paragraphs. An initial Vp and Vs propagation model is constructed. Ray tracing, travel time analysis and iterative tomographic velocity determination are then performed to obtain an updated Vp and Vs propagation model. Non-linear properties and Vp/Vs velocity ratio in the segments 1509-1 to 1509-k in the toroid intersection volume can then be extracted and mapped to the spatial locations on the toroidal intersection volume using ray tracing analysis applied to the updated Vp and Vs model using the workflow described in the above paragraphs. By repeating the measurements for all elevation angles of conical broadcast and vertical locations in the borehole of the transmitters, the images of non-linear properties and Vp/V s velocity ratio can be constructed for all mixing zones surrounding the borehole. The images for all scanned mixing zones can then be combined to yield a complete 3D image of non-linear properties and Vp/Vs velocity ratio using a suitable processing method known to those who are skilled in the art of seismic imaging such as, for example, the weighted stacking method, of all the images.

The system and method using low frequency conical acoustic broadcasts for a vertical well can work well when the rock volume has no azimuthal symmetry. However, if the rock volume has a very high degree of azimuthal symmetry for the propagating velocities Vp and Vs and therefore by implication their Vp/V s velocity ratio, such a system may encounter some difficulties in resolving the azimuthal variations in formation properties and generating 3D images. Referring to Fig. 15(b), the non-linear signals generated by non-linear interaction from segments 1509-1 to 1509-K in the toroid intersection volume arrive simultaneously at each receiver when there is total azimuthal symmetry. In this case, it would be difficult to separate the pulses originating from multiple segments 1509-1 to 1509-K in the toroidal intersection volume as these occur at the same time in the correlated signal M. This introduces ambiguity into the mapping exercise. The above difficulty can be avoided if a borehole or borehole system is designed to overcome the limitations due to the azimuthal symmetry. In the following paragraphs, a non-limiting measurement system and method utilizing various wellbore configurations designed to achieve this goal will be described.

Figures 16a and 16b show an example of non-collinear mixing arrangement with receiver array 1603 located in a straight section of the borehole and centers of transmitters 1601 and 1602 located on the receiver array's extended axis 1625 . For the case of linear phased arrays, the transmitters broadcast acoustic energy in cones 1604 and 1605 with axes aligned with the transmitter arrays. The cones and their intersection locus 1607 are illustrated in figures 16a and 16b, in accordance with aspects of the present disclosure. As shown in figure 16a, the intersection locus may be considered as a series of adjacent intersection volume segments, 1608-k, where k=l to K, four of which are identified as 1608a, 1608b, 1608c and 1608d. From geometric analysis it can be demonstrated that if two intersection volume segments have the same convergence angles, their corresponding time of flight differences (from lower and upper transmitter to intersection volume segment) tl - t2 are different. Two locations on opposite sides of the intersection locus corresponding to the mixing zones, with similar convergence angles and Vp/Vs velocity ratio, may therefore be individually activated by controlling the frequency ratio d and the time delay δ between coded broadcasts from the two transmitters. Similarly, two mixing zones with the same time of flight difference tl — t2 have different convergence angles Θ, so that an acoustic signal arriving simultaneously at two mixing segments could only honor the selection rule requirements for a Vp/Vs velocity ratio, frequency ratio d and convergence angle Θ at one of the two locations. As it can be appreciated from the above paragraphs, the intersection of two conical acoustic broadcasts defines an intersection volume in a form of a toroid. For the purpose of the present disclosure, a toroid is defined as an "annular" shape that is generated by revolving a plane geometrical shape such as a polygonal shape, circle, ellipse or other shape to define a closed volume. In one embodiment, the toroid can be a ring torus or O-ring where the plane geometrical shape or cross-sectional shape is a circle that is revolved around an axis. In another embodiment, the toroid can be defined as a polygonal shape that is revolved elliptically to form a closed volume. The toroid volume can be segmented in to a plurality of intersection volume segments. Depending upon non-linear selection rules and a judicious selection of various parameters including a start time difference and frequency ratio between the broadcast acoustic signals, one a more intersection volume segments can be "activated" to provide one or more mixing zones where the two conical broadcast signals interact non-linearly within the mixing zone(s) to generate a signal that is representative of the non-linear properties of the rock formation in that zone.

The coding methods discussed in the above paragraphs, the measurement and processing protocols described in the above paragraphs and the imaging methods discussed in the above paragraphs can be applied to the measurement system illustrated in Figures 16a and 16b. The applications follow the general approach described in the above paragraphs. Unlike the straight borehole configuration of figure 15, there is no issue of simultaneous arrivals in the space of (t, δ, d). Thus, measurement system described in figure 16 using a borehole with an intentionally bent trajectory would be more robust for general 3D imaging of non-linear properties and Vp/Vs velocity ratio. A by-product of the imaging method is the Vp and Vs propagation model generated by tomographic velocity inversion of the intermediate step of the imaging process.

It should be noted that the absence of azimuthal symmetry in the configuration of the curved or angled borehole trajectory and tilted conical acoustic wave broadcast as shown on figure 16 can also be achieved on a smaller scale by offsetting the axes of smaller transmitter arrays within a straight borehole, or by a similar configuration on a wireline or pipe conveyed logging tool, etc.. Figures 17a and 17b show two additional examples of transmitter-receiver arrangements within various boreholes configurations of axial orientation and cone angle that can be used for the creation of the 3D images of Vp/Vs velocity ratio and non-linear properties of rock formation surrounding the borehole. Figure 17a shows an example single well arrangement with a cranked rat hole where there is a complete intersection of lower cone within the upper cone. For example, as shown in Figure 17a, upper transmitter 1702 is arranged to produce a wider cone of acoustic energy 1705 than the narrower cone 1704 produced by lower transmitter 1701, such that there is a complete intersection of the lower cone within the upper cone. Both transmitters 1701 and 1702 are positioned in rat hole 1706, while receiver 1703 is positioned in main vertical borehole 1700. As in the configurations of Figures 15 and 16a-16d, transmitters 1701 and 1702 can be arranged as an array, for example a linear array, of acoustic point sources in a borehole. Acoustic energy generated by transmitters 1701 and 1702 interacts in the non-linear material in the intersection zone and acoustic energy is received at receiver 1703 in accordance with the selection rules, as discussed above. The non-linear mixing zone includes an intersection locus 1709 between cones 1704 and 1705 spanning a closest point 1707 to the borehole 1700 and a farthest point 1708 from the borehole 1700.

Figure 17b shows another example of single well arrangement with a more severely cranked rat hole with a lower transmitter emitting energy near perpendicular to the axis of borehole 1700. In these configurations, the intersection between the cones is a hyperbola, not an ellipse as in the above examples. Nonetheless, similar measurements, coding, data processing and imaging protocols remain applicable. In addition, because the intersection between the cones is not a closed curve, depth of investigation is determined by source strength, receiver sensitivity, and the effectiveness of the signal processing algorithms.

One difference between a configuration resulting in a closed elliptical intersection locus with a configuration resulting in a parabolic open curve locus is that in a closed elliptical intersection locus a frequency ratio scan may start at a low frequency ratio f2/fl (see, for example, Fig. 10b) that corresponds to a convergence angle higher than that at the nearest point. The frequency ratio f2/fl can then be increased, for example, until the nearest point on the ellipse is activated. The scan can be continued until the farthest intersection point is reached. By increasing further frequency ratio f2/fl no signal would be received as no zone is activated. On the other hand, if the intersection is an open curve, the frequency ratio f2/fl can be scanned starting at the nearest point and continue until the returning signals to the receiver from either side of the hyperbola become undetectable.

Figure 18a shows an example of a vertical well and sidetrack with receivers in the vertical part of the well. Figure 18b shows another example of a vertical pilot hole and horizontal sidetrack with receivers in the sidetrack in accordance with aspects of the present disclosure. The main vertical borehole 1800 includes transmitter 1801 arranged to produce a vertical cone of acoustic energy 1804, while sidetrack borehole 1806 includes transmitter 1802 arranged to produce acoustic cone 1805. Receiver 1803 can be arranged in either main vertical borehole 1800 as in Figure 18a, sidetrack borehole 1806 or both main borehole 1800 and sidetrack borehole 1806 depending on the particular application used. As in the examples of Figures 17a to 17c, transmitters 1801 and 1802 and receiver 1803 can include an array of transmitters and receivers, respectively. Transmitter 1802 in sidetrack borehole 1806 can be located above or below transmitter 1801 in main vertical borehole 1800.

Although the configurations of Figures 17a to and 17b, 18a and 18b differ somewhat in complexity from the drilling and operation perspectives, there are many benefits to using different configurations including the ability to perform deeper remote sensing from the borehole by maximizing the distance to the intersection of the two cones. Moreover, transmitters generating acoustic energy in a direction near perpendicular to the borehole can provide more power and more angular resolution for the scans defined by smaller cone angles. It should be noted that the above Figures are only example configurations of the use of a multitude of possible borehole configurations. As it can be appreciated, there are many other borehole configurations that also allow for the placement of the two primary acoustic arrays at different azimuth angles and elevation angles.

The techniques described above, a combination of invoking the selection rules and signal coding, may be used to scan and image a volume defined by two intersecting cones formed by the acoustic energy from the two transmitters, some distance away from the borehole. The discussion above described the P+P SV interaction. However, other permitted interactions can also be used in a similar fashion with examples shown in Table 1.

As can be appreciated from the above paragraphs, various methods of investigating a rock formation can be implemented with a system for survey planning, data acquisition and storage, and image processing and interpretation. Figure 19 depicts a system for survey planning, data acquisition and storage, and image processing and interpretation, according to an embodiment of the present invention. In one embodiment, the system can generate 3D images of non-linearity properties, Vp/V s velocity ratio and propagation compression and shear velocity of a cylindrical volume of rock formation, centered on the wellbore. In one embodiment, the 3D images of non-linearity and Vp/V s velocity ratio can include images of non-linearity and Vp/V s velocity ratio extending out to a radius of investigation of several hundred meters, for example.

In one embodiment, the system may be considered as a suite of hardware and/or software or module sub-systems. Sub-systems 1900 to 1904 are used for survey planning and execution, downhole tool conveyance, coding and broadcast of the transmitted acoustic waves, and recording and detection of the non-linear signal. Sub-systems 1901 and 1905 are used for post-survey image processing. Sub-system 1900 is used for survey design, data design acquisition, control and recording. Subsystem 1901 is used for preliminary non-linearity and velocity imaging. Sub-system 1902 is used for receiver and sensor control and transmission. Sub-system 1902 is configured to emit controlled acoustic broadcasts and receive acoustic energy. Subsystem 1903 is used for generating a broadcast signal. Sub-system 1904 is used for non-linear signal detection. Sub-system 1905 is used for imaging non-linearity and imaging velocity.

In one embodiment, sub-system 1900 includes modeling for acquisition module 1907, data acquisition controller 1908, data pre-processing and enhancement module 1930 and data storage device 1921. In one embodiment, sub-system 1901 includes initial velocity model module 1906 and module 1931 for preliminary generation of three dimensional images including propagation compressional and shear velocity images, Vp/Vs velocity ratio images and images of non-linearity which are associated with amplitudes of the measured signal originating from the non-linear interaction at the mixing zone. In one embodiment, sub-system 1902 includes tool deployment and transport module 1909, tool mechanical controller 1910, azimuth and elevation controller 1917 for controlling the azimuth and elevation angles of the first acoustic source (S I), and azimuth and elevation controller 1918 for controlling the azimuth and elevation angles of the second acoustic source (S2). In one embodiment, sub-system 1903 includes coded signal generator 191 1, frequency multiplier and time delay module 1912, signal amplifier 1913 for generating the signal sent to first acoustic source (S I), signal amplifier 1914 for generating the signal sent to second acoustic source (S2). In one embodiment, sub-system 1904 includes receiver module(s) 1922 for receiving signals u3(t) from receivers Rl, R2, Rn, non-linear signal enhancement module(s) 1924 for enhancing the signal received by the signal receiver module(s) 1922, a template signal generator module 1927 for generating a template signal u_s(t), and a signal correlation module 1928 for correlating the received signal w^t) with the template signal u_s(t), as described in the above paragraphs. In one embodiment, sub-system 1905 includes data pre-processing and enhancement module 1932, velocity model iteration module 1933, imaging iteration module 1934, output image module 1936 for velocity ratio images and/or non-linearity images, and output velocity module 1935 for outputting determined values of velocities Vp, velocity Vs and/or Vp/V s velocity ratio.

As it can be appreciated, the term module is used herein to encompass a hardware device, a software program, or both. For example, image iteration module can be a piece of hardware that is configured to perform the iteration or a software program that can be run on a computer to perform the iteration, or includes both a piece of hardware and software application.

In operation, prior to data acquisition, logs of compressional and shear slowness and information about formation lateral continuity can be used to build an initial layered earth model extending laterally away from the well bore using module 1906. The initial velocity model from module 1906 is used by a forward modeling acquisition module 1907 to provide a plan for data acquisition that an operator uses to program data acquisition controller 1908 using operator inputs.

After sub-system 1902 is deployed into the borehole via tool deployment and transport module 1909 and clamped via tool mechanical controller 1910 to the borehole wall if needed, input commands can be sent to the coded signal generator sub-system 1903 where coded pulse sequences are generated by coded signal generator module 1911. The coded signals generated by the coded signal generator module 191 1 are adjusted in frequency and delayed by the frequency multiplier and time delay module 1912 such that the signal amplifiers 1913 and 1914 provide signals to transmitters or sources SI 1915 and S2 1916 to broadcast the signals in a delayed fashion so that the signals from S I and S2 arrive simultaneously at the mixing zone, as depicted, for example, on figure 12a.

The geometry of the broadcasts including the elevation and azimuth angles of the acoustic broadcasts from sources SI 1915 and S2 1916 is controlled by commands from the data acquisition controller 1908 which sends angle control commands to downhole azimuth and elevation angle controllers 1917 and 1918. All data pertaining to tool configuration and broadcast geometry collected from tool deployment 1909, tool mechanical controller 1910, azimuth and elevation controllers 1917 and 1918 indicated herein as record acquisition parameters 1919 are recorded data storage device 1921. Similarly, the coding scheme from coded signal generator sub-system 1903 indicated herein as record broadcast information at 1920 is also recorded at data storage device 1921.

In one embodiment, signals from signal modules 1922 recorded at the receiver or receivers 1923 (e.g., each of the receivers may have for example hydrophones, 3 component geophones or both Rl to Rn) can be processed by non-linear signal enhancement module 1924 to enhance the content of non-linear origin and reduce or substantially suppress signals of linear interaction origin or potential noise. Raw signals indicated at 1925 from signal receive modules 1922 and enhanced signals indicated at 1926 are stored in data storage device 1921. A template signal u_s(t) generated by template signal generator 1927, which can be derived from signals ui(t) and u₂(t) generated by signal generators 1912 and frequency multiplier and time delay module 1912, is correlated with received signals indicated by 1925 at 1928, as described for example in figures 12 to 14. The correlation of template signal u_s(t) and received or detected signal ¾(ί) is used to extract via correlated signal output module 1929 a correlated signal M(t, d, δ, z3). The correlated signal is also stored in data storage device 1921. In one embodiment, the data acquisition process implemented using components or modules 1909 to 1929 may be repeated with different beam geometry or at a multiple locations (z3) within the well. In one embodiment, data already recorded may be used to guide changes in acquisition parameters. For example, stored data within storage device 1921 may be further enhanced, for example, by hodogram analysis within preprocessing and hodogram processing module 1930. The data stored in storage device 1921 can be further used in conjunction with initial velocity model from initial velocity model module 1906 to create a set of images of non-linear properties and/or Vp/Vs ratios with imaging module 1931. These may be used to refine modeling by acquisition design module 1907 and/or acquisition parameters controlled by data acquisition controller 1908.

After the survey is complete, further data processing may be implemented via data pre-processing and enhancement module 1932. In one embodiment, hodogram analysis can be conducted to pre-condition the data for final analysis. The initial velocity model in module 1906 and images from module 1931 can be used as a starting point for further iteration of the velocity model through velocity model iteration module 1933 and imaging iteration of non-linear properties and/or Vp/Vs ratio through imaging iteration between velocity modeling module 1933 and imaging module 1934. The final result of the iteration is an optimized velocity model output through output module 1935 and images including Vp/Vs ratio images and/or nonlinear properties images output through output module 1936.

In one embodiment, the above implementation of the system is suitable for imaging out to relatively large distances from the wellbore (for example out to several hundred meters) and includes components designed to maximize signal to noise ratio and the detection of tenuous signals from a complex petrophysical, stratigraphic and structural context. Near wellbore applications probing a smaller volume with less variation in properties and stronger returned signals could potentially dispense with some aspects or portions or modules of the system concerned with signal detection and velocity model iteration. Likewise, in a less demanding acquisition environment, some aspects of the hardware could be simplified, for example clamped three- component geophones could be replaced with non-directional hydrophone receivers mounted on a centralized sonde. Furthermore, as it can be appreciated, although the system is described above as linking the acquisition data portion of the system to the imaging portion of the system, the acquisition of data portion can be accomplished separately from the imaging of Vp/Vs ratio and/or non-linear properties. Indeed, the acquired data can be accomplished by a first entity and the data stored in data storage device 1921. The acquired data in data storage device 1921 can then be transferred to a second entity, which can be the same or different from the first entity, the second entity can employ the imaging sub-system or imaging method described in the above paragraphs to obtain the Vp/Vs and/or non-linear properties images.

Furthermore, although each module is described in the above paragraph as having a specific functionality, as it can be appreciated any functionality in one or more modules can be moved to any other one or more modules. For example, some or all functionality in the coded signal generator sub-system 1903 can be moved to the non-linear signal detection subsystem 1904.

In addition, it must be appreciated that the term processor is used herein to encompass one or more processors. The one or more processors can be configured to implement the methods or portions of the methods described herein. The one or more processors can be located in one or more computers such as, for example, in a distributed computing environment. In some embodiments, programs for performing methods in accordance with embodiments of the invention can be embodied as program products in a computer such as a personal computer or server or in a distributed computing environment comprising a plurality of computers. Where reference is made to a processor that term should be understood to encompass any of these computing arrangements. The computer may include, for example, a desktop computer, a laptop computer, a handheld computing device. The computer program products may include a computer readable medium or storage medium or media having instructions stored thereon used to program a computer to perform the methods described above. Examples of suitable storage medium or media include any type of disk including floppy disks, optical disks, DVDs, CD ROMs, magnetic optical disks, RAMs, EPROMs, EEPROMs, magnetic or optical cards, hard disk, flash card (e.g., a USB flash card), PCMCIA memory card, smart card, or other media. Alternatively, a portion or the whole computer program product can be downloaded from a remote computer or server via a network such as the internet, an ATM network, a wide area network (WAN) or a local area network.

Although the invention has been described in detail for the purpose of illustration based on what is currently considered to be the most practical and preferred embodiments, it is to be understood that such detail is solely for that purpose and that the invention is not limited to the disclosed embodiments, but, on the contrary, is intended to cover modifications and equivalent arrangements that are within the spirit and scope of the appended claims. As a further example, it is to be understood that the present invention contemplates that, to the extent possible, one or more features of any embodiment can be combined with one or more features of any other embodiment.

Claims

WHAT IS CLAIMED IS:

1. A system for investigating rock formations outside a borehole, the system comprising:

a first acoustic source configured to generate a first acoustic compressional wave at a first frequency;

a second acoustic source configured to generate a second acoustic

compressional wave at a second frequency, wherein the first and the second acoustic sources are arranged within a localized area of the borehole and wherein the first and the second acoustic waves intersect in an intersection volume outside the borehole; and

a receiver arranged in the borehole and configured to receive a third shear acoustic wave at a third frequency, the third shear acoustic wave returning to the borehole due to a non-linear mixing process in a non-linear mixing zone within the intersection volume, wherein the third frequency is equal to a difference between the first frequency and the second frequency.

2. The system according to claim 1, wherein the first acoustic compressional wave and the second acoustic compressional wave are acoustic beams.

3. The system according to claim 1, wherein the first acoustic source, the second acoustic source or both are arranged to include a linear array of transducers.

4. The system according to claim 1, wherein the receiver is arranged to be clamped to a wall of the borehole and includes a three component receiver array.

5. The system according to claim 1, further comprising a controller arranged to steer in azimuth and elevation the first acoustic wave at the first frequency, the second acoustic wave at the second frequency, or both, in a given direction.

6. The system according to claim 1, further comprising a processor arranged to process machine-readable instructions, that when executed cause the processor to analyze the returning third acoustic wave and determine the location of the mixing area in the rock formation where non-linear interactions resulted in the returned energy in the borehole.

7. The system according to claim 1, wherein an axis of the first acoustic source is arranged at an angle relative to an axis of the second acoustic source.

8. The system according to claim 1, wherein the ratio of the second frequency to the first frequency is selected to maximize a mixing coefficient of the non-linear mixing process within the non-linear mixing zone.

9. The system according to claim 1, further comprising an encoder configured to encode the first acoustic wave, the second acoustic wave or both, with a time or frequency varying attribute or both.

10. The system according to claim 1, further comprising a storage device configured to record data corresponding to the received third shear acoustic wave, and a processor, arranged to process machine-readable instructions, that when executed cause the processor to estimate a compressional and shear velocity ratio according to selection rules governing the frequency ratio and convergence angle relationship, or non-linear properties of the rock formation around the borehole using the recorded data from the storage device and by repeating the generating of the first acoustic wave at the first frequency, the generating of the second acoustic wave at the second frequency, the receiving of the third acoustic wave returning to the borehole and identifying a location of the intersection volume based on an arrangement of the first and the second acoustic sources and a plurality of orientations of the first acoustic wave and the second acoustic wave.

1 1. The system according to claim 1 , further comprising a storage device configured to record data corresponding to the received third shear acoustic wave, and a processor, arranged to process machine-readable instructions, that when executed cause the processor to calculate a three-dimensional image of non-linear properties of the rock formation or a three-dimensional image of compressional and shear velocity ratio, using the recorded data from the storage device by repeating the generating of the first acoustic wave at the first frequency, the generating of the second acoustic wave at the second frequency, the receiving of the third acoustic wave returning to the borehole and identifying a location of the intersection volume based on an arrangement of the first and the second acoustic sources and a plurality of orientations of the first acoustic wave and the second acoustic wave.

12. A method for investigating rock formations outside a borehole, the method comprising:

generating a first compressional acoustic wave at a first frequency by a first acoustic source;

generating a second compressional acoustic wave at a second frequency by a second acoustic source, wherein the first and the second acoustic sources are arranged within a localized area of the borehole and wherein the first and the second acoustic waves intersect in an intersection volume outside the borehole; and

receiving a third shear acoustic wave at a third frequency, the third shear acoustic wave returning to the borehole due to a non-linear mixing process in a nonlinear mixing zone within the intersection volume at a receiver arranged in the borehole, wherein the third frequency is equal to a difference between the first frequency and the second frequency.

13. The method according to claim 12, wherein generating the first acoustic compressional wave and generating the second acoustic compressional wave comprises generating a first acoustic beam and generating a second acoustic beam.

14. The method according to claim 12, further comprising arranging the first acoustic source, the second acoustic source or both as a linear array of transducers.

15. The method according to claim 12, further comprising clamping the receiver to a wall of the borehole, the receiver including a three component receiver array.

16. The method according to claim 12, further comprising steering the acoustic wave at the first frequency, the acoustic wave at the second frequency, or both in a given direction by orienting in azimuth or in elevation, or both.

17. The method according to claim 12, further comprising analyzing the returning energy and determining the location of the mixing area in the rock formation where non-linear interactions resulted in the returned third acoustic wave in the borehole.

18. The method according to claim 12, further comprising arranging the axial orientation of the first acoustic source at an angle relative to an axial orientation of the second acoustic source.

19. The method according to claim 12, further comprising encoding, using an encoder, the first acoustic wave, the second acoustic wave or both with a time or frequency varying attribute, or both.

20. The method according to claim 12, further comprising calculating a location of the intersection zone based on an arrangement of the first and the second acoustic sources and a direction of the returned energy.

21. The method according to claim 20, further comprising estimating a compressional velocity to shear velocity ratio or non-linear properties of rock formation around the borehole or both, using data recorded by repeating the generating and receiving at a plurality of orientations of the first acoustic wave and the second acoustic wave.

22. The method according to claim 12, further comprising determining a mixing coefficient of the non-linear mixing process within the non-linear mixing zone and selecting a ratio of the second frequency and the first frequency so as to maximize the mixing coefficient.

23. A system for investigating rock formations outside a borehole using coded signals from acoustic sources, the system comprising:

a first acoustic source configured to generate a first acoustic signal comprising a first plurality of pulses arranged as a time sequence, the first plurality of pulses being separated in time, each pulse comprising a first modulated signal at a central frequency, wherein central frequencies of two consecutive pulses are different;

a second acoustic source configured to generate a second acoustic signal comprising a second plurality of pulses arranged as a time sequence, the second plurality of pulses being separated in time, wherein a separation in time between centers of two consecutive pulses is the same as a separation in time between centers of two corresponding pulses in the first plurality of pulses, wherein a start time difference is provided between a start time of a broadcast of the second plurality of pulses and a start time of a broadcast of the first plurality of pulses, wherein each pulse comprises a second modulated signal, and a central frequency of the second modulated signal within each pulse in the second plurality of pulses is a selected fraction of the central frequency of the first modulated signal within each

corresponding pulse in the first plurality of pulses;

wherein the first acoustic source and the second acoustic source are arranged within the borehole and are controllable such that trajectories of the first and the second acoustic signals intersect in an intersection volume outside the borehole; a receiver arranged within the borehole, the receiver being configured to detect a received signal returning to the borehole, the received signal including a signal being generated by a non-linear mixing process from the first acoustic signal and the second acoustic signal in a non-linear mixing zone within the intersection volume; and

a processor configured to perform data processing on the received signal to extract the signal generated by the non-linear mixing process over noise or over signals generated by a linear interaction process, or both.

24. The system according to claim 23, wherein a separation in time between centers of two consecutive pulses in the first plurality of pulses is greater than a time duration of each pulse.

25. The system according to claim 23, wherein the processor is configured to extract signals generated by the non-linear mixing process using waveform recognition.

26. The system according to claim 23, wherein the processor is configured to extract signals generated by the non-linear mixing process using band pass filtering based on forecast properties of the non-linear signal according to selection rules of non-linear mixing in the region where the non-linear mixing process occurs.

27. The system according to claim 23, wherein the processor is configured to extract signals generated by the non-linear mixing process by correlating the detected signal with a template signal designed according to forecast properties of the signal according to the selection rules of non-linear mixing in the region where the nonlinear mixing process occurs.

28. The system according to claim 23, wherein the signal generated by a nonlinear mixing process comprises a third plurality of pulses arranged in a time sequence and separated in time, wherein a separation in time between centers of two consecutive pulses is the same as the separation in time between centers of two consecutive pulses in the first plurality of pulses,

wherein each pulse in the third plurality of pulses comprises a third modulated signal having a third central frequency that is equal to a difference between the first central frequency and the second central frequency in each corresponding pulse of the first plurality of pulses and the second plurality of pulses,

wherein an arrival time at the receiver of each pulse of the third plurality of pulses is time delayed relative to a generation of a corresponding pulse of the first plurality of pulses by a total of the travel time from the first acoustic source to a center of the mixing zone and the travel time from the center of the mixing zone to the receiver.

29. The system according to claim 23, wherein the processor is configured to control the first source and the second source to repeat the generating of the first signal, the generating of the second signal, to control the receiver to repeat the detecting of the received signal returning to the borehole and repeat the performing of data processing on the received signal to extract the signal generated by the non-linear mixing process over noise or over signals generated by a linear interaction process, or both, for a range of start time differences and a range of frequency fractions d to obtain signals generated at a plurality of mixing zones.

30. The system according to claim 23, wherein the signal generated by non-linear interaction at the mixing zone is generated and detected only if the start time difference between the second plurality of pulses of the second acoustic source and the first plurality of pulses of the first acoustic source is equal to the difference in travel time between the travel time from the first acoustic source to the center of the mixing zone and the travel time between the second acoustic source to the center of the mixing zone.

31. The system according to claim 23, wherein increasing a number of the first plurality of pulses and increasing a number of the second plurality of pulses enhances discrimination of the signal generated by the non-linear mixing process over noise or over signals generated by a linear interaction process, or both.

32. The system according to claim 23, wherein the first signal comprises a sum of a plurality of pulses, wherein each pulse has a signal amplitude equal to the product of an envelope function and a modulated signal function.

33. The system according to claim 23, wherein the first signal ui(t) is in the form

-∑Eln(t - «)* exp (ί 2π* fi * (t—Tn ))* exp( ¾

n

where:

Eln(t-T_n) corresponds to an envelope of an nth pulse in the first plurality of pulses, the envelope having a positive or negative amplitude,

exp (i2jt* f_n* (t-T_n)) corresponds to the first modulated signal within the nth pulse in the first plurality of pulses,

f_n is a central frequency of the modulated signal within the nth pulse, t is a signal time, T_n is a time when the nth pulse is generated, and ζ_η is a phase of each pulse n and εχρ(ίζ_η) is a phase term of each pulse n.

34. The system according to claim 23, wherein the second signal comprises a sum of a plurality of pulses, wherein each pulse has a signal amplitude equal to the product of an envelope function and a modulated signal function.

35. The system according to claim 23, wherein the second signal ¾ (t-δ) is in the form: u₂ (t - S) =∑ E2«(t - (Tn + S))* exp(z 2π * d * /» * (t - (Tn + S))) * exp(

n where:

E2n(t-(T_n+5)) corresponds to an envelope of an nth pulse in the second plurality of pulses,

exp (i2jt* d* f_n* (t-(T_n+5))) corresponds to the second modulated signal within the nth pulse in the second plurality of pulses,

d * fn is a frequency of the modulated signal within the nth pulse in the second plurality of pulses,

d is a selected frequency ratio,

t is the signal time,

(T_n + δ) is a time when the nth pulse in the second plurality of pulses is broadcast,

δ is the start time difference between the start time of the broadcast of the second acoustic signal and the start time of the broadcast of the first acoustic signal, and ζ_η is a phase of each pulse n and εχρ(ίζ_η) is a phase term of each pulse n.

36. The system according to claim 23, wherein the signal generated by a nonlinear mixing process is proportional to a sum over a number of a third plurality of pulses of a product of an envelope function for each pulse by a modulated signal function for each pulse.

37. The system according to claim 23, wherein the signal ¾(ί) being generated by a non-linear mixing process is in the form:

ii3( -∑ E3n(t - (T_n + T)) * exp(J2;r * (l - d) * f_n * (t - (T_n + T))) * exp(i£) n where:

E3n(t - (T_n + T)) is an envelope of a pulse in a third plurality of pulses of the signal u3 generated by non-linear mixing of the first signal and the second signal,

exp(i2 * (l-d) * f_n* (t-(T_n+T))) corresponds to a third modulated signal within the nth pulse in the third plurality of pulses,

(l-d) * fn is a frequency of the modulated signal within the nth pulse in the third plurality of pulses which is a difference between the central frequency fn of the first modulated signal within the nth pulse of the first plurality of pulses and the central frequency d* fh of the second modulated signal within the nth pulse of the second plurality of pulses,

d is a selected frequency ratio,

t is the signal time,

(T_n + T) is a time when the nth pulse in the third plurality of pulses arrives at the receiver,

T is the total travel time from the first acoustic source to a center of the mixing zone and from the center of the mixing zone to the receiver, and ζ_η is a phase of each pulse n and εχρ(ίζ_η) is a phase term of each pulse n.

38. The system according to claim 23, wherein the processor is configured to extract signals generated by the non-linear mixing process by correlating the detected signal with a template signal wherein the template signal comprises a sum of a plurality of pulses wherein each pulse has a signal amplitude equal to the product of an envelope function and a modulated signal function.

39. The system according to claim 23, wherein the processor is configured to extract signals generated by the non-linear mixing process by correlating the detected signal with a template signal, wherein the template u_s is in the form: u,(t)

- T_n εχρ(ί2π * g (f_n) * (t - T_n )) * exp(i £)

n

where:

W_n(t-T_n) is an amplitude envelope of an nth pulse in the template signal, exp (ϊ2π* g(f_n) * (t-T_n)) corresponds to modulated signal within the nth pulse in the plurality of pulses of the template signal,

g(f_n) is a selected function of the modulated signal within the nth pulse in the plurality of pulses of the template,

t is a propagation time,

T_n is a time when the nth pulse in the plurality of pulses is simulated to arrive at the receiver, and ζ_η is a phase of each pulse n and εχρ(ίζ_η) is a phase term of each pulse n.

40. The system according to claim 23,

wherein the processor is configured to extract signals generated by the non- linear mixing process by correlating the detected signal with a template signal designed according to forecast properties of the signal according to the selection rules of non-linear interaction of acoustic waves, wherein the processor is further configured to repeat for a plurality of start time differences δ between the second plurality of pulses and the first plurality of pulses and for a plurality of frequency ratios d between the central frequency of the modulated signal within each pulse in the second plurality of pulses and the central frequency of the modulated signal with each corresponding pulse in the first plurality of pulses, and the correlation between the template signal and the detected signal to yield, for each start time difference δ and for each frequency ratio d, a correlated signal M(t, δ, d) which contains a non-linear interaction signal.

41. The system according to claim 40, wherein the processor is further configured to search all values of arrival time t, start time difference δ and frequency ratio d in the signal M(t, δ, d) to determine a value of arrival time TNL, start time difference 6NL and frequency ratio dm, for the presence of a correlated band-limited spike signal corresponding to the occurrence of a signal generated by the non-linear mixing process.

42. The system according to claim 41, wherein the determined start time difference 5NL, where the band-limited spike signal generated by the non-linear process occurs in the correlated signal M(t, δ, d), is substantially equal to a time difference between a travel time from the first acoustic source to a center of the mixing zone and a travel time from the second acoustic source to the center of the mixing zone.

43. The system according to claim 41, wherein the determined arrival time T L, where the band-limited spike signal generated by the non-linear mixing process occurs in the correlated signal M (t, δ, d), is substantially equal to a total travel time of the travel time from the first acoustic source to a center of the mixing zone and the travel time from the center of the mixing zone to the receiver.

44. A method for investigating rock formations outside a borehole, the method comprising: generating, by a first acoustic source, a first acoustic signal comprising a first plurality of pulses arranged as a time sequence, the first plurality of pulses being separated in time, each pulse comprising a first modulated signal at a central frequency, wherein central frequencies of two consecutive pulses are different;

generating, by a second acoustic source, a second acoustic signal comprising a second plurality of pulses arranged as a time sequence, the second plurality of pulses being separated in time, wherein a separation in time between centers of two consecutive pulses is the same as a separation in time between centers of two corresponding pulses in the first plurality of pulses, wherein a start time difference is provided between a start time of a broadcast of the second plurality of pulses and a start time of a broadcast of the first plurality of pulses, wherein each pulse comprises a second modulated signal and a central frequency of the second modulated signal within each pulse in the second plurality of pulses is a selected fraction of the central frequency of the first modulated signal within each corresponding pulse in the first plurality of pulses;

wherein the first acoustic source and the second acoustic source are arranged within the borehole and are controllable such that trajectories of the first and the second acoustic signals intersect in an intersection volume outside the borehole; receiving, by a receiver arranged within the borehole, a detected signal returning to the borehole, the detected signal including a signal being generated by a non-linear mixing process from the first acoustic signal and the second acoustic signal in a non-linear mixing zone within the intersection volume; and

performing, by a processor, data processing on the received signal to extract the signal generated by the non-linear mixing process over noise or over signals generated by a linear interaction process, or both.

45. The method according to claim 44, wherein generating the first acoustic signal comprises generating the first acoustic signal such that a separation in time between the centers of two consecutive pulses in the first plurality of pulses is greater than a time duration of each pulse.

46. The method according to claim 44, further comprising extracting, by the processor, signals generated by the non-linear mixing process using waveform recognition.

47. The method according to claim 44, further comprising extracting, by the processor, signals generated by the non-linear mixing process using band pass filtering based on forecast properties of the non-linear signal according to the selection rules of non-linear mixing in the region where the non-linear mixing process occurs.

48. The method according to claim 44, further comprising extracting, by the processor, signals generated by the non-linear mixing process by correlating the detected signal with a template signal designed according to forecast properties of the signal according to selection rules of non-linear mixing in the region where the non- linear mixing process occurs.

49. The method according to claim 44, wherein the signal generated by the nonlinear mixing process comprises a third plurality of pulses arranged in a time sequence and separated in time, wherein a separation in time between centers of two consecutive pulses is the same as the separation in time between centers of two consecutive pulses in the first plurality of pulses,

wherein each pulse in the third plurality of pulses comprises a third modulated signal having a third central frequency that is equal to a difference between the first central frequency and the second central frequency in each corresponding pulse of the first plurality of pulses and the second plurality of pulses,

wherein an arrival time at the receiver of each pulse of the third plurality of pulses is time delayed relative to a generation of a corresponding pulse of the first plurality of pulses by a total of the travel time from the first acoustic source to a center of the mixing zone and the travel time from the center of the mixing zone to the receiver.

50. The method according to claim 44, further comprising repeating the generating of the first signal, the generating of the second signal, to repeat the detecting by the receiver of the received signal returning to the borehole and repeating the performing of data processing on the received signal to extract the signal generated by the non-linear mixing process over noise or over signals generated by a linear interaction process, or both, for a range of start time differences and a range of frequency fractions d to obtain signals generated at a plurality of mixing zones.

51. The method according to claim 44, further comprising detecting the signal generated by non-linear interaction at the mixing zone only if the start time difference between the second plurality of pulses of the second acoustic source and the first plurality of pulses of the first acoustic source is equal to the difference in travel time between the travel time from the first acoustic source to the center of the mixing zone and the travel time between the second acoustic source to the center of the mixing zone.

52. The method according to claim 44, further comprising increasing a number of the first plurality of pulses and increasing a number of the second plurality of pulses to enhance discrimination of the signal generated by the non-linear mixing process over noise or over signals generated by a linear interaction process, or both.

53. The method according to claim 44, wherein generating the first signal comprises generating a first signal that is a sum of a plurality of pulses, wherein each pulse has a signal amplitude equal to the product of an envelope function and a modulated signal function.

54. The method according to claim 44, wherein generating the first signal comprises generating a first signal ui (t) is in the form

^M i( -∑Eln(t - «)* exp (ί 2π* fi * (t—Tn ))* exp( ¾

n where:

Eln(t-T_n) corresponds to an envelope of an nth pulse in the first plurality of pulses, the envelope having a positive or negative amplitude, exp (i2jt* f_n* (t-T_n)) corresponds to the first modulated signal within the nth pulse in the first plurality of pulses,

f_n is a central frequency of the modulated signal within the nth pulse,

t is a signal time,

T_n is a time when the nth pulse is generated, and ζ_η is a phase of each pulse n and εχρ(ίζ_η) is a phase term of each pulse n.

55. The method according to claim 44, wherein generating the second signal

comprises generating a second signal that is a sum of a plurality of pulses, wherein

each pulse has a signal amplitude equal to the product of an envelope function and a modulated signal function.

56. The method according to claim 44, wherein generating the second signal

comprises generating a second signal ¾ (t-δ) is in the form: u₂ (t - d) =∑ E2«(t - (Τη + δ)) * exp(z 2π * ά * β * (ί - (Τη + S))) * exp(

n where:

E2n(t-(T_n+5)) corresponds to an envelope of an nth pulse in the second plurality of pulses,

exp (ί2π* d* f_n* (t-(T_n+5))) corresponds to the second modulated signal within the nth pulse in the second plurality of pulses,

d * fn is a frequency of the modulated signal within the nth pulse in the second plurality of pulses,

d is a selected frequency ratio,

t is the signal time,

(T_n + δ) is a time when the nth pulse in the second plurality of pulses is broadcast, δ is the start time difference between the start time of the broadcast of the second acoustic signal and the start time of the broadcast of the first acoustic signal, and ζ_η is a phase of each pulse n and εχρ(ίζ_η) is a phase term of each pulse n.

57. The method according to claim 44, wherein the signal generated by the nonlinear mixing process comprises a sum of a plurality of pulses, wherein each pulse has a signal amplitude equal to the product of an envelope function and a modulated signal function.

58. The method according to claim 44, wherein the signal ¾(ί) being generated by a non-linear mixing process is in the form:

ii3( -∑ E3n(t - (T_n + T))* ε^χρ(ί2π * (l - d) * f_n * (t - (T_n + T))) * exp(i£) n where:

E3n(t - (T_n + T)) is an envelope of a pulse in a third plurality of pulses of the signal u3 generated by non-linear mixing of the first signal and the second signal,

exp(i2 * (l-d) * f_n* (t-(T_n+T))) corresponds to a third modulated signal within the nth pulse in the third plurality of pulses,

(l-d) * fh is a frequency of the modulated signal within the nth pulse in the third plurality of pulses which is a difference between the central frequency fn of the first modulated signal within the nth pulse of the first plurality of pulses and the central frequency d* fh of the second modulated signal within the nth pulse of the second plurality of pulses,

d is a selected frequency ratio,

t is the signal time, (T_n + T) is a time when the nth pulse in the third plurality of pulses arrives at the receiver,

T is the total travel time from the first acoustic source to a center of the mixing zone and from the center of the mixing zone to the receiver, and ζ_η is a phase of each pulse n and εχρ(ίζ_η) is a phase term of each pulse n.

59. The method according to claim 44, further comprising extracting, by the processor, signals generated by the non-linear mixing process by correlating the detected signal with a template signal wherein the template signal comprises a sum of a plurality of pulses, wherein each pulse has a signal amplitude equal to the product of an envelope function and a modulated signal function.

60. The method according to claim 44, further comprising extracting, by the processor signals generated by the non-linear mixing process by correlating the detected signal with a template signal, wherein the template u_s is in the form: u,(t) = X Wn{t - T_n εχρ(ί2π * g (f_n) * (t - T_n )) * exp(^)

n

Where:

W_n(t-T_n) is an amplitude envelope of an nth pulse in the template signal, exp (ί2π* g(f_n) * (t-T_n)) corresponds to modulated signal within the nth pulse in the plurality of pulses of the template signal,

g(f_n) is a selected function of the modulated signal within the nth pulse in the plurality of pulses of the template,

t is a propagation time,

T_n is a time when the nth pulse in the plurality of pulses is simulated to arrive at the receiver, and ζ_η is a phase of each pulse n and εχρ(ίζ_η) is a phase term of each pulse n.

61. The method according to claim 44, further comprising:

extracting, by the processor, signals generated by the non-linear mixing process, by correlating the detected signal with a template signal designed according to forecast properties of the signal according to the selection rules of non-linear interaction of acoustic waves, and repeating for a plurality of start time differences δ between the second plurality of pulses and the first plurality of pulses and for a plurality of frequency ratios d between the central frequency of the modulated signal within each pulse in the second plurality of pulses and the central frequency of the modulated signal with each corresponding pulse in the first plurality of pulses, the correlating between the template signal and the detected signal to yield, for each start time difference δ and for each frequency ratio d, a correlated signal M(t, δ, d) which contains a non-linear interaction signal.

62. The method according to claim 61, further comprising searching, by the processor, all values of arrival time t, start time difference δ and frequency ratio d in the signal M(t, δ, d) to determine a value of arrival time T_NL, start time difference 6NL and frequency ratio dm, for the presence of a correlated band-limited spike signal corresponding to the occurrence of a signal generated by the non-linear mixing process.

63. The method according to claim 62, wherein the determined start time difference 5NL where the band-limited spike signal generated by the non-linear mixing process occurs in the correlated signal M(t, δ, d) is substantially equal to a time difference between a travel time from the first acoustic source to a center of the mixing zone and a travel time from the second acoustic source to the center of the mixing zone.

64. The method according to claim 62, wherein the determined arrival time T L where the band-limited spike signal generated by the non-linear mixing process occurs in the correlated signal M(t, δ, d) is substantially equal to a total travel time of the travel time from the first acoustic source to a center of the mixing zone and the travel time from the center of the mixing zone to the receiver.

65. A system for investigating rock formation outside a borehole, the system comprising:

a first acoustic source configured to generate a first acoustic signal at a first frequency;

a second acoustic source configured to generate a second acoustic signal at a second frequency;

wherein the first acoustic source and the second acoustic source are arranged and positioned within the borehole and are controllable such that the first and the second acoustic signals intersect in an intersection volume outside the borehole and a start time difference is provided between the second acoustic signal and the first acoustic signal;

a receiver arranged within the borehole, the receiver being configured to receive a detected signal returning to the borehole having a frequency equal to a difference between the first frequency and the second frequency, the detected signal being generated by a non-linear mixing process from the first acoustic signal and the second acoustic signal in a non-linear mixing zone within the intersection volume; and

a recording system configured to record the detected signal and store the detected signal in a storage device and to record measurement parameters including a position of the first acoustic source, a position of the second acoustic source, a position of the receiver, elevation angle and azimuth angle of the first acoustic signal and elevation angle and azimuth angle of the second acoustic signal.

66. The system according to claim 65, further comprising a processor configured to scan a range of start time differences between the start time of a broadcast of the first acoustic signal and the start time of a broadcast of the second acoustic signal so that the range of start time differences includes a time difference for which the first acoustic signal and the second acoustic signal arrive substantially simultaneously at the mixing zone.

67. The system according to claim 65, further comprising hardware or software or both, configured to control a phase of the first acoustic signal or a phase of the second acoustic signal, or both, such that received signals from the interaction of the first acoustic signal and the second acoustic signal have opposite polarity.

68. The system according to claim 67, wherein the received signals from the interaction of the first acoustic signal and the second acoustic signal having opposite polarities are combined to enhance the signal generated by the non-linear interaction over noise or signals originating from a linear interaction or both.

69. The system according to claim 65, wherein at least a portion of the detected signal generated by the non-linear interaction is proportional to a product of amplitudes of the first acoustic signal and the second acoustic signal.

70. The system according to claim 65, further comprising a time variant band-pass filtering device configured to apply a time-variant frequency band-pass filter to the detected signal so as to keep a narrow frequency bandwidth around an expected frequency bandwidth of the signal generated by non-linear interaction, the narrow frequency bandwidth being selected around the difference between the first frequency and the second frequency.

71. The system according to claim 65, wherein:

the first acoustic signal comprises a first plurality of pulses arranged as a time sequence, the first plurality of pulses being separated in time, each pulse comprising a first modulated signal at a central frequency, wherein central frequencies of two consecutive pulses are different;

the second acoustic signal comprises a second plurality of pulses arranged as a time sequence, the second plurality of pulses being separated in time, wherein a separation in time between centers of two consecutive pulses is the same as a separation in time between centers of two corresponding pulses in the first plurality of pulses, wherein a start time difference is provided between a start time of a broadcast of the second plurality of pulses and a start time of a broadcast of the first plurality of pulses, wherein each pulse comprises a second modulated signal, and wherein a central frequency of the second modulated signal within each pulse in the second plurality of pulses is a selected fraction of the central frequency of the first modulated signal within each corresponding pulse in the first plurality of pulses; and

the start time difference is controlled such that the first acoustic signal and the second acoustic signal intersect in the mixing zone.

72. The system according to claim 71, wherein a separation in time between centers of two consecutive pulses in the first plurality of pulses is greater than a time duration of each pulse.

73. The system according to claim 71, further comprising a processor configured to correlate the detected signal with a template signal comprising a plurality of pulses, the plurality of pulses being arranged in a time sequence and separated in time, wherein a separation in time between centers of two consecutive pulses is the same as the separation in time between centers of two consecutive pulses in the first plurality of pulses, wherein each pulse in the plurality of pulses comprises a modulated signal having a central frequency that is equal to a difference between a central frequency of the first modulated signal and a central frequency of the second modulated signal in each corresponding pulse of the first plurality of pulses and the second plurality of pulses, to obtain a correlated signal containing the signal generated by non-linear interaction in the non-linear mixing zone.

74. The system according claim 73, further comprising a plurality of receivers configured to receive a plurality of detected signals, wherein the processor is configured to apply multi-dimensional filtering, or time moveout analysis and stacking or both, to the plurality of detected signals that are received by the plurality of receivers.

75. The system of claim 73, wherein the receiver comprises a hydrophone or multi-component detectors, or both, and wherein the detected signal comprises a pressure signal detected by the hydrophone or multi-component signals detected by the multi-component detectors or both.

76. The system of claim 75, wherein the processor is further configured to perform hodogram analysis on the detected signal.

77. The system of claim 76, wherein the processor is further configured to apply a hodogram analysis to the three detected component signals corresponding to vertical shear SV, horizontal shear SH and compressional P signals.

78. The system of claim 77, wherein the processor is further configured to apply multi-dimensional filtering, or time moveout analysis and stacking, or both, to the three detected component signals that are received by the three detectors, or to apply multi-dimensional filtering, or time moveout analysis and stacking, or both, to the pressure signal detected by the hydrophone.

79. The system of claim 77, further comprising a time variant band-pass filtering device configured to apply a time-variant frequency band pass filter to the detected signal so as to keep a narrow frequency bandwidth around an expected frequency bandwidth of the signal generated by non-linear interaction, the narrow frequency bandwidth being selected around the difference between the first frequency and the second frequency.

80. The system of claim 79, wherein the processor is further configured to control the first acoustic source, the second acoustic source and the receiver, and repeat a plurality of times the generating of the first acoustic signal, the generating of the second acoustic signal, the receiving of the detected signal, the applying of the time variant frequency band pass filter, the correlating of the detected signal with the template signal, the applying of the hodogram analysis, and the combining of signals to improve signal to noise ratio.

81. The system of claim 79, wherein the processor is further configured to vary the start time difference and repeat a plurality of times the generating of the first acoustic signal, the generating of the second acoustic signal, the receiving of the detected signal, the applying of the time variant frequency band pass filter, the correlating of the detected signal with the template signal, and the applying of the hodogram analysis for each of three detected signals to separate different propagation signal modes including compressional mode P, vertical shear mode SV and horizontal shear mode SH.

82. The system of claim 79, wherein the processor is further configured to vary the start time difference and repeat a plurality of times the generating of the first signal, the generating of the second signal, the receiving of the detected signal, the applying of the time variant frequency band pass filter, and the correlating of the signal detected by a hydrophone with the template signal.

83. The system according to claim 79, wherein the processor is further configured to control the first acoustic source, the second acoustic source and the receiver to scan a plurality of azimuth and elevation angles of the first acoustic signal, a plurality of azimuth and elevation angles of the second acoustic signal, or a position of the receiver, and repeating a plurality of times the generating of the first acoustic signal, the generating of the second acoustic signal, the receiving of the detected signal, the applying of the time-variant frequency band pass filter to obtain a filtered signal, the correlating of the filtered signal with the template signal, and the applying of the hodogram analysis on the three detected component signals to separate the component signals to different propagation signal modes including compressional mode P, vertical shear mode SV and horizontal shear mode SH.

84. The system according to claim 83, wherein the processor is further configured to repeat varying the start time difference between the start time of a broadcast of the first acoustic signal and the start time of a broadcast of the second acoustic signal, to repeat varying the frequency ratio between the first frequency and the second frequency, and to repeat scanning the plurality of azimuth and elevation angles of the first acoustic signal, the plurality of azimuth and elevation angles of the second acoustic signal, and the position of the receiver to obtain a plurality of signal measurements.

85. The system according to claim 84, further comprising a storage device in communication with the processor, the storage device being configured to store the obtained plurality of signal measurements.

86. A method of investigating rock formation outside a borehole comprising: generating, by a first acoustic source, a first signal at a first frequency;

generating, by a second acoustic source, a second signal at a second frequency;

wherein the first acoustic source and the second acoustic source are arranged and positioned within the borehole and are controllable such that the first and the second acoustic signals intersect in an intersection volume outside the borehole and a start time difference is provided between the second acoustic signal and the first acoustic signal;

receiving, by a receiver, a detected signal returning to the borehole having a frequency equal to a difference between the first frequency and the second frequency, the detected signal being generated by a non-linear mixing process from the first acoustic signal and the second acoustic signal in a non-linear mixing zone within the intersection volume; and

recording the detected signal and storing the detected signal in a storage device and recording measurement parameters including a position of the first acoustic source, a position of the second acoustic source, a position of the receiver, elevation angle and azimuth angle of the first acoustic signal and elevation angle and azimuth angle of the second acoustic signal.

87. The method according to claim 86, further comprising scanning a range of start time differences between the start time of a broadcast of the first acoustic signal and the start time of a broadcast of the second acoustic signal so that the range of start time differences includes a time difference for which the first acoustic signal and the second acoustic signal arrive substantially simultaneously at the center of the mixing zone.

88. The method according to claim 86, further comprising controlling a phase of the first acoustic signal or a phase of the second acoustic signal, or both such that received signals from the interaction of the first acoustic signal and the second acoustic signal have opposite polarity.

89. The method according to claim 88, further comprising combining received signals from the interaction of the first acoustic signal and the second acoustic signal having opposite polarities to enhance the signal generated by the non-linear interaction over noise or signals originating from linear interaction, or both.

90. The method according to claim 86, wherein at least a portion of the detected signal generated by the non-linear interaction is proportional to a product of amplitudes of the first acoustic signal and the second acoustic signal.

91. The method according to claim 86, further comprising applying, by a time variant band-pass filter, a time-variant frequency band-pass filter to the detected signal so as to keep a narrow frequency bandwidth around an expected frequency bandwidth of the signal generated by non-linear interaction, the narrow frequency bandwidth being selected around the difference between the first frequency and the second frequency.

92. The method according to claim 86, wherein:

generating the first acoustic signal comprises generating a first plurality of pulses arranged as a time sequence, the first plurality of pulses being separated in time, each pulse comprising a first modulated signal at a central frequency, wherein central frequencies of two consecutive pulses are different;

generating the second acoustic signal comprises generating a second plurality of pulses arranged as a time sequence, the second plurality of pulses being separated in time, wherein a separation in time between centers of two consecutive pulses is the same as a separation in time between centers of two corresponding pulses in the first plurality of pulses, wherein a start time difference is provided between a start time of a broadcast of the second plurality of pulses and a start time of a broadcast of the first plurality of pulses, wherein each pulse comprises a second modulated signal, and wherein a central frequency of the second modulated signal within each pulse in the second plurality of pulses is a selected fraction of the central frequency of the first modulated signal within each corresponding pulse in the first plurality of pulses; and the start time difference is controlled such that the first acoustic signal and the second acoustic signal intersect in the mixing zone .

93. The method according to claim 92, wherein a separation in time between centers of two consecutive pulses in the first plurality of pulses is greater than a time duration of each pulse.

94. The method according to claim 92, further comprising

correlating, by a processor, the detected signal with a template signal comprising a plurality of pulses, the plurality of pulses being arranged in a time sequence and separated in time, wherein a separation in time between centers of two consecutive pulses is the same as the separation in time between centers of two consecutive pulses in the first plurality of pulses, wherein each pulse in the plurality of pulses comprises a modulated signal having a central frequency that is equal to a difference between a central frequency of the first modulated signal and a central frequency of the second modulated signal in each corresponding pulse of the first plurality of pulses and the second plurality of pulses, to obtain a correlated signal containing the signal generated by non-linear interaction in the non-linear mixing zone.

95. The method according claim 94, further comprising receiving a plurality of detected signals at a plurality of receivers and applying, by the processor, multidimensional filtering, or time moveout analysis and stacking, or both, to the plurality of detected signals that are received by the plurality of receivers.

96. The method according to claim 94, further comprising detecting a pressure signal by a hydrophone or detecting multi-component signals by multi-component detectors, or both.

97. The method according to claim 94, further comprising detecting each of three components signals of the received signal by each of three detectors of the receiver.

98. The method according to claim 97, further comprising performing, by the processor, hodogram analysis on the detected signal.

99. The method according to claim 98, wherein performing hodogram analysis comprises applying a hodogram analysis to the three detected component signals corresponding to vertical shear SV, horizontal shear SH and compressional P signals.

100. The method according to claim 99, further comprising applying multidimensional filtering, or time moveout analysis and stacking, or both, to the three detected component signals that are received by the three detectors, or applying multidimensional filtering, or time moveout analysis and stacking, or both, to the pressure signal detected by the hydrophone.

101. The method according to claim 99, applying a time variant frequency band pass filter to the detected signal so as to keep a narrow frequency bandwidth around an expected frequency bandwidth of the signal generated by non-linear interaction, the narrow frequency bandwidth being selected around the difference between the first frequency and the second frequency.

102. The method according claim 101, further comprising repeating a plurality of times the generating of the first acoustic signal, the generating of the second acoustic signal, the receiving of the detected signal, the applying of the time variant frequency band pass filter, the correlating of the detected signal with the template signal, the applying of the hodogram analysis, and the combining of signals to improve signal to noise ratio.

103. The method according claim 101 further comprising varying the start time difference and repeating a plurality of times the generating of the first acoustic signal, the generating of the second acoustic signal, the receiving of the detected signal, the applying of the time variant frequency band pass filter, the correlating of the detected signal with the template signal, and the applying of the hodogram analysis for each of three detected signals to separate different propagation signal modes including compressional mode P, vertical shear mode SV and horizontal shear mode SH.

104. The method according to claim 103, varying the start time difference and repeating a plurality of times the generating of the first signal, the generating of the second signal, the receiving of the detected signal , and the applying of the time variant frequency band pass filter, and the correlating of the signal detected by the hydrophone with the template signal.

105. The method of claim 97, further comprising scanning a plurality of azimuth and elevation angles of the first acoustic signal, a plurality of azimuth and elevation angles of the second acoustic signal, or a position of the receiver, and repeating a plurality of times the generating of the first acoustic signal, the generating of the second acoustic signal, the receiving of the detected signal, the applying of the time variant frequency band pass filter to obtain a filtered signal, the correlating of the filtered signal with the template signal, and the applying of the hodogram analysis on the three detected component signals to separate the component signals to different propagation signal modes including compressional mode P, vertical shear mode SV and horizontal shear mode SH.

106. The method according to claim 105, further comprising repeating the varying the start time difference between the start time of a broadcast of the first acoustic signal and the start time of a broadcast of the second acoustic signal, and repeating the varying the frequency ratio between the first frequency and the second frequency, and repeating the scanning the plurality of azimuth and elevation angles of the first acoustic signal, the plurality of azimuth and elevation angles of the second acoustic signal, and the position of the receiver to obtain a plurality of signal measurements.

107. The method according to claim 106, further comprising storing in a storage device in communication with the processor, the obtained plurality of signal measurements.

108. A method of generating a three-dimensional image of a rock formation comprising:

generating a first acoustic signal by a first acoustic source, wherein the first acoustic signal comprises a first plurality of pulses arranged in a time sequence, the first plurality of pulses being separated in time, each pulse comprising a first modulated signal at a central frequency, wherein central frequencies of two consecutive pulses are different;

generating a second acoustic signal by a second acoustic source, wherein the second acoustic signal comprises a second plurality of pulses arranged as a time sequence, the second plurality of pulses being separated in time, wherein a separation in time between centers of two consecutive pulses is the same as a separation in time between centers of two corresponding pulses in the first plurality of pulses, wherein a start time difference is provided between a start time of a broadcast of the second plurality of pulses and a start time of a broadcast of the first plurality of pulses, wherein each pulse comprises a modulated signal and wherein a central frequency of the modulated signal within each pulse in the second plurality of pulses is a selected fraction of the central frequency of the modulated signal within each corresponding pulse in the first plurality of pulses; wherein the first acoustic source and the second acoustic source are arranged within the borehole and are controllable such that trajectories of the first and the second acoustic signals intersect in an intersection volume outside the borehole;

receiving, by a receiver arranged within the borehole, a detected signal returning to the borehole, the detected signal including a signal being generated by a non-linear mixing process from the first acoustic signal and the second acoustic signal in a non-linear mixing zone within the intersection volume;

performing, by a processor, data processing on the received signal to extract the signal generated by the non-linear mixing process over noise or over signals generated by a linear interaction process, or both; and

generating a three dimensional image of a propagation compressional velocity, a propagation of shear velocity, a ratio of compressional velocity and shear velocity, or non-linear properties of a rock formation, or any combination of two or more thereof based on the signal generated by the non-linear mixing process.

109. The method according to claim 108, wherein generating the first acoustic signal comprises generating the plurality of pulses such that a separation in time between centers of two consecutive pulses in the first plurality of pulses is greater than a time duration of each pulse.

1 10. The method according to claim 108, wherein generating the three dimensional image of the propagation compressional velocity, the propagation of shear velocity, the ratio of compressional velocity and shear velocity, or the non-linear properties of the rock formation, or any combination of two or more thereof is performed by using Kirchhoff imaging, beam imaging or wave equation imaging.

1 11. The method according to claim 110, further comprising determining a value of the propagation compressional velocity or a value of the shear velocity or both using tomographic velocity inversion or full wave form inversion or by iterative imaging in combination with tomographic velocity inversion.

1 12. The method according to claim 108, wherein generating the three dimensional image of the propagation compressional velocity, the propagation of shear velocity, the ratio of compressional velocity to shear velocity, or the non-linear properties of a rock formation, or any combination of two or more thereof is performed using an analysis or mapping method comprising:

estimating an initial propagating velocity model using well logs in a borehole and assumptions about lateral continuity of a rock formation away from the borehole; calculating a first travel time of the first acoustic signal from the first acoustic source to a center of the mixing zone;

calculating a second travel time of the second acoustic signal from the second acoustic source to the center of the mixing zone;

calculating a third travel time between the center of the mixing zone and a receiver configured to receive a detected signal returning to the borehole having a frequency equal to a difference between the first frequency and the second frequency, the detected signal being generated by a non-linear mixing process from the first acoustic signal and the second acoustic signal in the mixing zone;

calculating an arrival time Tp corresponding to a total travel time of a signal returning to the borehole by adding the first travel time and the third travel time, and calculating a start time difference δρ between the first travel time and the second travel time;

extracting from the signal received at the receiver signals generated by the non-linear mixing process to yield for each start time difference δ and for each frequency ratio d, a correlated signal M(t, δ, d), which is a function of arrival time t, the start time difference δ and the frequency ratio d, and which contains non-linear interaction signals;

searching all values of arrival time t, start time difference δ, and frequency ratio d in the correlated signal M(t, δ, d) to determine a value of arrival time TNL, start time difference 6NL and frequency ratio dNL for which a band-limited spike signal generated by the non-linear mixing process occurs; comparing the calculated arrival time Tp with a measured arrival time T L and comparing the calculated start time difference δρ with a measured start time difference 6NL_;

determining whether the calculated arrival time Tp is different from the measured arrival time T L and the calculated start time difference δρ is different from the measured start time difference 6NL and, if so, updating the propagating velocity model using a difference between the measured and the calculated arrival times and a difference between the calculated and the measured start time differences; and

repeating iteratively the above steps until the calculated arrival time corresponding to the total travel time is substantially equal to a measured arrival time, or the calculated start time difference is substantially equal to a measured start time difference, or both, for a range of values of the position of the first acoustic source, the position of the second acoustic source, the position of the receiver, the azimuth angle of the first acoustic signal, the azimuth angle of the second acoustic signal, the elevation angle of the first acoustic signal, or the elevation angle of the second acoustic signal, or any combination of two or more thereof, to obtain a propagating velocity Vp and propagating velocity Vs that results in the calculated arrival time being substantially equal to the measured arrival time or the calculated time difference being substantially equal to the measured start time difference, or both; and generating a three dimensional image of a propagation compressional velocity, a propagation of shear velocity, a ratio of compressional velocity to shear velocity, or non-linear properties, or any combination of two or more thereof.

1 13. The method according to claim 108, wherein performing, by a processor, data processing on the received signal to extract the signal generated by the non-linear mixing process comprises extracting signals generated by the non-linear mixing process by correlating the detected signal with a template signal designed according to forecast properties of the signal according to the selection rules for non-collinear mixing in a non-linear medium, and repeating for a plurality of start time differences δ between the second plurality of pulses and the first plurality of pulses and for a plurality of frequency ratios d between the central frequency of the modulated signal within each pulse in the second plurality of pulses and the central frequency of the modulated signal within each corresponding pulse in the first plurality of pulses, the correlating between the template signal and the detected signal to yield for each start time difference δ and each frequency ratio d a correlated signal M(t, δ, d) which contains a non-linear interaction signal.

1 14. The method according to claim 112, further comprising calculating, based on a propagation velocity model, trajectories of the first acoustic signal emitted by a first acoustic source and the second acoustic signal emitted by the second acoustic source using the position of the first acoustic source, the elevation angle and the azimuth angle of the first acoustic signal and the position of the second acoustic source, the elevation angle and the azimuth angle of the second acoustic signal.

1 15. The method according to claim 114, further comprising determining a location of a mixing zone where the first acoustic signal and the second acoustic signal interact non-linearly based on the position of the first acoustic source, the elevation angle and the azimuth angle of the first acoustic signal and the position of the second acoustic source, the elevation angle and the azimuth angle of the second acoustic signal.

1 16. The method according to claim 1 15, further comprising calculating spatial coordinates of the mixing zone where the two acoustic signals interact non-linearly, calculating a convergence angle between the first acoustic signal from the first acoustic source and the second acoustic signal from the second acoustic source based on the elevation angle and the azimuth angle of the first acoustic signal, and the position of the first acoustic source, and the elevation angle and the azimuth angle of the second acoustic signal, and a position of the second acoustic source.

1 17. The method according to claim 116, further comprising mapping the correlated signal M(TNL, 5NL, d L) to the spatial coordinates of the mixing zone corresponding to the non-linear properties of the rock formation at the mixing zone.

1 18. The method according to claim 1 17, further comprising calculating the velocity ratio between the compressional velocity and the shear velocity at the mixing zone from the frequency ratio dNL and the convergence angle using a selection rule relationship between velocity ratio of compressional velocity and shear velocity, frequency ratio and angle of convergence and mapping the calculated velocity ratio to the spatial coordinates of the mixing zone.

1 19. The method of claim 1 18, further comprising repeating mapping the correlated signal M(T_NL, 6NL, CINL) to the spatial coordinates of the mixing zone, or mapping the calculated velocity ratio to the spatial coordinates of the mixing zone, or both, for a range of values of the position of the first acoustic source, the position of the second acoustic source, the position of the receiver, the azimuth angle of the first acoustic signal, the azimuth angle of the second acoustic signal, the elevation angle of the first acoustic signal, or the elevation angle of the second acoustic signal, or any combination of two or more thereof, to obtain a three dimensional image of a strength of the non-linear properties or a three dimensional image of the velocity ratio, or both, by combining the measured values through repeated measurements.

120. The method according to claim 108, further comprising determining a value of the propagation compressional velocity, a value of the propagation of shear velocity, or a value of the ratio of compressional velocity and shear velocity, or any combination of two or more thereof.

121. The method according to claim 108, wherein determining the value of propagation compressional velocity, the value of propagation of shear velocity, the value of ratio of compressional velocity and shear velocity, or any combination of two or more thereof is performed using a mapping method comprising:

estimating an initial propagating velocity model using well logs in a borehole and assumptions about lateral continuity of a the rock formation away from the borehole; calculating a first travel time of the first acoustic signal from the first acoustic source to a center of the mixing zone;

calculating a second travel time of the second acoustic signal from the second acoustic source to the center of the mixing zone;

calculating a third travel time between the center of the mixing zone and a receiver configured to receive a detected signal returning to the borehole having a frequency equal to a difference between the first frequency and the second frequency, the detected signal being generated by a non-linear mixing process from the first acoustic signal and the second acoustic signal in a mixing zone;

calculating an arrival time Tp corresponding to a total travel time of a signal returning to the borehole by adding the first travel time and the third travel time, and calculating a start time difference δρ between the first travel time and the second travel time;

extracting from the signal received at the receiver, signals generated by the non-linear mixing process to yield for each start time difference δ and for each frequency ratio d, a correlated signal M(t, δ, d), which is a function of arrival time t, the start time difference δ and the frequency ratio d, and which contains non-linear interaction signals;

searching all values of arrival time t, start time difference δ, and frequency ratio d in the correlated signal M(t, δ, d) to determine a value of arrival time TNL, start time difference 6NL and frequency ratio dNL for which a signal generated by the nonlinear mixing process occurs;

comparing the calculated arrival time Tp with a measured arrival time T L and comparing the calculated start time difference δρ with a measured start time difference 6NL;

determining whether the calculated arrival time Tp is different from the measured arrival time T L and the calculated time difference δρ is different from the measured start time difference 6NL and, if so, updating the propagating velocity model using a difference between the measured and the calculated arrival times and a difference between the calculated and the measured start time differences; and repeating iteratively the above steps until the calculated arrival time corresponding to the total travel time is substantially equal to a measured arrival time, or the calculated time difference is substantially equal to a measured start time difference, or both, for a range of values of the position of the first acoustic source, the position of the second acoustic source, the position of the receiver, the azimuth angle of the first acoustic signal, the azimuth angle of the second acoustic signal, the elevation angle of the first acoustic signal, or the elevation angle of the second acoustic signal, or any combination of two or more thereof, to obtain a propagating velocity Vp and propagating velocity Vs that result in the calculated arrival time being substantially equal to the measured arrival time or the calculated start time difference being substantially equal to the measured start time difference, or both.

122. The method according to claim 121, further comprising calculating based on a propagation velocity model trajectories of the first acoustic signal having the first frequency emitted by the first acoustic source and the second acoustic signal having the second frequency emitted by the second acoustic source using a position of the first acoustic source, the elevation angle and the azimuth angle of the first acoustic signal and the position of the second acoustic source, the elevation angle and the azimuth angle of the second acoustic signal, and the position of the receiver.

123. The method according to claim 122, further comprising determining a location of a mixing zone where the first acoustic signal and the second acoustic signal interact non-linearly based on the position of the first acoustic source, the elevation angle and the azimuth angle of the first acoustic signal and the position of the second acoustic source, the elevation angle and the azimuth angle of the second acoustic signal.

124. The method according to claim 123, further comprising calculating spatial coordinates of the mixing zone where the two acoustic signals interact non-linearly, calculating a convergence angle between the first acoustic signal from the first acoustic source and the second acoustic signal from the second acoustic source based on the elevation angle and the azimuth angle of the first acoustic signal and the position of the first acoustic source, and the elevation angle and the azimuth angle of the second acoustic signal and the position of the second acoustic source.

125. The method according to claim 124, further comprising calculating the velocity ratio between the compressional velocity and the shear velocity at the mixing zone from the frequency ratio dNL and the convergence angle using a selection rule relationship between velocity ratio of compressional velocity and shear velocity, frequency ratio and angle of convergence and mapping the calculated velocity ratio to the spatial coordinates of the mixing zone.

126. The method of claim 125, further comprising repeating mapping the calculated velocity ratio to the spatial coordinate of the mixing zone for a range of values of the position of the first acoustic source, the position of the second acoustic source, the position of the receiver, the azimuth angle of the first acoustic signal, the azimuth angle of the second acoustic signal, the elevation angle of the first acoustic signal, or the elevation angle of the second acoustic signal, or any combination of two or more thereof, to obtain the value of the velocity ratio for a volume of the rock formation around the borehole.

127. A method of generating a three-dimensional image of non-linear properties of a rock formation, the method comprising:

reading, by a processor, measurement parameters including a position of a first acoustic source configured to generate a first acoustic signal at a first frequency, a position of a second acoustic source configured to generate a second acoustic signal at a second frequency, a position of a receiver configured to receive a detected signal returning to a borehole from a mixing zone within the rock formation at the intersection of the first acoustic signal and the second acoustic signal, and elevation and azimuth angles of the first and second acoustic signals;

reading, by the processor, a correlated signal M(t, δ, d), which is a function of arrival time t, start time difference δ and frequency ratio d between the first frequency and the second frequency, the correlated signal M(t, δ, d) containing signals generated by a non-linear mixing process from a first acoustic signal and a second acoustic signal at a mixing zone within the rock formation;

reading, by the processor, an initial propagation compressional and shear velocity model;

calculating, by the processor, a first travel time of the first acoustic signal from the first acoustic source to a center of the mixing zone;

calculating, by the processor, a second travel time of the second acoustic signal from the second acoustic source to the center of the mixing zone;

calculating, by the processor, a third travel time between the center of the mixing zone and the receiver configured to receive the detected signal returning to the borehole having a frequency equal to a difference between the first frequency and the second frequency, the detected signal being generated by a non-linear mixing process from the first acoustic signal and the second acoustic signal in the mixing zone;

calculating, by the processor, an arrival time Tp corresponding to a total travel time of a signal returning to the borehole by adding the first travel time and the third travel time, and calculating a start time difference δ_ρ between the first travel time and the second travel time;

searching all values of arrival time t, start time difference δ, and the frequency ratio d in the correlated signal M(t, δ, d) to determine for which value of arrival TNL, start time difference 6NL and frequency ratio dNL a band-limited spike signal

M(TNL, 8NL, d L) generated by the non-linear mixing process occurs; and

calculating, by the processor, trajectories of acoustic propagation using an elevation angle of the first acoustic signal, and a position of the first acoustic source, and an elevation angle of the second acoustic signal, and a position of the second acoustic source based on the propagation compressional and shear velocity model; mapping the band-limited spike signal M(TNL, 8NL, d L) generated by the nonlinear mixing process to spatial coordinates of an interaction mixing zone where the first acoustic signal and the second acoustic signal intersect; and

repeating the mapping of the band-limited spike signal M(T_NL,6NL,dNL) to the spatial coordinates of the interaction mixing zone for a plurality of values of the position of the first acoustic source, the position of the second acoustic source, the position of the receiver, an elevation angle of the first acoustic signal, an azimuth angle of the first acoustic signal, an elevation angle of the second acoustic signal, or an azimuth angle of the second acoustic signal, or any combination of two or more thereof, and generating a three dimensional image of a strength of the non-linear mixing process.

128. The method according to claim 127, further comprising determining a location of the mixing zone where the first acoustic signal and the second acoustic signal intersect using the position of the first acoustic source, the elevation angle and the azimuth angle of the first acoustic signal and the position of the second acoustic source, the elevation angle and the azimuth angle of the second acoustic signal.

129. A method of generating a three-dimensional image of a propagation compressional velocity, a propagation of shear velocity, a ratio of compressional velocity to shear velocity of a rock formation, or any combination of two or more thereof, the method comprising:

reading, by a processor, measurement parameters including a location of a first acoustic source configured to generate a first acoustic signal at a first frequency, a position of second acoustic source configured to generate a second acoustic signal at a second frequency, a position of a receiver configured to receive a detected signal returning to a borehole from a mixing zone within the rock formation at the intersection of the first acoustic signal and the second acoustic signal, and elevation and azimuth angles of the first and second acoustic signals;

reading, by the processor, a correlated signal M(t, δ, d), which is a function of arrival time t, start time difference δ and frequency ratio d between the first frequency and the second frequency, the correlated signal M(t, δ, d) containing signals generated by a non-linear mixing process from a first acoustic signal and a second acoustic signal at a mixing zone within the rock formation;

reading, by the processor, an initial propagation compressional and shear velocity model; estimating, by the processor, an initial propagating velocity model using well logs in a borehole and assumptions about lateral continuity of the rock formation away from the borehole;

calculating, by the processor, a first travel time of the first acoustic signal from the first acoustic source to a center of the mixing zone;

calculating, by the processor, a second travel time of the second acoustic signal from the second acoustic source to the center of the mixing zone;

calculating, by the processor, a third travel time between the center of the mixing zone and the receiver configured to receive the detected signal returning to the borehole having a frequency equal to a difference between the first frequency and the second frequency, the detected signal being generated by a non-linear mixing process from the first acoustic signal and the second acoustic signal in mixing zone;

calculating, by the processor, an arrival time Tp corresponding to a total travel time of a signal returning to the borehole by adding the first travel time and the third travel time, and calculating a start time difference δ_ρ between the first travel time and the second travel time;

searching all values of arrival time t, start time difference δ, and the frequency ratio d in the correlated signal M(t, δ, d) to determine for which value of arrival time T_NL, start time difference 6NL and frequency ratio dm, a band-limited spike signal M(TNL, 5NL, d L) generated by the non-linear mixing process occurs;

comparing the calculated arrival time Tp with a measured arrival time T L or comparing the calculated start time difference δ_ρ with a measured start time difference updating the propagating velocity model by tomographic velocity inversion or full wave form inversion using a difference between the measured and the calculated arrival times if the calculated arrival time is different from the measured arrival time, or using a difference between the calculated and the measured start time differences if the calculated start time difference is different from the measured start time difference, or both; repeating iteratively the above steps until the calculated arrival time is substantially equal to the measured arrival time or the calculated start time difference is substantially equal to the measured start time difference, or both, for a range of values of the position of the first acoustic source, the position of the second acoustic source, the position of the receiver, the azimuth angle of the first acoustic signal, the azimuth angle of the second acoustic signal, the elevation angle of the first acoustic signal, or the elevation angle of the second acoustic signal, or any combination of two or more thereof, to obtain a propagating velocity Vp and propagating velocity Vs that results in the calculated arrival time being substantially equal to the measured arrival time or the calculated start time difference being substantially equal to the measured start time difference, or both; and

generating a three dimensional image of a propagation compressional velocity, a propagation of shear velocity, a ratio of compressional velocity to shear velocity, or any combination of two or more thereof.

130. The method according to claim 129, further comprising determining a value of the propagating compressional velocity, a value of the propagating shear velocity, or a value of the velocity ratio of the compressional velocity and the shear velocity, or any combination of two or more thereof.

131. The method according to claim 129, further comprising calculating spatial coordinates of the mixing zone where the two acoustic signals intersect, calculating a convergence angle between the first acoustic signal from the first acoustic source and the second acoustic signal from the second acoustic source based on the elevation angle and the azimuth angle of the first acoustic signal, and the position of the first acoustic source, and the elevation angle and the azimuth angle of the second acoustic signal, and the position of the second acoustic source.

132. The method according to claim 131, further comprising calculating the velocity ratio of compressional velocity and shear velocity at the mixing zone from the measured frequency ratio dNL and the calculated the convergence angle using a selection rule relationship between velocity ratio of compressional velocity and shear velocity, frequency ratio and angle of convergence.

133. The method of claim 132, further comprising repeating mapping the calculated velocity ratio to the spatial coordinate of the mixing zone for a range of values of the position of the first acoustic source, the position of the second acoustic source, the position of the receiver, the azimuth angle of the first acoustic signal, the azimuth angle of the second acoustic signal, the elevation angle of the first acoustic signal, or the elevation angle of the second acoustic signal, or any combination of two or more thereof to obtain a three dimensional image of the velocity ratio.

134. A system for generating a three-dimensional image of a rock formation comprising:

a first acoustic source configured to generate a first acoustic signal, wherein the first acoustic signal comprises a first plurality of pulses arranged in a time sequence, the first plurality of pulses being separated in time, each pulse comprising a first modulated signal at a central frequency, wherein central frequencies of two consecutive pulses are different;

a second acoustic source configured to generate a second acoustic signal by a second acoustic source, wherein the second acoustic signal comprises a second plurality of pulses arranged as a time sequence, the second plurality of pulses being separated in time, wherein a separation in time between centers of two consecutive pulses is the same as a separation in time between centers of two corresponding pulses in the first plurality of pulses, wherein a start time difference is provided between a start time of a broadcast of the second plurality of pulses and a start time of a broadcast of the first plurality of pulses, wherein each pulse comprises a modulated signal, and wherein a central frequency of the modulated signal within each pulse in the second plurality of pulses is a selected fraction of the central frequency of the modulated signal within each corresponding pulse in the first plurality of pulses; wherein the first acoustic source and the second acoustic source are arranged within the borehole and are controllable such that trajectories of the first and the second acoustic signals intersect in an intersection volume outside the borehole; a receiver arranged within the borehole, the receiver being configured to detect a received signal returning to the borehole, the received signal including a signal being generated by a non-linear mixing process from the first acoustic signal and the second acoustic signal in a non-linear mixing zone within the intersection volume; a first processor configured to perform data processing on the received signal to extract the signal generated by the non-linear mixing process over noise or over signals generated by a linear interaction process, or both; and

a second processor configured to generate a three dimensional image of a propagation compressional velocity, a propagation of shear velocity, a ratio of compressional velocity and shear velocity, or non-linear properties of a rock formation, or any combination of two or more thereof based on the signal generated by the non-linear mixing process.

135. The system according to claim 134, wherein a separation in time between centers of two consecutive pulses in the first plurality of pulses is greater than a time duration of each pulse.

136. The system according to claim 134, wherein the second processor is configured to generate the three dimensional image of the propagation compressional velocity, the propagation of shear velocity, the ratio of compressional velocity and shear velocity, or the non-linear properties of the rock formation, or any combination of two or more thereof by using Kirchhoff imaging, beam imaging or wave equation imaging.

137. The system according to claim 136, wherein the second processor is configured to determine a value of the propagation compressional velocity or a value of the shear velocity or both using tomographic velocity inversion or full wave form inversion or by iterative imaging in combination with tomographic velocity inversion.

138. The system according to claim 134, wherein the second processor is configured to:

estimate an initial propagating velocity model using well logs in a borehole and assumptions about lateral continuity of a rock formation away from the borehole; calculate a first travel time of the first acoustic signal from the first acoustic source to a center of the mixing zone;

calculate a second travel time of the second acoustic signal from the second acoustic source to the center of the mixing zone;

calculate a third travel time between the center of the mixing zone and a receiver configured to receive a detected signal returning to the borehole having a frequency equal to a difference between the first frequency and the second frequency, the detected signal being generated by a non-linear mixing process from the first acoustic signal and the second acoustic signal in the mixing zone;

calculate an arrival time Tp corresponding to a total travel time of a signal returning to the borehole by adding the first travel time and the third travel time, and calculating a start time difference δ_ρ between the first travel time and the second travel time;

extract from the signal received at the receiver, signals generated by the nonlinear mixing process to yield for each start time difference δ and for each frequency ratio a correlated signal M(t, δ, d), which is a function of arrival time t, the start time difference δ and the frequency ratio d, and which contains non-linear interaction signals;

search all values of arrival time t, start time difference δ, and frequency ratio d in the correlated signal M(t, δ, d) to determine a value of arrival time TNL, start time difference 6 L and frequency ratio dNL for which a band-limited spike signal generated by the non-linear mixing process occurs;

compare the calculated arrival time Tp with a measured arrival time TNL and compare the calculated start time difference δρ with a measured start time difference

¾STL; determine whether the calculated arrival time Tp is different from the measured arrival time TNL and the calculated start time difference δρ is different from the measured start time difference 6NL and, if so, update the propagating velocity model using a difference between the measured and the calculated arrival times and a difference between the calculated and the measured start time differences; and

repeat iteratively the above steps until the calculated arrival time corresponding to the total travel time is substantially equal to a measured arrival time or the calculated start time difference is substantially equal to a measured start time difference, or both, for a range of values of the position of the first acoustic source, the position of the second acoustic source, the position of the receiver, the azimuth angle of the first acoustic signal, the azimuth angle of the second acoustic signal, the elevation angle of the first acoustic signal, or the elevation angle of the second acoustic signal, or any combination of two or more thereof, to obtain a propagating velocity Vp and propagating velocity Vs that results in the calculated arrival time being substantially equal to the measured arrival time or the calculated time difference being substantially equal to the measured start time difference, or both; and generate a three dimensional image of a propagation compressional velocity, a propagation of shear velocity, a ratio of compressional velocity to shear velocity, or non-linear properties, or any combination of two or more thereof.

139. The system according to claim 134, wherein the first processor is configured to extract signals generated by the non-linear mixing process by correlating the detected signal with a template signal designed according to forecast properties of the signal according to the selection rules for non-collinear mixing in a non-linear medium, and repeating for a plurality of start time differences δ between the second plurality of pulses and the first plurality of pulses and for a plurality of frequency ratios d between the central frequency of the modulated signal within each pulse in the second plurality of pulses and the central frequency of the modulated signal within each corresponding pulse in the first plurality of pulses, the correlating between the template signal and the detected signal, to yield for each start time difference δ and each frequency ratio d a correlated signal M(t, δ, d) which contains a non-linear interaction signal.

140. The system according to claim 138, wherein the second processor is configured to calculate, based on a propagation velocity model, trajectories of the first acoustic signal emitted by a first acoustic source and the second acoustic signal emitted by the second acoustic source using the position of the first acoustic source, the elevation angle and the azimuth angle of the first acoustic signal and the position of the second acoustic source, the elevation angle and the azimuth angle of the second acoustic signal.

141. The system according to claim 140, wherein the second processor is configured to determine a location of a mixing zone where the first acoustic signal and the second acoustic signal interact non-linearly based on the position of the first acoustic source, the elevation angle and the azimuth angle of the first acoustic signal and the position of the second acoustic source, the elevation angle and the azimuth angle of the second acoustic signal.

142. The system according to claim 141 , wherein the second processor is configured to calculate spatial coordinates of the mixing zone where the two acoustic signals interact non-linearly, calculating a convergence angle between the first acoustic signal from the first acoustic source and the second acoustic signal from the second acoustic source based on the elevation angle and the azimuth angle of the first acoustic signal, and the position of the first acoustic source, and the elevation angle and the azimuth angle of the second acoustic signal, and a position of the second acoustic source.

143. The system according to claim 142, wherein the second processor is configured to map the correlated signal M(T_NL, 5NL, CINL) to the spatial coordinates of the mixing zone corresponding to the non-linear properties of the rock formation at the mixing zone.

144. The system according to claim 143, wherein the second processor is configured to calculate the velocity ratio between the compressional velocity and the shear velocity at the mixing zone from the frequency ratio dNL and the convergence angle using a selection rule relationship between velocity ratio of compressional velocity and shear velocity, frequency ratio and angle of convergence and mapping the calculated velocity ratio to the spatial coordinates of the mixing zone.

145. The system of claim 144, wherein the second processor is configured to map the correlated signal M(T_NL, 6NL, C NL) to the spatial coordinates of the mixing zone or map the calculated velocity ratio to the spatial coordinate of the mixing zone, or both, for a range of values of the position of the first acoustic source, the position of the second acoustic source, the position of the receiver, the azimuth angle of the first acoustic signal, the azimuth angle of the second acoustic signal, the elevation angle of the first acoustic signal, or the elevation angle of the second acoustic signal, or any combination of two or more thereof, to obtain a three dimensional image of a strength of the non-linear properties or a three dimensional image of the velocity ratio, or both, by combining the measured values through repeated measurements.

146. The system according to claim 134, wherein the second processor is configured to determine a value of the propagation compressional velocity, a value of the propagation of shear velocity, or a value of the ratio of compressional velocity and shear velocity, or any combination of two or more thereof.

147. The system according to claim 146, wherein the second processor is configured to:

estimate an initial propagating velocity model using well logs in a borehole and assumptions about lateral continuity of a the rock formation away from the borehole;

calculate a first travel time of the first acoustic signal from the first acoustic source to a center of the mixing zone; calculate a second travel time of the second acoustic signal from the second acoustic source to the center of the mixing zone;

calculate a third travel time between the center of the mixing zone and a receiver configured to receive a detected signal returning to the borehole having a frequency equal to a difference between the first frequency and the second frequency, the detected signal being generated by a non-linear mixing process from the first acoustic signal and the second acoustic signal in a mixing zone;

calculate an arrival time Tp corresponding to a total travel time of a signal returning to the borehole by adding the first travel time and the third travel time, and calculating a start time difference δρ between the first travel time and the second travel time;

extract from the signal received at the receiver, signals generated by the nonlinear mixing process, to yield for each start time difference δ and for each frequency ratio d, a correlated signal M(t, δ, d) which is a function of arrival time t, the start time difference δ and the frequency ratio d, and which contains non-linear interaction signals;

search all values of arrival time t, start time difference δ, and frequency ratio d in the correlated signal M(t, δ, d) to determine a value of arrival time TNL, start time difference 6NL and frequency ratio dNL for which a signal generated by the non-linear mixing process occurs;

compare the calculated arrival time Tp with a measured arrival time T L and compare the calculated start time difference δ_ρ with a measured start time difference determine whether the calculated arrival time Tp is different from the measured arrival time T L and the calculated time difference δρ is different from the measured start time difference 6NL and, if so, update the propagating velocity model using a difference between the measured and the calculated arrival times and a difference between the calculated and the measured start time differences; and

repeat iteratively the above steps until the calculated arrival time corresponding to the total travel time is substantially equal to a measured arrival time, or the calculated time difference is substantially equal to a measured start time difference, or both, for a range of values of the position of the first acoustic source, the position of the second acoustic source, the position of the receiver, the azimuth angle of the first acoustic signal, the azimuth angle of the second acoustic signal, the elevation angle of the first acoustic signal, or the elevation angle of the second acoustic signal, or any combination of two or more thereof, to obtain a propagating velocity Vp and propagating velocity Vs that result in the calculated arrival time being substantially equal to the measured arrival time, or the calculated start time difference being substantially equal to the measured start time difference, or both.

148. The system according to claim 147, wherein the second processor is configured to calculate, based on a propagation velocity model, trajectories of the first acoustic signal having the first frequency emitted by the first acoustic source and the second acoustic signal having the second frequency emitted by the second acoustic source using the position of the first acoustic source, the elevation angle and the azimuth angle of the first acoustic signal and the position of the second acoustic source, the elevation angle and the azimuth angle of the second acoustic signal, and the position of the receiver.

149. The system according to claim 148, wherein the second processor is configured to calculate spatial coordinates of the mixing zone where the two acoustic signals interact non-linearly, calculate a convergence angle between the first acoustic signal from the first acoustic source and the second acoustic signal from the second acoustic source based on the elevation angle and the azimuth angle of the first acoustic signal, and the position of the first acoustic source, and the elevation angle and the azimuth angle of the second acoustic signal, and the position of the second acoustic source.

150. The system according to claim 149, wherein the second processor is configured to calculate the velocity ratio between the compressional velocity and the shear velocity at the mixing zone from the frequency ratio dNL and the convergence angle using a selection rule relationship between velocity ratio of compressional velocity and shear velocity, frequency ratio and angle of convergence and to map the calculated velocity ratio to the spatial coordinates of the mixing zone.

151. The system of claim 150, wherein the second processor is configured to repeat mapping the calculated velocity ratio to the spatial coordinates of the mixing zone for a range of values of the position of the first acoustic source, the position of the second acoustic source, the position of the receiver, the azimuth angle of the first acoustic signal, the azimuth angle of the second acoustic signal, the elevation angle of the first acoustic signal, or the elevation angle of the second acoustic signal, or any combination of two or more thereof, to obtain the value of the velocity ratio for a volume of the rock formation around the borehole.

152. A system for generating a three-dimensional image of non-linear properties of a rock formation, comprising a processor configured to:

read measurement parameters including a position of a first acoustic source configured to generate a first acoustic signal at a first frequency, a position of a second acoustic source configured to generate a second acoustic signal at a second frequency, a position of a receiver configured to receive a detected signal returning to a borehole from a mixing zone within the rock formation at the intersection of the first acoustic signal and the second acoustic signal, and elevation and azimuth angles of the first and second acoustic signals;

read a correlated signal M(t, δ, d), which is a function of arrival time t, start time difference δ and frequency ratio d between the first frequency and the second frequency, the correlated signal M(t, δ, d) containing signals generated by a non-linear mixing process from a first acoustic signal and a second acoustic signal at a mixing zone within the rock formation;

read an initial propagation compressional and shear velocity model;

calculate a first travel time of the first acoustic signal from the first acoustic source to a center of the mixing zone; calculate a second travel time of the second acoustic signal from the second acoustic source to the center of the mixing zone;

calculate a third travel time between the center of the mixing zone and the receiver configured to receive the detected signal returning to the borehole having a frequency equal to a difference between the first frequency and the second frequency, the detected signal being generated by a non-linear mixing process from the first acoustic signal and the second acoustic signal in the mixing zone;

calculate an arrival time corresponding to a total travel time Tp of a signal returning to the borehole by adding the first travel time and the third travel time, and calculating a start time difference δ_ρ between the first travel time and the second travel time;

search all values of arrival time t, start time difference δ, and the frequency ratio d in the correlated signal M(t, δ, d) to determine for which value of arrival TNL, start time difference 6NL and frequency ratio dm, a band-limited spike signal M(TNL, 8NL, d L) generated by the non-linear mixing process occurs; and

calculate trajectories of acoustic propagation using an elevation angle of the first acoustic signal, and a position of the first acoustic source, and an elevation angle of the second acoustic signal, and a position of the second acoustic source based on the propagation compressional and shear velocity model;

map the band-limited spike signal M(T L,5NL,d L) generated by the non-linear mixing process to spatial coordinates of an interaction mixing zone where the first acoustic signal and the second acoustic signal intersect; and

repeat the mapping of the band-limited spike signal M(T_NL,6NL,dNL) to the spatial coordinates of the interaction mixing zone for a plurality of values of the position of the first acoustic source, the position of the second acoustic source, the position of the receiver, an elevation angle of the first acoustic signal, an azimuth angle of the first acoustic signal, an elevation angle of the second acoustic signal, or an azimuth angle of the second acoustic signal, or any combination of two or more thereof, and generating a three dimensional image of a strength of the non-linear mixing process.

153. The system according to claim 152, further comprising determining a location of the mixing zone where the first acoustic signal and the second acoustic signal intersect using the position of the first acoustic source, the elevation angle and the azimuth angle of the first acoustic signal and the position of the second acoustic source, the elevation angle and the azimuth angle of the second acoustic signal.

154. A system for generating a three-dimensional image of a propagation compressional velocity, a propagation of shear velocity, a ratio of compressional velocity to shear velocity of a rock formation, or any combination of two or more thereof, the system comprising a processor configured to:

read measurement parameters including a location of a first acoustic source configured to generate a first acoustic signal at a first frequency, a position of second acoustic source configured to generate a second acoustic signal at a second frequency, a position of a receiver configured to receive a detected signal returning to a borehole from a mixing zone within the rock formation at the intersection of the first acoustic signal and the second acoustic signal, and elevation and azimuth angles of the first and second acoustic signals;

read a correlated signal M(t, δ, d), which is a function of arrival time t, start time difference δ and frequency ratio d between the first frequency and the second frequency, the correlated signal M(t, δ, d) containing signals generated by a non-linear mixing process from a first acoustic signal and a second acoustic signal at a mixing zone within the rock formation;

read an initial propagation compressional and shear velocity model;

estimate an initial propagating velocity model using well logs in a borehole and assumptions about lateral continuity of the rock formation away from the borehole;

calculate a first travel time of the first acoustic signal from the first acoustic source to a center of the mixing zone;

calculate a second travel time of the second acoustic signal from the second acoustic source to the center of the mixing zone; calculate a third travel time between the center of the mixing zone and the receiver configured to receive the detected signal returning to the borehole having a frequency equal to a difference between the first frequency and the second frequency, the detected signal being generated by a non-linear mixing process from the first acoustic signal and the second acoustic signal in mixing zone;

calculate an arrival time Tp corresponding to a total travel time of a signal returning to the borehole by adding the first travel time and the third travel time, and calculating a start time difference δ_ρ between the first travel time and the second travel time;

search all values of arrival time t, start time difference δ, and the frequency ratio d in the correlated signal M(t, δ, d) to determine for which value of arrival time T_NL, start time difference 6NL and frequency ratio dm, a band-limited spike signal M(TNL, 8NL, d L) generated by the non-linear mixing process occurs;

compare the calculated arrival time Tp with a measured arrival time TNL or comparing the calculated start time difference δ_ρ with a measured start time difference update the propagating velocity model by tomographic velocity inversion or full wave form inversion using a difference between the measured and the calculated arrival times if the calculated arrival time is different from the measured arrival time, or using a difference between the calculated and the measured start time differences if the calculated start time difference is different from the measured start time difference, or both;

repeat iteratively the above steps until the calculated arrival time is substantially equal to the measured arrival time or the calculated start time difference is substantially equal to the measured start time difference, or both, for a range of values of the position of the first acoustic source, the position of the second acoustic source, the position of the receiver, the azimuth angle of the first acoustic signal, the azimuth angle of the second acoustic signal, the elevation angle of the first acoustic signal, or the elevation angle of the second acoustic signal, or any combination of two or more thereof, to obtain a propagating velocity Vp and propagating velocity Vs that results in the calculated arrival time being substantially equal to the measured arrival time or the calculated start time difference being substantially equal to the measured start time difference, or both; and

generate a three dimensional image of a propagation compressional velocity, a propagation of shear velocity, a ratio of compressional velocity to shear velocity, or any combination of two or more thereof.

155. The system according to claim 154, wherein the processor is configured to determine a value of the propagating compressional velocity, a value of the propagating shear velocity, or a value of the velocity ratio of the compressional velocity and the shear velocity, or any combination of two or more thereof.

156. The system according to claim 154, wherein the processor is configured to calculate spatial coordinates of the mixing zone where the two acoustic signals intersect, calculating a convergence angle between the first acoustic signal from the first acoustic source and the second acoustic signal from the second acoustic source based on the elevation angle and the azimuth angle of the first acoustic signal, and the position of the first acoustic source, and the elevation angle and the azimuth angle of the second acoustic signal, and the position of the second acoustic source.

157. The system according to claim 156, wherein the processor is configured to calculate the velocity ratio of compressional velocity and shear velocity at the mixing zone from the measured frequency ratio dNL and the calculated the convergence angle using a selection rule relationship between velocity ratio of compressional velocity and shear velocity, frequency ratio and angle of convergence.

158. The method of claim 157, further comprising repeating mapping the calculated velocity ratio to the spatial coordinate of the mixing zone for a range of values of the position of the first acoustic source, the position of the second acoustic source, the position of the receiver, the azimuth angle of the first acoustic signal, the azimuth angle of the second acoustic signal, the elevation angle of the first acoustic signal, or the elevation angle of the second acoustic signal, or any combination of two or more thereof to obtain a three dimensional image of the velocity ratio.

159. A method of investigating rock formations outside a borehole, the method comprising:

generating, by a first acoustic source disposed within a borehole, a first conical acoustic signal, the first conical acoustic signal being broadcast at a first frequency; generating, by a second acoustic source disposed within the borehole, a second conical acoustic signal, the second conical acoustic signal being broadcast at a second frequency different from the first frequency, the first frequency and the second frequency being in the frequency range between approximately 500 Hz and 500 kHz, configuring the first acoustic source and the second acoustic source so that the first conical acoustic signal and the second conical acoustic signal intersect in a desired intersection volume outside the borehole; and

receiving, by a receiver disposed within the borehole, a detected signal returning to the borehole having a frequency equal to a difference between the first frequency and the second frequency, the detected signal being generated by a nonlinear mixing process from the first conical acoustic signal and the second conical acoustic signal in a non-linear mixing zone within the intersection volume;

recording the detected signal and storing the detected signal in a storage device; and recording, in the storage device, measurement parameters including a position of the first acoustic source, a position of the second acoustic source, a position of the receiver, and a cone angle and a cone axis orientation of the first signal and a cone angle and a cone axis orientation of second conical acoustic signal.

160. The method according to claim 159, wherein generating the first conical acoustic signal and generating the second conical acoustic signal comprise generating the first conical acoustic signal and the second conical acoustic signal while the first acoustic source and the second acoustic are located in a single borehole.

161. The method according to claim 159, wherein the borehole is not straight or the borehole includes sidetracks such that the first conical acoustic signal and the second conical acoustic signal have axes of symmetry that are not collinear.

162. The method according to claim 159, wherein generating the first conical acoustic signal comprises generating a first acoustic conical signal comprising a first plurality of pulses arranged in a coded time sequence, the first plurality of pulses being separated in time, each pulse comprising a first modulated signal at a central frequency, wherein central frequencies of two consecutive pulses are different.

163. The method according to claim 160, wherein a separation in time between centers of two consecutive pulses in the first plurality of pulses is greater than a time duration of each pulse.

164. The method according to claim 159, wherein generating the second conical acoustic signal comprises generating a second plurality of pulses arranged in a coded time sequence, the second plurality of pulses being separated in time, wherein a separation in time between centers of two consecutive pulses is the same as a separation in time between centers of two corresponding pulses in the first plurality of pulses, wherein a start time difference δ is provided between a start time of the broadcast of the second plurality of pulses and a start time of the broadcast of the first plurality of pulses, wherein each pulse comprises a modulated signal, and a central frequency of the modulated signal within each pulse in the second plurality of pulses is a selected fraction frequency ratio d of the central frequency of the modulated signal within each corresponding pulse in the first plurality of pulses.

165. The method according to claim 164, further comprising repeating the generating of the first acoustic signal, the generating of the second acoustic signal and the receiving of the detected signal for a plurality of start time differences δ between the first acoustic signal and the second acoustic signal and for a plurality of frequency ratios d between the first frequency and the second frequency.

166. The method according claim 159, wherein generating the first conical acoustic signal and generating the second conical acoustic signal comprises controlling a phase of the first conical acoustic signal or a phase of the second conical acoustic signal, or both such that the broadcast signals have opposite polarity.

167. The method according to claim 166, further comprising adding the signal detected at the receiver, generated by the mixing process from the first conical acoustic signal and the second conical acoustic signal, to the signal detected at the receiver, generated by the mixing process from the first conical acoustic signal having opposite polarity and the second acoustic signal having opposite polarity, in order to enhance non-linear interaction signals and minimize linear interaction signals and noise.

168. The method according claim 159, further comprising receiving a plurality of detected signals by a plurality of detectors and applying multi-dimensional filtering, or time moveout analysis and stacking, or both, to the plurality of detected signals that are received by the plurality of detectors.

169. The method according to claim 168, wherein receiving the detected signal by the plurality of detectors comprises receiving component detected signals of each of the detected signals by a plurality of component detectors of each of the plurality of detectors, and applying multi-dimensional filtering, or time moveout analysis and stacking, or both, to the plurality of detected component signals that are received on the plurality of component detectors of each of the plurality of receivers.

170. The method according to claim 159, further comprising repeating the generating of the first conical acoustic signal, the generating of the second conical acoustic signal and the receiving of the detected signal for different positions of the first acoustic source, different positions of the second acoustic source, different positions of the receiver, broadcast cone angles of the first acoustic source, or broadcast cone angles of the second acoustic source, or any combination of two or more thereof.

171. The method according to claim 159, further comprising extracting and enhancing the signal generated by the non-linear mixing process and suppressing signals generated by a linear interaction process or noise, or both, from the recorded detected signal generated by non-linear interaction in the rock formation from the first and second acoustic sources.

172. The method according to claim 159, wherein generating the first conical acoustic signal comprises generating a first plurality of pulses arranged in a coded time sequence, the first plurality of pulses being separated in time, each pulse comprising a first modulated signal at a central frequency, wherein central frequencies of two consecutive pulses are different,

wherein generating the second conical acoustic signal comprises generating a second plurality of pulses arranged in a coded time sequence, the second plurality of pulses being separated in time, wherein a separation in time between centers of two consecutive pulses is the same as a separation in time between centers of two corresponding pulses in the first plurality of pulses, wherein a start time difference δ is provided between a start time of the broadcast of the second plurality of pulses and a start time of the broadcast of the first plurality of pulses, wherein each pulse comprises a second modulated signal, and a central frequency of the second modulated signal within each pulse in the second plurality of pulses is a selected fraction frequency ratio d of the central frequency of the first modulated signal within each corresponding pulse in the first plurality of pulses, and

reading the recorded data and correlating, using a processor, the detected signal with a template signal comprising a plurality of pulses arranged in a time sequence, wherein a separation in time between centers of two consecutive pulses is substantially the same as a separation in time between centers of two consecutive pulses in the first plurality of pulses, wherein each pulse in the plurality of pulses comprises a modulated signal having a central frequency that is equal to a difference between a central frequency of the first modulated signal and a central frequency of the second modulated signal in each corresponding pulse of the first plurality of pulses and the second plurality of pulses.

173. The method according to claim 172, wherein a separation in time between centers of two consecutive pulses in the first plurality of pulses is greater than a time duration of each pulse.

174. The method according to claim 172, repeating for a plurality of start time differences δ and for a plurality of frequency ratios d between the first frequency and the second frequency, the correlating between the template signal and the detected signal to yield a correlated signal M(t, δ, d) which contains non-linear interaction signals.

175. The method according to claim 174, further comprising applying a time- variant frequency band pass signal to the detected signal so as to keep a narrow frequency band around an expected frequency bandwidth of the detected signal, the narrow frequency bandwidth being selected around the difference between the first frequency and the second frequency of each pulse.

176. The method according to claim 172, further comprising performing, using a processor, a hodogram analysis of the detected signal.

177. The method of claim 159, further comprising scanning a frequency ratio of the first frequency to the second frequency or scanning a range of start time difference, or both, so as to control a position of the mixing zone within the intersection volume in accordance with selection rules for non-linear interaction of acoustic signals within the rock formation.

178. The method of claim 159, further comprising: positioning the first acoustic source in a first axial orientation within the borehole and positioning the second acoustic source in a second axial orientation within the borehole, the first axial orientation and the second axial orientation being non-coaxial orientations or non- parallel.

179. The method of claim 178, further comprising controlling the frequency ratio of the first frequency and the second frequency and the start time difference to select a specific mixing zone within the intersection volume located at a locus of the intersection of the first conical acoustic signal and the second conical acoustic signal.

180. The method of claim 179, further comprising discriminating between mixing zones located on the locus of the intersection of the first conical broadcast and the second conical broadcast by a total of a travel time between the first acoustic source to the center of the mixing zone and from the center of the mixing zone to the receiver, or a convergence angle between the first conical acoustic signal and the second conical acoustic signal to the mixing zone according to non-linear interaction selection rule relationship between convergence angle and frequency ratio of the first and second acoustic signals, or both.

181. The method of claim 159, further comprising positioning the first acoustic source within a first borehole and positioning the second acoustic source within a second borehole and positioning the receiver within either the first borehole or the second borehole, or a third borehole.

182. The method of claim 181, further comprising selecting the cone angles or cone axis orientation or both such that the first conical acoustic signal and the second conical acoustic signal intersect at a targeted region in the rock formation surrounding the borehole.

183. A method of determining a value of propagation compressional or shear velocity in a volume of rock formation surrounding a borehole, or both, by travel time tomographic inversion, the method comprising:

estimating an initial propagating velocity model using well logs in a borehole and assumptions about lateral continuity of a rock formation away from the borehole; calculating, by a processor, spatial coordinates of mixing zones from positions of a first acoustic source broadcasting a first conical acoustic signal and a second acoustic source broadcasting a second conical acoustic signal and the elevation angles of the first and second conical conical acoustic signals;

calculating, by the processor, total travel times from the first acoustic source to the centers of the mixing zones and from the centers of the mixing zones to the receiver;

calculating, by the processor, differences in travel times from the first acoustic source to centers of the mixing zones and from second acoustic source to the centers of the mixing zones;

reading, by the processor, a correlated signal M(t, δ, d), which is a function of arrival time t, the start time difference δ between the second acoustic signal and the first acoustic signal and the frequency ratio d between a frequency of the first acoustic signal and a frequency of the second acoustic signal, and which contains non-linear interaction signals;

scanning, by the processor, a range of values of arrival time t, start time difference δ and frequency ratio d in the signal M(t, δ, d) to determine a plurality of values for the arrival time T_NL, a plurality of values for the start time difference 5_NL, a plurality of values for the frequency ratio dNL for which the measured arrival time T_NL, the start time difference 5_NL and the frequency ratio dm, correspond to a band-limited spike signal M(TN_L, 8N_L, d L) generated by the non-linear mixing process;

comparing, by the processor, the calculated arrival times with measured arrival times T_NL or comparing the calculated start time differences with measured start time differences 5N_L, or both;

updating the propagating velocity model by tomographic velocity inversion or full wave form inversion using the differences between the measured and the calculated arrival times and the differences between the calculated and the measured start time differences for all receiver, first source, second source and elevation of the first conical acoustic signal and elevation of the second conical acoustic signal, if the calculated arrival times are different from the measured arrival times and the calculated start time difference are different from the measured start time differences, repeating iteratively the above steps for a plurality of detected signals obtained for a plurality of positions of the first acoustic source, a plurality of positions of the second acoustic source, a plurality of elevation angles of the first conical acoustic signal and a plurality of elevation angles of the second conical acoustic signal until the calculated arrival time corresponding to the total travel time is substantially equal to a measured arrival time or the calculated start time difference is substantially equal to a measured start time difference, or both, to obtain a propagating velocity Vp and propagating velocity Vs that results in the calculated arrival times being substantially equal to the measured arrival times and the calculated start time differences being substantially equal to the measured start time differences.

184. The method of claim 183, further comprising calculating trajectories of acoustic propagation and travel times, based on the propagating velocity model, using ray tracing or other acoustic numerical modeling techniques.

185. The method according to claim 183, further comprising generating a three dimensional image of non-linear properties or a three-dimensional image of velocity ratio of compressional velocity and shear velocity, or any combination of two or more thereof from the recorded signal data by mapping, Kirchhoff imaging, beam imaging or wave equation imaging or determining propagation compressional or shear velocity or both, using tomographic velocity inversion or full wave form inversion or iterative imaging in combination with tomographic velocity inversion.

186. A method of generating a three dimensional image of non-linear properties, the method comprising:

reading, by a processor, measurement parameters including a position of the first acoustic source configured to generate a first conical acoustic signal at a first frequency, a position of second acoustic source configured to generate a second conical acoustic signal at a second frequency, a position of a receiver configured to receive a detected signal returning to a borehole from a mixing zone within the rock formation at the intersection of the first conical acoustic signal and the second conical acoustic signal, and cone angle and conical axis orientation of the first conical acoustic signal and the second conical acoustic signal;

reading, by the processor, a correlated signal M(t, δ, d) which contains nonlinear interaction signals, which is a function of arrival time t, the start time difference δ between a start time of a broadcast of the first conical acoustic signal and a start time of a broadcast of second conical acoustic signal and the frequency ratio d between the first frequency and the second frequency;

reading, by the processor, an initial propagation compressional and shear velocity model;

calculating, by the processor, based on the propagation compressional and shear velocity model, trajectories of acoustic propagation, using a position of the first acoustic source and a position of the second acoustic source, and a cone angle and axial orientation of the first conical acoustic signal and a cone angle and axial orientation of the second conical acoustic signal;

calculating spatial coordinates of an intersection volume zone where the first conical acoustic signal and the second conical acoustic signal intersect;

searching, by the processor, the correlated signal M(t, δ, d) in a range of arrival time t, start time difference δ and frequency ratio d to determine a plurality of values for the arrival time TNL, a plurality of values for the start time difference 5NL, a plurality of values for the frequency ratio dNL for which the measured arrival time TNL, the start time difference 5NL and the frequency ratio dm, correspond to a band-limited spike signal M(TNL, 8NL, d L) generated by the non-linear mixing process; and

mapping the band-limited spike signals M^NL^NL^NL) generated by the nonlinear mixing process to spatial coordinates of the interaction mixing zones where the first acoustic conical broadcast signal and the second acoustic conical broadcast signal intersect; and

repeating the mapping of the band-limited spike signals Μ(Τ^,δ^ ^) to the spatial coordinates of the interaction mixing zones for a plurality of values of the position of the first acoustic source, the position of the second acoustic source, the position of the receiver, a cone angle and axial orientation of the first conical acoustic signal or a cone angle and axial orientation of the second conical acoustic signal, or any combination of two or more thereof, and combining mapped correlated signals to obtain a three dimensional image of a strength of the non-linear mixing process.

187. A method of generating a three dimensional image of a velocity ratio of compressional velocity and shear velocity, the method comprising:

reading, by a processor, measurement parameters including a position of the first acoustic source configured to generate a first conical acoustic signal at a first frequency, a position of second acoustic source configured to generate a second conical acoustic signal at a second frequency, a position of a receiver configured to receive a detected signal returning to a borehole from a mixing zone within the rock formation at the intersection of the first conical acoustic signal and the second conical acoustic signal, and cone angle and conical axis orientation of the first conical acoustic signal and the second conical acoustic signal;

reading, by the processor, a correlated signal M(t, δ, d) which contains non- linear interaction signals, which is a function of arrival time t, the start time difference δ between a start time of a broadcast of the first conical acoustic signal and a start time of a broadcast of second conical acoustic signal and the frequency ratio d between the first frequency and the second frequency;

reading, by the processor, an initial propagation compressional and shear velocity model;

calculating, by the processor, based on the propagation compressional and shear velocity model, trajectories of acoustic propagation, using a position of the first acoustic source and a position of the second acoustic source, and a cone angle and axial orientation of the first conical acoustic signal and a cone angle and axial orientation of the second conical acoustic signal;

calculating spatial coordinates of an intersection volume zone where the first conical acoustic signal and the second conical acoustic signal intersect;

calculating a convergence angle between the first conical acoustic signal and the second conical acoustic signal based on the cone angle and the axial orientation of the first conical acoustic signal and the cone angle and the axial orientation of the second conical acoustic signal, and a position of the first acoustic source, and a position of the receiver;

searching, by the processor, the correlated signal M(t, δ, d) in a range of values of arrival time t, start time difference δ and frequency ratio d to determine the plurality of values for the arrival time TNL, the plurality of values for the start time difference 5NL, the plurality of values for the frequency ratio dNL for which the measured arrival time TNL, the start time difference 5NL and the frequency ratio dNL correspond to a band-limited spike signal M(T_NL,6NL,dNL) generated by the non-linear mixing process;

calculating the velocity ratios using the frequency ratios dNL resulting in a detected non-linear signal and calculating the convergence angle using a selection rule relationship between the frequency ratio and the angle of convergence;

mapping the calculated velocity ratios to the spatial coordinates of the mixing zones;

repeating the mapping of the velocity ratios to the spatial coordinates of the mixing zone for a plurality of values of the position of the first acoustic source, the position of the second acoustic source, the position of the receiver, the cone angle and the axial orientation of the first conical acoustic signal or the cone angle and the axial orientation of the second conical acoustic signal, or any combination of two or more thereof, and combining mapped correlated signals to obtain a three dimensional image of the velocity ratio.

188. A system for investigating rock formations outside a borehole, the system comprising:

a first acoustic source disposed within a borehole, the first acoustic source being configured to generate a first conical acoustic signal, the first conical acoustic signal being broadcast at a first frequency;

a second acoustic source disposed within the borehole, the second acoustic source being configured to generate a second conical acoustic signal, the second conical acoustic signal being broadcast at a second frequency different from the first frequency, the first frequency and the second frequency being in the frequency range between approximately 500 Hz and 500 kHz,

wherein the first acoustic source and the second acoustic source are configured such that the first conical acoustic signal and the second conical acoustic signal intersect in a desired intersection volume outside the borehole; and

a receiver disposed within the borehole, the receiver being configured to receive a detected signal returning to the borehole having a frequency equal to a difference between the first frequency and the second frequency, the detected signal being generated by a non-linear mixing process from the first conical acoustic signal and the second conical acoustic signal in a non-linear mixing zone within the intersection volume;

a storage device configured to record the detected signal and store the detected signal and record and store measurement parameters including a position of the first acoustic source, a position of the second acoustic source, a position of the receiver, and a cone angle and a cone axis orientation of the first conical acoustic signal, and a cone angle and a cone axis orientation of second conical acoustic signal.

189. The system according to claim 188, wherein the first acoustic source generating first conical acoustic signal generating and the second acoustic source generating the second conical acoustic signal are located in a single borehole.

190. The system according to claim 188, wherein the borehole is not straight or the borehole includes sidetracks such that the first conical acoustic signal and the second conical acoustic signal have axes of symmetry that are not collinear.

191. The system according to claim 188, wherein the first conical acoustic signal comprises a first plurality of pulses arranged in a coded time sequence, the first plurality of pulses being separated in time, each pulse comprising a first modulated signal at a central frequency, wherein central frequencies of two consecutive pulses are different.

192. The system according to claim 191, wherein a separation in time between centers of two consecutive pulses in the first plurality of pulses is greater than a time duration of each pulse.

193. The system according to claim 188, wherein the second conical acoustic signal comprises a second plurality of pulses arranged in a coded time sequence, the second plurality of pulses being separated in time, wherein a separation in time between centers of two consecutive pulses is the same as a separation in time between centers of two corresponding pulses in the first plurality of pulses, wherein a start time difference δ is provided between a start time of the broadcast of the second plurality of pulses and a start time of the broadcast of the first plurality of pulses, wherein each pulse comprises a modulated signal, and a central frequency of the modulated signal within each pulse in the second plurality of pulses is a selected fraction frequency ratio d of the central frequency of the modulated signal within each corresponding pulse in the first plurality of pulses.

194. The system according to claim 193, further comprising a processor configured to control the first acoustic source and the second acoustic source and the receiver to repeat a generating of the first acoustic signal, a generating of the second acoustic signal and a receiving of the detected signal for a plurality of start time differences δ between the first acoustic signal and the second acoustic signal and for a plurality of frequency ratios d between the first frequency and the second frequency.

195. The system according claim 188, further comprising a phase controller configured to control a phase of the first conical acoustic signal or a phase of the second conical acoustic signal, or both such that the first conical acoustic signal and the second conical acoustic signal have opposite polarity.

196. The system according to claim 195, further comprising a processor configured to add the signal detected at the receiver, generated by the mixing process from the first conical acoustic signal and the second conical acoustic signal, to the signal detected at the receiver, generated by the mixing process from the first conical acoustic signal having opposite polarity and the second acoustic signal having opposite polarity, in order to enhance non-linear interaction signals and minimize linear interaction signals and noise.

197. The system according claim 188, further comprising a plurality of detectors configured to receive a plurality of detected signals and applying multi-dimensional filtering, or time moveout analysis and stacking, or both, to the plurality of detected signals that are received by the plurality of detectors.

198. The system according to claim 197, wherein the detector comprises a plurality of component detectors wherein the detected signal comprises component detected signals by the component detectors and further comprising a processor configured to apply multi-dimensional filtering, or time moveout analysis and stacking, or both, to the plurality of detected component signals that are received on the plurality of component detectors.

199. The system according to claim 188, further comprising a processor configured to control the first acoustic source, the second acoustic source and the receiver to repeat the generating of the first conical acoustic signal, the generating of the second conical acoustic signal and the receiving of the detected signal for different positions of the first acoustic source, different positions of the second acoustic source, different positions of the receiver, broadcast cone angles of the first acoustic source, or broadcast cone angles of the second acoustic source, or any combination of two or more thereof.

200. The system according to claim 188, further comprising a processor configured to extract and enhance the signal generated by the non-linear mixing process and suppress signals generated by a linear interaction process or noise, or both, from the recorded detected signal generated by non-linear interaction in the rock formation from the first and second acoustic sources.

201. The system according to claim 188, wherein the first conical acoustic signal comprises a first plurality of pulses arranged in a coded time sequence, the first plurality of pulses being separated in time, each pulse comprising a first modulated signal at a central frequency, wherein central frequencies of two consecutive pulses are different,

wherein the second conical acoustic signal comprises a second plurality of pulses arranged in a coded time sequence, the second plurality of pulses being separated in time, wherein a separation in time between centers of two consecutive pulses is the same as a separation in time between centers of two corresponding pulses in the first plurality of pulses, wherein a start time difference δ is provided between a start time of the broadcast of the second plurality of pulses and a start time of the broadcast of the first plurality of pulses, wherein each pulse comprises a second modulated signal, and a central frequency of the second modulated signal within each pulse in the second plurality of pulses is a selected fraction frequency ratio d of the central frequency of the first modulated signal within each corresponding pulse in the first plurality of pulses, and

further comprising a processor configured to read the recorded data and correlate the detected signal with a template signal comprising a plurality of pulses arranged in a time sequence, wherein a separation in time between centers of two consecutive pulses is substantially the same as a separation in time between centers of two consecutive pulses in the first plurality of pulses, wherein each pulse in the plurality of pulses comprises a modulated signal having a central frequency that is equal to a difference between a central frequency of the first modulated signal and a central frequency of the second modulated signal in each corresponding pulse of the first plurality of pulses and the second plurality of pulses.

202. The system according to claim 201, wherein a separation in time between centers of two consecutive pulses in the first plurality of pulses is greater than a time duration of each pulse.

203. The system according to claim 201, further comprising a processor configured to repeat for a plurality of start time differences δ and for a plurality of frequency ratios d between the first frequency and the second frequency, the correlation between the template signal and the detected signal to yield a correlated signal M(t, δ, d) which contains non-linear interaction signals.

204. The system according to claim 203, further comprising a processor configured to apply a time-variant frequency band pass signal to the detected signal so as to keep a narrow frequency band around an expected frequency bandwidth of the detected signal, the narrow frequency bandwidth being selected around the difference between the first frequency and the second frequency of each pulse.

205. The system according to claim 201, further comprising a processor configured to perform a hodogram analysis of the detected signal.

206. The system of claim 188, further comprising a processor configured to scan a frequency ratio of the first frequency to the second frequency, or scan a range of start time difference, or both, so as to control a position of the mixing zone within the intersection volume in accordance with selection rules for non-linear interaction of acoustic signals within the rock formation.

207. The system of claim 188, wherein the first acoustic source is positioned in a first axial orientation within the borehole and the second acoustic source is positioned in a second axial orientation within the borehole, the first axial orientation and the second axial orientation being non-coaxial orientations or non-parallel.

208. The system of claim 207, further comprising a processor configured to control the frequency ratio of the first frequency and the second frequency and a start time difference to select a specific mixing zone within the intersection volume located at a locus of the intersection of the first conical acoustic signal and the second conical acoustic signal.

209. The system of claim 208, further comprising a processor configured to discriminate between mixing zones located on the locus of the intersection of the first conical broadcast and the second conical broadcast by a total of a travel time between the first acoustic source to the center of the mixing zone and from the center of the mixing zone to the receiver, or a convergence angle between the first conical acoustic signal and the second conical acoustic signal to the mixing zone according to nonlinear interaction selection rule relationship between convergence angle and frequency ratio of the first and second acoustic signals, or both.

210. The system of claim 188, wherein the first acoustic source is positioned within a first borehole and the second acoustic source is positioned within a second borehole and the receiver is positioned within either the first borehole or the second borehole, or a third borehole.

211. The system of claim 210, wherein the cone angles or cone axis orientation or both are selected such that the first conical acoustic signal and the second conical acoustic signal intersect at a targeted region in the rock formation surrounding the borehole.

212. A system for determining the value of propagation compressional or shear velocity in a volume of rock formation surrounding a borehole, or both, by travel time tomographic inversion, the system comprising a processor configured to:

estimate an initial propagating velocity model using well logs in a borehole and assumptions about lateral continuity of a rock formation away from the borehole; calculate spatial coordinates of mixing zones from the positions of a first acoustic source broadcasting a first acoustic signal in a form of a first conical signal broadcast and a second acoustic source broadcasting a second acoustic signal in a form of a second conical signal broadcast and the elevation angles of the first and second conical broadcasts;

calculate total travel times from the first acoustic source to the centers of the mixing zones and from the centers of the mixing zones to the receiver;

calculate the differences in travel times from the first acoustic source to the centers of the mixing zones and from second acoustic source to the centers of the mixing zones; read a correlated signal M(t, δ, d), which is a function of arrival time t, the start time difference δ between the second acoustic signal and the first acoustic signal and the frequency ratio d between a frequency of the first acoustic signal and a frequency of the second acoustic signal, and which contains non-linear interaction signals;

scan all values of arrival time t, start time difference δ and frequency ratio d in the signal M(t, δ, d) to determine the pluralities of values for the arrival time T L , the start time difference 5NL , the frequency ratio d^ and M(T_NL,6NL,dNL) for which the measured arrival time TNL, the start time difference 5NL and the frequency ratio dNL correspond to a band-limited spike signal M(T_NL,6NL,dNL) generated by the non-linear mixing process; ;

compare, by the processor, the calculated arrival times with measured arrival times T L or comparing the calculated start time differences with measured start time differences 5NL, or both;

update the propagating velocity model by tomographic velocity inversion or full wave form inversion using the differences between the measured and the calculated arrival times and the differences between the calculated and the measured start time differences for all receiver, first source, second source and broadcast elevation configurations, if the calculated arrival times are different from the measured arrival times and the calculated start time difference are different from the measured start time differences,

repeat iteratively the above steps for a plurality of detected signals obtained for a plurality of positions of the first acoustic source, a plurality of positions of the second acoustic source, a plurality of broadcast elevation angles of the first acoustic signal and the second acoustic signal until the calculated arrival time corresponding to the total travel time is substantially equal to a measured arrival time or the calculated start time difference is substantially equal to a measured start time difference, or both, to obtain a propagating velocity Vp and propagating velocity Vs that results in the calculated arrival times being substantially equal to the measured arrival times and the calculated start time differences being substantially equal to the measured start time differences.

213. The system of claim 212, further comprising a processor configured to calculate trajectories of acoustic propagation and travel times, based on the propagating velocity model, using ray tracing or other acoustic numerical modeling techniques.

214. The system according to claim 212, further comprising a processor configured to generate a three dimensional image of non-linear properties or a three dimensional image of velocity ratio of compressional velocity and shear velocity, or any combination of two or more thereof from the recorded signal data by mapping,

Kirchhoff imaging, beam imaging or wave equation imaging or determining propagation compressional or shear velocity or both, using tomographic velocity inversion or full wave form inversion or iterative imaging in combination with tomographic velocity inversion.

215. A system for generating a three dimensional image of non-linear properties, the system comprising a processor configured to:

read measurement parameters including a position of the first acoustic source configured to generate a first conical acoustic signal at a first frequency, a position of second acoustic source configured to generate a second conical acoustic signal at a second frequency, a position of a receiver configured to receive a detected signal returning to a borehole from a mixing zone within the rock formation at the intersection of the first conical acoustic signal and the second conical acoustic signal, and cone angle and conical axis orientation of the first conical acoustic signal and the second conical acoustic signal;

read a correlated signal M(t, δ, d) which contains non-linear interaction signals, which is a function of arrival time t, the start time difference δ between a start time of a broadcast of the first conical acoustic signal and a start time of a broadcast of second conical acoustic signal and the frequency ratio d between the first frequency and the second frequency;

read an initial propagation compressional and shear velocity model; calculate, based on the propagation compressional and shear velocity model, trajectories of acoustic propagation, using a position of the first acoustic source and a position of the second acoustic source, and a cone angle and axial orientation of the first conical acoustic signal and a cone angle and axial orientation of the second conical acoustic signal; calculate spatial coordinates of an intersection volume zone where the first conical acoustic signal and the second conical acoustic signal intersect; search the correlated signal M(t, δ, d) in a range of arrival time t, start time difference δ and frequency ratio d to determine a plurality of values for the arrival time T_NL, a plurality of values for the start time difference 5_NL, a plurality of values for the frequency ratio dNL for which the measured arrival time T_NL, the start time difference 5N_L and the frequency ratio dNL correspond to a band-limited spike signal M(TN_L, 8N_L, d L) generated by the non-linear mixing process;

map the band-limited spike signals M(TN_L, 8N_L, dNL) generated by the non- linear mixing process to spatial coordinates of the interaction mixing zones where the first acoustic conical broadcast signal and the second acoustic conical broadcast signal intersect; and repeat the mapping of the band-limited spike signals M(TNL,5NL,dNL) to the spatial coordinates of the interaction mixing zones for a plurality of values of the position of the first acoustic source, the position of the second acoustic source, the position of the receiver, a cone angle and axial orientation of the first conical acoustic signal or a cone angle and axial orientation of the second conical acoustic signal, or any combination of two or more thereof, and combining mapped correlated signals to obtain a three dimensional image of a strength of the non-linear mixing process.

216. A system for generating a three dimensional image of velocity ratio of compressional velocity and shear velocity, the system comprising a processor configured to:

read measurement parameters including a position of the first acoustic source configured to generate a first conical acoustic signal at a first frequency, a position of second acoustic source configured to generate a second conical acoustic signal at a second frequency, a position of a receiver configured to receive a detected signal returning to a borehole from a mixing zone within the rock formation at the intersection of the first conical acoustic signal and the second conical acoustic signal, and cone angle and conical axis orientation of the first conical acoustic signal and the second conical acoustic signal;

read a correlated signal M(t, δ, d) which contains non-linear interaction signals, which is a function of arrival time t, the start time difference δ between a start time of a broadcast of the first conical acoustic signal and a start time of a broadcast of second conical acoustic signal and the frequency ratio d between the first frequency and the second frequency;

read an initial propagation compressional and shear velocity model;

calculate based on the propagation compressional and shear velocity model, trajectories of acoustic propagation, using a position of the first acoustic source and a position of the second acoustic source, and a cone angle and axial orientation of the first conical acoustic signal and a cone angle and axial orientation of the second conical acoustic signal;

calculate spatial coordinates of an intersection volume zone where the first conical acoustic signal and the second conical acoustic signal intersect;

calculate a convergence angle between the first conical acoustic signal and the second conical acoustic signal based on the cone angle and the axial orientation of the first conical acoustic signal and the cone angle and the axial orientation of the second conical acoustic signal, and a position of the first acoustic source, and a position of the receiver;

search the correlated signal M(t, δ, d) in a range of values of arrival time t, start time difference δ and frequency ratio d to determine the plurality of values for the arrival time TNL, the plurality of values for the start time difference 5NL, the plurality of values for the frequency ratio dNL for which the measured arrival time TNL, the start time difference 5NL and the frequency ratio dNL correspond to a band-limited spike signal M(TNL, 8NL, dNL) generated by the non-linear mixing process; ; calculate the velocity ratios using the frequency ratios CINL resulting in a detected non-linear signal and the convergence angle using a selection rule relationship between the frequency ratio and the angle of convergence;

map the calculated velocity ratios to the spatial coordinates of the mixing zones;

repeating the mapping of the velocity ratios to the spatial coordinates of the mixing zone for a plurality of values of the position of the first acoustic source, the position of the second acoustic source, the position of the receiver, the cone angle and the axial orientation of the first conical acoustic signal or the cone angle and the axial orientation of the second conical acoustic signal, or any combination of two or more thereof, and combining mapped correlated signals to obtain a three dimensional image of the velocity ratio.

217. A system for investigating non-linear properties of a rock formation around a borehole, comprising:

a first sub-system configured to perform data acquisition, control and recording of data;

a second subsystem in communication with the first sub-system and configured to perform non-linearity and velocity preliminary imaging;

a third subsystem in communication with the first subsystem and configured to emit controlled acoustic broadcasts and receive acoustic energy;

a fourth subsystem in communication with the first subsystem and the third subsystem and configured to generate a source signal directed towards the rock formation; and

a fifth subsystem in communication with the third subsystem and the fourth subsystem and configured to perform detection of signals representative of the nonlinear properties of the rock formation.

218. The system according to claim 217, wherein the first subsystem comprises an acquisition modeling module configured to receive inputs from the second subsystem and operator inputs.

219. The system according to claim 218, wherein the first subsystem comprises a data acquisition controller configured to receive operator inputs and inputs from an acquisition design module.

220. The system according to claim 219, wherein the first subsystem comprises a data pre-processing and enhancement module and a data storage device, the data preprocessing and enhancement module configured to read data from the data storage device and input data into the second subsystem.

221. The system according to claim 217, wherein the second subsystem comprises an initial velocity model module configured to provide an initial velocity model derived from logs of compressional and shear slowness and information about formation lateral continuity to the first subsystem.

222. The system according to claim 221, wherein the second subsystem comprises an image initializing module configured to initialize volumes of propagation compressional and shear velocity values, velocity ratio images and images of non- linearity which are associated with amplitudes of a measured signal generated by the non-linear interaction at a mixing zone.

223. The system according to claim 217, wherein the third subsystem comprises a first acoustic source, a second acoustic source, one or more receivers, a tool deployment and transport module, a tool mechanical controller, a directional controller for controlling an orientation of the first acoustic signal, and a directional controller for controlling an orientation of the second acoustic signal.

224. The system according to claim 223, wherein the third system is configured to receive from the first subsystem commands to control the directional controllers within the third subsystem to control azimuth and elevation angles of the first acoustic signal, the second acoustic signal, or both.

225. The system according to claim 223, wherein acquisition record parameters pertaining to tool configuration and broadcast geometry are collected from the third subsystem and are recorded in a data storage device within the first subsystem.

226. The system according to claim 217, wherein the fourth subsystem comprises a signal generator configured to generate a first signal, a signal generator configured to generate a second acoustic signal, a frequency multiplier and start time difference module, and a coded signal generator.

227. The system according to claim 226, wherein the fourth subsystem is configured to receive input commands from the first subsystem.

228. The system according to claim 227, wherein the fourth subsystem is configured to receive input commands from the first subsystem so as to generate pulse sequences using the coded signal generator module.

229. The system according to claim 228, wherein the coded signal generated by the coded signal generator module is input into a first signal control and amplifier module to input a first signal at a first frequency to the first acoustic source, wherein the coded signal generated by the coded signal generator module is input into frequency multiplier and start time difference module and then input into a second signal control and amplifier module to provide a second signal to the second acoustic source having a start time difference relative to the first signal and having a second frequency which is a selected fraction of the first frequency.

230. The system according to claim 229, wherein the start time difference and the frequency fraction and are stored as record broadcast information in a data storage device within the first sub-system.

231. The system according to claim 217, wherein the fifth subsystem comprises a receiver module configured to receive a signal from one or more receivers, a nonlinear signal enhancement module configured to enhance the signal received by the signal receiver module, a template signal generator module configured to generate a template signal, and a signal correlation module configured to correlate the received signal with the template signal.

232. The system according to claim 231, wherein signals from the receiver module are processed by the non-linear signal enhancement module to enhance a content of non-linear origin and reduce noise, wherein the signals from the receiver module and enhanced signals are stored in a data storage device within the first sub-system.

233. The system according to claim 231, wherein the template generator is configured to generate a template signal, wherein the correlation module is configured to correlate the template signal with the received signal so as to extract a correlated signal the correlated signal identifying a signal generated at an intersection of a first signal from a first acoustic source and a second signal from a second acoustic in a non-linear mixing zone within the rock formation around the borehole, where t corresponds to signal time, d to a frequency ratio of the first frequency and the second frequency, δ is a time delay between the first signal and the second signal.

234. The system according to claim 233, wherein the correlated signal is stored in a data storage device within the first subsystem.

235. The system according to claim 234, wherein a data acquisition controller within the first subsystem is configured to repeat recording within the storage device of the first subsystem and repeat extracting the correlated data from a plurality of values of the position of the first acoustic source, the position of the second acoustic source, the position of the receiver, an elevation angle of the first acoustic signal, or an elevation angle of the second acoustic signal, or any combination of two or more thereof.

236. The system according to claim 235, further comprising a hodogram analysis module configured to process stored data within the storage device to further enhance the correlated data.

237. The system according to claim 217, further comprising a sixth subsystem in communication with the second subsystem and configured to perform imaging of non-linear properties of the rock formation and determining acoustic signal velocity within the rock formation.

238. The system according to claim 237, wherein the sixth subsystem comprises sub-system a data pre-processing and enhancement module, a velocity model iteration module, an imaging iteration module, an output image module for velocity ratio images and/or non-linearity images, and an output velocity module for outputting determined values of compressional velocity Vp, shear velocity Vs, or velocity ratio between the compressional velocity and the shear velocity, or any combination of two or more thereof.