JP5919869B2 - Sound source position estimation device - Google Patents

Description

  The present invention relates to a sound source position estimation device, and more particularly to a device that receives a direct wave and a multipath wave of a pulsed sound wave radiated from a sound source in water and estimates the position of the sound source from the arrival time difference therebetween.

  In underwater sound wave propagation, in addition to the sound waves (direct waves) that arrive directly at the receiver from the sound source, there are multipath waves that arrive after being reflected at the sea surface or the sea floor. Arrives at different times. As a technique for estimating the position of a sound source by receiving a pulsed sound wave (pulse wave) from a sound source in water, there is one disclosed in Patent Document 1 below.

  In the conventional sound source position estimation device disclosed in Patent Document 1, the difference in arrival time between a direct wave and a multipath wave is detected based on a signal representing a sound wave (arrival wave) received by an acoustic sensor (receiver). Output as arrival time difference observation value. On the other hand, the arrival time of the sound wave is estimated by performing sound ray calculation of the sound wave from each of the plurality of virtual sound sources based on the underwater sound velocity profile, and the estimation result is stored as an arrival time difference table value (theoretical value). deep. Then, a cost function value reflecting the degree of inconsistency is calculated from the observed arrival time difference value and the table value, and the virtual sound source position at which the cost function is minimized is estimated as the sound source position of the estimation target (target). .

  When the above-described conventional apparatus is used, the value of the cost function on the line connecting the sound source and the receiver becomes a small value over a considerably wide range before and after the position of the sound source, so that the sound source position can be accurately estimated. There was a problem that it was difficult. That is, since the arrival time difference includes an observation error, there is a high possibility that the estimated position is greatly shifted at a place where the value of the cost function is small. Therefore, there is a problem that the observation error is relatively large.

  Therefore, even if there is an error in the arrival time difference observation value, in order to reduce the estimation error of the sound source position, the arrival time difference table is separately provided for the multipath wave whose arrival elevation angle is upward and the multipath wave whose arrival elevation angle is downward. Value is memorized, it is judged whether the incoming elevation angle is upward or downward for each received multipath wave, and the discrepancy between the observed arrival time difference and the table value is determined as a multipath wave having an upward elevation angle. And a multipath wave having a downward elevation angle are separately calculated, and the sound source position is estimated using the sum as a cost function (Patent Document 2).

JP 2006-194627 A JP 2011-149864 A

Finn B.H. Jensen, William A. Kuperman, Michael B. et al. Porterand Henrik Schmidt, “Computational Ocean Acoustic”, Chp. 3, Alp Press

  Non-patent document 1 will be described later.

However, even with the method disclosed in Patent Document 2, the sound source position may not be accurately estimated. There are the following two factors.
One is an error in the correspondence between the observation value of the arrival time difference and the table. This can be reduced by using the method of Patent Document 2, compared with the method of Patent Document 1, but it is not sufficient. The other is an error in the estimation result due to erroneous detection or detection omission of the incoming wave.

  An object of the present invention is to enable accurate estimation of a sound source position even if there is an erroneous detection or omission of an incoming wave.

The sound source estimation apparatus of the present invention is
In a sound source estimation device that estimates the position of a sound source based on an incoming wave received by a receiver group composed of a plurality of receivers ,
Detection value of arrival time difference between arrival waves in a series of arrival wave trains or average value of a plurality of measurement results of the arrival time difference, detection value of arrival elevation angle of each arrival wave, or measurement results of the arrival elevation angle multiple times An arrival wave analysis unit that outputs the standard deviation of the arrival time difference of each arrival wave and the standard deviation of the arrival elevation angle of each arrival wave;
A plurality of pulse wave trains that are in an order opposite to the arrival order of the incoming waves in the series of incoming wave trains, and in which the time difference between successive pulse waves is the same as the arrival time difference of the incoming waves or the average value thereof. Estimate the propagation path of each pulse wave by sound ray calculation when the pulse wave is sequentially transmitted from the receiver group at the same elevation angle as the arrival angle of the corresponding incoming wave or the average value thereof, A back propagation estimator for estimating the likelihood that each pulse wave exists at each position on and around the estimated propagation path at time;
Integrated likelihood calculation for obtaining the sum of all the pulse waves of the likelihood that each pulse wave of the pulse wave train exists at each position at each time after the transmission of the first pulse wave of the pulse wave train And
A convergence time determination unit that obtains a distance between likelihood distributions at the respective times estimated for the plurality of pulse waves in the back propagation estimation unit, and estimates a time at which the distance is minimum as a convergence time; ,
A convergence position determining unit configured to estimate, as a sound source position, a position at which the sum of the likelihoods calculated by the integrated likelihood calculating unit is maximum at the convergence time estimated by the convergence time determining unit. To do.

  According to the present invention, since the sound source position is estimated based on the likelihood distribution taking into account erroneous detection and detection omission, the error is increased even if there is an erroneous detection or detection omission of the propagation path arrival wave. The sound source position can be accurately estimated without any problem.

It is a block diagram which shows the sound source position estimation apparatus of Embodiment 1 of this invention. It is a figure which shows an example of the propagation path from the sound source of sound source position WS. FIG. 3 is a diagram showing an arrangement of receivers in the first embodiment. It is a graph which shows an example of the envelope of the absolute value of the sample value in a receiver. (A) And (b) is a figure which shows an example of the timing of reception in a pair of receiver. It is a block diagram which shows the structural example of the back propagation estimation part of FIG. (A) And (b) is a figure which shows the incoming wave train detected by the back propagation estimation part of FIG. 1, and the pulse wave train which reversed the order produced | generated by the order inversion part of FIG. It is a figure which shows the local coordinate for sound ray calculation. It is a figure which shows the example of the likelihood distribution of the arrival position about two pulse waves. It is a figure which shows the structural example of the convergence time determination part of FIG. It is a figure which shows the example of the sound source position assumed for simulation, and the position of a receiver group. 12 is a graph showing the relationship between the back propagation time and the sum of the inter-distribution distances obtained in the case of the arrangement of FIG. (A)-(c) is a figure which shows the change of the likelihood distribution of an arrival position with respect to the time passage from transmission of a pulse wave. FIG. 14 is an enlarged view showing a portion of the estimated sound source position in FIG.

Embodiment 1 FIG.
FIG. 1 shows a sound source position estimation apparatus according to Embodiment 1 of the present invention. The illustrated sound source position estimation device receives a direct wave and a multipath wave of a pulsed sound wave (pulse wave) radiated from a sound source, and estimates the position of the sound source based on the arrival time difference between them. The estimated position is represented by depth and horizontal distance from the receiver. Note that information indicating the direction in which the sound source is located is obtained separately.

  The illustrated sound source position estimation apparatus includes first and second receivers 201a and 201b, first and second arrival time determination units 202a and 202b, an incoming wave analysis unit 203, a back propagation estimation unit 205, , An integrated likelihood calculation unit 206, a convergence time determination unit 207, a convergence position determination unit 208, and a database 209.

The database 209 stores information indicating the sound velocity profile 209a, the seabed topography 209b, and the receiver position 209c.
As the information indicating the sound velocity profile 209a, it is desirable to use “the sound velocity actually measured in the field” or “the sound velocity profile calculated from the temperature distribution and the salinity concentration by an empirical formula”. When these cannot be obtained, the value of a commercially available database can be used, or the sound speed can be processed as a constant value (for example, a constant value of 1500 m / s which is usually used) is used.

As the information indicating the seabed topography 209b, an actual measurement value or a value of a database that is generally available on the market is used.
As information indicating the receiver position 209c, information indicating the depth at which the receivers 201a and 201b are set and the position in the horizontal direction is input.

The propagation path from the sound source WS is assumed as shown in FIG. In FIG. 2, only five propagation paths WP1 to WP5 are shown for simplicity. The reason why the propagation path (sound ray: flow of sound energy) is not linear is that the speed of sound is not constant due to non-uniformity such as water temperature and salinity.
A wave that reaches the receiving point WR without being reflected once is called a direct wave, and a wave that has reflected at the sea surface or the bottom of the sea once and reached the receiver is called a multipath wave.

  When there is no other sound source between the sound source and the receiver group, the sound wave propagation is expressed by a wave equation in free space. The solution of the wave equation in free space is represented by the Green function. One major feature of the Green function is reciprocity. As a characteristic related to the present invention, since reciprocity is established, the order of the direct wave and the multipath wave that have come out of the sound source and reached the receiver group is reversed (reversed in time) in the reverse direction from the receiver group. When propagating, the sound waves emitted from the receiver group converge at the sound source position.

  This is also true for the sound ray calculation that approximates the wave equation, observing the arrival time and arrival angle at the time of reception, reversing the order of the arrival waves, and reversing at the same angle (the same elevation angle and direction) as the arrival angle Assuming that pulse waves were sent from the receiver group in the direction of the sound source, the arrival position of each transmitted pulse wave at each time point is estimated by sound ray calculation. The arrival position of each pulse wave changes with time. However, all pulse waves should reach (converge) the sound source position at the same time. In the present invention, the sound source position is estimated using this property.

  The first receiver 201a and the second receiver 201b are for receiving an incoming wave, and are provided at different positions on the same vertical line in order to obtain the elevation angle of the incoming wave. A receiver group 201 is constituted by the receivers 201a and 201b.

For example, as shown in FIG. 3, the first and second receivers 201a and 201b are arranged with a gap DWR. The interval DWR is set to a value equal to or greater than the product of the maximum sound speed assumed at the reception position and the minimum time difference measurable by the receivers 201a and 201b.
Since the first and second receivers 201a and 201b are disposed on the upper side and the lower side, respectively, they may be referred to as “upper receiver” and “lower receiver”.

  The outputs of the upper and lower receivers 201a and 201b are input to the first and second arrival time determination units 202a and 202b, respectively.

  The first and second receivers 201a and 201b output data representing sample values obtained by sampling the sound waves received by them, but the absolute values of the sample values obtained by sampling the sound waves The envelopes of those arranged on the time axis are, for example, as shown in FIG. Note that an envelope of the square value of the sample value may be used instead of the absolute value of the sample value.

  The first and second arrival time determination units 202a and 202b determine arrival times (arrival times) of incoming waves based on data strings from the upper and lower receivers 201a and 201b, respectively. For the determination of the arrival time, a method based on replica correlation, a method based on rising edge detection, or the like can be used as in the past. These methods are described, for example, as part of the operation of the arrival time difference measurement unit in Patent Document 1 described above.

5A and 5B, the arrival times of the pulse waves at the upper receiver 201a and the lower receiver 201b determined or detected by the first and second arrival time determination units 202a and 202b. (Detection timing) is shown. 5 (a) is a detection timing in the upper receivers 201a, _TA 1 are listed in the ascending order, TA _2, shows ..., FIG. 5 (b), the detection timing of the lower of the receivers 201b, early in turn _TB 1, _{TB 2,} shown in ....

  Note that any method can be applied to the method of observing the arrival time difference and the arrival angle of the incoming wave as long as the desired information can be observed. In addition to using the information related to the reception timing between the two receivers. In addition, it is possible to configure a receiver array and perform phasing processing, and obtain desired information from the result.

  The outputs of the first and second arrival time determination units 202 a and 202 b are input to the arrival wave analysis unit 203.

  The arrival wave analysis unit 203 receives the outputs of the arrival time determination units 202a and 202b, and calculates an arrival time difference and an arrival elevation angle for each of a plurality of pulse waves (direct wave and multipath wave) constituting a series of pulse waves. To do. The arrival time difference means a time difference from reception (arrival) of the first pulse wave in a series of pulse waves to reception (arrival) of each pulse wave. For the first pulse wave, if the arrival time difference is treated as zero, the subsequent processing can be performed uniformly. In calculating the arrival time difference and the arrival angle, the following processing is performed.

  First, although the reception time of the same pulse wave is different between the upper receiver 201a and the lower receiver 201b, the difference between the reception times, the interval between the receivers 201a and 201b, and the receiver group The elevation angle is calculated based on the speed of sound at the position where. In calculating the arrival time difference, an intermediate value between the reception times of the two receivers 201 a and 201 b is used as the reception time of the receiver group 201.

Based on only a series of pulse waves (based on reception of the first pulse wave and the subsequent pulse wave of a series of pulse waves), the arrival time difference and the arrival angle can be obtained. There may be measurement errors due to missed detection of incoming waves and erroneous detection of incoming waves. Therefore, in the present embodiment, an average value of arrival time differences is taken based on a result of repeating a series of reception of pulse waves a plurality of times, and the average value is used as an “arrival time difference”. Similarly, the average value of the reception results of a plurality of times is used for the arrival angle.
At the same time, based on the measurement results of a plurality of times, the standard deviation _i σ _t of the arrival time difference as an index indicating the degree of variation in the arrival time difference observed for each arrival wave, and the arrival elevation angle observed for each arrival wave A standard deviation _i σ _θ of the arrival angle is obtained as an index representing the degree of variation.

  It should be noted that the position of the sound source or the receiver group may change during the measurement results of a plurality of times, but the change is sufficiently small. In other words, a plurality of measurements are performed in a sufficiently short time with respect to changes in the positions of the sound source and the receiver group.

  Information indicating the arrival time difference (average value) and the arrival elevation angle (average value) for the series of arrival waves obtained by the arrival wave analysis unit 203 is supplied to the back propagation estimation unit 205.

The back-propagation estimation unit 205 calculates the time difference between the incoming and outgoing pulse waves in the order opposite to the arrival order of the incoming waves in the series of incoming wave sequences obtained by the incoming wave analysis unit 203, and the arrival time difference ( When a plurality of pulse waves constituting the pulse wave train having the same average value) are sequentially transmitted from the receiver group 201 at the same transmission elevation angle as the arrival elevation angle (average value) of the corresponding incoming wave, The propagation path of each pulse wave is estimated by sound ray calculation, and the likelihood that each pulse wave exists at each position on and around the estimated propagation path at each time is estimated. Specifically, the estimated position on the propagation path is the center of likelihood, and each pulse wave exists at each position around the center of likelihood based on the standard deviations _i σ _t and _i σ _θ. The likelihood to be obtained is obtained, and the relationship between each position and the likelihood is output as a likelihood distribution. In this case, it is assumed that the likelihood distribution is a normal distribution.

  In estimating the propagation path, information indicating the sound velocity profile 209a, the seafloor topography 209b, and the receiver position 209c stored in the database 209 is also used.

  For example, as illustrated in FIG. 6, the back propagation estimation unit 205 includes an order inversion unit 51, first to Nth propagation calculation units 52-1 to 52-N, and first to Nth likelihood calculation units. 53-1 to 53-N and first to Nth variance-covariance calculation units 54-1 to 54-N.

The order reversing unit 51 is a pulse wave obtained by reversing the order of a series of pulse waves (arrival waves) based on the information about the arrival time difference and the arrival elevation angle input from the arrival wave analysis unit 203 to the back propagation estimation unit 205. Generate information that represents a column.
That is, when the information supplied from the incoming wave analysis unit 203 receives the first, second,... Nth pulse waves in order, the order is reversed, and the Nth, N−1th,. Information representing a pulse wave train rearranged in the order of waves is generated. In this case, the interval between the i-th pulse wave (i is any one of 1 to N-1) and the i + 1-th pulse wave is the same as the interval in the received incoming wave train. As will be described later, information on these pulse waves is used for estimating a propagation path when the pulse waves are virtually transmitted.

Specifically, the arrival time difference (the difference between the reception time of the first arrival wave and the reception time of the i-th arrival wave) for the i-th arrival wave (i-th arrival pulse wave) of the arrival wave train is defined as Fi. If the arrival time difference for the Nth (last) arrival wave is FN, the time difference (time difference from the first pulse wave) Gi of the i-th pulse wave corresponding to the i-th arrival wave is
Gi = FN-Fi (1)
And This time difference Gi is used to determine the timing of virtual transmission for calculating a propagation path, which will be described later, and is called a transmission time difference.
As a result of such processing, the i-th pulse wave in the pulse wave train corresponding to the i-th incoming wave in the incoming wave train is positioned i-th when counted from the back in the pulse wave train. Become.

For example, when the number of arriving waves constituting the arriving wave train is 5, and the arriving waves FW1 to FW5 are received at the timing shown in FIG. 7A, the order is reversed and adjacent to each other. Information indicating a pulse wave train composed of the pulse waves GW5 to GW1 as shown in FIG. 7B is generated with the original time interval between the pulse waves.
From equation (1),
G5 = F5-F5 = 0
G4 = F5-F4
G3 = F5-F3
G2 = F5-F2
G1 = F5-F1 = F5
There is a relationship.

  The order reversing unit 51 generates information on the transmission time difference between the pulse waves constituting the pulse wave train as described above, and also generates information indicating the elevation angle of each pulse wave. In this case, the incoming elevation angle for each incoming wave of each incoming wave train is generated as it is as the outgoing elevation angle for the corresponding pulse wave of the pulse train.

  The output of the order inversion unit 51 (information on the transmission time difference and the transmission elevation angle for each pulse wave of the pulse wave train) is supplied to the first to Nth propagation calculation units 52-1 to 52-N, respectively. That is, information about the i-th pulse wave is supplied to the i-th propagation calculation unit 52-i.

  The first to Nth propagation calculation units 52-1 to 52-N respectively include information on the first to Nth pulse waves generated by the order inversion unit 51 and the sound velocity profile 209a stored in the database 209. Each of the first to Nth pulse waves is transmitted at a virtual transmission timing determined by the transmission time difference, with the information indicating the seafloor topography 209b and the receiver position 209c being input. The propagation path of is calculated.

Each of the first to N-th propagation calculation units 52-1 to 52-N, that is, the i-th propagation calculation unit 52-i (i is any one of 1 to N) is supplied from the order inversion unit 51. The information about the arrival time difference and the arrival angle for the i-th pulse wave, and the i-th pulse wave (pulsed sound wave) at the timing corresponding to the transmission time difference Gi for the pulse wave, Estimate the propagation path (sound ray) of the pulse wave when it is virtually transmitted (launched) in the direction of the sound source at the elevation angle (in the same direction as the incoming wave), and each time (the first pulse wave transmitted) The arrival position of the pulse wave is estimated at the elapsed time after the transmission). Note that the “elapsed time after the transmission of the first transmitted pulse wave” may be simply referred to as the back propagation time.
The “timing corresponding to the transmission time difference Gi” is related to the first to Nth propagation calculation units 52-1 to 52-N. That is, the transmission timing set for the i-th pulse wave in the i-th propagation calculation unit 52-i is based on the transmission timing set for the N-th pulse wave in the N-th propagation calculation unit 52-N. Is also after the time difference Gi.

As a propagation path estimation (sound ray calculation) method, for example, the method described in Non-Patent Document 1 can be used. Specifically, the depth is represented by z, the direction from the receiver group 201 toward the sound source is represented by r, and calculation is performed using the following equations (2a) to (2d).

In the expressions (2a) to (2d), C is the speed of sound.
ζ and ξ are based on the sending elevation angle θ ζ = sin θ / C (0)
ξ = cos θ / C (0)
Is required. Here, C (0) is the speed of sound at the position of the receiver group.
s is a parameter and represents the length of the propagation path along the propagation path.

The time required for each pulse wave to reach each position on the propagation path is calculated by the following equation (3).

  The integral term in equation (3) represents the propagation time along the propagation path. That is, since this integral term is an integral with respect to the parameter s, it is obtained by integrating the reciprocal of the sound velocity along the propagation path. Further, from the equations (2a) to (2d), the path length (the length along the propagation path) from the start point of the sound wave propagation (the receiver group to which the pulse wave is transmitted) can be obtained.

The i-th propagation calculation unit 52-i uses r _ih and z _ih representing the position on the propagation path for the corresponding i-th pulse wave in place of r and z in the above formula, for an appropriate time. The arrival position of each pulse wave (r _ih (t), z _ih (t)) and the received wave at every time, that is, every time a predetermined time elapses after the head pulse wave (Nth pulse wave) is transmitted. The path length s _i (t) from the instrument group 201 is obtained by numerical calculation and is output to the corresponding i-th likelihood calculation unit 53-i and i-th variance-covariance calculation unit 54-i.

On the other hand, from the incoming wave analysis unit 203, information indicating the standard deviation _i σ _t of the arrival time difference and the standard deviation _i σ _θ of the incoming elevation angle for the i-th pulse wave observed by the receiver group 201 is the likelihood calculation. To the unit 53-i and the i-th covariance calculation unit 54-i.

  As described above, in the i-th propagation calculation unit 52-i, each time of the i-th pulse wave is determined based on the arrival time difference (average value) and the arrival elevation angle (average value) of the corresponding i-th arrival wave. The arrival position (the center of likelihood) is calculated by the incoming wave analysis unit 203. The result of analysis by the incoming wave analysis unit 203 may be due to erroneous detection or detection omission of the pulse wave, and is received by the receiver group 201. There are variations in the arrival time difference and the arrival angle of the arrival wave. Taking these into consideration, the arrival position at each time is not only the position on the propagation path estimated by the propagation calculation unit 52-i, but also its surroundings. It may exist in the position. Therefore, the likelihood that each pulse wave exists (has reached) at each time is considered to be distributed around the estimated position on the propagation path as well, and at these positions. A likelihood distribution in which a pulse wave exists is obtained.

Specifically, the likelihood of the arrival position is calculated using the standard deviation _i σ _t of the arrival time difference and the standard deviation _i σ _θ of the arrival angle obtained by the arrival wave analysis unit 203. The likelihood obtained in this way is the arrival position (the likelihood of the likelihood) at each time obtained by sound ray calculation based on the arrival time difference (average value) and arrival elevation angle (average value) of the corresponding i-th arrival wave. (Center). In the following description, it is assumed that the distribution is a normal distribution as described above.

The likelihood that the arrival position of the i-th pulse wave at the time t (elapsed time from the transmission time of the first pulse wave) is the position (r _i (t), z _i (t)) is expressed by the following equation (4). Defined by
_{L i (t) = p i} (r i (t), z i (t) | Path i) p (Path i) (4)

Here, p (Path _i ) is a prior probability that the i-th arrival wave corresponding to the i-th pulse wave is generated (the i-th arrival wave reaches the receiver group).
Normally, it is considered that information such as the directivity of the sound source cannot be obtained from the receiver group side, and therefore the prior probabilities for each incoming wave should be equal based on the principle of Laplace's insufficient reason.

In addition, the conditional probability should be configured in accordance with the information that makes the probability density function clear.
When there is no information for clarifying the probability density function, the arrival time difference (average value) and the arrival elevation angle (with s as the ray path direction and η as the direction perpendicular to the ray path as shown in FIG. Local coordinates (u _i (t), v _i ) having the origin (r _ih (t), z _ih (t)) obtained by sound ray calculation using the transmission time difference and the elevation angle corresponding to the average value) as the origin The likelihood that the position represented by (t)) (position (position represented by r _i (t), z _i (t)) is the arrival position of the i-th transmitted wave is expressed by the following equation (5). The likelihood calculation unit 53-i calculates the likelihood L _i (t) by performing an operation according to the following equation.

In equation (5),
s _i (t) is the length of the propagation path to the above position (r _ih (t), z _ih (t)).
_i σ _t is the arrival time difference standard deviation for the _i- th pulse wave, and _i σ _θ is the arrival elevation standard deviation (expressed in radians) for the i-th pulse wave.
Further, in the derivation of the equation (5), since the relative size, not the absolute value of the likelihood, is a problem in the subsequent processing, the processing is performed with p (Path _i ) = 1 in the equation (4). doing.

The likelihood distribution represented by the equation (5) is conceptually shown in FIG.
In FIG. 9, two pulse waves (the first arrival wave and the transmission pulse wave corresponding to the second arrival wave) reach each position at time t (elapsed time from the transmission of the first pulse wave) t. Likelihood to do is shown by the contour line. That is, the likelihood p ₁ (r ₁ (t), z ₁ (t) | Path ₁ ) that the arrival position at the time t of the transmitted pulse wave corresponding to the first incoming wave is each position (r, z). And the likelihood p ₂ (r ₂ (t), z ₂ (t) | Path ₂ ) that the arrival position of the transmitted pulse wave corresponding to the second incoming wave is each position (r, z) at time t. Is shown.

  As described above, a position shifted by a value corresponding to the standard deviation with respect to a position (center of likelihood) obtained based on the average value of the arrival time difference and the average value of the arrival angle is indicated by a chain line. As indicated by the chain line, the position corresponding to the standard deviation is elliptical.

The above processing is performed in each of the first to Nth likelihood calculation units 53-1 to 53-N, and the likelihoods L _i (t (t) obtained by the respective likelihood calculation units 53-1 to 53-N ) Is output to the integrated likelihood calculation unit 206 and the convergence time determination unit 207.

In the variance-covariance calculation unit 54-i, the variance and covariance are obtained by the following equations (6A) to (6C).

The dispersion and covariance for the i-th pulse wave represented by the formulas (6A) to (6C) are collectively represented by the symbol VC _i .
The above processing is performed in each of the first to Nth variance-covariance calculation units 54-1 to 54-N, and the variance and covariance obtained by the respective variance-covariance calculation units 54-1 to 54-N VC _i is output to the convergence time calculation unit 207.

The integrated likelihood calculation unit 206 inputs the _{first to} Nth likelihoods L ₁ (t) to L _N (t) output from the first to Nth likelihood calculation units 53-1 to 53-N. Then, the likelihoods are integrated, and the integration result is output to the convergence position determination unit 208 as the integrated likelihood Ls (t).
The integrated likelihood calculating unit 206 performs the above integration by adding the _{first to} Nth likelihoods L ₁ (t) to L _N (t) at each time t according to the following equation (7).

式（７）の計算は、各時刻ｔについてそれぞれ行われ、その結果尤度Ｌｓ（ｔ）の時系列情報
Ｌｓ（ｔ１）、Ｌｓ（ｔ２）、Ｌｓ（ｔ３）、…
が得られる。

The calculation of Expression (7) is performed for each time t, and as a result, time series information Ls (t1), Ls (t2), Ls (t3),...
Is obtained.

The convergence time determination unit 207 obtains the distance between the distributions of the likelihoods L ₁ (t) to L _N (t) estimated by the likelihood calculation units 53-1 to 53-N of the back propagation estimation unit 205, The time at which the distance is minimum is estimated.
Specifically, the convergence time determination unit 207 obtains the distance between the likelihood distributions for the two pulse waves constituting each of all combinations obtained by extracting two from the plurality of pulse waves, and The sum for all combinations is obtained and output as the distance between the likelihood distributions for a plurality of pulse waves.

For example, as shown in FIG. 10, the convergence time determination unit 207 includes a plurality of inter-distribution distance calculation units 71-(1, 2) to 71 (N−1, N) and a distance addition unit 72. Between the total number U of the inter-distribution distance calculation units 71- (1,2) to 71 (N-1, N) and the number N of pulse waves,
U = (N−1) · N / 2
There is a relationship.

Each inter-distribution distance calculation unit 71- (i, j) (i is any one from 1 to (N-1), j is any one from (i + 1) to N, and therefore j ≠ i) Likelihood L _i (t) output from the likelihood calculator 53-i, likelihood Lj (t) output from the j-th likelihood calculator 53-j, and i-th covariance calculation and variances and covariances VC _i is output from the section 54-i, receives the variance and covariance VC _j outputted from the variance-covariance calculation unit 54-j of the j, the likelihood distribution of delivery wave of the i And the inter-distribution distance d (i, j) of the likelihood distribution of the j-th transmitted wave.

ここで、

は第ｉの送出波の尤度の中心位置、

は第ｊの送出波の尤度の中心位置、
Σ_ｉは第ｉの送出波の尤度の分散共分散行列であり、

で表され、
Σ_ｊは第ｊの送出波の尤度の分散共分散行列であり、

で表される。
また、^ｔ（μ_ｉ−μ_ｊ）は（μ_ｉ−μ_ｊ）の転置ベクトルである。 The inter-distribution distance d (i, j) is obtained by the equation (8).

here,

Is the central position of the likelihood of the i-th transmitted wave,

Is the central position of the likelihood of the jth transmitted wave,
Σ _i is the variance-covariance matrix of the likelihood of the i-th transmitted wave,

Represented by
Σ _j is the variance-covariance matrix of the likelihood of the j-th transmitted wave,

It is represented by
Also, ^t (μ _i −μ _j ) is a transposed vector of (μ _i −μ _j ).

  The value obtained by equation (8) is generally called the Mahalanobis distance. FIG. 9 conceptually illustrates this Mahalanobis distance.

In FIG. 9, the Mahalanobis distance between _two likelihood distributions L ₁ (t) and L ₂ (t) at a specific time t is indicated by a symbol d (1, 2).

  Each of the inter-distribution distance calculation units 71- (i, j) (i = 1 to (N−1), j = (i + 1) to N) is a distance between the input likelihood distributions (transmitted wave distribution). Distance) d (i, j) is calculated, and the calculation result is output to the distance adding unit 72.

The distance adding unit 72 calculates the inter-distribution distances d (1, 2) to d (N−1, N) obtained by the inter-distribution distance calculating units 71 (1, 2) to 71- (N−1, N). As an input, the convergence time is calculated and the calculation result is output. That is, the distance adding unit 72 receives the inter-distribution distances d (1, 2) to d (N−1, N) as input, and obtains the sum of the inter-distribution distances using the equation (11).

The total distance Ds (t) between the distributions is a function of the reverse propagation time (elapsed time from the transmission of the first transmission wave) t. Using the sum Ds (t) of the inter-distribution distance, the convergence time t _c is obtained by the equation (12).
t _c = t such that Ds (t) → minimum (12)
That is, the time that Ds (t) is minimum the convergence time t _c.
Convergence time determination unit 207 outputs information representing the convergence time t _c determined.

Convergence position determination unit 208 receives the likelihood function Ls (t) from the integrated likelihood calculating unit 206 receives the convergence time _{t c} from the convergence time determination unit 207, in the likelihood function Ls (t) at time _{t c} (In the distribution of likelihood Ls (t) for different values of (r, z)), the value of (r, z) that maximizes the value of likelihood Ls (t) is obtained and obtained (r , Z) is output as a position representing the convergence position.

As shown in FIG. 11, the receiver group WR is at a distance r = 0 m and a depth z = 70 m, the sound source WS is at a distance r = 150 m and a depth z = 10 m, the sea floor is 75 m flat and the sound velocity distribution is constant. FIG. 12 shows the relationship between the back propagation time t and the sum Ds (t) of the inter-distribution distances D _i (t) when the simulation is performed at 1500 m / s.
In the case of the positional relationship of FIG. 11, it is correct to determine that the sound ray path has converged at 113 msec from the start of propagation of the transmitted wave.

  As shown in FIG. 12, the sum of the inter-distribution distances is the smallest at the time when the sound rays (propagation paths) of all the outgoing waves are closest to each other (most converged). It can be seen that the most converged time can be correctly determined based on the sum.

As described above, the convergence position determination unit 208 receives the integrated likelihood Ls (t) output from the integrated likelihood calculation unit 206 and the convergence time t _c output from the convergence time determination unit 207 as input, and determines the convergence position. Detect and output information indicating the convergence position.

  13A to 13C, the range corresponding to the standard deviation of the arrival time difference and the standard deviation of the arrival angle is represented by an ellipse, and shows how the range changes with the back propagation time. FIG. 14 shows an enlarged portion of the estimated sound source position in FIG. FIGS. 13A to 13C and FIG. 14 also show simulation results when it is assumed that the positional relationship between the sound source and the receiver group is as shown in FIG.

In the above description, it is assumed that an incoming wave sequence including a number of incoming waves equal to the number of propagation calculation units 52-1 to 52-N provided in the back propagation estimation unit 205 is obtained. When the number n of arriving waves constituting the arriving wave train obtained by the wave analyzing unit 203 is smaller than the number N of the propagation calculating units 52-1 to 52-N, the propagation calculating units 52-1 to 52- The propagation path is calculated using only n of N, and the likelihood is calculated using only n of the likelihood calculating units 53-1 to 53-N. Similarly, the convergence time determination unit 207, based on the n likelihoods L ₁ (t) to L _n (t), the inter-distribution distance calculation units 71- (i, j) to 71- (N−1, The distribution distances d (1,2) to d (n-1) using only the distribution distance calculation units 71- (i, j) to 71- (n-1, n) which are a part of N). , N). Further, the order inversion unit 51 inverts the first to n-th incoming waves to generate first to n-th pulse waves. In this case, the transmission timing is obtained by the following equation (13) instead of the above equation (1).
Gi = Fn-Fi (13)
Here, Fn is the arrival time difference of the nth arrival wave.

  When the method according to the above embodiment is used, the arrival time difference-observed value table used in the conventional apparatus is not used. Therefore, it is not necessary to associate the observed arrival wave with the arrival time difference-observed value table, and a ranging error due to an association error that has occurred in the past cannot occur.

  In the above embodiment, a receiver group consisting of receivers arranged above and below is used to receive sound waves from a sound source. However, a receiving array may be used instead. Any configuration that observes the elevation angle of arrival is acceptable.

  Moreover, it is good also as providing not only one receiver group but several receiver groups. For example, if a plurality of receiver groups are arranged at different depths and information is collected at each position, the sound source estimation accuracy can be improved.

  In addition, error variance is obtained by actual measurement. However, when the number of samples in the same situation is small and the number of samples sufficient to calculate the standard deviation cannot be secured, the SN ratio at the time of reception is measured, and Cramer-Rao's The value obtained from the inequality can be used as the standard deviation. In this case, the arrival wave analyzing unit 203 is configured to obtain the standard deviation as described above.

  51 order reversal unit, 52-1 to 52-N propagation calculation unit, 53-1 to 53-N likelihood calculation unit, 71- (1,2) to 71- (n-1, n) inter-distribution distance calculation unit 72 distance adder, 201a, 201b receiver, 202a, 202b arrival time determination unit, 203 arrival wave analysis unit, 205 back propagation estimation unit, 206 integrated likelihood calculation unit, 207 convergence time determination unit, 208 convergence position determination Department.

Claims (7)

In a sound source estimation device that estimates the position of a sound source based on an incoming wave received by a receiver group composed of a plurality of receivers ,
Detection value of arrival time difference between arrival waves in a series of arrival wave trains or average value of a plurality of measurement results of the arrival time difference, detection value of arrival elevation angle of each arrival wave, or measurement results of the arrival elevation angle multiple times An arrival wave analysis unit that outputs the standard deviation of the arrival time difference of each arrival wave and the standard deviation of the arrival elevation angle of each arrival wave;
A plurality of pulse wave trains that are in an order opposite to the arrival order of the incoming waves in the series of incoming wave trains, and in which the time difference between successive pulse waves is the same as the arrival time difference of the incoming waves or the average value thereof. Estimate the propagation path of each pulse wave by sound ray calculation when the pulse wave is sequentially transmitted from the receiver group at the same elevation angle as the arrival angle of the corresponding incoming wave or the average value thereof, A back propagation estimator for estimating the likelihood that each pulse wave exists at each position on and around the estimated propagation path at time;
Integrated likelihood calculation for obtaining the sum of all the pulse waves of the likelihood that each pulse wave of the pulse wave train exists at each position at each time after the transmission of the first pulse wave of the pulse wave train And
A convergence time determination unit that obtains a distance between likelihood distributions at the respective times estimated for the plurality of pulse waves in the back propagation estimation unit, and estimates a time at which the distance is minimum as a convergence time; ,
A sound source position estimation apparatus comprising: a convergence position determination unit configured to estimate, as a sound source position, a position at which the total of the likelihoods calculated by the integrated likelihood calculation unit is maximum at the convergence time estimated by the convergence time determination unit. . The back propagation estimator is
Using the estimated position on the propagation path as the center of likelihood, based on the standard deviation, the likelihood that each pulse wave exists at each of the center of the likelihood and the surrounding positions is obtained, The sound source position estimation apparatus according to claim 1, wherein a relationship with likelihood is output as a likelihood distribution. The convergence time determination unit
Obtaining a distance between likelihood distributions of two pulse waves constituting each of all combinations obtained by extracting two from the plurality of pulse waves, obtaining a sum of all the combinations of the distances, The sound source position estimation apparatus according to claim 1, wherein the sound source position estimation apparatus outputs a distance between likelihood distributions of a plurality of pulse waves. The back propagation estimator is
Pulse waves in which the time difference between successive pulse waves is the same as the arrival time difference of the arrival waves or the average value thereof in the order opposite to the arrival order of the first to n-th arrival waves in the series of arrival wave trains An order reversing unit that generates information indicating the timing of virtual transmission of the first to n-th pulse waves constituting the sequence and information indicating the arrival angle of the pulse waves;
When the information about the 1st to n-th pulse waves generated by the order reversing unit is input, and the 1st to n-th pulse waves are respectively transmitted at the virtual transmission timing, respectively. First to n-th propagation calculation units for calculating propagation paths of pulse waves of
First to n-th likelihood calculators for calculating likelihoods at respective positions based on calculation results of propagation paths for the first to n-th pulse waves by the first to n-th propagation calculators; Have
The i-th propagation calculation unit (i is any one of 1 to n) transmits the i-th pulse wave at the virtual transmission timing corresponding to the arrival time difference of the corresponding i-th arrival wave. Estimating the propagation path of the i-th pulse wave when transmitted from the receiver group at the same elevation angle as the arrival elevation angle or the average value of the arrival waves by sound ray calculation,
The i-th likelihood calculation unit uses the position on the propagation path of the i-th pulse wave estimated by the i-th propagation calculation unit as the center of likelihood, and the standard deviation of the arrival time difference of the i-th arrival wave And the likelihood that the i-th pulse wave exists at each of the center of the likelihood of the i-th pulse wave and its surroundings based on the standard deviation of the angle of elevation of the i-th incoming wave. 2. The sound source position estimating apparatus according to claim 1, wherein the sound source position estimating apparatus outputs the relationship between each position and likelihood as a likelihood distribution for the i-th pulse wave. The convergence time determination unit
A plurality of inter-distribution distance calculation units for obtaining distances between likelihood distributions at each time for all combinations obtained by extracting two from the first to n-th pulse waves;
A distance addition unit for obtaining a total of the inter-distribution distances obtained by the plurality of inter-distribution distance calculation units,
Each inter-distribution distance calculation unit calculates the likelihood distribution for the i-th pulse wave constituting the corresponding combination and the likelihood for the j-th pulse wave (j is 1 to n, where j ≠ i). The sound source position estimation apparatus according to claim 4, wherein a distance of degree distribution is obtained. The sound source position estimation apparatus according to claim 2 or 4, wherein the likelihood of each position is obtained assuming that the likelihood distribution is a normal distribution. The incoming wave analysis unit
Obtain the standard deviation of the arrival time difference of each arrival wave based on the measurement results of the arrival time difference between arrival waves in a series of arrival wave sequences,
The sound source according to any one of claims 1 to 6, wherein a standard deviation of the arrival angle of each incoming wave is obtained based on a plurality of measurement results of the elevation angle of each incoming wave in a series of incoming wave trains. Position estimation device.

Images