GB2298920A - Method of estimating positions of seismic elements in a marine seismic array - Google Patents

Abstract

Estimates of positions (m(tn)) of elements within a seismic array are formed by fitting theoretical inter-element distances to measure inter-element distance (dk(tn)) for each shot (n) of a seismic survey. The fitting uses a constrained lp norm. The constraint requires each estimate of an element's position (m(tn)) to be related to that estimated for the previous shot (m(o)(tn)). For each shot, the estimate is formed iteratively within a loop comprising steps 3 to 8 of Figure 7.

Description

METHOD OF ESTIMATING POSITIONS OF SEISMIC ELEMENTS IN A MARINE SEISMIC SURVEY.

The present invention relates to a method of estimating the position of seismic elements, i.e. guns and streamers, within a marine seismic survey arrangement.

Undertaking a marine seismic survey requires the towing of a spatially extended array of energy sources, i.e. guns and hydrophones, behind one or more survey vessels. The positions of each of the energy sources and each of the hydrophones needs to be known in order to properly carry out processing of the seismic data collected by such an array. However, the position of each array element is not rigidly fixed. The guns and hydrophones are connected to the survey vessel or vessels by cables.

The cables are flexible and consequently cannot hold each element in a rigidly fixed position. Furthermore the array is subjected to the motion of the waves, to currents and, where an element extends above the sea surface, to the force of the wind. All these factors cause the movement of elements of the array away from their nominal positions. Analysis of the collected seismic data cannot proceed until the positions of the element within the array have been established.

It is known to determine the positions of the elements within the array by estimating the position of each element, calculating acoustic ranging data and other ranging data for a plurality of inter-element distances, comparing the measured and calculated ranges, bearings and positions and modifying the estimates of the positions in accordance with a least squares fitting procedure.

A problem arises in that the least squares fitting procedure is non-robust in the presence of poor data. The acoustic ranging data can be of very variable quality. The marine environment can be a noisy place, e.g.

wave action on the surveying equipment and noise from the survey vessels, can degrade the signal-to-noise ratio of the acoustic ranging measurements to such an extent that some of the data are unsuitable for use within a least squares fitting procedure. It is known to pre-process the data such that data which are obviously bad are removed. Gaps in the remaining data are then interpolated and the data are smoothed. The least squares fit is then applied to the smoothed data. Such a method is used in the TRINAV (Trade Mark) system used by Geco-Prakla.

The pre-processing has become very sophisticated, and human intervention is not as frequently required as was the case with the earlier pre-processing techniques. However, human intervention may still occasionally be required and this limits the turn-around time for processing the navigation data. The removal of obviously bad data may result in the filtering procedure becoming ill-conditioned.

Bruce R. Harvey, "Survey network adjustments by the L1 method" Aust, J.

Ceod. Photogram, Surv. No.59, December, 1993, pp 39-52, describes the use of an L, method to make adjustments in site survey data prior to processing the data using a least squares fitting procedure.

According to a first aspect of the present invention, there is provided a method for estimating the positions of a plurality of elements within a marine seismic surveying arrangement, the method comprising the steps of fitting estimates of the element positions to measured inter-element distances using a first function which forms a deviation between a theoretical inter-element distance based on the estimates of the element positions and measured inter-element distances, the deviation being raised to a first power less than two, and a second function which urges the estimate of an element's position towards a position based on a preceding estimate of the respective element's positions.

According to a second aspect of the present invention, there is provided a method of estimating the positions of a plurality of elements within a marine seismic surveying arrangement, the method comprising fitting estimates of the element positions to measured inter-element distances using a function which forms a deviation between a theoretical interelement distance based on estimates of the element positions and measured inter-element distances, the method further comprising the steps forming a series of simultaneous equations representing the estimates of the element positions and solving the simultaneous equations iteratively using a conjugate gradient method.

Preferred embodiments of the present invention are defined in the appended claims.

It is thus possible to fit to data using a constrained Qp norm. The constraint is applied by requiring the solution for the current shot to be related to the solution for the previous shot. The use of Qp norm i.e. a procedure similar to least squares but with a power less than two, decreases the sensitivity of the fit to bad data.

The present invention will be further described, by way of example, with reference to the accompanying drawings, in which: Figure 1 is a schematic diagram illustrating the positions of nodes in a marine seismic survey arrangement; Figure 2 is a diagram illustrating the measurements taken between nodes in Figure 1; Figures 3a to 3g show raw inter-node measurement data; Figure 4 shows an estimate of the locus of a single node using the raw data; Figure 5 shows the data of Figure 3 after pre-processing as used in the prior art; Figure 6 shows the locus of the node calculated using the pre-processed data exemplified in Figure 5; Figure 7 is a flow chart of processing steps constituting an embodiment of the present invention; and Figure 8 shows the locus of the node calculated using a method according to the present invention.

Figure 1 schematically illustrates the positions of a plurality of nodes at the front of a marine seismic survey using two survey vessels B1 and B2.

The positions are plotted relative to an arbitrary origin in the network.

Translational and rotational movement of the network as a whole is determined by other measurements and not range data, and is not considered herein. Each survey vessel is arranged to take range and/or bearings to other nodes within the survey from the positions within the vessel. For vessel B1, the forward measuring position is denoted as b1t1; whereas the rear measuring position is denoted as bird. Four guns G1, G2, G3 and G4 and four streamers S1, S2, S3 and S4 are towed behind the vessels B1 and B2. Positioning devices for the guns and streamers are indicated by glhl, girt1, g2hl, g2t1, slhl, sits, etc.The position of each hydrophone within a streamer is located with reference to the front portion of the streamer by reference to the known length of interhydrophone sections within the streamer and direction measurement using compasses within a streamer. The positioning of Figure 2 schematically shows the internode distance measurement for the network shown in Figure 1.

The arrangement shown in Figure 1 and 2 comprises twenty four nodes located on the ships, floats, guns and streamers. In total ninety one internode distances were measured for each shot of the seismic survey.

The spatial extent of the network is approximately 400m square. The vertical positions, i.e. depths and heights, of the nodes are assumed to be fixed and known. The typical motion of a node with respect to the network is + 10m. The translational motion a node with respect to the network between shots is normally less than lm and the positioning accuracy required of the node with respect to the network is about 2m.

The absolute position ot the ships and the translational and rotational motion of the network as a whole is determined by data other than the ranging data and is not of concern here.

Figures 3a to 3g show raw data for some of the internode distance measurements. The figures show a typical cross-section of the data, ranging from good to bad. In each of the figures, the horizontal scale covers 250 shots, the vertical axis covers a range of +5m about the mean distance, and the figures have been plotted to show the deviation of each internode distance about its mean value, d.

If the raw data are Immediately supplied to a fitting scheme using a least squares fitting procedure, the results are unsatisfactory and are shown in Figure 4 in which the position of one node s2t2 is plotted, each dot representing the node position at a shot and every 10th shot being numbered. The graph is only plotted over a limited range of positions and some of the measurements fall outside the limits of the graph.

The prior art approach is to filter and interpolate the internode distance measurements, so as to obtain the results shown in Figures 5an59, and then to supply these data to a least squares fitting procedure. Figure 6 shows the result of this manipulation on the position data for the node s2t2. It should be noted that this pre-processing must inherently affect the result obtained during the fitting.

The method of the present invention also solves the positions of the nodes as a function of time, but it does so without the need to preprocess the data.

During the surveying a number of shots are made at times tn. For a network having j nodes the position of each node at each shot can be represented bv xl(tn) where j is an integer from 1 to j. K internode distance measurements are made for each shot.

For simplicity, the method will only be described with reference to internode distance measurements, but fully determined positions and bearings may also be incorporated by suitable modification of some of the x,(tn) data.

The internode distances dk for a kth distance measurement (where k= 1 to K) between the ith and ith nodes are corrected for the vertical displacement (assumed known from the geometry of the network) between the nodes, by applying the following operation:

where d'k is the corrected internode distance, dk is the measured distance, zl is the vertical position of the ith node, and z is the vertical position of the jth node.

The theoretical internode distances D'k are calculated from:

where x, and x are the x components of theoretical node positions x and Xj, and y, and y, are y components of theoretical node positions yi and yi.

The fitting procedure finds node locations that minimise the sum of the objective function #d + (t)m for each shot t,, where

(3 and 4) where #d is the data norm objective function, rn is the regularisation norm objective function, #d@ is the standard deviation of the kth internode distance measurement, #m@ is the standard deviation for the jth theoretical internode distance measurement, and # is a regularisation parameter.

The unknown model coordinates xj and y, are expressed as a single vector m, where m2ll = x1 for j = 1 to J m21 = y1 for j = 1 to The locations m(0) are an estimate of the solution. Various estimates are possible, such as the nominal positions of the nodes, the previous location of each node, a linear extrapolation of node position based on a constant velocity model or a high order extrapolation method. The previous location has proven to be a satisfactory choice, such that mj(0)(tn) = mj(tn-1) (5) The regularisation parameter A can be chosen for the objective solution (prescribed data fit, d=1), 1), the subjective solution (prescribed model fit, #m= = 1), the stochastic solution (#= = 1), or the corner of the trade-off curve, i.e.

# # @ # # m = # # # # If the standard deviations Gdk and #mk are appropriate, the stochastic solution and the 'corner' solutions should be similar.

The function #d + ##m is minimised by linearising the equations and solving iteratively.

In this equation P, Y, Q and Z are diagonal matrices, and bk( ) = dk' - Dk'( ) for k = 1 to K, m( )2j-1 = X@( ) for j = 1 to J, m( )2@ = y@( ) for j = 1 to J, Ak@( ) = (X@( )-X@( ))/Dk'( ) if l = 2i-1 and {i, j} # k, = (Y( )-Y@( ))/Dk'( ) if l = 2i and {i, j} # k, = -(x@( )-X@( ))/Dk'( ) if l = 2j-1 and {i, j} # k, = -(Y@( )-Y@( ))/Dk'( ) if l = 2j and {i, j} # k = 0 otherwise

for k = I to K 1 Z = #m@#2J for j = 1 to 2J

for k = 1 to K

for j = 1 to 2J.

and p and U are iteration indices.

The term pu) YA belongs to the data norm, whereas the term belongs to he regularisation norm. The expressions for P(@) are summed for Q = 1 to 2J. The cut-off factors Adk and Am, are selected to produce a Huber filter (least-squares weighting) for small deviations.

The least squares solution of these equations minimises In'l, which is equivalent to minimising the objective function (#d + ##m) However, P and Q depend on the solution (except when p=2 or q=2), so the equations must be solved iteratively, adjusting P and Q at each iteration.

The iterative, reweighed least squares solution is equivalent to the required (p and (p solution.

If q#2, the Q reweighting can be modified to apply the q-norm to the distance rather than the components (combining the odd and even components ot m).

The order of computation and iteration is shown in Figure 7. Values of the measured inter-node distances and standards deviations for the Nth shot are read in at step 1. An estimate of m(o)(tn) is formed at step 2. If the shot involved is the first shot, the nominal positions of the nodes may be used. For subsequent shots, the positions at the (tn-1) shot may be used. Values for AW and b are calculated at step 3. Initial values P(o) and Q("', where Pkk(o) =Q@(o)= 1 , are set at step 4.

Estimates of m are formed at step 5 by iterative solution of equation 6.

Updating of the values of P''') and Q'9' by further iteration over U is performed at step 6 which forms an iteration loop with step 5. Steps 5 and 6 can in fact be performed together. The 'least squares' solution at step 5 can be efficiently performed using the CG (methods of conjugate gradients for solving linear systems) described by M. Hestenes and E. Stiefel, NBS J. Research, 49, 409-435, 1952, (which is incorporated herein by reference) since the matrix is large (2J + K) x 2J but sparse (4K + 2J non-zero terms).To describe Hetsenes and Stiefel's method, we simplify the notations of equation (6) to produce: b' = Am' + n'. (7) Starting with m'(0)=0, s(0)=b', r(0)=P(o)=A'Ts(0) (8) we iteratively solve (e= = 1,2...) where # within parenthesis represents a number of an iteration and the first iteration (#= = 1) starts with equation (11)::

p(@)=r(@@@)+ss(@)p(@ (10) @@=A'pCG(@) (11)

m'(@) = m'(@-1)+α(#)pCG(@) (13) - (14) r(#)=A'Ts(#). (15) The method, and other variants of it, are efficient because the sparse matrix A' is only used in the products (11) and (15) which can take full advantage of the sparseness to utilise fast matrix manipulation techniques. Normally, the iterations over c are continued until r(#) converges to approximately zero. However, this method is itself iterative and it has been shown by J. A. Scales, A. Gersztenkorn and S.Treitel, 1988, "Fast (p solution of large, sparse linear systems: application to seismic travel time tomography", J. Comp. Phys, 75, 314-333, that iterative reweighing can be performed before the CG algorithm converges, eg. after JIS iterations. Step 7 examines the gradient of the trade-off curve and the value of k can be refined in a regularization parameter loop containing steps 4 to 7. However, this loop may be avoided by setting # to a constant value, such as unity. Step 8 performs a check for convergence. If sufficient convergence has not been achieved, a further iteration in is performed and control is passed back to step 3.If sufficient convergence has been achieved, the values of m (tn) are stored at step 9 and control is returned to step 1 where data for the next shot is read in.

Figure 8 shows the calculated path of the node s2t2 using all of the data, without using any pre-processing and using the values p=1 and q=2.

The motion is not as smooth as that shown in Figure 6 as 'ne input data had not been pre-processed by excision of 'bad' data and smoothing, but the results are very similar. The region 20 which does not occur in data processed by the prior art, relates to the first few (10 or so) position calculations where the effects of using the nominal positions of the nodes were still noticeable.

The use of the Qp norm (i.e. p less than two) and the regularisation norm to constrain the solutions allows the method to work well in the presence of bad data and is robust against ill-conditioning. The processing techniques, such as the CG method, are not crucial but do allow efficient solution of the physical problem. The Huber filter enhances the stability of the processing.

In the example, p= 1 was chosen and fixed. However, a time varying p may be used depending on the data quality. The smoothness of the solution can be improved by adding extra constraints. i.e. extra norms to the objective function and extra augmenting equations.

It is thus possible to provide an efficient method of estimating the position of equipment within a marine seismic survey. The method does not rely on pre-processing of the data and hence can avoid the use of the sophisticated and time consuming manipulation of data to perform the pre-processing. The method has been described with reference to the positions of elements at the front of a seismic survey arrangement.

However the method can also be used to establish the relative positions of elements at the rear of a seismic survey arrangement, or to establish the positions of elements within the entire survey arrangement, i.e. the front network, the streamers and the tail network.

Claims (15)

1. A method of estimating the positions of a plurality of elements within a marine seismic surveying arrangement, the method comprising the steps of fitting estimates of the element positions to measured interelement distances using a first function which forms a deviation between a theoretical inter-element distance based on the estimates of the element positions and measured inter-element distances, the deviation being raised to a first power less than two, and a second function which urges the estimate of an element's position toward a position based on a preceding estimate of the respective element's position.

2. A method as claimed in Claim 1, in which the first function forms a data norm by dividing each deviation by a respective standard deviation to form a weighted deviation, raising the weighted deviations to the first power and summing the results, and dividing the sum by the number of deviations.

3. A method as claimed in Claim 1 or 2, in which the first power is greater than or equal to unity.

4. A method as claimed in any one of the preceding claims, in which the second function forms a regularization norm by forming a second deviation between the estimate ot an element's position and a further estimate of the respective element's position, the further estimate being based on at least one preceding estimate of the respective element's position.

5. A method as claimed in Claim 4, in which the further estimate is the previous estimate of the position of the respective element.

6. A method as claimed in Claim 4, in which the further estimate is a constant velocity extrapolation of the respective element's trajectory.

7. A method as claimed in Claim 4, in which the further estimate is an extrapolation of the respective element's trajectory based on at least velocity and acceleration of the element.

8. A method as claimed in any one of Claims 4 to 7, in which the second function further comprises forming a second weighted deviation for each element by dividing the second deviation for each element by a respective standard deviation, forming a second sum of the weighted deviations raised to a second power and dividing the second sum by the number of second deviations.

9. A method as claimed in Claim 8, in which the second power is greater than or equal to unity and less than or equal to two.

10. A method as claimed in any one of the preceding claims, in which the second function is weighted by a regularization parameter.

11. A method in which the fitting is achieved by minimising the sum of the first and second functions iteratively.

12. A method as claimed in Claim 11, in which the conjugate gradient method is used.

13. A method as claimed in any one of the preceding claims, in which the first power is a variable.

14. A method of estimating the positions of a plurality of elements within a marine seismic surveying arrangement, the method comprising fitting estimates of the element positions to measured inter-element distances using a function which forms a deviation between a theoretical inter-element distance based on estimates of the element positions and measured inter-element distances, the method further comprising the steps forming a series ot simultaneous equations representing the estimates of the element positions and solving the simultaneous equations iteratively using a conjugate gradient method.

15. A method as claimed in Claim 14, in which the function forms the modulus of the difference between a theoretical inter-element distance and the corresponding measured inter-element distance and further applies a constraint which urges each estimate of an element position towards a further estimate of the element's position based on a preceding estimate of the respective element's position.