Method and system for microseismic event wavefront estimation

Abstract

In the field of seismic data analysis, the problem of estimating features of a microseismic event waveform, detected at the sensors of a geophone array in the presence of background seismic energy and sensor noise, is considered. Embodiments of the present disclosure provide a method and system for estimating wavefront arrival times and polarizations across the sensors of an array, while simultaneously extracting the associated event sensor waveforms from the seismic background. The resulting waveform estimates are estimated in such a way that estimation errors are expected to be noise-like, and free of any systematic changes in the magnitude or phase of the arriving wavefronts.

Claims

1. A method of analyzing a microseismic event, comprising: receiving a plurality of traces associated with the microseismic event that were respectively generated by a plurality of sensors disposed underground; processing the plurality of traces to obtain a respective plurality of sets of samples that are time-aligned in accordance with a respective estimate of arrival time of a wavefront associated with the microseismic event at each of the plurality of sensors; and decomposing the plurality of sets of samples into a respective plurality of sets of one or more basis waveform shapes and corresponding sensor weights to form an estimate of the wavefront at the sensors while suppressing background noise associated with the traces, wherein the plurality of traces are associated with a microseismic survey of a subsurface region, and wherein the estimate of the wavefront associated with the microseismic event assists in analyzing the microseismic survey.

2. A method according to claim 1, wherein each of the plurality of traces includes two or more trace components, and wherein processing includes combining together the two or more trace components to form a respective single trace for each of the plurality of sensors from the two or more components.

3. A method according to claim 2, wherein combining is performed in accordance with an estimate of an orientation between the respective sensor and the microseismic event.

4. A method according to claim 2, wherein combining is performed in accordance with an estimate of an orientation at which SNR of the single trace is maximized.

5. A method according to claim 1, wherein processing includes windowing each of the plurality of traces around the estimated arrival time of the wavefront at the respective sensor to obtain the respective set of samples.

6. A method according to claim 1, further comprising correlating the estimate of the wavefront with the plurality of traces to re-estimate the respective arrival times for each of the sensors.

7. A method according to claim 6, further comprising estimating the respective arrival times.

8. A method of analyzing a microseismic event, comprising: receiving a plurality of traces associated with the microseismic event that were respectively generated by a plurality of sensors disposed underground; processing the plurality of traces to obtain a respective plurality of sets of samples that are time-aligned in accordance with a respective estimate of arrival time of a wavefront associated with the microseismic event at each of the plurality of sensors; and decomposing the plurality of sets of samples into a respective plurality of sets of one or more basis waveform shapes and corresponding sensor weights; and forming cleaned versions of the plurality of traces using some of the respective plurality of sets of sensor weights, wherein the plurality of traces are associated with a microseismic survey of a subsurface region, and wherein the cleaned versions of the plurality of traces assists in analyzing the microseismic survey.

9. A method according to claim 8, wherein forming includes: projecting the plurality of traces into a sensor basis in accordance with the plurality of corresponding weights.

10. A method according to claim 9, wherein forming further includes: time shifting the projected plurality of traces in accordance with the estimated respective arrival times for the plurality of sensors.

11. A system for analyzing a microseismic event, comprising: an input for receiving a plurality of traces associated with the microseismic event that were respectively generated by a plurality of sensors disposed underground; an alignment block that is adapted to process the plurality of traces to obtain a respective plurality of sets of samples that are time-aligned in accordance with a respective estimate of arrival time of a wavefront associated with the microseismic event at each of the plurality of sensors; and a decomposition block that is adapted to decompose the plurality of sets of samples into a respective plurality of sets of one or more basis waveform shapes and corresponding sensor weights to form an estimate of the wavefront at the sensors while suppressing background noise associated with the traces, wherein the plurality of traces are associated with a microseismic survey of a subsurface region, and wherein the estimate of the wavefront associated with the microseismic event assists in analyzing the microseismic survey.

12. A system according to claim 11, wherein each of the plurality of traces includes two or more trace components, and wherein the alignment block is further adapted to combine together the two or more trace components to form a respective single trace for each of the plurality of sensors from the two or more components.

13. A system according to claim 12, wherein combining is performed in accordance with an estimate of an orientation between the respective sensor and the microseismic event.

14. A system according to claim 12, wherein combining is performed in accordance with an estimate of an orientation at which SNR of the single trace is maximized.

15. A system according to claim 11, wherein the alignment block is further adapted to window each of the plurality of traces around the estimated arrival time of the wavefront at the respective sensor to obtain the respective set of samples.

16. A system according to claim 11, further comprising a correlation block that is adapted to correlate the estimate of the wavefront with the plurality of traces to re-estimate the respective arrival times for each of the sensors.

17. A method according to claim 1, wherein the processing and decomposing are performed by a computer.

18. A method according to claim 1, further comprising: presenting the estimated wavefront by forming a new plurality of traces for the respective plurality of sensors with the background noise suppressed.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) These and other aspects and features of the present embodiments will become apparent to those ordinarily skilled in the art upon review of the following description of specific embodiments in conjunction with the accompanying figures, wherein:

(2) FIG. 1 is a block diagram illustrating an example application of embodiments of the invention;

(3) FIG. 2 is a plot illustrating raw traces from a sensor array such as that shown in FIG. 1;

(4) FIG. 3 is a plot illustrating traces such as those shown in FIG. 2 that have been bandpass filtered according to embodiments;

(5) FIG. 4 is a functional block diagram illustrating an example system and method according to embodiments;

(6) FIG. 5 is a plot illustrating time aligned, oriented and windowed versions of traces such as those shown in FIG. 3 according to embodiments;

(7) FIG. 6 is a plot illustrating normalized singular values that are obtained from a matrix formed from trace data such as that shown in FIG. 5;

(8) FIG. 7 are plots illustrating various orders of waveform estimates generated from processed trace data according to embodiments;

(9) FIG. 8 are the plots from FIG. 7 shown as overlays for each of the orders;

(10) FIGS. 9 and 10 are plots illustrating clean traces that are generated using different orders according to embodiments; and

(11) FIG. 11 is a functional block diagram illustrating another example system and method according to embodiments.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

(12) The present embodiments will now be described in detail with reference to the drawings, which are provided as illustrative examples of the embodiments so as to enable those skilled in the art to practice the embodiments and alternatives apparent to those skilled in the art. Notably, the figures and examples below are not meant to limit the scope of the present embodiments to a single embodiment, but other embodiments are possible by way of interchange of some or all of the described or illustrated elements. Moreover, where certain elements of the present embodiments can be partially or fully implemented using known components, only those portions of such known components that are necessary for an understanding of the present embodiments will be described, and detailed descriptions of other portions of such known components will be omitted so as not to obscure the present embodiments. Embodiments described as being implemented in software should not be limited thereto, but can include embodiments implemented in hardware, or combinations of software and hardware, and vice-versa, as will be apparent to those skilled in the art, unless otherwise specified herein. In the present specification, an embodiment showing a singular component should not be considered limiting; rather, the present disclosure is intended to encompass other embodiments including a plurality of the same component, and vice-versa, unless explicitly stated otherwise herein. Moreover, applicants do not intend for any term in the specification or claims to be ascribed an uncommon or special meaning unless explicitly set forth as such. Further, the present embodiments encompass present and future known equivalents to the known components referred to herein by way of illustration.

(13) According to certain general aspects, embodiments of the invention provide techniques for extracting event arrival waveforms and estimating sensor wavefront features such as arrival times and polarizations, which techniques can operate in the presence of significant seismic background energy and which do not result in any systematic distortion of the arrival waveforms.

(14) FIG. 1 illustrates an example application of principles of the invention in a down-hole array of sensors.

(15) Referring to FIG. 1, consider a source 102 of transient microseismic energy, such as energy resulting from a fracture event, located at an underground position relative to a set of M sensors (e.g. geophones) 104 at respective locations, for example different separated locations along a cable suspended in an underground bore hole. As shown, the orientation or polarization of a wavefront from source 102 arriving at the mth sensor 104 is .sub.m. It should be noted that a sensor 104 can be comprised of a pod of sensor components, such as sensor components respectively configured to detect energy along three orthogonal axes, in which case the orientation would also have three components.

(16) According to certain aspects, the present inventors recognize that the wavefronts arriving at the sensors 104 from source 102 are often corrupted by significant seismic background energy and sensor noise, which obscures important event wavefront features. Processing of the sensor traces such as bandpass filtering or Wiener filtering is commonly performed to suppress frequency bands in which the event wavefront signal-to-noise ratio is poor (note that here noise refers to both the seismic background and sensor noise); however, in many cases such processing does not remove a sufficient amount of the seismic background to effectively understand details of the arriving wavefront. Other non-linear processing techniques, such as adaptive filtering methods, suppress noise effectively, but can alter the arriving waveforms in doing so.

(17) According to certain aspects, therefore, embodiments provide a seismic trace processing method which can remove a significant amount of seismic background during an event wavefront arrival, while preserving important features of the arriving waveforms.

(18) As in FIG. 1, consider an array of M sensors, each receiving signals g.sub.m(t), measured along one or more sensor components as a function of time t. The measured signal at sensor m is assumed to be an event waveform s.sub.m(t.sub.m) arriving at time .sub.m with polarization .sub.m in additive noise n.sub.m(t) consisting of seismic background and sensor noise,
g.sub.m(t)=s.sub.m(t.sub.m).sub.m.sup.T+n.sub.m(t)(1)
where .sup.T denotes a matrix transpose. The goal is to estimate the arriving event waveforms s.sub.m(t) at each sensor m from the measured seismic traces g.sub.m(t); m=1 to M. It should be noted that for P-wave arrivals, s.sub.m(t) and .sub.m will have one dimension or column, while for S-wave arrivals, s.sub.m(t) and .sub.m can have two dimensions or columns. For clarity in the following, s.sub.m(t) and .sub.m are treated as having one dimension or column, representing P-waves or S-waves polarized along a single direction. It will be apparent to those skilled in the art how to extend the invention described herein to the case of s.sub.m(t) and .sub.m having two dimensions or columns.

(19) An example set of raw measured traces is shown in FIG. 2. In this example, there are M=20 three-component sensors (i.e. each sensor m generates a set of three traces 202), and a P-wave arrival is present in a background of seismic energy. It should be noted that the set of trace data can be received as a computer readable file in a standard or proprietary format such as SEG-D or SEG-Y, perhaps depending on the sensor manufacturer, and having a defined time duration such as a half-second or ten seconds. However, the principles of the invention are not limited to these examples. A conventional approach is to bandpass filter the trace data according to the arrival spectrum and seismic background, as seen in FIG. 3. Here, the applied bandpass filtering has been adjusted for the arriving wavefront so that the background is suppressed at frequencies where the SNR is less than 0 dB. Note that the arrivals 304 of the wavefront at each sensor is more prominent in traces 302 after bandpass filtering. While the bandpass filtering suppresses some of the seismic background and the arriving wavefronts are more clearly visible, important waveform features remain obscured.

(20) A functional block diagram of an example method and system for implementing aspects of embodiments of the invention is shown in FIG. 4. It should be known that FIG. 4 can be implemented using any combination of hardware and software using techniques that are known to those skilled in the art after being taught by the present disclosure. For example, an example system can be implemented using an x86 computer system loaded with software that implements the functionality described in connection with FIG. 4 and having a conventional graphics system and display for generating plots such as those provided in the figures, as well as a conventional file system for reading and saving files containing trace data. Many alternatives as possible.

(21) Returning to FIG. 4, as shown in this example, the trace data that is received and processed is conventional bandpass filtered trace data such as that described above in FIG. 3. The filtered trace data is time aligned in block 402 according to an estimate of the arrival moveout so as to have the event wavefront arrive at all sensors simultaneously. More particularly, the respective arrival time .sub.m for each sensor m is estimated with respect to a common time reference for each of the filtered traces. These initial estimates of .sub.m can be performed using various techniques, such as those described in the co-pending applications. Then the trace data for each sensor is respectively shifted in time with respect to the common time reference by the respective time .sub.m.

(22) Additionally in block 402, if the sensors contain multiple components, the respective traces from these components are combined into a single trace for each sensor. In one example, using an estimate for each respective polarization .sub.m, the data for each sensor is pointed in the source direction, so as to combine the multiple traces for each sensor into a single component and single trace. The initial estimate for .sub.m (i.e. estimate of source location relative to the position of the sensor in the array) can be determined in various ways, such as those described in the co-pending applications.

(23) In another possible example, a single trace for each sensor is generated by pointing each sensor in the direction at which it achieves its maximum SNR. It should be noted that SNR may be estimated by examining the signal level before and during the event arrival, and that the direction of maximum SNR may be found by a search over arrival directions. It should be further noted that the direction at which the maximum SNR is achieved is not necessarily determined by the source direction. This is because while the source energy is maximized when the sensor is oriented according to the source direction, the maximum is rather broad. By contrast, energy from a noise source is minimized when the sensor is oriented perpendicular to the noise source polarization, but the minimum is relatively narrow. As a result, the SNR maximizing orientation will often be more aligned relative to a dominant noise source than according to the source direction.

(24) Finally, in block 402, the set of time-aligned, oriented single traces is time-windowed to select the duration of trace data over which the wavefront appears at each sensor. An example of how block 402 processes the traces of FIG. 3 to be time-aligned, oriented, and windowed as described above is shown in FIG. 5.

(25) Denote by h.sub.m(t) the time-aligned, oriented, and windowed sensor signals as produced by block 402, and evaluated at sample time t,
h.sub.m(t)=w(t).Math.g.sub.m(t+.sub.m).sub.m; m=1,2, . . . M(2)
where w(t) represents the time window, which is zero outside the T-sample long duration of the aligned arrivals. Denote by H.sub.m the column of T samples of h.sub.m(t),
H.sub.m=[h.sub.m(1)h.sub.m(2) . . . h.sub.m(T)].sup.T(3)
and by H the TM matrix of time-aligned, oriented, and windowed sensor signals,
H=[H.sub.1H.sub.2 . . . H.sub.M](4)
The columns of the arrival waveform matrix H consist of sampled event waveforms s.sub.m(t) in a background of noise n.sub.m(t+.sub.m).

(26) The event waveforms are expected to be similar sensor to sensor, and should be accurately represented by the weighted sum of a relatively small number of basis shapes. The noise, on the other hand, is expected to vary in an uncorrelated manner from sensor to sensor. In this way, embodiments decompose the arrival waveform matrix into signal and noise subspaces (for example, via singular value decomposition), and estimate the event waveforms at each sensor as that portion of the arrival waveform matrix H lying in the subspace associated with its few larger singular values.

(27) Denoting by S the matrix of arrival waveforms,
S=[s.sub.1s.sub.2 . . . s.sub.M](5)
where s.sub.m is the vector of sensor m event waveform samples s.sub.m(t); t=1, 2, . . . T, embodiments model the arrival waveforms as
=U.sub.LW.sub.L(6)
where U.sub.L is an L-column matrix of T-long orthonormal waveform basis shapes, and W.sub.L is an M-column matrix of L-element weights. The basis waveform shapes and weights (i.e. the matrices U and V) can be found via singular value decomposition of H as performed in block 404. In one non-limiting example, singular value decomposition is performed using the MATLAB function call svd(H, econ). The outputs of this function call using the matrix H formed using the outputs of block 402 can be expressed as:
UDV.sup.T=H(7)
where U is the TM matrix of orthonormal basis shapes, D is an MM diagonal matrix of singular values d.sub.l; l=1, 2, . . . M, arranged largest to smallest, and V is an MM orthonormal matrix of normalized weights.

(28) In forming the model waveform basis shapes U and sensor weights W, it is desired to take into consideration the relative sensor SNRs, in particular, it is desired to give more weight to those sensors having larger SNRs. A weighted decomposition may be formed as follows. Denote by the sensor noise covariance matrix, e.g. a diagonal matrix of sensor noise variances. The matrix H.sup.1 then weighs the columns of H by the respective sensor inverse variances. Performing a singular value decomposition on H.sup.1 gives (after multiplying through by ),
U.sub.D.sub.V.sub..sup.T=H,
with the orthonormal U.sub. being the waveform basis shapes, and D.sub.V.sub..sup.T being the sensor weights.

(29) Using these matrices, in block 406, to estimate the sensor waveforms during the arrival, the order L is selected, and Equation (6) is used to construct the estimated waveforms, with U.sub.L being the first L columns of U, and W.sub.L being the first L rows of the product DV. (Note that the singular values D could be apportioned in any manner between U and V.sup.T in computing the waveform estimates.) In some sense, reconstructs the event arrival waveforms without using the noise-like components.

(30) As an example, FIG. 6 shows the singular values of the sensor signal matrix H produced using samples of the waveforms shown in FIG. 3. As can be seen, there are a handful of larger singular values 602, with the rest being similarly small. The few larger singular values 602 are associated with the arriving event waveforms; the others are associated with noise.

(31) FIG. 7 shows the estimated sensor arrival waveforms formed from matrices U and W as described above using model orders L=1, 2, 5, 10 and 20 as a waterfall plot, and FIG. 8 shows the same estimated waveforms overlaid for each model order. Note that as the order decreases from the total number of sensors, L=M (which exactly reproduces the arrival waveform matrix H) to L=1, the noise energy is monotonically reduced, with the order 1 estimate .sub.1 having roughly 1/M of the noise energy in the original traces H. Also note that with a model order of L=1, the estimated waveform shape (but not its amplitude) is unchanging from sensor to sensor. As the model order is increased, differences in waveform shape sensor-to-sensor become better modeled. Accordingly, to pick the model order, it is preferred to select the smallest order L that appears to capture the arrival waveform changes across the array. This can be done by comparing the estimated waveforms and original sensor signals H. It can also be achieved by finding the order after which the singular value decreases significantly, i.e. the order L such that d.sub.L/d.sub.L+1 is relatively large.

(32) Returning to the example embodiment shown in FIG. 4, while estimates of the arriving waveforms as generated in block 406 are useful in determining arrival parameters, additional or alternative processing can be performed using the matrices U, D and V output from block 404. For example, embodiments generate cleaned versions of the original trace data that has been similarly processed to preserve event arrival waveforms while reducing background noise. Note that the singular value decomposition above produces a time basis U and a sensor basis V. To clean an entire seismic record containing the wavefront arrival, the entire record may be projected onto the sensor basis.

(33) First, in block 408, the first L columns of matrix V are obtained from the output of block 404.

(34) Next, in block 408, embodiments generate time-aligned, oriented traces (essentially h.sub.m(t), but not windowed),
k.sub.m(t)=g.sub.m(t+.sub.m).sub.m; m=1,2, . . . M(8)
and form the matrix K using k.sub.m(t) as its t, mth row-column element. In this way, the entire duration of the traces is processed, not just the time window used to find V.

(35) Then, in block 408, embodiments project the traces onto the first L columns of V, denoted V.sub.L to form the estimated traces , using the following equation:
=KV.sub.LV.sub.L.sup.T(9)
Note that the estimated traces will match the estimated arrival waveforms in the arrival time window.

(36) Finally in block 410, the estimated waveforms are time shifted back to their original moveout using .sub.m, and each sensor signal is rotated according to its polarization using .sub.m to generate multicomponent cleaned traces. As an example, processed versions of the raw trace data appearing in FIG. 2 are shown in FIG. 9, where the processed traces have been generated according to Equation (9) using V as constructed from the event arrival data H and with order L=2 (i.e. V.sub.2). Note the P-wave arrival 902 is clearly visible, and that the seismic background 904 is suppressed.

(37) Trace estimates made using L=5 are shown in FIG. 10. Compared with the L=2 estimate, arrival waveform changes 1002 across the array are more closely tracked using the L=5 estimate; however, the seismic background 1004 is less suppressed.

(38) One question concerns the nature of the waveform estimation errors in or . An expectation is that the estimation errors will be additive and noise-like as the estimates are constructed additively, and their errors are expected to result from including unwanted portions of additive noise in the signal basis or including noise subspace components by choosing an order L which is too large. Therefore, this process is not expected introduce systematic changes in the phase or magnitude of the arriving wavefronts. This seems to be the case, as evidenced by the similarity between the cleaned and bandpass filtered arrivals as the wavefront traverses the array in.

(39) According to certain additional aspects, the estimated arrival waveforms can be used to generate improved event arrival-time and polarization estimates for each sensor in embodiments of the invention. FIG. 11 is a functional block diagram illustrating an example method and system according to these embodiments.

(40) First, in block 1102, the estimated arrival waveforms are generated as described above (e.g. in connection with blocks 402, 404 and 406).

(41) Next, in block 1104, event arrival polarizations may be estimated using the projection of the estimated arrival waveforms S onto the corresponding sensor's multi-component trace data:

(42) $\begin{array}{cc}{}_m=\left[{g_m\left(1+{}_m\right)}^T\mathrm{.Math.}{g_m\left(T={}_m\right)}^T\right]{\hat{s}}_m& \left(10\right)\\ {}_m={}_m/.Math.{}_m.Math.& \left(11\right)\end{array}$
where .sub.m is the estimated arrival waveform at the mth sensor.

(43) Similarly, in block 1104, arrival times can be estimated by cross-correlating the estimated arrival waveforms with the corresponding sensor's trace data, preferably oriented according to the arrival polarization or the direction of maximum SNR, and finding the correlation lag at which the cross correlation is maximized. Note that these cross-correlations could also be performed using multi-component input trace data and the multi-component signal estimates made using the method described below directly, without reducing the data to a 1-dimensional signal per sensor.

(44) The case of planar motion (two-dimensional s.sub.m(t) and m) may be accommodated by first solving for the first column of .sub.m1, then subtracting that direction .sub.m1.sup.T from g.sub.m, and finally repeating the 1-D process with that first dimension of motion removed.

(45) Improved event arrival waveform estimates may be found iteratively by first using the process above to estimate the arriving waveforms in block 1102, and then using those waveforms to form new arrival-time and polarization estimates in block 1104. As shown in FIG. 11, the new arrival times and polarizations are then used to generate a new set of arrival waveform estimates, and so on, for example, until the process converges, or for a given number of iterations.

(46) According to still further aspects, another embodiment of the invention forms the arrival waveform matrix without first combining multi-component sensor signals to produce a reduced-dimension signal for each sensor such as a single, directionally oriented trace for each sensor. Denote by .sub.m(t) the time-aligned and windowed component signals of sensor m, evaluated at sample time t,
.sub.m(t)=w(t).Math.g.sub.m(t+.sub.m); m=1,2, . . . M(12)
where w(t) represents the time window, which is zero outside the T-sample long duration of the aligned arrivals. For a collection of sensors each having N components, denote by the TNM matrix of time-aligned and windowed multi-component sensor signals,

(47) $\begin{array}{cc}\hat{H}=\left[\begin{array}{cccc}{\tilde{h}}_1\left(1\right)& {\tilde{h}}_2\left(1\right)& \mathrm{.Math.}& {\tilde{h}}_M\left(1\right)\\ \mathrm{.Math.}& \mathrm{.Math.}& \mathrm{.Math.}& \mathrm{.Math.}\\ {\tilde{h}}_1\left(T\right)& {\tilde{h}}_2\left(T\right)& \mathrm{.Math.}& {\tilde{h}}_M\left(T\right)\end{array}\right]& \left(13\right)\end{array}$
where the multi-component arrival waveform matrix is N-times as wide as the oriented arrival waveform matrix H.

(48) To estimate the event arrival waveforms, the multi-axis waveform matrix is split into signal and subspace matrices. As above, a singular value decomposition can be used,
{circumflex over ({dot over (U)})}{tilde over (V)}.sup.T={hacek over (H)}(14)
with {circumflex over ({dot over (U)})}, and {tilde over (V)} being the orthonormal basis shapes, diagonal matrix of singular values (arranged largest to smallest), and orthonormal matrix of normalized weights. As above, the multi-component arriving wavefront signal may be estimated using a given number of singular values L,
={circumflex over ({dot over (U)})}.sub.L{acute over (W)}.sub.L(15)
where {circumflex over ({dot over (U)})}.sup.L is the matrix made from the first L columns of {circumflex over ({dot over (U)})}, and {acute over (W)}.sub.L is the first L rows of {circumflex over (V)}.sup.T. Similarly, any duration of trace data may be projected onto the columns of {circumflex over (V)}.sub.L (as above) to estimate the arrival signal from an entire record.

(49) Note that the N columns of {acute over (W)}.sub.L associated with a given sensor are likely to have a rank which is the lesser of N and L. Accordingly, when L>1, the different waveform shapes specified by {circumflex over ({dot over (U)})}.sub.L will be added in different proportions to estimate the arrival waveform along each component of a sensor. However, a P-wave arrival has motion along only the arrival direction; an S-wave arrival can involve motion only in the plane perpendicular to the arrival direction. As a result, it is preferred to replace the columns of {acute over (W)}.sub.L associated with each sensor by either a rank-one or rank-two approximation, depending on the arrival type, for instance found by singular value decomposition of each sensor's LN portion of {acute over (W)}. Note that the row space of this rank-one or rank-two approximation can be used as an estimate of the arrival polarization. Finally, note that this procedure may be applied in the case that different sensors within the array have varying number of components.

(50) Although the present embodiments have been particularly described with reference to preferred ones thereof, it should be readily apparent to those of ordinary skill in the art that changes and modifications in the form and details may be made without departing from the spirit and scope of the present disclosure. It is intended that the appended claims encompass such changes and modifications.