Marine Streamer Position Determination System

The present invention relates to marine seismic exploration, and in particular to a method and system for determining the position of a towed marine seismic streamer.

In marine seismic exploration, an exploration vessel tows a seismic streamer which houses a plurality of pressure sensitive detectors commonly referred to as hydrophones. Impulsive sources such as explosives or air guns, which may also be towed by the exploration vessel, are fired so as to impart propagating energy to the body of water, and ultimately to the underlying crustal material of the earth. As this energy propagates in a generally downward direction, portions of the energy are reflected by subsurface discontinuities and ultimately detected as pressure variations by the hydrophones. The corresponding electrical signals are coupled by the streamer to recording apparatus aboard the exploration vessel for subsequent use in interpretation of the subsurface crustal structure. For such interpretations to be meaningful, however, it is necessary to know with some precision the location of the exploration vessel at the time the data are required. A number of sophisticated navigational systems have been developed for providing this information. One such system, for example, is the GEONAV* seismic marine navigation and location system manufactured by Texas Instruments Incorporated of Dallas, Texas.

While such systems are capable of providing reliable measurements of the location of the exploration vessel, there still remain ambiguities in knowledge of the position of the remote parts of the streamer relative to the track followed by the vessel. It is rare, for example, that the streamer trails directly along track. While the streamer is attached to the vessel at a point above the water surface, the majority of the streamer, which may be a mile or greater in length, is typically maintained at a nominal depth below the surface of the water by a plurality of depth controllers located along the length of the streamer. One such depth controller, for example, is the subject of U.S. Pat. No. 3,372,666 to Buford M. Baker, and assigned to the assignee of the present application. One reason for failure of the streamer to trail along the track of the vessel is the fact that the cross-track current velocity at the depth of the streamer often differs from the cross-track current velocity operating on the hull of the vessel itself. Other factors contributing to this problem are boat crab angle and along track current gradients.

Previously, estimates of the streamer position have been made using a radar reflector or transponder located on a surfaced tail buoy which is attached to the end of the streamer remote from the exploration vessel. This permits estimation of the location of the remote end of the streamer, but with inaccuracies which are too great to be compatible with the resolution provided by modern interpretive techniques available for the seismic data itself. Moreover, the radar location technique assumes a straight line path for the streamer between the stern of the vessel and the tail buoy.

It is, therefore, an object of the present invention to provide an improved method and system for the determination of the location of a towed marine seismic streamer.

It is another object to provide a seismic streamer position determining method and system which employs measurements which are conveniently made on the seismic streamer.

It is a further object of the invention to provide a method and system for determining the position of a marine seismic streamer wherein the complexity of the method and system is related to the accuracy of the results desired.

It is yet another object of the invention to provide apparatus for determining the orientation of that portion of the seismic streamer adjacent to the towing vessel.

In a paper entitled "Three-dimensional Boundary Value Problems for Flexible Cables", John W. Bedenbender, Second Annual Off-Shore Technology Conference, Houston, Texas, April 22-24, 1970, Paper No. OTC1281, there is developed a mathematical model for the equilibrium configuration of a marine seismic cable, in the form of six differential equations. The paper further presents a technique for solving the set of differential equations. This technique, however, requires knowledge of a large number of parameters, some of which are poorly known, if at all, in a towed seismic streamer situation. However, using the methods presented in the paper for the study of simulated cases, we have observed that the lateral displacement of the streamer is closely dependent upon the yaw angle of the cable in an earth-relative coordinate system (measured in the horizontal plane), relative to the track traversed by the vessel.

Accordingly, in one aspect of the invention, there is employed a unique apparatus for coupling the marine seismic streamer to the exploration vessel. The body of the apparatus may be attached fixedly at the stern of the vessel. The marine seismic streamer is coupled to the body of the apparatus by a universal joint. In the preferred embodiment, the apparatus includes two synchro units, which cooperate with the universal joint to provide signals representative of the yaw and pitch of the portion of the seismic streamer adjacent the vessel. These measurements of yaw and pitch are given in a ship-relative coordinate system. The apparatus further includes pendulum units, which provide signals representative of the roll and pitch of the exploration vessel. These signals are combined with the signals provided by the synchro units to yield measures of the cable yaw and pitch in an earth-relative coordinate system. As will be disclosed in greater detail subsequently, these measurements of yaw and pitch are used to develop models of varying degrees of complexity to represent the position of the streamer.

In the case of one of the models, there is employed additionally a magnetic compass, located at a point on the streamer remote from the exploration vessel. The signal provided by the magnetic compass, after correction for known magnetic variation in the area, is representative of the lateral orientation of the streamer at the compass location. The model comprises a two-line approximation of the position of the streamer. The first of the lines, passing through the stern of the exploration vessel, is constructed tangent to the measured yaw of the streamer at the vessel. The second of the lines is constructed parallel to the corrected compass indication and intersects the first line at a point estimated by an empirical relation.

Other models which may in some cases provide more accurate estimates of the location of the streamer will be discussed hereinbelow.

A more complete understanding of the invention, as well as of other features and advantages thereof, may be had by a consideration of the following detailed description in connection with the attached drawings wherein:

FIG. 1 is a plan view of a vessel towing a marine seismic streamer.

FIG. 2 illustrates the yaw and pitch angles of the streamer.

FIG. 3 shows a fixture for connecting the streamer to the vessel and for measuring the yaw and pitch angles.

FIGS. 4a and 4b illustrates the connection of the streamer to the fixture and the mounting of a tensiometer in the streamer.

FIG. 5 is a block diagram of a stramer positioning system.

FIG. 6 is a flow chart for the generation of a position model in the preferred embodiment.

FIG. 7 illustrates the position model of the preferred embodiment.

FIGS. 8a through 8i a flow chart for the generation of a position model in an alternate embodiment.

In FIG. 1, there is shown diagramatically an exploration vessel towing a marine seismic streamer, the view representing that which would be seen looking down on the surface of the sea from directly above. The vessel 10 travels along a track, which is represented by dashed line 12. While in FIG. 1 the longitudinal axis of vessel 10 is seen to be aligned with track 12, this will in general not be the case, the vessel generally being at some crab angle with respect to the track along which it is traveling. The position shown for seismic streamer 14 may be thought of as the location of the projection of the seismic streamer on the surface of the sea. Shown disposed at regular intervals along the arc length of the streamer, are a plurality of depth controllers 16 of the type previously discussed. While seismic streamers are commonly constructed to be approximately neutrally buoyant, the depth controllers 16 ensure that the portions of the cable at which they are affixed will be maintained at a constant predetermined depth below the surface of the sea. Thus, it will be seen that while the portion of the streamer adjacent vessel 10 will be affixed to the vessel at a point above the surface of the sea, the portion of the streamer at the first depth controller 16 will be maintained at the aforementioned predetermined depth. At the tail end of streamer 14, that is the end most remote from vessel 10, the streamer is affixed to a tail buoy 18, which serves to maintain the tail end of streamer 14 at the surface of the sea.

Line 20 in FIG. 1 is constructed in the horizontal plane of an earth-relative coordinate system and is tangent to the projection of streamer 14 on the surface of the sea. The point of tangency is the point at which streamer 14 is attached to vessel 10. Line 20 along with the longitudinal axis of vessel 10 define an angle .theta., which will be referred to as the yaw of the cable at the vessel 10. With the attitude of vessel 10 shown in FIG. 1, the longitudinal axis of the vessel lies along track 12.

Yaw angle .theta., as well as pitch angle.phi., are further illustrated in FIG. 2 which shows a unit vector 22. Unit vector 22 is codirectional with the streamer at the point where the streamer is attached to vessel 10, this point being given by the origin of the coordinate system in FIG. 2. In the earth-relative coordinate system illustrated in FIG. 2, the x-axis is coincident with the longitudinal axis of vessel 10, the longitudinal axis of vessel 10 being assumed, for the moment, to be codirectional with the track 12 of the vessel. The z-axis, which also lies in the horizontal plane, is the transverse axis of the vessel, while the y-axis is the vertical axis. The projection of unit vector 22 on the horizontal plane is illustrated by dashed vector 24, while the projection of unit vector 22 on the x-axis is represented by that portion of the x-axis terminating at point 26. It will again be seen from FIG. 2 that yaw angle .theta. is the angle between the horizontal projection of the streamer at the vessel and the ship's longitudinal axis. Pitch angle .phi. is the angle between the unit vector, representing the streamer at the point where it meets the vessel, and the horizontal projection of this unit vector.

An apparatus for the measurement of yaw angle .theta. and pitch angle .phi. at the stern of vessel 10 is illustrated in FIG. 3, generally at 30. The apparatus includes a support structure 32, which is fixedly attached to the deck at the stern of vessel 10, by means of mounting bolts 34. Bracket 36 is coupled to support structure 32 by means of universal joint 38. Cable 40 is coupled to bracket 36 by means of pulley 42. As seen in FIG. 4a, cable 40 is coupled at its other extremity to an exchange tool 70, which in turn supports the end of seismic streamer 14, adjacent vessel 10. Since universal joint 38 permits two degrees of freedom between bracket 36 and support structure 32, it will be seen from the foregoing that the face 44 of bracket 36 will be maintained at all times normal to the unit vector which represents the orientation of that portion of seismic streamer 14, adjacent vessel 10. Accordingly, the orientation of bracket 36, relative to support structure 32, is seen to be directly representative of the orientation of seismic streamer 14, relative to vessel 10.

Rotation of bracket 36, relative to support structure 32, about an axis parallel to the vertical axis of vessel 10, is sensed by synchro unit 46. This angular rotation is coupled from universal joint 38 to synchro unit 46 by means of flexible shaft 48. As a result, the electrical signal provided by synchro unit 46 on electrical cable 50, is representative of the yaw angle .theta.of streamer 14, relative to a ship's coordinate system.

In a similar manner, rotation of bracket 36, relative to support structure 32, about a transverse axis, is coupled by a second flexible shaft (not shown) to synchro unit 52, which provides on electrical cable 54 a signal representative of the pitch .phi. of seismic streamer 14. Again, the signal on electrical cable 54 represents the pitch .phi. measured in a ship-relative coordinate system. The synchros in units 46 and 52 may each be of a type commercially available as Model No. 1802880 from Sperry Marine Systems.

In the preferred emboidment, it will be desired that the yaw angle .theta. and pitch angle .phi.be expressed in an earth-relative coordinate system. Accordingly, there is rigidly attached to support structure 32, a pendulum unit 56. Pendulum unit 56 contains a first pendulous potentiometer which is constrained in one dimension, so as to be capable of rotating freely only about the longitudinal access of vessel 10. It will be seen, therefore, that this first potentiometer is sensitive to roll maneuvers of vessel 10. Pendulum unit 56 contains a second pendulous potentiometer, which is also constrained in one dimension, so as to be able to rotate freely only about the transverse axis of vessel 10. Accordingly, this second potentiometer is sensitive to pitch maneuvers of vessel 10. The two electrical signals provided by these two potentiometers are available on electrical connector 58 for further processing. A device satisfactory for use as the pendulous potentiometer in unit 56 is commercially available from Humphrey Inc. of San Diego, California as Model No. CP17-0601-1.

Again, referring to FIG. 2, it will be seen that the x, y, and z components of unit vector 22 will have the values shown in equation (1), where .theta. is the cable yaw angle and .phi. is the cable pitch angle, both in an earth-relative coordinate system. The electrical signals appearing on electrical cables 50 and 54, respectively, however, are representative of the cable yaw and pitch angles relative to the ship's coordinate system, whose longitudinal, vertical, and transverse coordinate axes will be labeled 1, v, and t, respectively. The 1, v, and t components of unit vector 22 are given by equation (2), where .theta..sub.s and .phi..sub.s are the cable yaw and pitch angles relative to the ship's coordinate system, these angles being represented by the signals appearing on electrical cables 50 and 54. ##EQU1## Unit vectors u.sub.22 and u.sub.s are related by transformation matrix [T] as shown in equation (3).

The individual elements of transformation matrix [T] are given in equation (4). ##EQU2## In equation (4), the quantities R and P are the ship roll and pitch angles, respectively, these angles being represented by the signals appearing on electrical connector 58. The signals appearing on electrical connectors 50, 54 and 58 are combined in accordance with equation (3), so as to provide the x, y, and z components of unit vector 22 as shown in equation (1).

The individual components of the cable unit vector u.sub.22, given in equation (1), are combined as shown in equation (5) to expressly yield the cable yaw and pitch angles .theta. and .phi..

Where

U.sub.x is the x componeent of U.sub.22 ;

U.sub.y is the y component of U.sub.22 ; and

U.sub.z is the z component of U.sub.22.

To this point, the yaw angle .theta. is measured relative to the longitudinal axis of vessel 10. The angle .theta. is further combined with the heading angle of vessel 10 as provided by a gyrocompass, which comprises a portion of the ship's navigational system to re-express the yaw angle .theta. relative to north. It will be assumed from this point on that the computed yaw angle .theta. has been so corrected to be relative to north.

Fixture 30 cooperates with the coupling structure illustrated in FIG. 4a to couple the streamer 14 to the stern of vessel 10. With reference to FIG. 4a, cable 40 connects pulley 42 to exchange tool 70. Streamer 14 is rigidly affixed to exchange tool 70 at point 72. Emerging from streamer 14 at point 72 where it engages exchange tool 70 is a cable bundle 74. In a typical streamer arrangement, cable bundle 74 includes a pair of electrical conductors for each of the hydrophone groups disposed along streamer 14. In the present embodiment, cable bundle 74 will additionally include electrical conductors leading to tensiometers and a magnetic compass which are also located at points disposed along streamer 14. As is well known in the art, cable bundle 74 is coupled to signal recording and processing equipment located aboard vessel 10.

In the preferred embodiment, there is disposed a first tensiometer within streamer 14 at a point 74 near exchange tool 70. FIG. 4b is an exploded view of the interior of streamer 14 in the vicinity of the first tensiometer. A suitable tensiometer is a Transducer Load Cell Model No. WML2-251-10K-3926, commercially available from Transducers, Inc. of Whittier, California. In FIG. 4b, tensiometer 76 is seen to have a generally hexagonal cross-section. Located at either end of tensiometer 76, as expressly shown at 78, is a threaded recess by means of which tensiometer 76 is coupled to streamer 14. Four pins 80 are provided to permit external electrical connection to the Wheatstone bridge of tensiometer 76. A bracket 82 is provided to couple one end of streamer 14 to tensiometer 76. Bracket 82 includes a threaded stud 84 which engages with one of threaded recesses 78, so as to rigidly affix bracket 82 to tensiometer 76. The three stainless steel strain members 86, which comprise the tension bearing members in a conventional seismic streamer, are affixed to bracket 82 at points 88. Portions of bracket 82 are cut away as at 90 to permit the cable bundle of the streamer to pass through bracket 82, along side tensiometer 76, and on up the streamer 14 to vessel 10. Additional electrical leads connected to pins 80 also become part of the cable bundle leading to vessel 10. A second bracket similar to bracket 82 is provided to couple the forward end of streamer 14 to the forward end of tensiometer 76. As is well known in the art, the entire seismic streamer 14, including the portion comprising tensiometer 76, is encased within a sleeve or jacket which is typically extruded from a material such as polyvinylchloride.

It is common practice with marine seismic streamers to include in the portion of the streamer between the tow vessel and the first depth controller, a "stretch" section. The stretch section differs from the conventional streamer section in that the stainless steel strain members are replaced by nylon strain members which are relatively elastic. The presence of the stretch section reduces the probability of shock damage to the streamer and provides isolation of the hydrophones from mechanical noises induced by the tow vessel. In order to know the precise arc distance from the vessel to any point on the streamer beyond the stretch section, however, it is necessary to know the percentage stretch of the stretch section.

Tensiometer 76 provides a signal representative of the longitudinal tension acting on streamer 14 at the point where the streamer is coupled to the stern of vessel 10. In the preferred embodiment, a second similar tensiometer is located in streamer 14 in the vicinity of the depth controller 16 which is located nearest vessel 10. As discussed below in greater detail, the tension measurements provided by these two tensiometers may be employed to determine the degree of elastic expansion of that portion of streamer 14 intermediate to the two tensiometers.

A magnetic compass is disposed in streamer 14 at a point sufficiently remote from vessel 10 to ensure that distortions of the earth's magnetic field by the field of vessel 10 will not introduce unacceptable errors in the sensed magnetic "north". In the preferred embodiment, for example, the magnetic compass is located at a point about midway along the length of streamer 14. A suitable magnetic compass is Model 314 Heading Sensor, available from Digicourse, Inc. This particular model of magnetic compass has the desirable feature of being gimballed to be rotatable through a full 360.degree. about the longitudinal axis of streamer 14. As a result, the twists which commonly occur along the length of streamer 14, do not interfere with proper operation of the magnetic compass. Additional electrical conductors in the cable bundle couple the electrical signal provided by the magnetic compass to processing equipment aboard vessel 10.

The interconnection of the system components comprising the preferred embodiment of the invention is illustrated in the block diagram of FIG. 5. The various parameter sensing units employed in the system include tensiometers 100 and 102, magnetic compass 104, synchro units 46 and 52, and the two pendulous potentiometers in pendulim unit 56. The analog signals provided by tensiometers 100 and 102 are digitized in A/D converters 106 and 108, respectively. Similarly, the analog signals provided by the two potentiometers of pendulum unit 56 are digitized by A/D converters 116 and 118. The A/D converters may be any of a wide variety of commercially available units, one such suitable unit being the Model No. AD2003 panel meter with digital output produced by Analog Devices.

Magnetic compass 104, which as previously mentioned is a Model 314 Heading Sensor, is coupled to interface unit 110. In the preferred embodiment, interface unit 110 is a Model 251 Heading Sensor Interface Unit also produced by Digicourse, Inc. Interface unit 110 is coupled to block 120 wherein the pulse train provided by magnetic compass 104 and interface unit 110 is counted to obtain a digital indication of the compass reading. Unit 120 also serves to selectively switch power, under control of main processor 126, to magnetic compass 104 and interface unit 110. In a preferred method of operation, power to the magnetic compass is switched off at those times when seismic data are being acquired, so as to eliminate the possibility of noise contamination of the seismic data. The magnetic compass is powered up to provide a compass reading sample just prior to the initiation of a seismic disturbance and recording of the resulting seismic data.

The analog signals provided by synchro units 46 and 52 are digitized in S/D converters 112 and 114, respectively. S/D converters 112 and 114 may each be a Model M5000M36/3DP Synchro-to-Digital Converter manufactured and sold by AstroSystems, Inc. of Lake Success, New York.

The digital signals generated as discussed hereinabove are combined in main processor 126, so as to provide models of the streamer position. The main processor 126 may be a Model 980A computer, maufactured and sold by Texas Instruments Incorporated of Dallas, Texas. As illustrated by block 122, communication of data samples to main processor 126 is by means of a main processor I/O interface 122. The provision of the I/O interface to interrogate the various digital samples and communicate the digital samples to main processor 126 is well understood by those skilled in the art and requires no further elaboration here. By way of example, however, a portion of main processor I/O interface 122 has been expressly designated in FIG. 5 as I/O Data Module 124. This portion of the interface, which services the tensiometer and magnetic compass signals, comprises a Model 961648-2 I/O Data Module produced and sold by Texas Instruments Incorporated. Main processor I/O interface 122 periodically, i.e., once each second, acquires a data sample from each of the input units and passes these data samples to main processor 126 wherein the data samples are used to generate a model of the streamer position. The various data samples, as well as the model, are recorded on magnetic tape unit 128. Magnetic tape unit 128 may be a Model TI-979 magnetic tape unit produced and sold by Texas Instruments Incorporated.

The operations performed by main processor 126 in the preferred embodiment of the invention are illustrated by the flow diagram of FIG. 6. The processing used to obtain the streamer model for each set of samples from the various input units is initiated by entry at 140. At step 142, the streamer positioning system acquires the instantaneous latitude, longitude, and azimuth of vessel 10. In the preferred embodiment, the streamer position system functions in cooperation with the aforementioned GEONAV* seismic marine navigation and location system. In particular, main processor 126 is shared by the streamer positioning system and the marine navigation and location system. As a result, main processor 126 will contain at all times the most recent measurements of ship's latitude, longitude and azimuth, so that these quantities are always available for use by the streamer positioning system as at step 142. As will be discussed henceforth in greater detail, the ship's azimuth is used to permit computation of the streamer model in an earth-relative coordinate system. The invention, however, is not so limited and may be used in the case where ship's azimuth is not available to provide a streamer model expressed in a ship-relative coordinate system.

At step 144, main processor 126 acquires the most recent data sample from magnetic compass 104 which is positioned along streamer 14. At step 146, the raw compass reading is corrected to remove compass inaccuracies. Prior to locating compass 104 in streamer 114, the compass will have been tested to obtain an empirical curve expressing the true magnetic reading as a function of the raw compass reading. This empirical curve is preloaded into main processor 126 in the form of an approximation comprising eight linear segments. At step 146, main processor 126 applies the raw compass reading to this curve to determine the true magnetic reading. At step 148, the magnetic variation is removed from the compass reading to provide a heading expressed relative to true North. Magnetic variation is input to main processor 126 by an operator and may be updated periodically so that the magnetic variation employed by the main processor is appropriate for the area under exploration.

At step 149, the ship's azimuth is subtracted from the adjusted compass heading to provide the angle .theta.2 which expresses the orientation of the streamer relative to the ship's coordinate system.

At step 150, the processor acquires the most recent samples of streamer yaw and pitch, and of the ship's roll and pitch from fixture 30. At 152, the sines and cosines of these angles are determined, these trigonometric functions being employed as in equations (1) - (4) to obtain a streamer unit vector expressed in ship-relative coordinates. At step 154, the quantity D, which is defined in equation (4) is computed. At step 156, The processor determines the cable unit vector expressed in a ship-relative coordinate system as defined by equation (2). At step 158, the rotation represented by equation (3) is accomplished. Then, using the components of the cable unit vector expressed in an earth-relative coordinate system, the cable yaw angle .theta.1 and pitch angle .phi. are determined at step 160 in accordance with equation (5).

At step 162, the most recent samples provided by tensiometers 1 and 2 are acquired by main processor 126. At step 164, the difference between the tensions measured by tensiometers 1 and 2 are used to determine the percentage of stretch in the stretch section of the cable. It has been determined empirically that the percentage stretch, when expressed as a function of the difference in tension at the two ends of the stretch section, may be approximated by a quadratic function. The appropriate quadratic function is preloaded into main processor 126 and employed at step 164 in connection with the difference between the two measured tensions to yield the percentage stretch. At step 166, this percentage stretch is multiplied by the equilibrium length of the stretch section to provide the stretch length of the streamer. At step 168, this stretch length is added to the arc distances of the streamer head and tail, that is, the lengths along the unstressed streamer from the stern of vessel 10 to the head and tail, respectively, of the active portion of streamer 14. The streamer head in this context is that point along the streamer at which are located the hydrophones nearest vessel 10.

The two-line approximation comprising the streamer position model in the preferred embodiment of the invention is illustrated in FIG. 7. The first line 172 of the model passes through the stern of vessel 10 and is oriented at an angle .theta.1 with respect to the x or longitudinal axis of vessel 10. The second line 174 is at an angle .theta.2 with respect to the x axis. The point at which lines 172 and 174 intersect has ccordinates (XB, ZB), which are estimated at step 170 of FIG. 6 by the empirical relation given in equation (6).

Zb = xb (tan .theta.1)

where:

at step 176, the polar coordinates in a North-relative system of the head and tail of the streamer are computed, using their arc distances, the angles .theta.1 and .theta.2, the coordinates XB and ZB, and the ship azimuth. A transformation is also applied to the location of the breakpoint (XB, ZB) to express this breakpoint in the polar coordinate system. At step 178, the raw data, the streamer head and tail locations, the breakpoint location in polar coordinates, and the ship latitude, longitude and azimuth are output to the magnetic tape unit 128.

At step 180, the next sample from the magnetic compass 104 is read in preparation for the next computation. After this next sample from the compass is acquired, a compass "idle"command is sent to the data module in step 182. As explained previously, this places the magnetic compass in a powered down mode during acquisition of seismic data. Termination of the model generation occurs at exit point 184. The series of operations illustrated by the flow diagram of FIG. 6 is repeated at periodic intervals, typically once per second. In this way, there is placed on the magnetic tape the raw data as well as a periodically updated model of the instantaneous streamer position.

There will next be discussed an alternative method and system for providing a streamer position model using the outputs of the various sensors of the system illustrated in FIG. 5. In the aforementioned paper by John Bedenbender, there are presented six differential equations of equilibrium for a seismic streamer. These equations are reproduced here as equation (7) - (12).

Where:

x, y and z are rectangular coordinates,

T is cable tension,

s is arc length along the streamer,

W is weight per unit length of the streamer in the fluid medium,

.theta. is the angle between the positive x axis and the vertical plane containing element ds,

.PHI. is the angle between ds and the x-z plane in the vertical plane of element ds,

F, h and G are hydrodynamic forces per unit length of streamer acting on ds in respective directions tangent to ds, normal to ds, and normal to vertical plane .PHI., and normal to ds in the vertical plane of .PHI..

In these equations, s the arc length along the streamer is the independent variable, while x, y, z, T, .PHI. and .theta. are the dependent variables. In many cases, including the present streamer positioning system, the initial values of some of the dependent variables and the hydrodynamic forces are not sufficiently known to permit a determinate numerical integration of the differential equations. The paper by Bedenbender presents a solution method which is applicalbe to some such problems. In the method, the unknown initial values are estimated from a general knowledge of the physical situation. Using these initial values, the differential equations are integrated along the streamer to a point at which certain of the dependent variables are known. The known values of dependent variables (boundary conditions) are compared with the values resulting from integration of the differential equations and the differences used in Newton-Raphson iteration technique to provide new estimates for the unknown initial values. The procedure is repeated iteratively using the new initial value estimates until the values produced by the integrated differential equations are within preselected distances of the known boundary conditions.

The aforementioned iterative solution technique has been extended for use with the input variables available in the seismic streamer positioning system. Integration of the differential equations begins at some point along the length of the streamer, for example, at the point nearest vessel 10 at which there is located one of the depth controllers 16 in FIG. 1. Since the depth controller maintains the cable at this point at a predetermined depth, the initial value of dependent variable y is known. Depenent position variables x and z are arbitrarily assigned the value 0. Estimates of dependent variables T, .PHI. and .theta. are made, based on familiarity with conditions generally attending towing of the seismic streamer. Beginning with these initial conditions, the differential equations are numerically integrated along the streamer toward vessel 10 to produce solution values for the dependent variables at the stern of vessel 10. The solution values for dependent variables y, T and .PHI. are compared with known values y, T, and .PHI. at the stern of vessel 10; where y is the known elevation of the streamer at vessel 10; T is the tension measured by the tensiometer located at point 74 in FIG. 4a; and .PHI. is the streamer pitch angle measured with the aid of fixture 30. In the Newton-Raphson iteration technique, this comparison is used to provide new estimates for the initial values of T .theta. and .PHI. at the first depth controller. The iterative procedure is continued until the solution values for y, T and .PHI. converge to within preselected thresholds of the known values.

In the ordinary case, the hydrodynamic forces F, H and G stemming from the velocity of the water relative to the streamer will not be known precisely. The streamer position is strongly influenced by cross-tract water velocity which is known to sometimes vary with depth. In performing the aforementioned iteration on initial conditions at the first depth controller, an initial cross-track velocity gradient is assumed. An outer Newton-Raphson iteration loop is performed to obtain better estimates of the cross-track velocity gradient. The dependent variable which is used in this outer Newton-raphson loop is the yaw angle .theta. at the stern of vessel 10. After the inner and outer loops have converged to provide satisfactory estimates of the three unknown initial conditions, as well as of the cross-track velocity gradient, the final integration of the differential equations using these improved estimates constitutes the model of streamer position for that position of the streamer over which the integration is performed.

Using the new estimates of initial conditions at the first depth controller as known values for the next segment of the streamer, that is, the segment between the first and second depth controllers, the procedure may be repeated to obtain refined estimates of the initial conditions at the location of the second depth controller. In this way, the solution may be stopped down the streamer, so as to obtain a model of the streamer position over its entire length.

Alternatively, as will be discussed hereinbelow, a hyperbolic curve may be fitted to the solution for that portion of the streamer between vessel 10 and the first depth controller.

A FORTRAN listing for implementing the iterative solution is appended hereto as Appendix A. The variables which appear in the FORTRAN program are defined in Table I. FIGS. 8a-i comprise a flowchart corresponding to the FORTRAN program.

With reference to FIG. 8 at step 200, various input quantities are listed. These include the input stream velocity and cable drag coefficients which serve to define the hydrodynamic forces occurring in the differential equations. Various parameters defining the resolution desired in the integration procedure are also established. Finally, estimates of the initial values for tension, yaw angle, and dip angle at the first depth controller, as well as known boundary values at vessel 10 for tension, yaw angle, dip angle and elevation are input. At step 201, a first estimate of cross-track velocity gradient VS1 is made and the variable VSLOPE is set equal to VS1. Step 202 initiates a counter IFLAG for the outer iteration loop and sets the initial Y value at the first depth controller to the predetermined depth maintained by this depth controller. At step 203, the inner loop iteration counter IIT is initiated. At 204, counter JJ is given an initial value of 1. The conditional branch represented by step 205 branches to point 505 when the value of JJ is equal to one. In step 206, the first set of initial values for tension, yaw angle, and dip angle at the first depth controller is assigned.

Using these initial values at step 207, the differential equations are numerically integrated along the streamer from the location of the first depth controller to the stern of the vessel 10. Any of the various methods of numercial integration may be used at this point, one suitable method being the Adams-Moulton prediction corrector method with a Runge-Kutta starter. This particular numercial integration method is implemented by subroutine RKAM for which the FORTRAN code is included in Appendix A.

The integration result in final values for the variables Y, T, and .phi. as shown at step 208. At step 209, the difference between these computed final values and the corresponding boundary conditions is determined. The flow then jumps from step 210 to point 510 of FIG. 8d from which, at step 211, the counter JJ is incremented by 1 and thereafter has the value 2. As a result, at 212, the answer to the test is "NO" and the flow loops back to point B of FIG. 8a. In this case, at step 205, the flow branches to point 506. The sequence of steps 213-217 is almost identical to the sequence of steps 206-210, the only difference being that in step 213, the initial value for tension TEE2 is set to a value slightly different than the value of TEE1 used at step 206. As a result, the error terms formed at step 216 correspond to a slightly different initial condition on tension than do the error terms computed at step 209. The flow than jumps from step 217 to point 510, the variable JJ is incremented again, flow loops back to point B, and ultimately branches to point 507. The initial values set at step 218 differ from those at step 206 in that the initial value for yaw angle is varied slightly, thereby leading to a third set of error terms at step 221. In a similar manner, after the flow loops back to point 508, the initial value on dip angle is varied slightly at step 223, thereby resulting in a fourth set of error terms at step 226. From step 227, flow jumps to point 510 and in this case at step 211, the counter JJ is incremented to have a value of 5. The test at step 212 is again failed and flow loops back to point B from which at step 205 the flow branches to point 509.

Since the flow at this point is passing through the inner boundary loop for the first time, the test at step 228 fails and a fifth set of initial values is established at step 229. It will be seen that the initial values set at step 229 are all slightly different than the values originally set at step 206. This results in a fifth set of error terms at step 232 and flow jumps from step 233 to point 510. Now, counter JJ is incremented to have a value of 6 at step 211 and the answer to test 212 is "YES" so that flow moves to step 234.

The various error terms XLMDA (I, J) resulting from the integrations performed in steps 206-233 are used in the aforementioned Newton-Raphson iterative technique to produce improved initial value estimates TEE3, THETA3, and PHI3, in accordance with matrix equation (13). ##EQU3##

The matrix elements A (I, J) are computed at steps 234 and 235 and are seen to be approximations to the partial derivatives of the error terms with respect to the three initial value quantities TEE, THETA, AND PHI. Thus, matrix equation (13) is implemented by steps 234-240 of the flow chart to yield the new initial value estimates TEE3, THETA3, and PHI3. The new set of initial value estimates is employed in steps 241-244 wherein the differential equations are again integrated to obtain a set of error terms SAVE16, SAVE26, and SAVE36 corresponding to the new set of initial values. These error terms, representing errors in the Y value, the tension value, and the cable dip angle value, are tested in steps 245-247 against preselected upper bounds EPSLN1, EPSLN2, and EPSLN3. If any of these bounds are exceeded by an error term, flow passes to step 248.

At step 248, in preparation for the next pass through the inner loop, the distances of TEE3 from TEE1 and TEE2 are determined. In steps 249 through 251, the quantity TEE2 for the next pass through the inner loop is set to the newly computed quantity TEE3. The quantity TEE1 for the next pass through the loop may or may not be modified in accordance with the result of the test at step 249. Similarly, in steps 252-259, the initial conditions on THETA and PHI for the next pass through the inner loop are chosen closest to the newly computed initial values.

At step 260, the difference between the new quantities TEE2 and TEE1 is tested against a minimum TMIN. If this difference becomes excessively small, it will lead to an unbounded result at step 234. Accordingly, if the answer to the test at step 260 is "YES", the new value TEE2 is adjusted at step 261 to be at least a minimum distance from TEE1.

Similarly, the initial condition sets for PHI and THETA are tested at steps 262 through 265 and adjusted to minimum differences if necessary. At step 266 the inner boundary loop counter IIT is incremented and is tested at step 267 against an iteration loop counter NITER. If the iteration limit is not exceeded at step 267, flow returns to point A for the next pass through the inner boundary loop. On the second and subsequent passes through the inner boundary loop, the only difference from the first pass through this loop stems from the fact that the integration and error evaluation performed in steps 229-233 will already have been performed in the last pass through the loop. Accordingly, at step 228 on these subsequent passes flow branches to point 5095 where these particular error terms are set equal to the values computed in the last pass through the inner boundary loop. The inner boundary loop is exited in one of two ways. If convergence is not ultimately achieved, the iteration limit will be exceeded at step 267 and flow proceeds to step 269 where the program is ended. In the normal mode of operation, prior to this event convergence will be achieved, and the tests at steps 245, 246, and 247 will all be satisfied thereby causing flow to branch to point 600. It will be recalled that the value of VSLOPE used in all of the estimations to this point had been set equal to an initial estimate VS1 in step 201. Also, the counter on the outer loop IFLAG was set equal to -1 at step 202. Accordingly, the test at step 270 will fail, but the test at step 271 will result in a "YES" answer so that flow moves to step 272. Here the quantity VSLOPE is set to a new estimate, VS2, in preparation for the next pass through the outer loop. The error in theta, THTER1, corresponding to the cross-track velocity gradient VS1 is also computed at step 272 and the outer loop counter incremented by one. Flow then returns from step 273 to point 1 to initiate the inner boundary loop for the new value of VSLOPE, that is VS2.

Once again, after convergence is achieved within the inner loop, flow will return to point 600. In this case, the counter IFLAG will have the value 0 and the test at step 270 will result in a "YES" answer. Flow moves to point 3 and at step 274 the error in theta, THTER2, corresponding to the cross-track velocity gradient VS2 is computed. The two velocity gradient estimates, as well as the resulting error terms, are combined at step 275 to yield an improved estimate of cross-track velocity gradient VS3. The quantity VSLOPE is set equal to VS3 and the counter incremented at step 275 and flow moves from step 276 back to point 1 for the next pass through the inner boundary loop.

Once again, after convergence is achieved within the inner boundary loop, flow returns to point 600. In this case, however, the counter IFLAG will have the value 1 and none of tests 270, 271, or 277 will be satisfied. Accordingly, flow moves to point 4. At step 278 the quantities THTER1 and THTER2 are set to the errors corresponding to velocity gradients VS2 and VS3, respectively. If at step 279 the error term THTER2 is found to be less than a predetermined bound value, RPD, the program is terminated at step 280, since convergence will have been achieved.

If convergence is not detected at step 279, then in steps 281-284, the two initial values closest to VS3 from the set VS1, VS2, and VS3 are selected. These two values are used in the Newton- Raphson iteration formula at step 285 to compute an improved initial value estimate VS3. The outer loop counter is incremented at step 386 and flow returns from step 287 back to point 1 for a pass through the inner loop using the newly estimated value of cross-track velocity gradient. The outer loop iteration on cross-track velocity gradient continues either until convergence is detected at step 279 or until five iterations have been completed without achieving convergence, and the program is terminated at step 288. In the normal case where convergence is achieved, the values determined for the initial conditions on streamer tension, yaw angle .theta., and dip angle .phi. and for the cross-track velocity gradient determined in the last pass through the outer loop are the correct values for these quantities. The last integration performed in the last pass through the outer loop will have been performed using these quantities and accordingly, the values for the location quantities x, y, and z determined in the integration define the location model for the streamer.

Once having obtained the solution model for that portion of the streamer between the stern of vessel 10 and the first depth controller, the procedure may be repeated to obtain a model for that portion of the streamer between the first and second depth controllers. The values for the parameters at the first depth controller obtained from the first model are used as boundary values for applying the Newton-Raphson technique to the portion of the streamer between the first and second depth controllers. After a model is developed for that portion of the streamer between the first and second depth controllers, it is used to provide boundary values at the second depth controller for that portion of the streamer between the second and third depth controllers. In this way, the solution proceeds down the streamer until the most remote section of the streamer is modeled.

Another method for generating a model of streamer position involves solving the differential equations as set forth hereinabove only for that portion of the streamer located between vessel 10 and the first depth controller. We have discovered that the position of the projection of the seismic streamer at the water surface at any instant may be well represented by a hyperbola. The x,z coordinates of various sample points along that portion of the streamer between vessel 10 and the first depth controller are determined from the model provided by integrating the differential equations. A hyperbola is then determined which best fits these sample point locations in the least mean square error (LMSE) sense. The hyperbola may then be used as a model of streamer position both for that portion of the streamer between vessel 10 and the first depth controller, as well as for other portions of the streamer.

The general form of the hyperbola to be determined is given in equation (14) where the quantities Z(I) and X(I) are the known coordinates of the Ith sample point along that portion of the streamer between vessel 10 and the first depth controller.

where:

the quantities C(1) and C(2) are parameters of the hyperbola which are to be determined. In addition to these two unknown quantities, it is also true that the origin of the hyperbola which best fits the data samples in the LMSE sense does not occur at the stern of vessel 10, that is at the origin of the x,z coordinate system. As a result, it is necessary to also determine the value of D, the z axis offset of the best fitting hyperbola and CC, the x axis offset of the best fitting hyperbola. The quantities D and CC then, when properly determined, will define an origin vector which extends from the origin of the x,y coordinate system to the origin of the best fitting hyperbola.

To begin the determination of this hyperbola, values are assumed for the quantities D and CC. At this point for simplificity equation (14) may be re-expressed as shown in equation (15).

where:

the values of C(1) and C(2) which define the LMSE hyperbola for the known sample values are then determined in accordance with matrix equation (16). ##EQU4## Once having determined the values of these coefficients, equation (17) may be used to determine the mean square error (MSE) between the hyperbola and the known sampled values. ##EQU5## This hyperbola will be the true LMSE hyperbola, however, only if the origin vector for the hyperbola defined by the assumed values for D and CC is the true origin vector. This of course will rarely be the case. In order to seek a better origin vector, the values of D and CC are varied, and the procedure just described is repeated leading to a new mean square error as defined by equation (17). This iterative process is repeated until values for D and CC are found which lead to a mean square error which is less than a preselected upper bound or until a preselected number of iterations has been performed.

Considering this iterative process in greater detail, after the mean square error has been determined for the hyperbola having the first assumed origin vector, a new random sample for the origin vector is determined. The new origin vector has the same length as the originally assumed origin vector but has a random orientation. The LMSE hyperbola for this new origin vector is then determined, and its mean square error determined in accordance with equation (17). This generation of a new randomly oriented origin vector and subsequent determination of the corresponding mean square error is repeated until an origin vector is found which leads to a mean square error less than that resulting from previous origin vectors. The orientation of this new origin vector is assumed, for the moment, to be the correct orientation.

Now, using this new orientation, a golden section search routine is employed to determine the length of the origin vector which leads to the least mean square error for origin vectors having the new orientation. After determining this new length, the orientation of the vector is again allowed to vary in a random fashion in search of an origin vector with the new length which leads to an even lower mean square error. If some new orientation does indeed provide a lower mean square error, then the golden section search routine is again used to find the LMSE length for this new orientation. The procedure is continued alternating between varying the length and varying the orientation of the origin vector until an origin vector is found which results in a mean square error lower than a preselected upper bound. The hyperbola leading to this final mean square error is the model of the streamer position.

A set of FORTRAN instructions which performs the above-described procedure is included as Appendix B. In addition to the main program, this instruction set contains a function routine called FUNC. This routine determines the C coefficients of the LMSE hyperbola as defined by equation (16) and further determines the mean square error for the hyperbola as expressed in equation (17). The listing also includes subroutine EXPMN2 which inverts a 2 by 2 matrix as is required in equation (16). The instruction set further includes subroutine GSRCH which performs the golden section search routine used to determine the LMSE length for an origin vector with fixed orientation. Finally, the set includes a function routine called YMIN1 which is called by subroutine GSRCH to determine which of a set of four variables has the minimum value.

As set forth above, the procedure provides a hyperbola which accurately represents the location of that portion of the streamer for which integration of the equilibrium differential equations has produced sample locations. Since the location of the more remote portions of the streamer is strongly influenced by the location of that portion of the streamer between vessel 10 and the first depth controller, the hyperbola also serves to accurately represent the location of the remote portions of the streamer.

Whereas there has been disclosed the preferred embodiment of the invention along with several variations, there may be suggested to those skilled in the art other embodiments which do not depart from the spirit and scope of the invention as set forth in the following claims.