POSITIONING OF SEISMIC EQUIPMENT IN A TOWED MARINE ARRAY

Abstract

A positioning system (100) for marine seismic surveying, comprising a towing vessel (110), a source array (120) and a receiver array (130) with several streamers (131). Each streamer (131) comprises at least three birds (134) and positioning sensors (134-137) wherein several seismic receivers (132) are placed between each pair of birds (134). The system (100) comprises a dynamic model wherein each streamer (131) is represented by a fitted B-spline curve and each bird (134) is associated with a constant velocity and a constant acceleration; and a Kalman filter using the dynamic model and observations from the positioning sensors (134-137) to provide a geodetic position of each seismic receiver (132) with better accuracy than the dynamic model alone and observations alone within a time interval t equal to or less than a minimum time between shots determined by the source array (120).

Claims

1-10. (canceled)

11. A positioning system for marine seismic surveying, comprising: a towing vessel; a source array; and a receiver array with several streamers, each streamer comprising at least three birds and positioning sensors, wherein several seismic receivers are placed between each pair of birds, wherein a) a dynamic model wherein each streamer is represented by a fitted B-spline curve and each bird is associated with a constant velocity and a constant acceleration; and b) a Kalman filter using the dynamic model and observations from the positioning sensors to provide a geodetic position of each seismic receiver with better accuracy than the dynamic model alone and observations alone within a time interval t equal to or less than a minimum time between shots determined by the source array.

12. The system according to claim 11, further comprising local support limited to at most four subsequent birds such that a change imposed anywhere in the receiver array affect no element in the system outside this range.

13. The system according to claim 11, wherein a position sensor further provides observations of velocity and/or acceleration.

14. The system according to claim 11, wherein the observations are expressed in geodetic (X, Y) and/or vessel bound (x, y, z) coordinates.

15. The system according to claim 11, further comprising a closed loop control system for controlling the physical components of the system, capable of receiving input from the Kalman filter and of changing sign on an input vector and presenting the resulting vector as response.

16. The system according to claim 15, wherein the input vector expresses a deviation from a desired position in geodetic coordinates (X, Y).

17. The system according to claim 15, wherein the input vector expresses a deviation from a desired streamer configuration in vessel bound coordinates (x, y, z).

18. The system according to claim 15, wherein the input vector is expressed in a tangent-normal-binormal frame local to a bird and/or positioning sensor.

19. The system according to claim 11, wherein a transformation between any pair of coordinate systems is performed by one matrix multiplication.

20. The system according to claim 11, further comprising linear extrapolation of the geodetic coordinates within the time interval t.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0024] The invention will be described by way of example and reference to the accompanying drawings, in which:

[0025] FIG. 1 illustrates a system according to the invention and

[0026] FIG. 2 illustrates a general setting for a spline fit to streamers.

DETAILED DESCRIPTION

[0027] The drawings are schematic and intended to illustrate the invention. Thus, they are not to scale, and numerous details known to one skilled in the art are omitted for clarity.

[0028] The general idea of the invention is to use a fast and effective dynamic model and Kalman filtering to improve estimates of geodetic coordinates. We distinguish between speed (scalar) and velocity (vector), and note that the velocity (speed and heading) of a towing vessel moving at 4-5 knots is approximately constant during each time interval t<5 s. If the local constant velocity {dot over (x)} is known at each geodetic position at the start of each time interval t, then even more precise geodetic position estimates may be extrapolated during t.

[0029] The notation throughout this description and the claims follow usual conventions: Vectors and matrices are in bold Latin letters with capital letters denoting matrices. Italics indicate that R.sup.3 and C(s) are not matrices. Mark notation, e.g. C(s), is used for space derivatives and dot notation, e.g. {dot over (x)}, is used for time derivatives. All scalars are real.

[0030] FIG. 1 illustrates a system 100 for marine seismic surveying comprising a seismic survey vessel 110 towing a source array 120 with several subarrays 121 and a receiver array 130 with several streamers 131. Arrow S.sub.q-1 illustrates a previous sailing line from top to bottom of the drawing and arrow 201 illustrates an area on both sides of the sailing line S.sub.q-1 covered by the previous pass. Similarly, arrow S.sub.q illustrate a present sailing line and a towing direction in the opposite direction, i.e. from bottom to top of the drawing. Arrow 202 indicates the stripe to be covered by the present pass. According to common practice, there is a small overlap 203 between the previous and present pass. The overlap is preferably small to achieve a smallest possible survey time, yet large enough to avoid or reduce infill.

[0031] The shaded area 204 illustrates an area with adequate coverage, and the irregular area 205 within the previous pass is an area with inadequate coverage, i.e. an area needing infill. The area 205 lack coverage for some unspecified reason, and some streamers 131 are steered laterally from the present sailing line S.sub.q in order to make an additional pass over the area 205.

[0032] In FIG. 1, each streamer 131 is connected to a lead-in cable 112 in a towing point 113. Paravanes 114, also known as deflectors and doors, provide a lateral force on a spreader cable 115 connected to the towing points 113 in order to maintain fixed distances between leading ends of the streamers 131. Referring to the introduction, we note that the fixed distances do not need to be equal as long as the streamers are maintained in a parallel configuration because the Langeland '156 patent belongs to free technique. However, if the spread at the tail end exceeds 5%, one may be accused for infringing on the '094 Eick patent, at least in US territorial waters. To be sure, it may be advisable to keep the leading ends at a fixed distance within 5%. No irony intended.

[0033] Several receivers 132 and birds 134 are arranged along the streamers. The birds 134 are known devices of any kind capable of imposing a lateral force on the streamer 131 and are used for steering, e.g. to compensate for lateral currents and adjust the towing depth.

[0034] A vessel-bound spatial Cartesian coordinate system x, y, z is illustrated by axes x, y with origin midship in vessel 110, y-axis parallel to and directed opposite the present sailing line or towing direction S.sub.q and x-axis horizontal and perpendicular to the y-axis. The z-axis (not shown) is conveniently perpendicular to the x and y-axes in a right handed Cartesian. The direction of the y-axis on FIG. 1 illustrates that the z-axis indicates depth. A navigation system might show the y-axis pointing in the towing direction, i.e. opposite to that shown in FIG. 1. Streamer bound coordinates are expressed in right handed Cartesian coordinates tangent, normal and binormal to the streamer, commonly called a TNB-frame. Similar to the z-axis, the B-axis is not shown in FIG. 1.

[0035] The area 205 is expressed in geodetic coordinates, for example in an Earth-bound Cartesian grid (X, Y). A graphical display similar to FIG. 1 merely requires a 2D map of the area. In particular, a point between a source 121 and a receiver 132 is sufficient to indicate a reflection in a graphical view similar to FIG. 1. The density of such points may indicate coverage. Similarly, underwater currents at the towing depth in geodetic coordinates could be displayed as a field of short arrows representing sea current velocities as known in the art.

[0036] For steering physical arrays, the system 100 needs a closed loop control system. In the present context, such a control system just changes sign on an input vector. For example, the proper response to a deviation from a desired position in geodetic resp. local coordinates is a displacement of equal size in the opposite direction expressed in input coordinates. Similarly, the response to an input force vector would be a response vector of equal size and opposite direction expressed in arbitrary input coordinates. The control system as such is known in the art and need no further explanation herein.

[0037] The Cartesian coordinates in FIG. 1 are chosen for convenience, and may be replaced with polar, spherical or other coordinates if desired. However, according to common practice calculations should be performed in a consistent coordinate system to prevent round-off errors from adding up. For example, transforming rotational coordinates roll , pitch and yaw to Cartesian coordinates (x, y, z) involves calculating sines and cosines which introduce round-off errors. As long as calculations are performed in one system, e.g. (x, y, z), the numerical errors caused by one transformation from (, , ) to (x, y, z) before calculation and another transformation back to (, , ) are insignificant. Similar considerations apply to representations that involves finite differences, e.g. Hermite splines and Jacobians.

[0038] All transformations between the systems in FIG. 1 comprise a translation and a rotation, which are examples of affine transformations. Other important affine transformations are scaling, reflection and projection, e.g. from 3D to a 2D computer screen. Cartesian coordinates form affine spaces, and affine transformations by definition preserves collinearity and ratio of length. Thus, a position vector and its time derivatives velocity and accelerations in Cartesian coordinates are invariant under affine transformations. In addition, 2D and 3D B-splines are affine invariant, and may be projected on a 2D computer screen without difficulty. For example, an affine transformation of a 3D B-spline curve is equivalent to transformation of its control points. After transformation of a few control points, the curve is easily recalculated in 2D for display. The following brief overview of affine transformations in matrix notation is included for reference.

[0039] Any sequence of affine transformations can be combined to one transformation in a higher dimensional space. For example, a 3D affine transformation can be expressed as:

y(x)=Ax+b(1)

where
x is a vector in an affine domain space, e.g. the (x, y, z) coordinates in FIG. 1,
y, b are vectors in an affine codomain (target domain), e.g. (T, N, B) or (X, Y, Z) in FIG. 1 and
A is a transformation matrix, e.g. a 33 rotation matrix.

[0040] The affine transformation (1) can be rewritten in augmented form in 4D:

[00001]$\begin{array}{cc}{y}_4=\left[\begin{array}{c}{y}_1\\ y_2\\ y_3\\ 1\end{array}\right]=M.Math..Math.x_4=\left[\begin{array}{cccc}{a}_{11}& a_{12}& a_{13}& b_1\\ a_{21}& a_{22}& a_{23}& b_2\\ a_{31}& a_{32}& a_{33}& b_3\\ 0& 0& 0& 1\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\\ 1\end{array}\right]& \left(2\right)\end{array}$

where a.sub.ij is a real scalar from row i and column j of the transformation matrix A.

[0041] Any sequence of affine augmented matrices can be multiplied to yield a new matrix of the same form. Moreover, M is invertible if and only if A is invertible, and is expressed by:

[00002]$\begin{array}{cc}{M}^{-1}=\left[\begin{array}{cc}{A}^{-1}& \left(-A^1.Math.b\right)\\ 0.Math..Math.\mathrm{.Math.}.Math..Math.0& 1\end{array}\right]& \left(3\right)\end{array}$

where A.sup.1=A.sup.T if A represents any rotation in R.sup.N, N=2, 3, . . . . Note the similarities with and differences from a 3D decomposition of rotation, e.g. yaw, pitch and roll expressed by a composite 3D rotation matrix R=R.sub.z,R.sub.yR.sub.x.

[0042] Returning to FIG. 1, the task is to determine the geodetic position of each source and receiver in the arrays 120 and 130 to within a few metres. The position at the end of a time interval t can be predicted as x.sub.1=x.sub.0+vt where x0 is the position at the start of the interval and v is a constant velocity. If many such steps are added, truncation errors add and the predicted sum diverge or drift from the true position. The general idea of a Kalman filter is to compensate this drift by measurements, but not so much that the output estimates change rapidly due to uncertain observations, e.g. GPS readings. In addition, a Kalman filter includes calculated responses, e.g. a force applied through a bird 134 to oppose feathering.

[0043] The simple Kalman filter is a recursive algorithm where each time step k includes a prediction phase essentially comprising the calculations, and an update phase for collecting observations (sensor measurements) and computing output estimates from the step. The output estimates are weighted sums of predictions and observations, and are passed on to the next time step (k+1).

[0044] The calculations in the prediction phase depend on a dynamic model of the system. We have already approximated velocity by a constant velocity v, and note that constant acceleration does not represent a loss of generality due to the mean value theorem over the time interval t. Thus, a simple dynamic model of the present system is x=vt+at.sup.2/2 with constant v and a applied to all elements in the system 100. Real embodiments may include, for example, a response to feathering forces acting on a length of the streamer. However, we know in advance that any model is imperfect: The present model includes truncation errors and round-off errors in the approximations of v and a. Best practice is to keep the dynamic model simple and leave as much as possible to unmodelled dynamics described below. If this approach fails, more features such as tension, curvature and/or rotational forces may be included in the dynamic model.

[0045] In the present example, only the paravanes 114, birds 134 and steerable tail buoys 133 can exert a force. In the prediction phase, all other elements are just shifted by x=vt, where v is the instantaneous towing velocity at the start of time interval t. Imagine that all sensor positions x in geodetic coordinates are stored in a large computer array. The positions would be stored in an ordered manner, for example as a look-up table where each row represent a streamer and element E[r, i] in row r represent the X or Y coordinate of sensor no i along streamer no. r. For numbers, 32b (4B) per dimension provides negligible round-off errors on the scale 1/20000. Using 32b numerical resolution, 50 000 geodetic positions in 2D (X, Y) would need about 400 kB of storage. Thus, several or all rows in the lookup table would fit nicely in the computer memory of an ordinary PC. Adding a constant vt to all elements in an array stored in fast hardware (computer memory) is handled by hardware during a few cycles in most CPUs and GPUs. Thus, we do not worry about the execution time of 10 000-30 000 additions in a loop, and remark that such straightforward parallelisation generally decreases execution times by orders of magnitude.

[0046] Next, we need a model for the streamers. The streamers are not mechanically coupled, so an appropriate dynamic model should represent each streamer as an individual curve in space. B-splines are preferred for their local support property and several other useful properties. B-splines are widely used in other fields of technology for similar reasons. For example, B-splines limit the number of pixels that must recalculated in a computer image with a moving foreground and static background.

[0047] FIG. 2 illustrates the general setting. We want to pass a smooth (C.sup.2-continuous) B-spline curve C(s) through a sequence of data points D.sub.i in (x, y, z) coordinates. Specifically, this means finding a sequence of control points P.sub.i in an enclosing polyline. Moreover, the polynomial degree of C(s) should be as low as possible for a best possible approximation of the streamer. This implies polynomials with degree j=3, i.e. cubic B-splines. The bolded notation C(s) means that the scalar parameter s generates separate curves x(s), y(s) and z(s), each of which is independent of polynomial degree. Once C(s) is found, the position of any seismic receiver 132 can be found by polynomial interpolation using de Boor's algorithm, which is known to be fast and numerically stable.

[0048] The seismic receivers 132 are commercially available devices used to sample a P-wavefield and typically contain hydrophones and/or other seismic sensors. Any suitable seismic receiver 132 may be included in a streamer 131. For later reference, we note that the tension in a 10-20 km long streamer may cause significant strain. However, the distance between two adjacent seismic receivers 132 will remain approximately constant during a survey. The two distances c.sub.k and c.sub.+1 on streamer 131b illustrate that distances between adjacent elements are approximately constant but not necessarily equal. Each inline distance c.sub.k is known in advance, e.g. measured during deployment.

[0049] Birds 134 and position sensors 135-137 are also included in a streamer. A bird 134 may include a position sensor, e.g. an acoustic transducer 137 for determining distances to nearby acoustic transducers 137. The dashed lines between streamers 131a and 131b in FIG. 2 illustrate distances between such acoustical transducers, not mechanical connections.

[0050] Other common position sensors shown in FIG. 2 include a reflector 135 mounted above the sea surface on a tail buoy 133 attached to the trailing end of streamer 131a. A ranging laser/interferometer system operating in the visible or infrared range determines the distance from vessel 110 to the reflector 135 as illustrated by a double-headed arrow. A laser ranging system has a limited range, so the reflector 135 might be mounted on a buoy closer to the towing vessel in an array with several km long streamers. As the reflector 135 bobs up and down and tilts back and forth due to waves on the sea surface, measurements of range and bearing will vary from one time interval to the next without a corresponding motion of the tail of streamer 131a. This is an example of observation noise and illustrate that a real variance in R.sup.3 may be an ellipsoid.

[0051] Streamer 131b is connected to a separate tail buoy or float 133 with a GPS antenna 136 above the sea surface. The associated data point has a different mean position and variance than those of the GPS antenna on vessel 110.

[0052] Other common and useful sensors not shown in FIG. 2 include, but are not limited to, inertial measurement units (IMUs) with MEMS-accelerometers, flow meters for measuring size and direction of underwater currents and stress or strain gauges mounted at any location where stress or strain is an issue, e.g. at the towing points 113. The skilled person is familiar with these and other suitable sensors.

[0053] The data points D.sub.i may be measured by any means. Referring to FIG. 1, the length of lead-in cables 112, positions of the paravanes 114 and known lengths in the spreader cable 115 may help to determine the towing point 113 for a streamer 131.

[0054] Only the last few data points D.sub.n2, D.sub.n1, D.sub.n are shown in FIG. 2. Do is the towing point 113 at the leading end and the last point D.sub.n corresponds to the tail end of the streamer. The sequence of data points D.sub.i defines n curve segments C.sub.i(s)[D.sub.i1, D.sub.i>; 1<i<n. C.sup.2-continuity means that the curve C(s) is required to be continuous and have continuous first and second derivatives C(s), C(s) at each data point D.sub.i.

[0055] In the following, we will use some generally known properties and formulas. An overview can be found online, e.g. in B-spline curves: Important properties [1] or Wikipedia.

[0056] For our fit, we require that P.sub.0=D.sub.0 and that P.sub.m=D.sub.n. For a cubic B-spline, this implies that m=n+2. The parameter value s.sub.i that corresponds to the joint at D.sub.i is called a knot, and the sequence of knots in ascending order is called the knot vector. According to common practice, the knot vector is normalised to [0, 1] to improve numerical stability due to the relatively high density of floating numbers in this interval. Referring to FIG. 2, a cumulative chord length parameterization seems reasonable. In recursive form:

[00003]$\begin{array}{cc}{s}_0=0.Math..Math.s_k=s_{k-1}+\frac{c_k}{L_{\mathrm{tot}}},1kn;& \left(4\right)\end{array}$

where
c.sub.k is a constant length between two adjacent sensors 132-137 on the streamer
L.sub.tot is the total length of the streamera measured value if strain is significant, and
n is the number of sensors 132-137 on the streamer

[0057] Note that this knot vector contains values s.sub.k for all sensors, while the curve fitting needs a much smaller knot vector with knots s.sub.i for the positioning sensors. Knot vectors representing s.sub.k resp. s.sub.i may conveniently be stored as rows in a look-up table similar to that for geodetic positions described above.

[0058] In order to clamp a B-spline of degree j to control points at its end, i.e. ensure that D.sub.0=P.sub.0 and D.sub.n=P.sub.m=P.sub.n+2, the first and last j+1 knots are repeated. Thus, for a cubic B-spline clamped at both ends s.sub.0=s.sub.1=s.sub.2=s.sub.3=0 and s.sub.n=s.sub.n+1=s.sub.n+2=s.sub.n+3=1 in the s.sub.i-knot vector.

[0059] The B-spline curve has the form

[00004]$\begin{array}{cc}C\left(s\right)=\underset{i=0}{\overset{n}{.Math.}}.Math.P_i.Math.N_{i,j}\left(s\right),s\left[0,1\right]& \left(5\right)\end{array}$

where N.sub.i,j(s) are B-spline basis functions defined by Cox-de Boor's recursion formula:

[00005]$\begin{array}{cc}{N}_{i,0}\left(s\right)=1.Math..Math.\mathrm{if}.Math..Math.s_{i}s<s_{i+1}.Math..Math.\mathrm{and}.Math..Math.0.Math..Math.\mathrm{otherwise}& \left(6.Math.a\right)\\ N_{i,j}\left(s\right)=\frac{s-s_i}{s_{i+j}-s_i}.Math.N_{i,j-1}\left(s\right)+\frac{s_{i+j+1}-s}{s_{i+j+1}-s_{i+1}}.Math.N_{i+1,j-1}\left(s\right)& \left(6.Math.b\right)\end{array}$

[0060] We do not assume that all knot spans are equal, so a closed form of (6a, 6b) would look complicated. More important, the recursive expressions (6a), (6b) execute faster than closed form expressions would. Due to the step function in (6a), the first term in No(s) is nonzero only on knot span i, and the second term is nonzero only on knot span i+1. Each subsequent iteration of (6b) adds a knot span until 4 knot spans are added for j=3. This is the local support mentioned above.

[0061] In the case of cubic splines, at most 3 curve intervals C.sub.i(s) are valid on any knot S.sub.i, so the fitting criterion is of the form:

N.sub.i,3(s.sub.i)P.sub.i+N.sub.i+2,3(s.sub.i)P.sub.i+2=D.sub.i(7)

[0062] Now, P.sub.0=D.sub.0 and P.sub.n+2=D.sub.n so equation (7) yields n1 independent equations in n+1 unknowns. We need two extra conditions for a unique solution, and set the second derivatives to zero at the ends, i.e. C(0)=0 and C(1)=0. This seems reasonable in view of the streamer to be modelled.

[0063] We need expressions for the spatial derivatives for later reference:

[00006]$\begin{array}{cc}{C}^{}\left(s\right)=\frac{\mathrm{dC}\left(s\right)}{\mathrm{ds}}=\underset{i=0}{\overset{n-1}{.Math.}}.Math.N_{i,2}\left(s\right).Math.Q_i;Q_i=\frac{3}{s_{i+4}-s_{i+1}}.Math.\left(P_{i+1}-P_i\right)& \left(8\right)\\ C^{}\left(s\right)=\frac{d^2.Math.C\left(s\right)}{{\mathrm{ds}}^2}=\underset{i=0}{\overset{n-2}{.Math.}}.Math.N_{i,1}\left(s\right).Math.R_i;R_i=\frac{6.Math.\left(P_{i+1}-P_i\right)}{\left(s_{i+4}-s_{i+1}\right).Math.\left(s_{i+3}-s_{i+1}\right)}& \left(9\right)\end{array}$

where
the knot vector for (8) is the original with one copy of so and one copy of s.sub.n removed and
the knot vector for (9) is the original with second copies of so and s.sub.n removed.

[0064] In our curve fitting of a cubic spline, we have s=s.sub.0=s.sub.1=0 corresponding to P.sub.0=D.sub.0 in (x, y, z) Similarly, P.sub.n+2=D.sub.n correspond to s.sub.n=s.sub.n+1=s.sub.n+2=1. Thus, the first and last three terms of (9) yields a.sub.1P.sub.1+a.sub.2P.sub.2=D.sub.0 and a.sub.nP.sub.n+a.sub.n+1P.sub.n+1=D.sub.n by simple algebra. Together with (7), these yield a system with n+1 equations in n+1 variables:

[00007]$\begin{array}{cc}A\left[\begin{array}{c}{P}_1\\ P_2\\ \mathrm{.Math.}\\ P_{n+1}\end{array}\right]=\left[\begin{array}{c}{D}_0\\ D_1\\ \mathrm{.Math.}\\ D_n\end{array}\right]& \left(10\right)\end{array}$

where A is a banded matrix with at most 3 non-zero elements per row.

[0065] Standard methods for solving (10) without numerically expensive inversion include LU-decomposition and QR-decomposition with forward substitution. In a variant of forward substitution, we could solve the upper left 33 matrix of A and compute each next control point P.sub.i from the previous 3 control points and the next D.sub.i.

[0066] The system (7) corresponds to a fat matrix, i.e. a matrix with more columns than rows (m>n), and the second derivatives need not be zero. In practice, minimal norms are widely used to provide the missing equations for a fat matrix system. For example, a least squares fit (minimal L.sub.2 norm) is appropriate in many applications. For B-splines, a minimal polyline (minimal L.sub.1 norm) might be an alternative due to the strong convex hull property of B-splines. In general, the QR decomposition involves an nn matrix R.sub.i padded with zeros in the last m-n rows. Another useful concept for numerical applications is pseudo inverse matrices A.sup. which are sufficiently inverse to make A.sup.AI. More precisely, a typical pseudo inverse matrix is computed by a recursion that stops when a residual is below some limit, for example a predefined tolerance of a few metres or a limit due to truncation or round-off errors. These and other algorithms from numerical linear algebra are available in software libraries and need no further explanation here.

[0067] The s.sub.k for any sensor 132-137 is found from a lookup table, cf eq. (4). C(s.sub.k) and de Boor's algorithm yields the position of the sensor in vessel-bound coordinates (x, y, z).

[0068] Next, we need a TNB-frame for birds 134, flow meters measuring sea current etc. The original Frenet-Serret formulas regard the kinematics of a particle moving along a continuous and differentiable parametric curve in R.sup.3 and define a TNB-frame with Cartesian unit vectors tangent, normal and binormal to the parametric curve. They are valid for any real scalar parameter s, not just the special case in equation (4). Again, we only need a few general results, and refer to textbooks or online articles for a thorough description.

[0069] Noting that C(s) is tangent to C(s), we express the TNB-frame in matrix notation:

[00008]$\begin{array}{cc}T\left(\stackrel{\_}{s}\right)=\left[\begin{array}{c}{T}_x\left(\stackrel{\_}{s}\right)\\ T_x\left(\stackrel{\_}{s}\right)\\ T_z\left(\stackrel{\_}{s}\right)\end{array}\right]=\frac{C^{}\left(\stackrel{\_}{s}\right)}{.Math.C^{}\left(\stackrel{\_}{s}\right).Math.}& \left(11\right)\\ N\left(\stackrel{\_}{s}\right)=\left[\begin{array}{c}{N}_x\left(\stackrel{\_}{s}\right)\\ N_y\left(\stackrel{\_}{s}\right)\\ N_z\left(\stackrel{\_}{s}\right)\end{array}\right]=\frac{C^{}\left(\stackrel{\_}{s}\right)\left[C^{}\left(\stackrel{\_}{s}\right)C^{}\left(\stackrel{\_}{s}\right)\right]}{.Math.C^{}\left(\stackrel{\_}{s}\right).Math..Math..Math.C^{}\left(\stackrel{\_}{s}\right)C^{}\left(\stackrel{\_}{s}\right).Math.}& \left(12\right)\\ B\left(\stackrel{\_}{s}\right)=TN& \left(13\right)\end{array}$

where C={square root over (C.sub.x.sup.2+C.sub.y.sup.2+C.sub.z.sup.2)} is the Euclidian norm, and a cross product in matrix form is

[00009]$\begin{array}{cc}c=ab=\left[\begin{array}{c}{c}_x\\ c_y\\ c_z\end{array}\right]=\left[\begin{array}{c}{a}_2.Math.b_3-a_3.Math.b_2\\ a_3.Math.b_1-a_1.Math.b_3\\ a_1.Math.b_2-a_2.Math.b_1\end{array}\right];.Math..Math.a=\left[a_x,a_y,a_z\right],b=\left[b_x,b_y,b_z\right]& \left(14\right)\end{array}$

[0070] A unit normal is undefined for a straight line and need special consideration: Let D.sub.i1, D.sub.i, D.sub.i+1 be corners in a triangle. If the height of the triangle is greater than a specified tolerance, compute N and B at D.sub.i from (12) and (13). Otherwise, set the values for N and B equal to those of D.sub.i1. Here, the maximum height of the triangle corresponds to half the tolerance specified for the (x, y, z) positions D.sub.i.

[0071] T and N are shown in FIG. 2. In right-handed coordinates, B would point downward or into the paper plane in FIG. 2 which is the preferred direction for seismic surveying applications. C(s) is efficiently calculated by de Boor's algorithm using the B-spline in equation (9). C(s) and its Euclidian norm from eq. (11) are reused in eq. (12). Similarly, de Boor's algorithm on the B-spline in eq. (9) yields C(s).

[0072] Derivation with respect to time does not change the direction of a vector, so measurements of position, orientation, velocity and acceleration in a TNB-frame can be transformed to any coordinates of choice using eq. (2) and pre-multiplication.

[0073] The model of streamers described so far may be enhanced by fitting B-splines in the crossline direction and/or by defining a polygonal, e.g. triangular, mesh over acoustic transducers 137. Spline fitting in 1D and 2D is available in readymade software, e.g. graphics libraries. Some 2D splines, e.g. thin plate splines, assume elastic properties in two spatial directions and are inappropriate for decoupled streamers.

[0074] Assume the framework above is used to define a desired streamer configuration in a GUI and that the computer model provides a list of desired data points d.sub.i. The Kalman filter maintains a similar list of estimated data points D.sub.i. Both lists are in vessel-bound coordinates (x, y, z) and are thus independent of towing velocity and associated time shifts in geodetic coordinates. In each time interval t, changing sign on D.sub.i and adding to d.sub.i yields a tentative response. As above, operations can be performed on entire rows in computer memory. All elements in the tentative response array smaller than a predefined constant tolerance are ignored, those remaining are proper responses to a deviation.

[0075] For each such deviation, equation (7) defines 3 control points that must be recomputed to fit the cubic B-spline C(s) through the deviant D.sub.i. No other part of the B-spline is affected due to the local support property. Thus, the curve fitting means solving a 33 system with known coefficients. QR-decomposition with forward substitution is effective also in the rare cases where several adjacent data points are out of line within t.

[0076] Once the three control points are recomputed, all sensor positions on the 4 affected knot spans must be recomputed using de Boor's algorithm on a sequence from the s.sub.k knot vector. The results are transformed to geodetic coordinates and replace elements in the geodetic lookup table. A few of these elements, e.g. birds 134 and flow meters may need additional computation for a new TNB-frame.

[0077] A GUI may conveniently have a local view for streamer configuration etc. in vessel-bound coordinates and a geodetic view for showing common midpoints and sea current fields in geodetic coordinates as described above. B-splines are affine invariant, so a graphics engine, i.e. graphics software and a graphics card, just need a few control points to display the curve. In addition, any such graphics engine combine affine transformations and thus enable zooming in on any part of the streamer array and viewing the section or the entire streamer array in perspective and from any angle in real time. The required input is just a few control points in vessel-bound coordinates. The local sea current may be displayed as a nicely coloured graded background. For this, the graphics engine typically uses Bzier curves, which are B-splines without knots between the clamping zeros at the start and end of the knot vector.

[0078] The Kalman filter may have several external degrees of freedom such as the geodetic and vessel-bonds positions above. The number of internal states and variables is a matter of design. In the following, we use a bird 134 as example, and add velocity and acceleration as internal states. Using the notation above and dot notation for time derivatives, we let:

[00010]$\begin{array}{cc}{x}_k=\left[\begin{array}{c}x\\ \stackrel{.}{x}\\ \stackrel{\mathrm{.Math.}}{x}\end{array}\right].Math..Math.F_k=\left[\begin{array}{ccc}1.Math.& .Math..Math.t& \frac{.Math..Math.t}{2}\\ 0& 1& .Math..Math.t\\ 0& 0& 1\end{array}\right].Math..Math.\mathrm{.Math.}.Math..Math.B_k=\left[\begin{array}{ccc}{b}_{11}& b_{12}& b_{13}\\ b_{21}& b_{22}& b_{23}\\ b_{31}& b_{32}& b_{33}\end{array}\right].Math..Math.\mathrm{.Math.}.Math..Math.u_k=\left[\begin{array}{c}0\\ 0\\ \mathrm{ca}\end{array}\right]& \left(15\right)\end{array}$

where x is in geodetic coordinates and B transforms from the bird's TNB-frame.

[0079] Let index min denote the state up to and including n. The prediction phase estimates a preliminary (a priori) state in step k from the output (a posteriori) estimates in step (k1):

x.sub.k|k-1=F.sub.k{circumflex over (x)}.sub.k-1|k-1+B.sub.ku.sub.k P.sub.k|k-1=F.sub.kP.sub.k-1|k-1F.sub.k.sup.T+Q.sub.k(16)

[0080] F.sub.k is the state transition model that brings the filter from step (k1).

[0081] B.sub.k is a control model and u.sub.k a control vector that express the response calculated for step k, e.g. based on an estimated local current and properties of the bird. The resulting force is transformed first to ka in vessel-bound coordinates, then to the relevant TNB-frame by ca=B.sup.1 ka. The two transformations are combined into one matrix multiplications by eq. (2). Paravanes, steerable tail buoys and other equipment that can apply a force have Bu-terms, whereas seismic receivers 132, acoustic transducers 137 etc. lack the Bu term.

[0082] The P's are prediction error covariance matrices defined by recursion. If little or nothing is known in advance, a common assumption is random walk leading to P.sub.0|0=I.sup.2, where I is the identity matrix and the system variance .sup.2 some scalar, e.g. 0.5. A properly designed filter converges, so the effort spent on estimating optimal values should be limited. The present KF will have performed 120 five second steps in ten minutes.

[0083] Q is the covariance of process noise w.sub.kN(0, Q). The process noise may contain unmodelled dynamics and other systematic faults that may cause the filter to diverge and that are difficult to distinguish from observation noise. Common practice is to compute the covariance matrix Q by statistical methods, e.g. an ALS technique. ALS is an acronym for auto-covariance least-squares and essentially uses data from several previous steps to look for systematic deviations. A comprehensive description can be found online or in the literature.

[0084] The update phase regards measurements, weighted sums and the output state. A first task is to find an observation residual from present and previous measurements:

{tilde over (y)}.sub.k=z.sub.kH.sub.k{circumflex over (x)}.sub.k|k-1(17)

[0085] The observation vector z.sub.k for a specific bird correspond to a data point D.sub.i in FIG. 2. However, the D.sub.i's were expressed in (x, y, z), whereas x and z are expressed in (X, Y). The observation model H simply picks previous estimates. For example, H=[1 0 0] picks a position estimate corresponding to a D.sub.i for a bird 134 or an acoustic transducer 137. If the bird also provides a velocity measurement, e.g. representing sea current, H would have an additional 1 in the velocity position, i.e. H=[1 1 0]. In an alternative version of the KF, the measurements from acoustic transducers 137 could be treated as a fourth degree of freedom and be picked by H=[0 0 0 1].

[0086] We also need a residual covariance S.sub.k=H.sub.kP.sub.k|k-1H.sub.k.sup.T+R.sub.k, where R is the covariance of observation noise and usually approximated by ALS.

[0087] The optimal Kalman gain K.sub.k=P.sub.k|k-1H.sub.k.sup.TS.sub.k.sup.1 is optimal in the sense that it minimises the least square difference between x and its estimate {circumflex over (x)}. Its elements are real scalars in the interval [0, 1]. The expressions for S and K may look complicated, but are easily implemented because H simply picks a previous variable as described above.

[0088] The output state or a posteriori estimates in step k are

{circumflex over (x)}.sub.k|k={circumflex over (x)}.sub.k|k-1+K.sub.k{tilde over (y)}.sub.k P.sub.k|k=(IK.sub.kH.sub.k)P.sub.k|k-1(18)

[0089] Large elements Kij in the Kalman gain K.sub.k, i.e. values in the interval <0.5, 1], puts more weight on observations and decrease the weight on prediction accordingly. The weighted sums in (18) are the values input to the next step (k+1).

[0090] In real embodiments, simple Kalman filters are numerically unstable if the values in the covariance matrices become small. The reason is that the covariance matrices P, Q, R by definition are positive definite, but small elements may cause negative round-off errors that render the matrices indefinite. The standard solution is to represent P, Q, R in square root form, preferably as UD-decompositions. Efficient algorithms for the Kalman prediction and update phases in square root form can be found in Thornton [2] and Bierman (1977) [3].

[0091] Adding velocity to the previous example with geodetic positions, the a posteriori position and velocity estimates from step (k1) can be used to extrapolate positions during step k. For a towing speed 2.6 m/s (5 knots), this linear approximation would improve geodetic estimates by about 1.3 m compared to the position output from step (k1). In a 2D-grid (X, Y) this extra degree of freedom would require an extra 400 kB of storage using the numbers above. Of course, the storage need will be different in a real embodiment.

[0092] During a turn between sailing lines S.sub.q-1, S.sub.q streamers are thrown out or separated to minimise operational cost and the risk for tangling. The shapes of the leftmost streamers in FIG. 1 resembles their shapes during such a throw out. A Kalman filter providing accurate position estimates reduces a need for safety margins due to uncertain position measurements. If velocity were added as a degree of freedom to such a KF, it would be possible to compare accurate crossline velocities to predetermined values optimised to minimise turning time. The added degrees of freedom would increase storage requirements, but the added processing power would only be marginally raised due to the local support property of B-splines.

[0093] The B-splines curves are expected to approximate physical streamers sufficiently well. Referring to FIG. 1, even during a right turn with maximum crossline speed at the tail buoy of the leftmost streamer, bird no. 3 from the end is unlikely to move much in t5 s.

[0094] The Kalman filter as such is cheap on storage and processing power. ALS techniques calculate autocorrelation on a series of previous estimates and are costlier. However, a rough estimate of a MB per degree of freedom and time lags of 1000 steps yield storage needs in the GB range. Thus, an inexpensive PC-class computer running reasonably effective common operating and filesystems, e.g. some Linux dialect with ext4, should be able to run the proposed KF in real-time. Alternatively, a lightweight real-time operating system installed on the PC-class computer would remove further dead cycles. A specialised file system may improve mapping from filesystem to memory compared to general purpose filesystems.

[0095] The above and other embodiments are within the scope of the present invention. The skilled person will recognize these and other applications and modifications of the invention defined in the appended claims.

REFERENCES

[0096] [1] http://www.cs.mtu.edu/shene/COURSES/cs3621/NOTES/spline/B-spline/bspline-curve-prop.html [0097] [2] Thornton, Catherine L: Triangular Covariance Factorizations for Kalman Filtering (PhD thesis), Nasa Technical Memorandum 33-798, 1976 [0098] [3] Bierman, G. J: Factorization Methods for Discrete Sequential Estimation, Academic Press, 1977