Wave-Signal-Processing System With a Two-Dimensional Recursive Digital Filter

My present invention relates to a digital system for processing wave signals propagating through a homogeneous medium. More particularly, my invention relates to a digital system employing two-dimensional recursive filters for detecting actual or reflective sources of shock waves, such as seismic tremors, in a homogeneous medium such as the soil.

Digital electronic systems for transmitting static or moving images or for processing acoustical recordings, radar echoes and seismic signals necessarily operate on signals defined in a two-dimensional domain, i.e. in time and space. Such signal-processing systems may include recursive filters, which considerably reduces the amount of computations while also providing a high application flexibility with respect to nonrecursive-filtering techniques. To augment processing efficiency, the pulse responses of the recursive filters are defined only in a half-plane of the space-time domain; the restriction of the pulse response does not preclude filter operation according to any specified frequency-response function. Although recursive filtering in a half-plane of the space-time domain is theoretically possible, conventional two-dimensional filters obtain good approximations only in a quadrant of the domain area. In the processing of images, constraining the signal response to a single quadrant achieves acceptable noise reduction or image restoration, but this technique is not suitable for the study of wave propagation through a homogeneous medium: a wavefront generated by a source in such a medium produces in a linear array of sensors a response defined by a hyperbola of the space-time domain, while a recursive filter with a pulse response defined only in a single quadrant is generally incapable of processing hyperbolically arrayed input data. In lieu of two-dimensional recursive filters effective in a half-plane, separate filtering has been proposed wherein the time variable is filtered recursively and the space variable is processed in a transversal filter. Such separate processing is likely to approximate the hyperbolic pulse response with a parabolic one, the symmetry of the response being preserved and times being roughly equivalent to those obtained with a two-dimensional recursive filter of like order. However, the simplification inherent in a parabolic approximation is too drastic to produce reliable results even by filters of a very high order.

The object of my present invention is to provide a reliable system for pinpointing, with good accuracy, the location of a source of shock waves such as the epicenter of a seismic tremor.

A processing system for determining the location of a shock-wave source in a homogeneous medium comprises, according to my present invention, a multiplicity of equispaced sensors disposed in a generally linear array for detecting magnitudes of a wave field generated by the source. Samplers operationally connected to the sensors extract pulses therefrom at successive sampling intervals equal to the quotient of the distance between consecutive sensors divided by the average wave velocity in the homogeneous medium. A two-dimensional phase-shifting filter connected to the samplers generates during each sampling interval, on the basis of sampling pulses extracted from the sensors, output signals approximating wave-field values at progressively increasing distances in depth from the sensor array. A two-dimensional attenuator at the output of the filter damps signals having a temporal frequency component less in magnitude than the product of an associated spatial frequency component times the aforementioned wave velocity, while a two-dimensional equalizer is linked to the attenuator for ensuring a global linear response thereof.

According to another feature of my present invention, a plurality of two-dimensional filters connected in cascade for estimating field values in incremental steps from the sensor array together approximate wave-field values at a distance from the sensor array equal to the sum of the incremental steps.

According to another feature of my invention, the filters are of the recursive digital type with multipliers and adders interconnected for implementing the following relations: ##EQU1## where p and q are integers respectively specifying the sampling intervals and identifying the sensors, functions ##EQU2## and f.sub.r (p,q) are input and output signals of any given filter, upper summation limit M represents the order of the respective filter, function v(p,q) is an intermediate processing signal generated by the respective filter in the p.sup.th sampling interval for the q.sup.th sensor, and real numbers a.sub.i,k are multiplication coefficients of the respective filter and satisfy the conditions a.sub.i,k =a.sub.k,i, a.sub.i,k =O for i+k odd, and a.sub.O,O =1.

According to yet another feature of my present invention, the attenuator and the equalizer are two-dimensional recursive filters, the attenuator including multipliers and adders for implementing the following relations: ##EQU3## where time and space coordinates p and q have the same significance as above, functions g(p,q) and g'(p,q) are input and output signals of the attenuator, upper summation limit N represents the order of the attenuator, function v'(p,q) is an intermediate processing signal generated by the attenuator in the p.sup.th sampling interval and associated with the q.sup.th sensor, and real numbers b.sub.i,k, c.sub.i,k are multiplicative coefficients satisfying the conditions b.sub.i,k =b.sub.k,i, c.sub.i,k =c.sub.k,i, b.sub.i,k =c.sub.i,k =0 for i+k odd, and c.sub.O,O =1.

Pursuant to another feature of my present invention, each filter includes a read/write memory storing, during each sampling interval and for assignment to a respective sensor, an intermediate signal generated during such sampling interval and M additional intermediate signals generated during M immediately preceding sampling intervals. The read/write memory also includes positions for storing initial and boundary data utilized in processing signals received during an initial sampling interval and signals assigned to peripherally disposed sensors. The adders and multipliers in each filter are connected to the read/write memory for receiving therefrom, during each sampling interval, intermediate signals generated during the M preceding sampling intervals and assigned to a respective sensor and intermediate signals generated during the preceding sampling interval and assigned to sensors adjacent this respective sensor.

These and other features of my present invention will now be described in detail, reference being made to the accompanying drawing in which:

FIG. 1A is a graph of a space-time domain of wave samples processed by a filter system according to my present invention;

FIG. 1B is a graph of a frequency domain of samples occurring in the space-time domain of FIG. 1A;

FIG. 2 is a block diagram showing a processing system according to my present invention;

FIG. 3 is a block diagram of a second-order pulse filter shown in FIG. 2;

FIG. 4 is a block diagram of a portion of the filter shown in FIG. 3;

FIG. 5 is a block diagram of a damping filter shown in FIG. 2;

FIG. 6 is a graph showing an input-data pattern for the filter system of FIG. 2; and

FIG. 7 is an output-data pattern for a series of pulse filters shown in FIG. 2.

According to my present invention, a multiplicity of sensors SE.sub.1, SE.sub.2, . . . SE.sub.n (FIG. 2) are separated by identical distances .DELTA.x along a line. As illustrated in FIG. 1A, a pulse u(p,q) extracted from a sensor SE.sub.q (q=1,2,..n) at time p .DELTA.t (p=0,1,2,..), where factor .DELTA.t is a constant sampling interval, is a function of a point (p,q) in a domain (t,x) having temporal and spatial coordinate axes t and x. A wavefront generated by a shock-wave source Q at a distance Z from sensor array SE.sub.1 -SE.sub.n (measured in a direction z transverse to line x) will give rise in the sensors to finite wave amplitudes located along a hyperbola A in the time-space domain (t, x), this hyperbola intersecting the t axis at a time coordinate T equal to the distance Z divided by an average wave velocity c in the medium between source Q and linear array SE.sub.1 -SE.sub.n. This medium, which in the case of a seismic wave will be the soil, is considered for present purposes to be substantially homogeneous. As described in detail hereinafter, a second-order filtering device according to my present invention estimates, at a time p.DELTA.t for a point (q,z) spaced a distance z from sensor SE.sub.q, a wave-field value u(p,q,z) at least partially in response to pulses u(p,q), u(p-1,q), u(p-2,q) extracted from sensor SE.sub.q at times p.DELTA.t, (p-1).DELTA.t, (p-2).DELTA.t and in response to pulses u(p-1,q-1), u(p-1,q+1) extracted from sensors SE.sub.q-1, SE.sub.q+1 at time (p-1).DELTA.t. A fourth-order filtering device according to my invention would operate on pulses extracted from a respective sensor SE.sub.q during a p.sup.th sampling interval and the four preceding sampling intervals and on pulses extracted from adjacent sensors SE.sub.q-1, SE.sub.q+1 during the three preceding sampling intervals (p-1,p-2,p-3) and on pulses extracted from two further sensors SE.sub.q-2, SE.sub.q+2 during the second interval (p-2) before the p.sup.th interval.

Pulses u(p,q) are extracted from sensors SE.sub.1 -SE.sub.n by respective samplers or signal-level extractors CA.sub.1, CA.sub.2, . . . CA.sub.n at intervals .DELTA.t=.DELTA.x/c. Such a sampling interval orients the asymptotes x.sub.1, x.sub.2 of hyperbola A in the t,x plane at a 45.degree. angle with respect to coordinate axes t,x, as shown in FIG. 1A. Asymptotes x.sub.1, x.sub.2 may be regarded as a second spatial coordinate system (x.sub.1, x.sub.2) related to the time/space domain (t,x) by a linear transformation, as described more fully hereinafter.

In FIG. 1B I have shown two pairs of coordinate axes .omega..sub.t, .omega..sub.x and .omega..sub.1, .omega..sub.2 defining two overlapping quadrants interrelated by a linear transformation and derived from quadrants (t,x) and (x.sub.1,x.sub.2) by respective Fourier transforms. As illustrated in FIG. 2, a system according to my present invention for determining the distance Z of shock-wave source Q from sensor array SE.sub.1 -SE.sub.n comprises a plurality of filters FF.sub.1, FF.sub.2, . . . FF.sub.R serving as phase shifters for frequencies in an area AA.sub.1 bounded by axes .omega..sub.1, .omega..sub.2 (FIG. 1B) and an attenuator filter FA which damps signals constituted by frequency pairs in adjoining octants forming an area AA.sub.2. Areas AA.sub.1 and AA.sub.2 are respectively defined by .vertline..omega..sub.t .vertline.>c.vertline..omega..sub.x .vertline.and .vertline..omega..sub.t .ltoreq.c.vertline..omega..sub.x .vertline..

Extractors CA.sub.1, CA.sub.2, . . . CA.sub.n transmit at intervals p.DELTA.t staggered pulses u(p,1), u(p,2), . . . u(p,n) from respective sensors SE.sub.1, SE.sub.2, . . . SE.sub.n to the cascaded filters FF.sub.1 -FF .sub.R, sensors SE.sub.1 -SE.sub.n being adapted to detect acoustic waves propagating through the intervening medium from shock-wave source Q. Filters FF.sub.1 -FF.sub.R are of the two-dimensional recursive digital type and have respective outputs successively approximating the wave field at progressively increasing distances from sensor array SE.sub.1 -SE.sub.n, as more fully described hereinafter. Filter FF.sub.R works into attenuator filter FA which is in turn linked to an equalizer filter EQ for ensuring a linear phase response thereof.

As shown in FIG. 3, a generic filter FF.sub.r (r=1,2, ..R) includes a digital adder S.sub.1 which emits on a lead 1, during the p.sup.th processing interval, a signal v(p,q) associated with sensor SE.sub.q. Signal v(p,q) is a spike fed to a binary multiplier ML.sub.4 for multiplication with a coefficient a.sub.2,2 and subsequent transmission to a digital adder S.sub.3 ; in addition to feeding the unit ML.sub.4, lead 1 supplies a time-delay unit T.sub.1 whose output lead 2 extends to a space-delay unit X.sub.1, to a space-advance unit X.sub.1.sup.-1, to a binary multiplier ML.sub.2 and to another time-delay unit T.sub.2. During the p.sup.th processing interval, lead 2 carries a pulse signal v(p-1,q) present on lead 1 during the (p-1).sup.th interval and delayed by a time .DELTA.t in unit T.sub.1. Upon receiving pulse v(p-1,q), space-delay unit X.sub.1 emits on a lead 4, extending to a digital adder S.sub.2, a pulse v(p-1) carried by lead 1 during the (p-1).sup.th processing interval and associated with the (q-1).sup.th sensor SE.sub.q-1. It is to be noted that, during the p.sup.th processing interval (p=0,1,2,..) for each sensor SE.sub.q (q=1,2,.. n), filter FF.sub.r processes an input signal f.sub.r-1 (p,q) to generate an output signal f.sub.r (p,q), pulse v(p,q) constituting an intermediate product. Thus, pulse v(p-1,q-1) is an intermediate signal in the processing of the output of the (q-1).sup.th sensor SE.sub.q-1.

During the p.sup.th processing interval, space-advance unit X.sub.1.sup.-1 receives pulse v(p-1,q) and transmits to adder S.sub.2 over a lead 5 a pulse v(p-1,q+1) whose magnitude is derived, at least in part, from the signal level extracted from sensor SE.sub.q +1, i.e. from a sensor adjacent sensor SE.sub.q in the array SE.sub.1 -SE.sub.n, during the (p-1).sup.th processing interval. The sum of signals v(p-1,q-1) and v(p-1,q+1) is transmitted from adder S.sub.2 to a binary multiplier ML.sub.1 for multiplication by a coefficient a.sub.0,2. Pulse v(p-1,q) is delayed by a time .DELTA.t in unit T.sub.2 ; this unit emits on a lead 3, working into adder S.sub.3 and a binary multiplier ML.sub.3, a pulse v(p-2,q) present on lead 1 two processing intervals prior to the p.sup.th interval. Transmitted to adder S.sub.1 from multipliers ML.sub.1, ML.sub.2, ML.sub.3 are, respectively, the product of coefficient a.sub.0,2 times the sum of pulses v(p-1,q-1), v(p-1,q+1), the product of a coefficient a.sub.1,1 times pulse v(p-1,q), and the product of a coefficient a.sub.2,2 times pulse v(p-2,q). Adder S.sub.1 subtracts the outputs of units ML.sub.1, ML.sub.2, ML.sub.3 from filter input signal f.sub.r-1 (p,q), thereby generating the processing intermediate v(p,q). The output signals of multipliers ML.sub.1, ML.sub.2, ML.sub.3 are also transmitted to adder S.sub.3 for summation with the product of coefficient a.sub.2,2 times pulse v(p,q) to produce filter output signal f.sub.r (p,q).

Time-delay units T.sub.1 and T.sub.2, space-delay unit X.sub.1 and space-advance unit X.sub.1.sup.-1 may be realized by a read/write memory RAM, as illustrated in FIG. 4. In addition to being tied to multiplier ML.sub.4, lead 1 works into memory RAM and into an address generator IN controlling the reading of the contents of that memory to a static register REG. Leads 2, 3 and 4,5 extend from register REG to multiplier ML.sub.2, to multiplier ML.sub.3 and adder S.sub.3, and to adder S.sub.2, respectively. Memory RAM includes a two-dimensional storage array having M rows and n columns, where M is the order of filter FF.sub.r and n is the number of sensors SE.sub.1 -SE.sub.n ; it also includes M storage positions for recording initial conditions of the system necessary for processing points on the contour of the discrete time-space domain (p,q), p=0,1,2 . . ., q=1,2 . . .n. Because the filter shown in FIG. 3 is of the second order, memory RAM has two rows.

Upon the completion of the p.sup.th processing cycle, i.e. at the end of the p.sup.th sampling interval, memory RAM has stored all the intermediate signals v(p,q) generated during this cycle, together with all the signals v(p-1,q) generated during the preceding or (p-1).sup.th cycle. Thus, at the onset of the p.sup.th cycle, signals v(p-1,q) and v(p-2,q) are stored in memory RAM for all the sensors SE.sub.q, q=1,2, . .n. Upon the transmission of a pulse f.sub.r (p,q) to adder S.sub.1, a signal on lead 1 induces address generator IN to read from memory RAM to register REG the signals v(p-1,q), v(p-2,q), v(p-1,q-1), v(p-1,q+1) which are then emitted on leads 2, 3, 4, 5, respectively. Signal v(p,q) is thereupon entered in memory RAM in the position theretofore occupied by signal v(p-2,q). Thus, in general, during the p.sup.th sampling interval an intermediate processing signal v(p,q) is generated for the q.sup.th sensor SE.sub.q (q=1,2,. .n) from the signals v(p-1,q) and v(p-2,q) derived during the preceding two cycles for the same sensor SE.sub.q and from the signals v(p-1,q-1) and v(p-1,q+1) produced during the preceding cycle for the two sensors SE.sub.q-1 and SE.sub.q+1 adjacent sensor SE.sub.q in the array SE.sub.1 -SE.sub.n. If q equals 1 or n, v(p-1,q-1) or v(p-1,q+1) take on pre-established values stored in memory RAM as initial boundary conditions. It is to be noted that the functions carried out by address generator IN, memory RAM and register REG (FIG. 4) are precisely those performed by time-delay units T.sub.1 and T.sub.2, space-delay unit X.sub.1 and space-advance unit X.sub.1.sup.-1 (FIG. 3).

Filter FF.sub.r implements the aforedescribed relationships: ##EQU4## where M=2 is the order of the filter, with a.sub.i,k =a.sub.k,i and with a.sub.i,k =0 if i+k is odd or zero. The process of deriving coefficients a.sub.i,k is described hereinafter with reference to FIGS. 1A, 1B and 2.

As illustrated in FIG. 5, attenuator FA comprises three adders S.sub.4, S.sub.5, S.sub.6, two time-delay units T.sub.3 and T.sub.4, seven binary multipliers ML.sub.5, ML.sub.6, ML.sub.7, ML.sub.8, ML.sub.9, ML.sub.10, ML.sub.11, a space-delay unit X.sub.2 and a space-advance unit X.sub.2.sup.-1. Adder S.sub.4 generates on a lead 11 an intermediate processing signal v'(p,q) from pulses emitted by multipliers ML.sub.9 -ML.sub.11 and from an attenuator input signal g(p,q) identical with an output signal f.sub.R (p,q) of filter FF.sub.R. Signal V'(p,q) is transmitted to multiplier ML.sub.5 for forming a product with a coefficient b.sub.0,0, this product being fed to adder S.sub.6. Besides being connected to multiplier ML.sub.5, lead 11 extends to unit T.sub.3 where signal v'(p,q) is delayed by a time .DELTA.t. Upon receiving signal v'(p,q), delay unit T.sub.3 emits on a lead 12 a pulse v'(p-1,q) which was carried by lead 11 during the preceding (p-1).sup.th processing cycle. Lead 12 conducts pulse v'(p-1,q) to multipliers ML.sub.7 and ML.sub.10, to time-delay unit T.sub.4, to space-delay unit X.sub.2 and to space-advance unit X.sub.2.sup.-1. Upon receiving pulse v'(p-1,q), unit T.sub.4 generates on a lead 13, working into multipliers ML.sub.8 and ML.sub.11, a signal v'(p-2,q) transmitted on lead 11 during the (p-2).sup.th processing cycle and delayed by a combined time 2.DELTA.t in units T.sub.3, T.sub.4. Pulses v'(p-1,q-1) and v'(p-1,q+1) are fed via leads 14 and 15 to adder S.sub.5 by space-delay unit X.sub.2 and space-advance unit X.sub.2.sup.-1,respectively, upon the reception thereby of pulse v'(p-1,q). Pulses v'(p-1,q-1) and v'(p-1,q+1) represent intermediate products in the processing of sampler output signals u(p-1,q-1) and u(p-1,q+1); adder S.sub.5 sums pulses v'(p-1,q-1) and v'(p-1,q+1) and feeds the results to multipliers ML.sub.6 and ML.sub.9 over a lead 16. Multipliers ML.sub.9, ML.sub.10 and ML.sub.11 form the products of a coefficient c.sub.0,2 times pulse sum v'(p-1,q-1)+v'(p-1,q+1), of a coefficient c.sub.1,1 times pulse v'(p-1,q), and of a coefficient c.sub.2,2 times signal v'(p-2,q), respectively, and transmit these products to adder S.sub.4 for summation with signal g(p,q) to produce signal v'(p,q). The sum v'(p-1,q-1)+v'(p-1,q+1) emitted by adder S.sub.5 on lead 16 is multiplied in unit ML.sub.6 by a weighting coefficient b.sub.0,2 and then transmitted to adder S.sub.6. Units ML.sub.7 and ML.sub.8 multiply pulses v'(p-1,q) and v'(p-2,q) by weighting coefficients b.sub.1,1 and b.sub.2,2 and feed the resulting products to adder S.sub.6. This adder sums the output signals of multipliers ML.sub.5 -ML.sub.8 to form the signal g'(p,q) which is transmitted to equalizer EQ. The latter may be a two-dimensional recursive digital filter having the structure shown in FIG. 3.

Attenuator FA may be designed as a read/write memory, as heretofore described with reference to FIG. 4, and implements the following relationships: ##EQU5## where b.sub.i,k =c.sub.i,k =0 if i+k is odd or zero and where b.sub.i,k =b.sub.k,i and c.sub.i,k =c.sub.k,i. The method of determining coefficients b.sub.i,k, C.sub.i,k is described hereinafter with reference to FIGS. 1A, 1B and 2.

Under the aforementioned assumption that the space between source Q (FIG. 1A) and sensor array SE.sub.1 -SE.sub.n is occupied by a homogeneous medium characterized by an average acoustical velocity c, and with disregard of multiple reflections in the medium, a pulse or wave packet generated by source Q at a distance Z from linear array SE.sub.1 -SE.sub.n creates a wavefront whose arrival at array SE.sub.1 -SE.sub.n is plotted along hyperbola A in the space-time domain (t,x) whose spatial dimension x and temporal dimenson t are discretized according to intersensor distance .DELTA.x and sampling interval .DELTA.t, respectively, as shown in FIG. 1A. Furthermore, because the problems of studying a wave field u(t,x,z) generated by source Q are reduced to linear ones under these assumptions, the principle of superposition may be applied, thereby allowing the utilization of the techniques of two-dimensional linear digital filtering, and particularly recursive filtering, to approximate the magnitude of the wave field u(t,x,z) at predetermined distances from the sensor array, SE.sub.1 -SE.sub.n.

Wave field u(t,x,z) satisfies the equation for a two-dimensional hyperbolic wave: ##EQU6## where variations along a third dimension y of Cartesian space (x,y,z) are neglected. The wave equation (3) may also be written in the form: ##EQU7## where function U(.omega..sub.t,.omega..sub.x,z) is the two-dimensional Fourier transform, in variables x and t, of the wave field u(t,x,z): ##EQU8##

The general solution of equations (3) and (4) in the frequency domain (.omega..sub.t,.omega..sub.x) of FIG. 1B with distance parameter z is: ##EQU9## Because waves traveling only in the negative direction of the z axis are to be considered (see FIG. 1A), coefficient A.sub.0 is set equal to 0 and general solution (6) takes on the form: ##STR1## where U(.omega..sub.t,.omega..sub.x, 0) is the Fourier transform of u(t,x,0) which is the boundary condition for wave equation (3). Wave field u(t,x,z) may be obtained from function U(.omega..sub.t,.omega..sub.x,z) by the inverse Fourier transform or by filtering the input data or boundary function u(t,x,0) with the processing system heretofore described with reference to FIGS. 2-5.

Let us define a frequency function H.sub.z (.omega..sub.t,.omega..sub.x) by the equation: ##EQU10##

Then a field value u(t,x,z) may be obtained from boundary or initial datum u(t,x,O) by convolving the same with a "migration" function h.sub.z (t,x):

where function h.sub.z (t,x) has dominant values along hyperbola A (FIG. 1A), has a predetermined spatial and termporal decay and is the inverse Fourier transform of frequency function H.sub.z (.omega..sub.t,.omega..sub.x). Equation (9) may be written: ##EQU11## or approximated by a discrete sum: ##EQU12## where .tau.=-.infin. . . . -3,-2,-1,0,1,2,3, . . . +.infin. and .sigma.=0,1,2,3 . . . .infin.. The wave-field value u(p,q,z+.DELTA.z) at a point in the homogeneous propagation medium spaced by a distance .DELTA.z from a point (p,q,z) along a perpendicular to linear array SE.sub.1 -SE.sub.n may be obtained from field value u(p,q,z) by the equation: ##EQU13## Thus, field values at a distance z+.DELTA.z may be estimated by summing field values u(.sigma.,.tau.,z) of a plane parallel to plane t,x. For an M.sup.th -order filter, equation (9c) may be rewritten in the form: ##EQU14## where f.sub.r (p,q)=u(p,q,z+.DELTA.z) and f.sub.r-1 (p,q)=u(p,q,z). Equation (9d) is fully equivalent to equations (1a), (1b) and, for a second-order filter (FIG. 3), may be reduced to: ##EQU15## where f.sub.r (p,q)=u(p,q,z+.DELTA.z) is the output of filter FF.sub.r during the p.sup.th processing interval for the q.sup.th sensor SE.sub.q, f.sub.r-1 (p,q)=u(p,q,z) is the output of filter FF.sub.r-1 during the p.sup.th interval for the q.sup.th sensor SE.sub.q, v(p-1,q-1), v(p-1,q), v(p-1,q+1) and v(p-2,q) are intermediate processing signals generated by the r.sup.th filter FF.sub.r, as heretofore described with reference to FIGS. 3 and 4, and a.sub.0,2, a.sub.1,1 and a.sub.2,2 are multiplier coefficients implemented by units ML.sub.1 -ML.sub.4 (FIG. 3).

As heretofore described, field u(t,x,z) may be obtained from function U(.omega..sub.t,.omega..sub.x,z) by applying the inverse Fourier transform. To determine the values of function U(.omega..sub.t,.omega..sub.x,z), a two-dimensional filter is needed whose phase and amplitude response is specified by frequency function H.sub.z (.omega..sub.t,.omega..sub.x) as per equation (8). Function H.sub.z (.omega..sub.t,.omega..sub.x) may be approximated within area AA.sub.1 (FIG. 1B) by a phase-shifting filter FF.sub.r having the frequency response: ##EQU16## and within area AA.sub.2 by an attenuator or fan-type filter FA with the frequency response: ##EQU17## A filter with the frequency response defined by equations (10) and (11) has a pulse response maximized along hyperbola A (FIG. 1A) which is defined by the equation:

Hyperbola A generally lies in a half-plane t.gtoreq.0 and has asymptotes t=.+-.(x/c). However, no convenient methods are available for designing two-dimensional recursive digital filters satisfying a given phase specification in the frequency domain. According to my present invention, it is possible to limit the response of a 2-D recursive filter to a single quadrant by normalizing the sampling intervals .DELTA.t and by rotating the coordinate axes t, x through a 45.degree. angle. Variables t and x are transformed according to the equations: ##EQU18## while sampling interval .DELTA.t is set equal to quotient .DELTA.x/c, as heretofore described with reference to FIG. 2. As illustrated in FIG. 1A, the new coordinate axes x.sub.1, x.sub.2 are superimposed upon the asymptotes of the hyperbola A whose equation (12) is transformed into:

Thus, the pulse response is confined to the quadrant (x.sub.1 .gtoreq.0 and x.sub.2 .gtoreq.0) and conventional techniques can be used to design the 2-D recursive phase filter. The sampling interval .DELTA.x.sub.1 =.DELTA.x.sub.2 is proportional to the real-time sampling interval .DELTA.t by the factor ##EQU19##

Upon the transformation of coordinates t and x into coordinates x.sub.1 and x.sub.2, equation (8) is converted into the form: ##EQU20## while equations (10) and (11) are transformed into:

and

A phase-shifting M.sup.th -order filter FF.sub.r having the frequency response specified by equation (16) must have a transfer function D(y.sub.1,y.sub.2) in complex variables y.sub.1 =exp(-j.omega..sub.1) and y.sub.2 =exp(-j.omega..sub.2): ##EQU21## where A.sup.-1 (y.sub.1,y.sub.2) and A(y.sub.1,y.sub.2) are polynomials of the type: ##EQU22## In order that transfer function D(y.sub.1,y.sub.2) optimally approximate frequency response P.sub.z (.omega..sub.1,.omega..sub.2), polynomial coefficients a.sub.i,k must satisfy the constraining equations:

and

Condition (21) ensures symmetry of filter response with respect to time variable t, while constraint (22) arises from the requirement that the filter coefficients match the rectangular sampling raster in the domain (x,t) of FIG. 1A after 45.degree. rotation. Coefficients a.sub.i,k having particular numerical values determined by an iterative method described hereinafter are the same as the identically labeled coefficients of equations (9d) and (9e) and the weighting coefficients multiplied with intermediate signals v(p,q) by units ML.sub.1 to ML.sub.4 (FIG. 3). Sample values for coefficients a.sub.i,k are listed for a fourth-order filter in the following Table 1.

According to my present invention, attenuating filter FA (FIG. 2) has a frequency response F.sub.z (.omega..sub.1,.omega..sub.2)--see equation (17)--approximated by a transfer function K(y.sub.1,y.sub.2) in complex variables y.sub.1, y.sub.2 : ##EQU23## where B(y.sub.1,y.sub.2) and C(y.sub.1,y.sub.2) are polynomials having the forms: ##EQU24## Coefficients b.sub.i,k and c.sub.i,k are the multiplicative factors formed into products with intermediate signals v'(p,q) by units ML.sub.5 to ML.sub.11 --see FIG. 5 and equation (2)--and satisfy symmetry conditions b.sub.i,k =b.sub.k,i, c.sub.i,k =c.sub.k,i for i=0,1,2,3, . . . and k=0,1,2,3, . . . and constraints b.sub.i,k =0, c.sub.i,k =0 for i+k odd. Sample values for coefficients b.sub.i,k, c.sub.i,k are given for a fourth-order attenuator in the following Tables 2 and 3, respectively.

Thus, filters FF.sub.r (FIGS. 2 and 3) approximate function P.sub.z (.omega..sub.1,.omega..sub.2) in a shaded triangular area shown in FIG. 1B, this area being bounded by the positive halves of axes .omega..sub.1,.omega..sub.2 and by the line .omega..sub.t =.pi.. It is to be noted that parameters .omega..sub.1,.omega..sub.2 represent coordinate axes rotated by 45.degree. from axes .omega..sub.t,.omega..sub.x of the domain (.omega..sub.t,.omega..sub.x, z) of Fourier-transform function U(.omega..sub.t,.omega..sub.x, z): ##EQU25## where upper limit P is the number of samples along time axis t (FIG. 1A) while 2Q+1 is equal to the number n of samples along the spatial dimension x. Wave numbers f.sub.x and wave frequencies f.sub.t detectable by the processing system of FIG. 2 range from zero to half the sensor density f.sub.sx and half the sampling rate f.sub.st, respectively, where the sensor density f.sub.sx is the reciprocal of the intersensor interval .DELTA.x and the sampling rate f.sub.st is the quotient of the seismic velocity c divided by distance .DELTA.x. Parameters .omega..sub.t,.omega..sub.x are related to wave frequency f.sub.t and wave number f.sub.x by the equations: ##EQU26## It is clear that the maximum detectable value of variables .omega..sub.t and .omega..sub.x is .pi.; thus, the line .omega..sub.t =.pi. is the lower limit of the shaded area in FIG. 1B which is the set of all parameter pairs .omega..sub.t, .omega..sub.x or .omega..sub.1,.omega..sub.2 for which a response specified by equation (10) or (16) must be given by a filter FF.sub.r (r=1,2,3, . . . R).

A method of designing a two-dimensional all-pass recursive filter FF.sub.r includes minimization of a mean squared error J(a .sub.i,k). Let frequency-response function P.sub.z (.omega..sub.1,.omega..sub.2) of equation (16) be separated into real and imaginary parts and evaluated at fixed parameter (.omega..sub.1m,.omega..sub.2m):

where .phi..sub.m =-z.sqroot.2.omega..sub.1m .omega..sub.2m. Two-dimensional all-pass phase filter FF.sub.r (r=1,2, . . . R) has a complex frequency response approximated by the transfer function D(y.sub.1m,y.sub.2m) at parameters .omega..sub.1m,.omega..sub.2m : ##EQU27## where .theta..sub.m =-M(.omega..sub.1m +.omega..sub.2m)+2<A(exp[-j.omega..sub.1m ], exp[-j.omega..sub.2m ]).

Mean-squared-error function J(a.sub.i,k) is formed from the square of the difference between equations (28) and (29), summed over the number m of predetermined parameter pairs (.omega..sub.1m,.omega..sub.2m): ##EQU28##

As described by C. K. Sanathanan and H. Tsukai in a paper entitled "Synthesis of Transfer Function From Frequency-Response Data" (International Journal of System Science, Vol. 5, No. 4, pp. 41-54, 1974), this is a typical nonlinear problem which can be converted into linear form by using a weighted error function J'(a.sub.i,k) in an iterative process: ##EQU29## where A.sub.m,l represents the value at parameters .omega..sub.1m,.omega..sub.2m of polynomial A(exp[-j.omega..sub.1m ],exp[-j.omega..sub.2m ]) whose coefficients are to be determined in the l.sup.th step of the iteration, complex number A.sub.m,l-1 similarly representing the value at parameters .omega..sub.1m,.omega..sub.2m of the polynomial A(exp[-j.omega..sub.1m ],exp[-j.omega..sub.2m ]) according to equation (20) whose coefficients have been determined at the (l-1).sup.th step of the iteration. The condition necessary for error function J(a.sub.i,k) to be minimized is satisfied by equating to zero the first-order partial derivatives of function J'(a .sub.i,k) with respect to the unknown coefficients a.sub.i,k at each step l of the iteration, which entails the solution of a linear system of equations:

where:

Matrix coefficients S.sub.s,w.sup.(l) are easily evaluated: ##EQU30## where: ##EQU31## and the constant term is given by:

Sum (34) is carried out on the fixed set of parameter pairs .omega..sub.1m,.omega..sub.2m in the approximation domain AA.sub.1 (FIG. 1B). Frequency-response function P.sub.z (.omega..sub.1m,.omega..sub.2m) of equation (28) must take into account the linear term due to the chosen order M of filter FF.sub.r.

I have used the method of error minimization to design recursive digital filters with linear phase, as described in my paper "On 2-D Recursive Digital Filter Design with Complex Frequency Specifications," published in Proceedings ICC 1977 (pages 9B-2-202 to 206, June 1977). It is to be noted that in designing all-pass filter FF.sub.r the coefficients a.sub.i,k chosen initially do not seriously affect the final approximation in the iteration. Convergence of the filter coefficients a.sub.i,k in the iteration corresponds to the minimization of the original unweighted error J(a.sub.i,k) as per equation (30). Furthermore, stability must be tested at each step and a stabilization strategy developed in order to guarantee a meaningful result. Owing to the limited order of the chosen two-dimensional one-quadrant all-pass recursive digital filters, these problems can be overcome satisfactorily with techniques described by T. S. Huang in the paper "Stability of Two-Dimensional Recursive Filters" published in IEEE Transactions on Audio Electroacoustics AU-20, No. 2, pages 158-163, June 1972.

Let us consider the problem of determining, according to equation (9), wave field u(t,x,z) at a distance z from sensor array SE.sub.1 -SE.sub.n. A phase function .phi..sub.z (.omega..sub.1,.omega..sub.2) to be approximated by .rho. all-pass filters FF.sub.1 -FF.sub..rho. (.rho..ltoreq.R) in the area AA.sub.1 of the rotated frequency domain (.omega..sub.1,.omega..sub.2) is: ##EQU32##

In order to minimize error function J(a.sub.i,k), it is necessary to connect in cascade filters FF.sub.1 -FF.sub..rho. each estimating a change in wave field u(t,x,z) over an incremental variation .DELTA.z in the distance from sensors SE.sub.1 -SE.sub.n. Filters FF.sub.1 -FF.sub..rho. each have a transfer function D.sub..DELTA.z (e.sup.-j.omega.1,e.sup.-j.omega.2) which together yield a transfer function D.sub.z (exp[-j.omega..sub.1 ],exp[-j.omega..sub.2 ]) according to the equation:

where function D.sub..DELTA.z (exp[-j.omega..sub.1 ],exp[-j.omega..sub.2 ]) approximates the phase response: ##EQU33## Taking the limit of the transfer function D.sub.z (exp[-j.omega..sub.1 ],exp[-j.omega..sub.2 ]) as distance increment .DELTA.z approaches 0 and filter number .rho. increases toward infinity, the product z=.rho..DELTA.z remaining constant, we obtain: ##EQU34## Thus, approximate operator D.sub..DELTA.z for the extension of the field through a short or incremental distance .DELTA.z is derived as a two-dimensional recursive filter.

Table 4 shows the performance index for a fourth-order filter of the phase operator D.sub..DELTA.z for decreasing values of the distance increment .DELTA.z.

Each term in this table has been determined by using a mask of coefficients a.sub.i,k for which the intermediate signal v(p,q) (FIG. 3) is obtained from a multiplicity of previously determined intermediate signals v(p-1,q-1), v(p-1,q), v(p-1,q+1) . . . whose space-time-domain points are located in a square with sides parallel to axes x.sub.1, x.sub.2, as shown in FIG. 1A. In the case of the fourth-order filter, only eight independent coefficients had to be optimized or derived via iteration, owing to the symmetry constraint a .sub.i,k =a.sub.k,i. Phase function .phi..sub..DELTA.z (.omega..sub.1,.omega..sub.2) was approximated and sum (34) carried out on a polar grid of 100 parameter pairs (.omega..sub.1m,.omega..sub.2m) in the first quadrant of the frequency domain (.omega..sub.1,.omega..sub.2) with the limit: ##EQU35## To determine the wave field u(t,x,z) at a distance z from sensor array SE.sub.1 -SE.sub.n it is necessary to utilize filters FF.sub.1 -FF.sub..rho. which approximate the field in small steps away from the sensor array. However, the squared error J(a.sub.i,k) of equation (30) increases approximately with the power of .rho.; the optimal choice for the elemental step .DELTA.z must be a compromise between the limit of equation (40) and the allowed maximum number of operations per sample.

A fourth-order digital filter was used to "migrate" the input data or response pattern u(t,x,O) shown in FIG. 6 and recorded at a surface z=O upon the explosion of a source Q at a distance z=20.DELTA.z from sensor SE.sub.0. Such a "migration" would have to produce an increasingly spiked pattern for the sensor SE.sub.0. The result of filtering the pattern u(t,x,O) shown in FIG. 6 in the opposite direction along axis t with 20 fourth-order all-pass phase operators or filters FF.sub.1 -FF.sub.20 connected in cascade is shown in FIG. 7. This allows the determination of the position of the source Q along array SE.sub.1 -SE.sub.n and the distance of the source from the array.