Various systems and methods of the present invention provide amorphous computing systems and methods for use thereof. In some cases, the amorphous computing systems include one or more amorphous hardware elements that are programmed under control of a computer processor such as, for example, an AMD™ or INTEL™ based personal computer. Portions of an algorithm can then be computed by the individual amorphous hardware elements that can be, as one example, an FPGA or some subset of hardware gates available on an FPGA. Methods can include a variety of computations including Wave Equations and Kirchhoff algorithms. Thus, in particular cases, the amorphous computing systems and methods for using such are applied to seismic imaging problems, as well as other problems.

Current U.S. Class:	702/17; 703/5
Current Intern'l Class:	G06F 19/00    (20060101)
Field of Search:	702/17,14,1,12,2 324/158.1 716/12 712/11 367/43 327/551 703/5,10

• Dp=8 *sqrt[(Sa*X_—distance_—spacing*X_—distance_—spacing)+(Sb* Y_—distance_—spacing*Y_—distance_—spacing), where Sa is the lateral distance in the x direction between the input trace and the image point and Sb is the lateral distance in the y direction between the input trace and the image point, and where X_—distance_—spacing is the distance between input trace rows in the x dimension and Y_—distance_—spacing is the distance between input trace rows in the y dimension.
• Kp=Dp*Dt/(Vrns*Time), where Vrms is the interpolated velocity provided from the aforementioned stage "A" of the processing pipeline, Time is the calculated time of the input trace from the aforementioned stage "B" of the processing pipeline that contributes to the image point, and Dt is the sample rate of the input trace typically measured in milliseconds.
  
  Filter Frequency=512/(2*Dt*Kp)
  
  Based on this disclosure, one of ordinary skill in the art will appreciate other functions that can be implemented to serve the function of F(Vrms, Sa, Sb). The following references provide additional information about anti-aliasing in relation to a Kirchhoff migration: Anti-aliased Kirchhoff 3-D migration, D. Lumley, John Claerbout, Stanford University, Stanford Exploration Project pp. 1282-1285; and Frequency Selective design of the Kirchhoff migration operator, Gray, S. H. 1992, Geophysical Prospecting, 40, pp. 565-571. The referenced portions of the foregoing publications are included in Appendix A attached hereto.

Based on the calculated Filter Frequency, a set of coefficients for a high frequency filter are accessed. These coefficients are pre-calculated by the program before any of the input traces are processed through the processing pipeline. The following function implemented in "C" code is exemplary of one method of calculating the coefficients for a Butterworth filter that gets a number of poles and assures that the poles are even. In the function, for each pole pair three coefficients are determined, and these coefficients along with the poles are stored.

void filter—trace( int npairs, int nstart, int nend, float *coeff—array, float

*trace—in, float *trace—out);

void bfdesign (float fpass, float apass, float fstop, float astop,int

*npoles, float

*f3db);

main( ) {

int i,npoles,npairs,jpair,filt—indx;

float fpass,apass,fstop,astop,f3db,coeff—array[50][30*3],

a,b1,b2,r,theta,scale,trace—in[1000],trace—out[1000];

for(i=0; i< 50; i++) { /* Assign freq, nyquist = .5 */

	fpass = .50 - (float)i * .01;
	apass = .99;
	fstop = MIN(fpass + .05,.499);
	astop = .05;
	bfdesign ( fpass, apass, fstop, astop, &npoles, &f3db);
	if( npoles % 2 > 0)
	npoles ++;
	printf("fpass %f fstop %f npoles %d f3db

%f\n",fpass,fstop,npoles,f3db);

	npairs = npoles / 2;
	coeff—array[i][0] = npairs;
	r = 2.0*tan(PI*fabs((double)f3db));
	for(jpair=0; jpair< npairs; jpair++){

	theta = PI*(2*jpair+1)/(2*npoles);
	scale = 4.0+4.0*r*sin(theta)+r*r;
	a = r*r/scale;
	b1 = (2.0*r*r-8.0)/scale;
	b2 = (4.0-4.0*r*sin(theta)+r*r)/scale;
	coeff—array[i][jpair*3+1] = a;
	coeff—array[i][jpair*3+2] = b1;
	coeff—array[i][jpair+3+3] = b2;} }

	filt—indx = 10;
	filter—trace( coeff—array[10][0], 600, 640, coeff—array[10],

trace—in, trace—out) }

A filter implemented as part of stage "C" of the processing pipeline is then used to filter a group of input trace samples that straddle the calculated Sample Number from stage "B" of the processing pipeline. In some cases, twenty samples straddling the Sample Number are used to interpolate. Thus, where twenty samples are used to interpolate the fractional Sample Number, ten samples on either side of the Sample Number are filtered resulting in twenty samples of the input trace filtered to the desired frequency. As one example, a Butterworth filter relying on the previously described coefficients can be implemented as part of stage "C" to perform the filtering. The following is a C-language implementation of an exemplary Butterworth filter that copies samples into the output trace for a given range, loops over the identified pole pairs using the previously calculated pole coefficients, and applies the filter. Based on the disclosure provided herein, one of ordinary skill in the art will appreciate a number of other filter techniques and/or implementations that can be used in accordance with the present invention. For example, an ordinary convolution could be used.

void filter—trace( int npairs, int nstart, int nend, float *coeff—array, float

*trace—in, float *trace—out)

{

int j,jpair;

float a,b1,b2,pjm2,pjm1,pj,qjm2,qjm1;

for (j=nstart; j<nend; j++)

trace—out[j] = trace—in[j];

for(jpair=0; jpair < npairs; jpair ++) {

	a = coeff—array[jpair*3+1];
	b1 = coeff—array[jpair*3+2];
	b2 = coeff—array[jpair*3+3];
	pjm1 = 0.0;
	pj = 0.0;
	qjm2 = 0.0;
	qjm1 = 0.0;
	for (j=nstart; j<nend; j++){

	pjm2 = pjm1;
	pjm1 = pj;
	pj = trace—out[j];
	trace—out[j] = a*(pj+2.0*pjm1+pjm2)-b1*qjm1-b2*qjm2;
	qjm2 = qjm1;
	qjm1 = trace—out[j];}}return;}

At stage "D" of the processing pipeline, the twenty filtered input trace samples from stage "C" are interpolated through application of a twenty coefficient sinc filter to determine the fractional distance between the calculated sample and the calculated sample+1. As with the preceding stage, the filter coefficients are pre-calculated by the program before any trace processing begins. From this, a set of coefficients are generated that will allow interpolating a sample to within one tenth of Dt, or the sample rate. To do this, ten sets of twenty coefficients are generated. The following is an exemplary function written in "C" code for producing the coefficients that are a least-squares-best approximation to the ideal obtained through solving a symmetric Toeplitz system for the sinc approximation.

	void mksinc—0( float d, int lsinc, float *sinc) {

	int j;
	double s[20],a[20],c[20],work[20],fmax;
	fmax = 0.066+0.265*log((double)1sinc);
	fmax = (fmax<1.0)?fmax:1.0;
	for (j=0; j<1sinc; j++) {

	a[j] = dsinc(fmax*j);
	c[j] = dsinc(fmax*(1sinc/2-j-1+d));}
	stoepd(1sinc,a,c,s,work);
	for (j=0;j<1sinc; j++)
	sinc[j] = s[j];}

The stoepd function called out in the preceding code segment is expressed below, again in "C" code. As previously discussed, these functions are merely exemplary and one of ordinary skill in the art will appreciate a number of other functions and/or implementations thereof that can be used in relation to the present invention.

	void stoepd (int n, double *r, double *g, double *f, double *a){

	int i,j;
	double v,e,c,w,bot;
	if(r[0] == 0.0) return;
	a[0] = 1.0;
	v = r[0];
	f[0] = g[0]/r[0];
	for (j=1; j<n; j++){
	/*-----------------------------------------*
	* solve Ra=v as in Claerbout, FGDP, p. 57 *
	*-----------------------------------------*/
	a[j] = 0.0;
	f[j] = 0.0;
	for (i=0,e=0.0; i<j; i++)

e += a[i]*r[j-i];

	c = e/v;
	v -= c*e;
	for (i=0; i<=j/2; i++) {

	bot = a[j-i]-c*a[i];
	a[i] -= c*a[j-i];
	a[j-i] = bot;}

	/*------------------------------------------------*
	* use a and v above to get f[i], i = 0,1,2,...,j *
	*------------------------------------------------*/

	for (i=0,w=0.0; i<j; i++)
	w += f[i]*r[j-i];
	c = (w-g[j])/v;
	for (i=0; i<=j; i++)

	f[i] -= c*a[j-i];}}

Of the ten sets of coefficients calculated, the closest to the calculated fractional distance between the calculated sample and the calculated sample+1 is used. A dot product is calculated using the twenty filtered samples and the twenty coefficients resulting in a single sample that represents the interpolated sample, between the calculated sample and the calculated sample+1. An exemplary implementation in "C" code of this function tailored for situations where the fractional distances are non-uniform is provided below.

float interp—samp—sinc( float *trace, int samp—base, float frac, int lsinc,

int nsamp) {

	int j,lsinc—h,indx;
	float table[20],sum;
	frac = frac * -1.0;
	mksinc—0( frac, lsinc, table);
	if( samp—base + 1 - lsinc—h >= 0 && samp—base + 1 + lsinc—h <

nsamp) {

	sum = 0;
	for(j=0; j< lsinc; j++)
	sum += table[j] * trace[samp—base+j+1-lsinc—h]; }

else {

	sum = 0;
	for(j=0; j< lsinc; j++) {
	indx = MAX(0,samp—base+j+1-lsinc—h);
	indx = MIN(indx, nsamp-1);
	sum += table[j] * trace[indx];}} return sum;}

Next, at stage "E" of the processing pipeline, the filtered and interpolated sample is added to the output trace at the image point. It should be noted that the preceding discussion has included examples of filters that one of ordinary skill in the art would know could be replaced with other filter types. Thus, for example, the discussed filters are twenty point filters, but could be replaced with a fifteen point or some other length of filter. Use of such a different filter would still be in accordance with the present invention, the scope of which is set forth in the claims below.

A finite difference wave equation is another example of an algorithm that can be implemented using systems and methods of the present invention. Such an equation is many times more computationally intensive than a Kirchhoff algorithm, and as a consequence it is seldom used in an uncompromised form. If it could be utilized without being compromised, it may produce images far superior to those of the Kirchhoff algorithm. In the prior art, the biggest compromises that are made to the wave equation involve limiting the frequency content of the result to less than twenty Hertz, and the maximum dips that are imaged are typically limited to less than forty-five degrees. It would be desirable to image frequencies higher than fifty Hertz and dips greater than seventy degrees. Furthermore, as now implemented, the data is migrated in the frequency domain which creates other compromises. In particular, to migrate in the frequency domain means that the data is regularized in some way to condition it for the multidimensional Fourier transform. This regularization degrades the input data somewhat.

Using systems and methods of the present invention, the wave equation can be implemented without the previous necessary compromises. The increased processing power offered in various embodiments of the present invention allows for implementation of the finite difference wave equation in the time/depth domain, rather than the frequency domain. Further, because of the implementation in the time/depth domain offset gathers can be output, and images utilizing higher frequencies and steeper dips can be developed.

The general process of implementing the wave equations using systems of the present invention is similar to that described in relation to the Kirchhoff algorithm. The principal difference is in the details of the processing pipeline. For the wave equation an exemplary processing pipeline includes the following stages expressed as loops:

	Outer loop over output z grid locations (depth)

Center loop over output x grid locations

Inner loop over output y grid locations

	Update the differencing equations
	Tri-diagonal matrix solver
	Thin lens adjustment ( vertical shift for changing
	velocity gradient)

	End of all three loops
	Update the output image traces.

Because the wave equation is a downward continuation function, the outer loop is the loop over the depth increments. With this looped implementation, many levels of parallelism can be achieved, and multiple processing pipelines can be implemented using amorphous hardware elements similar to that discussed in relation to the Kirchhoff algorithm.

Each of the processing pipelines images a different input trace into the set of output traces. Further, a great degree of parallelism can be achieved in the tri-diagonal matrix solver, so a single solver itself can be branched into many mini-pipelines. These parallel processing pipelines can be implemented in a number of amorphous hardware elements within a single FPGA, or across multiple FPGAs. An exemplary processing pipeline tailored to implementing the finite wave equation may include for each x, y and z: (A) update differencing equations, (B) tri-diagonal solver, (C) thin lense adjustment, and (D) sum the input trace's contribution into an output grid matrix. Each of these stages can include many clock cycles, and the overall processing pipeline can include hundreds of clock cycles.

In stage "A" a function F(Dx,Dy,Dz,Vint,Dt, Data) is set up to create a set of simultaneous equations. In the function, Dx, Dy and Dz are the distances between the grid cells in the respective x,y,z coordinate direction. Vint is the interval velocity model, Dt is again the time sampling of the input trace, and Data is the actual input trace data. Example code can be obtained from Claerbout, Imaging the Earth's Interior, p. 136. The referenced portion of the foregoing publication is included in Appendix A attached hereto. This set of simultaneous equations with boundary conditions form a tri-diagonal matrix.

In stage "B", the tri-diagonal matrix is solved producing a set of values that represents that input trace's contribution to the x,y,z image grid. These values are then filtered based on function F(Vint, Ds, Dr) to avoid aliasing similar to that discussed above in relation to the Kirchhoff algorithm. In stage "C" the thin lense term is applied to account for the lateral velocity gradient using a function, F(V(x,y,z), Dz). Such a function is more fully described on page 137, Claerbout, Imaging the Earth's Interior, p. 137. This adjusted value is then summed into the output grid matrix at stage "D".

In some cases, all processing is completed through programming the amorphous hardware elements a single time as previously described. In other cases, processing is completed in multiple programming stages where, for example, an amorphous hardware element is programmed to perform one function, and subsequently re-programmed to perform another function. Such an approach is described in relation to blocks 650-680 of FIG. 6. The results from one programmed computation can be received by the general purpose computer (block 650). The amorphous hardware elements can then be re-programmed to perform a different computation (block 660), and the results from that computation can also be received (block 670). The results from the various computations can then be assembled to create a final result (block 680).

Based on the disclosure provided herein, one of ordinary skill in the art will appreciate a variety of other methods and/or algorithms to which amorphous computing systems in accordance with the present invention can be applied. For example, some amorphous hardware elements can be programmed to perform one algorithm, while others are programmed to perform a different algorithm in parallel to the first algorithm. The results from the algorithm can then be applied to amorphous hardware elements programmed to implement yet another algorithm.

As will be appreciated by one of ordinary skill in the art, the degree of parallelism achievable is significant when compared to a general purpose microprocessor. While a general purpose processor offers some degree of parallelism through multiple pipe architectures, superscalar architectures, and the like, they are essentially serial devices when compared to the systems of the present invention. In some cases, multiple general purpose microprocessors can be arranged in parallel, however, this involves a great deal of inefficiency as much of the microprocessors goes unused, and memory sharing attributes associated with such parallel systems can be paralyzing.

As one example of the present invention, each input sample is interpolated, filtered, weighted and summed into five hundred thousand output samples. Each step of this process is unique for each output sample. On a general purpose microprocessor, about fifty clock cycles on a two Gigahertz chip are required to get one input sample through the interpolation process. To increase the throughput, about two hundred computers, each including two general purpose microprocessors are assembled to provide about four billion outputs per second.

In contrast, on the amorphous computing systems in accordance with the present invention where each amorphous computing element includes 0.5 Mbytes of memory, approximately five hundred input traces are loaded into the amorphous computing element, and outputs are sent to an eighteen Gigabyte memory communicably coupled to the amorphous computing elements. Each of the amorphous computing elements are programmed with a pipeline specifically tailored to perform interpolation, filter, weighting and summation of input traces into output traces. Using such an approach, each input trace can be formed into an output trace for each clock period of the amorphous computing element.

Thus, for example, for each output sample we start a pipeline on one of the input traces. Where both the input traces and output traces have five hundred samples each, the process proceeds as follows: (A) input trace 1, sums into output trace 1 sample one; (B) input trace 2, sums into output trace 1sample two; (C) input trace 3, sums into output trace 1 sample three . . . (D) input trace 400, sums into output trace 1 sample five hundred, (E) input trace 1, sums into output trace 2, sample one . . . (F) input trace 400 sums into output trace 400 sample five hundred. During the process, all of the pipelines programmed in the amorphous computing elements move to the next output sample on the next clock until all pipelines have summed all the input traces into their respective samples.

Thus, in this example, there are one thousand computational pipelines programmed in respective amorphous hardware elements within an FPGA. At startup, a few hundred clocks are expended setting the FPGA up for computation, after which one result per clock can be generated. The input traces sit in memory onboard the FPGA, and are summed into many Gbytes of output traces. As the input traces are summed into each output trace, the previous output traces are stored back in main memory, and the next output traces are fetched from main memory onto the amorphous hardware element. A simple double buffering scheme can be implemented to assure that the pipelines are always executing. After all five hundred input traces have been summed into the full eighteen Gbytes of output space, the next set of input traces are fetched, and the process continues.

Using this type of arrangement, one thousand results can be generated per clock on an amorphous hardware element operating at one hundred Megahertz. This results in one hundred billion results per second on a single FPGA. Thus, where ten FPGAs are programmed to operate in parallel, one 1 trillion results per second (or more where faster, larger, or more amorphous hardware elements are used) can be achieved. This is approximately 250 times the performance of a system created using the two hundred PC cluster described in the aforementioned example.

Turning to FIG. 8, processes associated with supporting varied trace lengths in accordance with some embodiments of the present invention are described. In contrast to a general purpose computer that typically uses some form of a "FOR" loop to implement a variable trace length, some embodiments of the present invention use FPGAs with preset trace lengths directly mapped to amorphous hardware elements on the device. Such FPGAs are not flexible enough to accommodate different length traces within the same direct mapping. Thus, if a function is implemented to accommodate one thousand samples, that map will work for only a trace length of one thousand samples. In some embodiments of the present invention, such a one thousand trace function can be extended to operate on trace lengths of one thousand or less by padding the trace with null data or zeroes out to the one thousand trace limit. Thus, for example, where a given trace is nine hundred samples, the samples from nine hundred one to one thousand are filled with null data to make the trace length compatible with the function. This approach does not allow for accommodating functions with a trace length greater than one thousand. Further, such an approach wastes area on an FPGA by implementing processing pipelines that operate on null data. Thus, in some embodiments, multiple functions of varying length are implemented. From this, a best fit function is selected allowing for accommodation of a particular trace length, while minimizing the amount of wasted area on the FPGA.

In one particular embodiment, the functions are VIVA™ based functions that are created and maintained as a reusable library. Each of the functions implement the same algorithm, but are each modified to accommodate a different trace length. Such trace lengths can be 500, 600, 700, 800, 900, 1000, 1100, 1200, etc. It should be recognized that any number of functions can be created depending upon the amount of space to be used in holding the functions, and the willingness to sacrifice space of an FPGA. Thus, in the extreme case, a function can be created to handle each increment in trace length such as, for example, 1000, 1001, 1002, etc. Each of the functions is then identified by a name that indicates the length of the function to be handled such as, for example, Kirchhoff_—Length1000. The various functions can then be compiled together, and at run time the general purpose computer responsible for programming the various FPGAs can select the appropriate function based in part on the desired trace length. In this way the direct mapping limitation associated with using amorphous hardware elements can be hidden from the end user by making a separate design for each increment of trace length that is to be supported.

FIG. 8a includes a flow diagram 800 that illustrates a process in accordance with embodiments of the present invention for producing functions capable of implementation on amorphous hardware elements operable to accommodate traces of differing lengths. Following flow diagram 800, a first trace length x is selected (block 810). This trace length x defines the shortest trace length function to be created. The function is then created for the trace length (block 820) and the function is named to indicate the trace length (block 830). It is then determined if the trace length is equal to or greater than a desired maximum trace length (block 840). If it is determined that the trace length is less than the desired maximum trace length (block 840), the trace length is incremented by a trace step y (block 850) and the processes of blocks 820, 830 and 840 are repeated for the incremented trace length. Alternatively, where the trace length is determined to be greater than or equal to the desired maximum trace length (block 840), all of the functions created to that point are compiled together (block 860), and the executable functions are stored to a database(block 870).

Turning to FIG. 8b, a flow diagram 801 illustrates a process that utilizes the previously described multi-length functions to implement a function algorithm applied to a particular trace length. Following flow diagram 801, an input trace and/or a survey of input traces is accessed (block 811). A length for a given input trace can be determined (block 821), or in the case where all input traces within a given survey are of the same length, a length for all of the traces can be determined. A database including the various multi-length functions is accessed (block 831), and one of the functions that is incrementally larger than the previously determined trace length is identified (block 841). Thus, for example, where the determined trace length is 1021 and the available functions are Function_—length900, Function_—length1000, Function_—length1100 and Function_—length1200, the selected function is Function_—length1100.

This identified function is accessed (block 851) and one or more amorphous hardware elements are programmed to perform the function (block 861). Input traces are then accessed (block 871), and the input traces are padded with null data such that the input trace length equals the length of the selected function. Thus, in the previous example, seventy-nine zeroes would be added to each input trace. These padded input traces are then fed to the amorphous hardware elements programmed to perform the selected function on the input traces (block 891).

The invention has now been described in detail for purposes of clarity and understanding. However, it will be appreciated that certain changes and modifications may be practiced within the scope of the appended claims. For example, it should be recognized that the present invention can be used in relation to a number of algorithms and/or variations of algorithms. For example, the present invention can be used in relation to the Kirchhoff pre-stack time and the Kirchhoff pre-stack depth, which are both variations on the same basic algorithm. The difference between the aforementioned variations is in the function that calculates the time of the input trace to contribute. In the Kirchhoff prestack-time it is calculated using the exemplary function previously described, and in the Kirchhoff pre-stack depth it uses a set of tables, that can be calculated using a program called a travel time generator (a ray tracer in one particular case). The program produces a disk file of times related to source/receiver combinations. The migration then reads these tables to get the time of the input trace, rather than uses the other portions of the aforementioned function.

Accordingly, it should be recognized that many other systems, functions, methods, and combinations thereof are possible in accordance with the present invention. Thus, although the invention is described with reference to specific embodiments and figures thereof, the embodiments and figures are merely illustrative, and not limiting of the invention. Rather, the scope of the invention is to be determined solely by the appended claims.

* * * * *

---