JP2012215490A - Sound source position estimation device - Google Patents

Abstract

PROBLEM TO BE SOLVED: To make an estimated error of a sound source position small even when there is an error in an observation value of arrival time difference.SOLUTION: Theoretical values of arrival time difference are separately stored for a sea level one time reflected wave and a sea bottom one time reflected wave (209M), degrees of noncoincidence between the arrival time difference and the theoretical values to a direct wave of the received sea level one time reflected wave and the received sea bottom one time reflected wave are separately calculated, the sum of them is used as a cost function (210), and a sound source position is estimated (211). Area filtering is performed based on context of arrival time observation values in first to third echo sounders (201, 202, 301) of each of the direct wave, the sea level one time reflected wave and the sea bottom one time reflected wave.

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 receiving a pulsed sound wave from a sound source in water and estimating the position of the sound source, there is one disclosed in Patent Document 1 below.

  FIG. 1 shows a conventional sound source position estimation apparatus disclosed in Patent Document 1. In the illustrated sound source position estimation apparatus, a signal representing a sound wave received by an acoustic sensor (receiver) is input to a multipath wave arrival time difference measuring device 102 via an input terminal 101. Multipath wave arrival time difference measuring device 102 detects the arrival time difference between the direct wave and the multipath wave, and outputs data indicating the detected arrival time difference observation value.

  On the other hand, the arrival time difference table creation unit 103 estimates the arrival time of the sound wave by calculating the sound ray of the sound wave from each virtual sound source based on the underwater sound velocity profile, and the estimated result is the theoretical value of arrival time difference (table value). As you remember.

  The cost function calculation unit 104 uses the observed multi-path wave arrival time difference obtained by the multi-path wave arrival time difference measuring unit 102 and the arrival time difference theoretical value obtained by the arrival time difference table creation unit 103 to reflect the degree of inconsistency between them. Calculate the value of the function.

A detailed block configuration of the cost function calculation unit 104 is shown in FIG.
The observed multipath wave arrival time difference obtained by the multipath wave arrival time difference measuring device 102 is input to the input terminal 111, and the theoretical value from the arrival time difference table creation unit 103 is input to the input terminal 112. The theoretical value-observed value difference calculation unit 113 calculates the difference between the observed arrival time difference value and the theoretical arrival time difference value.

Here, p is an index corresponding to the horizontal distance from the receiving position (receiving point) to the sound source position, and q is an index corresponding to the depth of the sound source position, and the nth obtained by the multipath wave arrival time difference measuring device 102. D _n is the arrival time difference for the arriving waves, and the n th arrival time difference theoretical value from the virtual sound source at the virtual sound source position (p, q) obtained by the arrival time difference table creation unit 103 is E _n (p, q). In this case, the theoretical value-observed value difference calculation unit 113 obtains the difference between the theoretical value and the observed value by the following equation (1).
ΔD _n (p, q) = E _n (p, q) −D _n (1)

  The sum-of-squares calculation unit 114 calculates the sum of all multipath waves (represented by the following expression (2)) as a cost function, although the result of the expression (1) is squared. Is output from the output terminal 115.

  The sound source position estimation unit 105 in FIG. 1 obtains coordinates (p, q) that minimize the cost function K (p, q) given by Equation (2), and the obtained coordinates (p, q) are the estimation targets ( The target sound source position is estimated, and the estimation result is output via the output terminal 106.

JP 2006-194627 A

  However, when the above-mentioned 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 is accurately estimated. There was a problem that I could not. That is, since the arrival time difference obtained by the multipath wave arrival time difference measuring device 102 includes an observation error, there is a high possibility that the estimated value is greatly shifted at a place where the value of the cost function is small, and thus the observation error is relatively large. was there.

As an example, consider the case where the sound speed is constant and the water depth Zb is constant as shown in FIG. In FIG. 3, Zs is the depth of the sound source, Zr is the depth of the receiving point, and r is the horizontal distance from the receiving point WR to the sound source WS.
FIG. 4 shows the result of calculating the cost function when there is no error in the arrival time difference observation values. The horizontal axis represents the horizontal distance r from the receiving point WR, and the vertical axis represents the depth (depth is greater in the downward direction) z. The curve in the figure is the contour of the cost function, and the position where the cost function is minimum is estimated as the sound source position.

  If there is no error in the arrival time difference observation value, the true value can be estimated. However, the cost function in the figure shows the cost along the direction of the line connecting the receiver and the sound source (direction of arrow DMC). Since the change in the function is small, if the observed value of the arrival time difference includes an error, the estimation result may deviate from the true sound source position along this direction. When empirically known observation errors were given and the simulation was tried 30 times under the conditions shown in FIG. 3, the standard error of the distance error was about 15% of the true distance value, which was a fairly large value. . Thus, there is a problem that the position estimation error of the sound source is large due to the characteristics of the cost function.

  Therefore, the present applicant stores the theoretical value of the arrival time difference separately for the multipath wave with the incoming elevation angle upward and the multipath wave with the incoming elevation angle downward, and arrives for each received multipath wave. Judgment is made on whether the elevation angle is upward or downward, and the discrepancy between the observed arrival time difference and the theoretical value is calculated separately for each of the multipath wave having an upward elevation angle and the multipath wave having a downward elevation angle. Thus, it has been proposed to estimate the sound source position using the sum as a cost function (Japanese Patent Application No. 2010-012186).

  However, even with this method, when the number of multipath waves that can be acquired is small, particularly when only a single reflected wave at the sea surface and a single reflected wave at the bottom of the sea are used, there is a problem that the estimation error cannot be sufficiently reduced. .

The sound source estimation apparatus of the present invention is
Reception of incoming waves at the first and second receivers arranged in a line on a common vertical line and a third receiver arranged at a position shifted from the vertical line at the reception position. In the sound source estimation device that estimates the position of the sound source based on the wave time,
For each of the plurality of virtual sound source positions, the theoretical value of the arrival time difference between the direct wave of the incoming waves at the receiving position and the once reflected wave at the sea surface, assuming that a virtual sound source exists at the virtual sound source position; And an arrival time difference theoretical value storage means for storing a theoretical value of an arrival time difference between the direct wave and the reflected light at the seabed once at the receiving position;
Based on the arrival wave data observed by the first and second receivers, an observed value of the arrival time difference between the direct wave from the target sound source, the sea surface reflected wave and the seabed reflected wave is calculated. Arrival time difference observation value calculation means,
Obtaining the inconsistency between the arrival time difference between the sea surface reflected wave and the direct wave and the theoretical value of the arrival time difference for the sea surface reflected wave as the first inconsistency;
The degree of inconsistency between the arrival time difference between the seabed once reflected wave and the direct wave and the theoretical value of the arrival time difference for the seabed once reflected wave is determined as a second inconsistency, and the first inconsistency and the A cost function calculating means for obtaining a cost function based on the sum of the second mismatch degree;
Sound source position estimating means for estimating the virtual sound source position at which the cost calculated by the cost function calculating means is the minimum as the position of the target sound source,
The cost function calculating means includes
Based on the anteroposterior relationship of the arrival time observation values in the first to third receivers of each of the direct wave, the sea surface once reflected wave and the seabed once reflected wave among the incoming waves, The sound source position is estimated based on the cost function for a position other than the position determined to be a sound source position among the virtual sound source positions.

  According to the present invention, it is possible to reduce the estimation error of the sound source position even when there is an error in the arrival time difference observation value and the number of multipath waves that can be acquired is small.

It is a block diagram which shows the conventional sound source position estimation apparatus. It is a block diagram which shows the detail of the cost function calculation part 104 of FIG. It is a figure which shows the propagation path of a sound wave in case sound velocity is constant and water depth is constant. It is a figure which shows an example of the contour line of the cost function calculated | required with the apparatus of FIG. It is a block diagram which shows the sound source position estimation apparatus by Embodiment 1 of this invention. 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)-(c) is a figure which shows an example of the timing of reception in the 1st thru | or 3rd receiver. It is a figure which shows an example of a virtual sound source position. It is a figure which shows an example of the propagation path from the sound source of virtual sound source position WS. FIG. 6 is a block diagram illustrating an example of an arrival time difference table creation unit 209 in FIG. 5. FIG. 6 is a block diagram illustrating an example of a cost function calculation unit 210 in FIG. 5. It is a block diagram which shows the structural example of the area filter production | generation part 305 of FIG. It is a block diagram which shows the structural example of the direct wave area | region filter production | generation part 311 of FIG. It is a figure which shows the characteristic of the area filter produced in the 1st direct wave area filter production part 331 of FIG. It is a figure which shows the characteristic of the area filter produced by the 2nd and 3rd direct wave area filter production parts 332 and 333 of FIG. It is a figure which shows the characteristic of the area filter produced in the sea surface reflected wave area filter production part 312 of FIG. It is a figure which shows the characteristic of the area filter produced in the seabed reflected wave area filter production part 313 of FIG. It is a figure which shows the characteristic of the area filter produced in the area filter production part 305 of FIG. It is a figure which shows an example of the contour line of the cost function calculated | required with the sound source position estimation apparatus of FIG. It is a figure which shows an example of the contour line of the cost function obtained when not applying area | region filtering with the apparatus of FIG. It is a figure which shows an example of the contour line of the cost function calculated | required with the apparatus of FIG. It is a block diagram which shows the modification of the cost function calculation part 210 of FIG. It is a block diagram which shows an example of the cost function calculation 210 used with the sound source position estimation apparatus of Embodiment 2 of this invention. It is a block diagram which shows an example of the time difference determination part 825 of FIG. (A)-(c) is a figure which shows an example of the timing of reception in the 1st thru | or 3rd receiver 201,202,301.

Embodiment 1 FIG.
FIG. 5 shows a sound source position estimation apparatus according to the first embodiment of the present invention.
The illustrated sound source position estimation device receives a direct wave and a multipath wave of a pulsed sound wave radiated from a sound source, and estimates the position of the sound source based on the arrival time difference between them. And the horizontal distance from the receiver. It should be noted that information indicating the direction in which the sound source is located (direction when viewed from above) is separately obtained.

The sound source position estimation apparatus according to Embodiment 1 includes first, second, and third receivers 201, 202, and 301.
As will be described later, the first receiver 201 and the second receiver 202 are also used to obtain information on whether the elevation angle of the incoming wave is upward or downward, as shown in FIG. Are provided at different positions on the same vertical line VL. With this arrangement, the elevation angle of the incoming wave is upward based on whether the incoming wave (signal) has arrived at the first receiver 201 or the second receiver 202 earlier. Or down.

Since the first, second, and third receivers 201, 202, and 301 are used to determine which receiver the same pulsed sound wave has reached first, as shown in FIG. DWRa, DWRb, and DWRc are opened. The interval DWRa 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 first and second receivers 201 and 202. The interval DWRb 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 first and third receivers 201 and 301. The interval DWRc 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 second and third receivers 202 and 301.
Since the first and second receivers 201 and 202 are disposed on the upper side and the lower side, respectively, they may be referred to as “upper receiver” and “lower receiver”.

  The third receiver 301 is disposed below the first receiver 201 and above the second receiver 202, and the first receiver 201 and the second receiver The position is shifted from the vertical line VL connecting 202, for example, the position shifted in the direction DWS of the sound source. As described above, since the direction of the sound source is separately detected, the third receiver 301 is rotated by rotating the third receiver 301 around the vertical line VL so as to be positioned in the detected direction. The device 301 can be located in the direction DWS of the sound source. In addition, in order to detect the direction of the sound source, when using a large number of receivers arranged along a ring or a cylindrical surface with a vertical line as the central axis or in a large number of directions, the above-described many received waves are used. Of these, a receiver located in the direction of the detected sound source may be used as the third receiver 301.

Note that the third receiver 301 does not have to be accurately displaced in the direction of the sound source, and the displacement only needs to have a positive or negative component in the direction of the sound source. In other words, it may be displaced in any other direction except the direction perpendicular to the direction of the sound source.
Hereinafter, in order to simplify the description, it is assumed that there is a deviation in the direction of the sound source. In the example shown in FIG. 6, the third receiver 301 passes through the midpoint of the first receiver 201 and the second receiver and is in a plane BL1 that is horizontal to the vertical line VL. It shall be.

  Furthermore, a straight line VL connecting the first receiver 201 and the second receiver 202, a straight line Lac connecting the first receiver 201 and the third receiver 301, and a second received wave. The condition that the straight lines Lbc connecting the device 202 and the third receiver 301 are not parallel to each other is satisfied. The significance of this will be apparent from the description made later with reference to FIGS.

Outputs of the first, second, and third receivers 201, 202, and 301 are input to first, second, and third arrival time determination units 205, 206, and 303 via input terminals 203, 204, and 302, respectively. Is input.
First, second, and third arrival time determination units 205, 206, and 303 determine or detect arrival times based on the outputs of the first, second, and third receivers 201, 202, and 301, respectively. .

  The first, second and third receivers 201, 202 and 301 output data representing sample values obtained by sampling the sound waves received by them, but are obtained by sampling the sound waves. The envelope of the absolute values of the sample values arranged on the time axis is, 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, second, and third arrival time determination units 205, 206, and 303 receive arrival waves based on data strings from the first, second, and third receivers 201, 202, and 301, respectively. Judge (arrival time). 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 multipath wave arrival time difference measuring device in Patent Document 1.

8A, 8B, and 8C, the first, second, and third received waves determined or detected by the first, second, and third arrival time determining units 205, 206, and 303 are shown. The arrival times (detection timings) of pulsed waves in the devices 201, 202 and 301 are shown. 8A shows the detection timing in the first receiver 201 in the order of TA ₁ , TA ₂ , TA ₃ ,... FIG. 8B shows the detection in the second receiver 202. The timings are indicated by TB ₁ , TB ₂ , TB ₃ ,... In order from the earliest. FIG. 8C shows the detection timings in the third receiver 301 by TC ₁ , TC ₂ , TC ₃ ,. Show.
If noise is erroneously detected as a signal or no detection omission occurs, the first detection timing is a direct wave detection timing, and the second and subsequent detection timings are multipath wave detection timings.

Information TA, TB, and TC indicating detection timings detected by the first, second, and third arrival time determination units 205, 206, and 303 are supplied to the cost function calculation unit 210.
Information TA and TB indicating the detection timings detected by the first and second arrival time determination units 205 and 206 are also supplied to the first and second multipath wave arrival time difference calculation units 207 and 208.

The first multipath wave arrival time difference calculation unit 207, based on the output of the first arrival time determination unit 205, arrives at the arrival time difference (first time) for each of the multipath waves received by the first receiver 201. The arrival time difference from the arrival wave is calculated. That is, the arrival time difference DU ₂ between the _first arrival wave PA ₁ and the second arrival wave PA ₂ (first multipath wave), the arrival time difference DU between the first arrival wave PA ₁ and the third arrival wave PA _{3. 3. In} general, the arrival time difference DU _n between the _first arrival wave PA ₁ and the nth arrival wave PA _n (n = 2, 3,...) Is calculated. The data string GDU (= DU ₂ , DU ₃ ,...) Indicating the arrival time difference DU _n calculated in this way is supplied to the cost function calculation unit 210.

Similarly, the second multipath wave arrival time difference calculation unit 208, based on the output of the second arrival time determination unit 206, receives the arrival time difference for each of the multipath waves received by the second receiver 202 ( The arrival time difference from the first arrival wave is calculated. That is, the arrival time difference DL ₂ between the _first arrival wave PB ₁ and the second arrival wave PB ₂ (first multipath wave), the arrival time difference DL between the first arrival wave PB ₁ and the third arrival wave PB _{3. 3. In} general, the arrival time difference DL _n between the _first arrival wave PB ₁ and the nth arrival wave PB _n (n = 2, 3,...) Is calculated. The data string GDL (= DL ₂ , DL ₃ ,...) Indicating the arrival time difference DL _n calculated in this way is supplied to the cost function calculation unit 210.

  The combination of the first arrival time determination unit 205 and the first multipath wave arrival time difference calculation unit 207 has the same function as the multipath wave arrival time difference measurement unit 102 in FIG. The combination of 206 and the second multipath wave arrival time difference calculation unit 208 also has the same function as that of the multipath wave arrival time difference measuring device 102 of FIG.

  As described above, the first multipath wave arrival time difference calculation unit 207 and the second multipath wave arrival time difference calculation unit 208 determine the arrival time (detection timing) TA and TB of the arrival wave observed at the receiving position. Based on this, an arrival time difference observation value calculating means for calculating an observation value of an arrival time difference between the first arrival wave from the sound source of the estimation target (target) and each of the second and subsequent arrival waves is configured, of which the first The multipath wave arrival time difference calculation unit 207 calculates an observation value of the arrival time difference between the first arrival wave and the second and subsequent arrival waves to the first receiver 201, and the second multipath wave arrival time The time difference calculation unit 208 calculates an observed value of the arrival time difference between the first incoming wave to the second receiver 202 and each of the second and subsequent incoming waves.

  The arrival time difference table creation unit 209 receives the virtual sound source setting position data WSP supplied via the input terminal 213, the sound velocity profile data SVP supplied via the input terminal 214, and the reception supplied via the input terminal. Based on the wave position data WRP, an arrival time difference table is created and stored in the storage unit 209M.

For example, the virtual sound source position when the arrival time difference table is created is set as shown in FIG.
In the example shown in FIG. 9, the search area is virtually divided into a mesh (lattice) shape, and each intersection is defined as a virtual sound source position (p, q) ((p, q) = (1, 1) to (P, Q). ). Note that p and q are not the values representing the distance r and the depth z, but are indices corresponding to the distance and the depth, and are represented by natural numbers 1 to P and 1 to Q.

  The mesh width (interval of virtual sound source positions) is set according to the estimation accuracy of a desired sound source position. As the mesh width is narrowed, measurement can be performed with higher accuracy. However, the amount of data and calculation increases, and real-time processing becomes difficult. Therefore, the mesh width is narrowed when high accuracy is required, and the mesh width is widened when high speed is important.

A propagation path from the virtual sound source position WS is assumed as shown in FIG. 10, for example. FIG. 10 shows only five propagation paths WP1 to WP5 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 SA or the seabed SB at least once and has reached the receiver is called a multipath wave. Of the multipath waves, one that is reflected once by the sea surface SA is called a sea surface once reflected wave, and one that is reflected once by the sea floor is called a sea bottom once reflected wave.

Based on such an assumption, the arrival time difference from each of the virtual sound source positions is estimated, an estimated value (theoretical value) is obtained as a theoretical time difference value, and stored in the storage unit 209M. The arrival time difference table creation unit 209 assumes that the pulsed sound wave is emitted at each of the virtual sound source positions, the estimated arrival time of the first arrival wave (direct wave) in the upper receiver 201, and the upper Estimated value or theoretical value EU _i (i = _i ) of the difference (arrival time difference) from the arrival time of the multipath wave estimated to arrive at the receiver 201 at the i-th (i = 1, 2,...) Upward elevation angle. 1 or ₂ ) is obtained for all the virtual sound sources, stored in the storage unit 209M, and similarly, the first receiver 202 in the lower receiver 202 has the _first or second group GEU (EU ₁ , EU ₂ ,... EU _I ). The estimated arrival time of one incoming wave (direct wave), and the arrival time of the multipath wave estimated to arrive at the jth (j = 1, 2,...) Estimated value or theoretical value EL _j (j = 1, 2) ,...), Or a sequence GEL (= EL ₁ , EL ₂ ,... EL _J ) is obtained for all virtual sound sources and stored in the storage unit 209M. In this way, the storage unit 209M plays a role as arrival time difference theoretical value storage means.

  Note that I is empirically estimated to be receivable by the receiver 201 and is a value equal to or greater than the maximum number of multipath waves whose elevation angle is upward, and J is empirically receivable by the receiver 202. The estimated angle of elevation is set to a value equal to or greater than the maximum number of multipath waves with downwards.

A sequence GEU of arrival time differences at the upper receiver 201 for each of all virtual sound source positions (p, q) ((p, q) = (1, 1) to (P, Q)), that is, EU ₁ (p, q), EU ₂ (p, q),... ((P, q) = (1, 1) to (P, Q))
Are obtained as table format data, and stored in the storage unit 209M of the arrival time difference table creation unit 209 as the upward elevation angle arrival time difference theoretical value.
Similarly, for each of the virtual sound source positions (p, q), a sequence of arrival time differences GEL at the lower receiver 202, that is, EL ₁ (p, q), EL ₂ (p, q), ... ((p, q) = (1,1) to (P, Q))
Are obtained as table format data and stored in the storage unit 209M of the arrival time difference table creation unit 209 as a downward elevation angle arrival time difference theoretical value.

  The cost function calculation unit 210 receives outputs (signals indicating arrival times) TA and TB of the first and second arrival time determination units 205 and 206 and calculates the first and second multipath wave arrival time difference. The arrival time difference data (observation values) DU and DL obtained by the units 207 and 208 and the arrival time difference table values (theoretical values) EU and EL created by the arrival time difference table creation unit 210 are input, and based on these A cost function is obtained, and the obtained cost function is supplied to the sound source position estimation unit 211.

  The sound source position estimation unit 211 uses the distance and depth of the virtual sound source that minimizes the cost function created by the cost function calculation unit 210 as the estimated value of the sound source position, and outputs the estimated value from the output terminal 212.

Hereinafter, a configuration example of the arrival time difference table creation unit 209 will be described in more detail with reference to FIG.
From the input terminal 213, the setting value of the index P indicating the distance r of the virtual sound source and the index Q indicating the depth z is input as the virtual sound source setting position data WSP, and the sound velocity profile data SVP is input from the input terminal 214. Receiver position data WRP is input from the input terminal 215.
The sound ray calculation processing unit 703 performs sound ray calculation based on these data. As a result of the sound ray calculation, an estimated value (theoretical value) of the arrival time of each incoming wave and an estimated value (theoretical value) of the arrival elevation angle are obtained.

  The elevation angle determination unit 704 determines whether the arrival elevation angle is upward or downward for each multipath wave from each virtual sound source position based on the calculation result in the sound ray calculation processing unit 703, and the determination result is received. This is supplied to the time difference calculation unit 705.

  The arrival time difference calculation unit 705 is based on the calculation result in the sound ray calculation processing unit 703, and further, on the multipath wave determined to be upward based on the determination result in the elevation angle determination unit 704, in the upper receiver 201. The arrival time difference (difference between the arrival times of the direct wave and the multipath wave) is calculated, and the calculated arrival time difference (estimated value) is written as a theoretical value in the up and down arrival time difference table 706 in the storage unit 209M. For the determined multipath wave, the arrival time difference (difference between the arrival time of the direct wave and the multipath wave) in the lower receiver 202 is calculated, and the calculated arrival time difference (estimated value) is stored in the storage unit 209M. Is written as a theoretical value in the up-down arrival time difference table 707.

As a result of performing such processing, the multipath arriving from the respective virtual sound source positions (p, q) and arriving at the upper receiver 201 from the upper side (at the upward elevation angle) in the uplift arrival time difference table 706. A set of data representing the arrival time difference for each of the waves is stored for each virtual sound source, and the rising and falling arrival time difference table 707 emits from each virtual sound source position (p, q) and receives the lower wave. A data set representing the arrival time difference for each of the multipath waves arriving from the lower side (at a downward elevation angle) to the device 202 is stored for each virtual sound source.
The data stored in these tables 706 and 707 are supplied to the cost function calculation unit 210 via terminals 708 and 709 as described later.

Next, a configuration example of the cost function calculation unit 210 will be described in more detail with reference to FIG.
The cost function calculation unit 210 includes an elevation angle determination unit 805, data extraction units 806 and 807, an incoming wave determination unit 810, an area filter creation unit 305, area filter application units 306 and 307, and an uplifting theoretical value -Observed value difference calculating unit 808, up and down theoretical value-observed value difference calculating unit 809, up and down sum of squares calculating unit 812, up and down sum of squares calculating unit 813, and sum calculating unit 814 .

  The input terminal 801 is supplied with the output TA of the first arrival time determination unit 205, the input terminal 802 is supplied with the output TB of the second arrival time determination unit 206, and the input terminal 811 has a third output. The output TC of the arrival time determination unit 303 is supplied, the output DU of the first multipath wave arrival time difference calculation unit 207 is supplied to the input terminal 803, and the second multipath wave arrival is supplied to the input terminal 804. The output DL of the time difference calculation unit 208 is supplied.

  The elevation angle determination unit 805 receives the outputs TA and TB of the first and second arrival time determination units 205 and 206 via the input terminals 801 and 802, respectively, and the incoming waves received by the receivers 201 and 202, respectively. For each of the above, it is determined whether the elevation angle is upward or downward, and a determination result YUL is output.

  In the example shown in FIGS. 8A and 8B, the first receiver 201 receives the first, second, and fourth incoming waves at an earlier timing, and the third, fifth, For the incoming wave, the lower receiver receives it at an earlier timing. From these facts, the first, second, and fourth incoming waves arrive from above, that is, the elevation angle is upward, and the third and fifth incoming waves arrive from below. It can be seen that the elevation angle is downward.

  Based on the determination result YUL in the elevation angle determination unit 805, the data extraction unit 806 uses the elevation angle determination unit 805 to determine the arrival time difference output from the first multipath wave arrival time difference calculation unit 207. Only data FU indicating the arrival time difference for a multipath wave whose elevation angle is determined to be upward is selected or extracted, and supplied to the elevation theoretical value-observation value difference calculation unit 808.

  Based on the determination result YUL in the elevation angle determination unit 805, the data extraction unit 807 includes the data DL representing the arrival time difference output from the second multipath wave arrival time difference calculation unit 208. Only the data FL indicating the arrival time difference for the multipath wave whose elevation angle is determined to be downward is selected or extracted and supplied to the elevation downward theoretical value-observed value difference calculation unit 809.

The arrival time difference in the upper receiver 201 for the multipath wave determined to be upward by the elevation angle determination unit 805 is represented by FU _a (a = 1, 2,...)
The arrival time difference at the lower receiver 202 for the multipath wave determined to be downward by the elevation angle determination unit 805 is represented by FL _b (b = 1, 2,...).
If each multipath propagation path is as shown in FIG. 10 (that is, as estimated), the first, second and fourth incoming waves (direct wave and first and third multipath waves) are upward. And the third and fifth incoming waves (second and fourth multipath waves) will have a downward elevation angle, and therefore
FU ₁ = DU ₂ ,
FU ₂ = DU ₄ ,
FL ₁ = DL ₃ ,
FL ₂ = DL ₅
It becomes.
The arrival time differences FU ₁ , FU ₂ ,... Are supplied to the uplifting theoretical value-observation value difference calculation unit 808, and the arrival time differences FL ₁ , FL ₂ ,. Is done.

  Note that since the difference in reception timing between the receivers 201 and 202 for the same incoming wave is extremely small, for example, if the difference in detection timing is equal to or less than a predetermined value, it may be determined that the same incoming wave has been received. . In addition, a sequence of sample values extending from and before and after the time point determined as the arrival time by the first arrival time determination unit 205, and a sample value extending from and before and after the time point determined as the arrival time by the second arrival time determination unit 206 When the correlation of the columns of the two is high, it may be determined that these arrival times are the same arrival wave detection timing.

  The incoming wave discrimination unit 810 receives information indicating arrival times by the first, second, and third arrival time determination units 205, 206, and 303 via the connection terminals 801, 802, and 811, and further by the elevation angle determination unit 805. As a result of the determination of the elevation angle, the arrival wave whose arrival time is detected by the first, second and third arrival time determination units 205, 206 and 303 is a direct wave, or a reflected wave once at the sea surface. It is determined whether there is a reflected wave once on the seabed, and the determination result DSA is supplied to the filter creation unit 305.

The area filter creation unit 305 receives the arrival times (first, second, and third receivers 201, 202, and 301) output from the first, second, and third arrival time determination units 205, 206, and 303. Detection timing) Information TA, TB and TC, virtual sound source setting position data WSP supplied via the input terminal 213, sound velocity profile data SVP supplied via the input terminal 214, and supply via the input terminal 215 Region filter information RFL is generated based on the received receiver position data WRP and the determination result DSA by the incoming wave determination unit 810 supplied via the input terminal 327,
The region filter application unit 306 performs region filtering on the theoretical value EU supplied from the rising-up arrival time difference table 706 via the input terminal 708 according to the generated region filter information RFL,
The region filter application unit 307 performs region filtering on the theoretical value EL supplied from the up-down arrival time difference table 707 via the input terminal 709 according to the generated region filter information RFL.

  The area filtering mentioned here means limiting the estimated position of the sound source. That is, by applying region filtering, the region that is not likely to be a sound source position (hereinafter referred to as “impossible region” or “ Is excluded from the sound source estimation result, or the virtual sound source position in the impossible area is not set as the sound source estimation result.

  For example, if the speed of sound in water is constant, the sound source for the direct wave is determined based on the order of arrival times in the first and second receivers 201 and 202. It can be seen whether it is located above or below the horizontal plane passing through the midpoint of the second receiver 202. Therefore, for example, if the arrival time (observation value) in the first receiver 201 is earlier than the arrival time (observation value) in the second receiver 202, an area below the horizontal plane is “impossible. If the arrival time at the second receiver 202 is earlier than the arrival time at the first receiver 201, the region above the horizontal plane is defined as the “impossible region”. Can be eliminated.

  Similarly, based on the “impossible region” based on the context of arrival times in the first and third receivers 201 and 301 and on the context of arrival times in the second and third receivers 202 and 301 “Impossible area” can be obtained.

In consideration of the fact that the sound speed in water is not constant, the boundary of the above "region" is determined by determining the propagation time from each virtual sound source position to each receiver, the position of the virtual sound source, and the sound speed. It is necessary to calculate based on the profile and the position of the receiver, and to make a determination based on the calculation result and the anteroposterior relationship of the arrival times at the respective receivers.
In addition, in consideration of detection errors and calculation errors of arrival times, it is desirable to slightly narrow the area to be excluded as an area that is not likely to be a sound source position.
A region other than the region that is not likely to be a sound source position is a region that may be a sound source position, and may be hereinafter referred to as a “possible region” or a “passing region”.
Information indicating whether each virtual sound source position is in the possible area or in the impossible area is binary information F (p, q) = true / for each of the virtual sound source positions (p, q). It takes the form of a data table composed of false.

  Filtering by the area filter application units 306 and 307 is a process (narrowing process) for preventing the virtual sound source position in the impossible area from being an estimation result, and cannot be performed prior to the calculation of the cost function. Exclude virtual sound source positions in the area (do not calculate cost function for excluded virtual sound source positions), or eliminate cost function for virtual sound source positions in impossible area after calculating cost function Thus, the virtual sound source position in the impossible area does not become the estimation result of the sound source position, and the area filter information created by the area filter creation unit 305 (the distinction between the possible area and the impossible area is determined). Information shown) Narrowing is performed based on RFL. The area filter creation unit 305 will be described in more detail later.

The uplifting theoretical value-observation value difference calculation unit 808 extracts the arrival time difference observation value DU supplied from the first multipath wave arrival time difference calculation unit 207 via the input terminal 803 by the data extraction unit 806. The elements FU ₁ , FU ₂ ,... Of the observed arrival time difference column GFU for the uplifted multipath wave, and the uplift supplied from the upturning arrival time difference table 706 via the input terminal 708. Upward arrival time difference theoretical value EU ₁ (p, q), EU ₂ (p, q),... ((P, q) = (1,1) to (P, Q))
Are subjected to region filtering by the filter application unit 306, that is,
EU ₁ (p, q), EU ₂ (p, q), (where (p, q) is (1,1) to (P, Q)) F (p, q) = true Satisfying)
Each of the elements is received, and the difference between the corresponding elements is calculated. In other words, the i-th element of the theoretical value of the rising-upward arrival time difference that has passed through the region filter application unit 306 is EU _i (p, q), and the observed value of the arrival time difference for the rising-up multipath wave is i. When the _ith element is FU _i , the calculation represented by the following equation is performed.
ΔFU _i (p, q) = EU _i (p, q) −FU _i (3)

The elevation-up-square sum calculation unit 812 obtains the sum (expressed by the following equation (4)) of all the multipath waves upward in the elevation angle although the calculation result of the equation (3) is squared.

  In Equation (4), I ′ is I ′ ≦ I, and when a multipath wave with an upward elevation angle is detected up to the I-th, I ′ = I and not detected until the I-th. In this case, the detected elevation angle represents the number of multipath waves having an upward direction.

Similarly, the up-down theoretical value-observed value difference calculation unit 809 is a data extraction unit 807 among the observed arrival time difference values supplied from the second multipath wave arrival time difference calculation unit 208 via the input terminal 802. .., Each element FL ₁ , FL ₂ ,... Of observation time difference row GFL for the up-and-down multipath wave extracted in step S, and a region filter from the up-down arrival time difference table 707 via the input terminal 709. Upward / downward arrival time difference values EL ₁ (p, q), EL ₂ (p, q),... ((P, q) = (1,1) to (P, Q) supplied via the application unit 306 ))
Are subjected to region filtering by the filter application unit 307, that is,
EL ₁ (p, q), EL ₂ (p, q), (where (p, q) is (1,1) to (P, Q)) F (p, q) = true Satisfying)
Each of the elements is received, and the difference between the corresponding elements is calculated. That is, the i-th element of the theoretical value of the rising and falling arrival time difference that has passed through the region filter application unit 307 is EL _i (p, q), and the observed value of the arrival time difference for the rising and falling multipath wave is j. If the _j- th element is EL _j , the calculation represented by the following equation is performed.
ΔFL _j (p, q) = EL _j (p, q) −FL _j (5)

The up-and-down square sum calculation unit 813 obtains the sum (represented by the following expression (6)) of all the multipath waves with the up-down angle downward, although the calculation result of Expression (5) is squared.

  In Equation (6), J ′ is J ′ ≦ J, and when a multipath wave with a downward elevation angle is detected up to the Jth, J ′ = J and not detected until the Jth. In this case, the detected elevation angle represents the number of multipath waves having a downward direction.

The sum calculation unit 814 performs the following calculation using the calculation results of the equations (4) and (6), and outputs the calculation result from the connection terminal 815.
K (p, q) = KU (p, q) + KL (p, q) (7)
K (p, q) in Equation (7) is a cost function used in the first embodiment.

  The sound source position estimation unit 211 obtains coordinates (p, q) at which the cost function given by equation (7) is minimized, and the obtained coordinates (p, q) represent the sound source position WS of the estimation target (target). It is estimated that.

  Both the sums of squares obtained by the equations (4) and (6) become smaller as the degree of coincidence between the observed value and the theoretical value (table value) increases. Therefore, the smaller the value of the cost function in the equation (7) given by the sum of the square sums of the equations (4) and (6), the larger the degree of coincidence (the degree of inconsistency) is small. Therefore, obtaining the coordinates (p, q) that minimize the cost function of Equation (7) means obtaining the coordinates (p, q) having the highest degree of coincidence (the smallest degree of inconsistency).

  As described above, the sound source position can be accurately estimated by obtaining the cost function and obtaining the coordinates (p, q) that minimize the cost function.

Next, the region filter creation unit 305 will be described.
FIG. 13 is a block diagram illustrating a configuration example of the region filter creation unit 305.
The illustrated region filter creation unit 305 includes a direct wave region filter creation unit 311, a sea surface reflected wave region filter creation unit 312, a seabed reflected wave region filter creation unit 313, a region filter combination unit 314, and an arrival time table creation unit. 315.

From the input terminal 321, the arrival time (information indicated) TA detected by the first arrival time determination unit 205 is input to the direct wave region filter creation unit 311 and the sea surface reflected wave region filter creation unit 312.
From the input terminal 322, the arrival time (information indicating) TB detected by the second arrival time determination unit 206 is input to the direct wave region filter creation unit 311 and the submarine reflected wave region filter creation unit 313.
From the input terminal 323, the arrival time (information shown) TC detected by the third arrival time determination unit 303 is generated by the direct wave region filter creation unit 311, the sea surface reflected wave region filter creation unit 312 and the sea bottom reflected wave region filter creation. Input to the unit 313.

From the input terminals 213, 214, and 215, virtual sound source setting position data WSP, sound velocity profile data SVP, and information indicating the positions of the first, second, and third receivers 201, 202, and 301 (receiver position data) ) WRP is input and input to the arrival time table creation unit 315.
The arrival time table creation unit 315 has a storage unit 315M, creates arrival time data as follows, and stores it in the storage unit 315M in the form of a table.

  The arrival time table creation unit 315 is configured so that the direct wave generated by the sound emitted from the virtual sound source position (p, q) is first based on the receiver position data WRP, the virtual sound source setting position data WSP, and the sound velocity profile SVP. , Times of arrival at the second and third receivers 201, 202, and 301 (relative times, for example, the time since the sound was generated by the virtual sound source, that is, the propagation time) SAd (p, q), SBd (p , Q) and SCd (p, q), and the time (relative to the first and third receivers 201 and 301) that the reflected waves of the sea surface from the sound emitted from the virtual sound source position (p, q) arrive at the first and third receivers 201 and 301, respectively. Time, for example, the time since the sound was generated by the virtual sound source, that is, the propagation time) SAa (p, q) and SCa (p, q), and the seabed once by the sound emitted from the virtual sound source position (p, q) The reflected waves are second and third, respectively. Times arriving at the receivers 202 and 301 (relative time, for example, time since sound is generated in the virtual sound source, that is, propagation time) SBb (p, q) and SCb (p, q) are used as all virtual sound sources. The calculation is performed for the positions (1, 1) to (P, Q), and the calculation result is stored in the form of a table in the storage unit 315M.

  Data of the arrival time table held in the storage unit 315M is read by the direct wave region filter creation unit 311, the sea surface reflected wave region filter creation unit 312 and the sea bottom reflected wave region filter creation unit 313, and each region filter information Used to create

  The direct wave region filter creation unit 311 is supplied from the first, second, and third arrival time determination units 205, 206, and 303 via the input terminals 321, 322, and 323, and includes the first, second, and third. Receiving information TA, TB and TC indicating the detection timing of the incoming waves in the receivers 201, 202 and 301, receiving the discrimination result DSA by the incoming wave discrimination unit 810 via the input terminal 327, and further receiving the arrival time table. The arrival time table data SAd (p, q), SBd (p, q), and SCd (p, q) are received from the creation unit 315, and the discrimination result DSA among the information TA, TB, and TC indicating the arrival detection timing is received. Based on this, information indicating the arrival detection timings TAd, TBd and TCd of the direct wave is identified or extracted, and the first, second and third receivers 201, 201 and 3 are identified. In 1, to create a direct wave region filter information RFd based on the context of the arrival time of the direct wave.

  The sea surface reflected wave region filter creation unit 312 creates information TA and TC indicating arrival wave detection timing and arrival time table creation in the first and third receivers 201 and 301 supplied from the input terminals 321 and 323. Based on the arrival time table data SAa (p, q) and SCa (p, q) from the unit 315 and the discrimination result DSA by the arrival wave discrimination unit 810, discrimination is made among the information TA and TC indicating the arrival detection timing. Based on the result DSA, the information TAa and TCa indicating the arrival detection timing of the once reflected sea surface are identified or extracted, and before and after the arrival time of the once reflected sea surface in the first and third receivers 201 and 301 Sea surface reflected wave region filter information RFa is created based on the relationship.

  The submarine reflected wave region filter creation unit 313 creates information TB and TC indicating arrival timing of arrival waves in the second and third receivers 202 and 301 supplied from the input terminals 322 and 323, and an arrival time table. Based on the data SBb (p, q) and SCb (p, q) of the arrival time table from the unit 315 and the discrimination result DSA by the arrival wave discrimination unit 810, the discrimination is performed among the information TB and TC indicating the arrival detection timing. Based on the result DSA, the information TBb and TCb indicating the detection timing of arrival of the once reflected seabed are identified or extracted, and the bottom reflected wave area filter for the once reflected seabeds in the second and third receivers 202 and 301 is obtained. Information RFb is created.

  The region filter combining unit 314 includes the direct wave region filter information RFd created by the direct wave region filter creating unit 311, the sea surface reflected wave region filter information RFa created by the sea surface reflected wave region filter creating unit 312, and the sea bottom reflected wave region filter. The submarine reflected wave region filter information RFb created by the creation unit 313 is combined to create region filter information RFL, which is output from the output terminal 328.

FIG. 14 shows a configuration example of the direct wave domain filter creation unit 311.
The illustrated direct wave domain filter creation unit 311 includes a first direct wave domain filter creation unit 331, a second direct wave domain filter creation unit 332, a third direct wave domain filter creation unit 333, and a direct wave domain filter combination. Part 334.

  The first direct wave region filter creation unit 331 is supplied via the input terminals 321 and 322, and includes information TA and TB indicating arrival wave detection timing in the first and second receivers 201 and 202, and arrival. The arrival time table data SAd (p, q) and SBd (p, q) from the time table creation unit 315 and the discrimination result DSA by the arrival wave discrimination unit 810 are input, and the first and second arrival time determination units First direct wave region filter information RFd1 is created based on the arrival times TAd and TBd of the direct waves detected at 205 and 206.

ここで、¬はブール代数における論理否定演算を表し、
Ｆ^＊ｄ１（ｐ，ｑ）は以下の式で定義されるブール値を表す。

ＳＡｄ（ｐ，ｑ）、ＳＢｄ（ｐ，ｑ）は到来時刻テーブル作成部３１５の記憶部３１５Ｍに記憶されている伝播時間の理論値である。 First, the first direct wave region filter creation unit 331 determines the front-rear relationship between the direct wave arrival time TAd in the first receiver 201 and the direct wave arrival time TBd in the second receiver 202. This determination is made by comparing the arrival time TAd with the arrival time TBd, that is, whether TAd <TBd. Then, based on the determination result, the calculation represented by the following equation (8) is performed, and information Fd1 (p, q) indicating whether or not each of the virtual sound source positions (p, q) is located within the possible area. q).

Where ¬ represents a logical negation in Boolean algebra,
F ^* d1 (p, q) represents a Boolean value defined by the following equation.

SAd (p, q) and SBd (p, q) are theoretical values of the propagation time stored in the storage unit 315M of the arrival time table creation unit 315.

  A set of the information Fd1 (p, q) (a set of all values (1, 1) to (P, Q) of the coordinates (p, q)) constitutes the first direct wave domain filter information RFd1. .

  If the sound speed in water is constant, the boundary of the region defined by the above equation (8) is as shown by the line BL1 in FIG. In FIG. 15, a line BL <b> 1 is a midpoint between the first receiver 201 and the second receiver 202 (on the vertical line VL, and intermediate between the first receiver 201 and the second receiver 202. A horizontal plane (a plane perpendicular to the vertical line VL) passing through the point) is shown.

With respect to the direct wave, as shown in FIGS. 8A and 8B, if the arrival time TA ₁ in the first receiver 201 is earlier than the arrival time TB ₁ in the second receiver 202, As shown in FIG. 15, the sound source WS is determined to be located above the horizontal plane BL1, and accordingly, the area below the horizontal plane BL1 is determined to be an area (impossible area) that cannot be the position of the sound source WS.

Concerning the direct wave, contrary to the example shown in FIGS. 8A and 8B, the arrival time TA _{1 at} the first receiver 201 is later than the arrival time TB _{1 at} the second receiver 202. If there is, the sound source WS is determined to be located below the horizontal plane BL1, and therefore the area above the horizontal plane BL1 is determined to be an impossible area.

  In FIG. 15, a line BL0 represents a plane that passes through the line VL connecting the first receiver 201 and the second receiver 202 and is perpendicular to the straight line that points to the direction DWS of the sound source. The direction of the sound source is already known by another method, and the third wave receiver 301 is arranged in the direction of the sound source, so that the side opposite to the third wave receiver 301 is centered on the boundary surface BL0. There can be no sound source. The line BL0 is also shown in FIGS. 16, 17, 18, and 19 described later, and has the same meaning. Note that the direction of the sound source is not represented by a specific angle, but is merely information indicating that it is within a certain angle range, so it is expressed as a “boundary surface” instead of a “boundary line”. The The same applies to the boundary represented by the above-described line BL1 and lines BL2 to BL5 described later.

  Information indicating an area to be excluded (impossible area) as described above, or information indicating an area that should not be excluded (possible area, that is, an area that may be the position of a sound source) The first direct wave region filter information RFd1 created by the direct wave region filter creation unit 331 is supplied to the region filter combination unit 334.

  The second direct wave region filter creation unit 332 is supplied via the input terminals 321 and 323, and includes information TA and TC indicating arrival timing of the incoming waves in the first and third receivers 201 and 301, and arrival. First and third arrival time determination units using the arrival time table data SAd (p, q) and SCd (p, q) from the time table creation unit 315 and the discrimination result DSA by the arrival wave discrimination unit 810 as inputs. Based on the direct wave arrival times TAd and TCd detected in 205 and 303, second direct wave region filter information RFd2 is created.

ここで、¬はブール代数における論理否定演算を表し、
Ｆ^＊ｄ２（ｐ，ｑ）は以下の式で定義されるブール値を表す。

ＳＡｄ（ｐ，ｑ）、ＳＣｄ（ｐ，ｑ）は到来時刻テーブル作成部３１５の記憶部３１５Ｍに記憶されている到来時刻の理論値である。 First, the second direct wave region filter creation unit 332 determines the front-to-back relationship between the arrival time TAd of the direct wave in the first receiver 201 and the arrival time TCd of the direct wave in the third receiver 301. This determination is made based on whether Ad <TCd. Then, based on the determination result, the calculation represented by the following expression (10) is performed, and information Fd2 (p, q) indicating whether or not each of the virtual sound source positions (p, q) is located within the possible area. q).

Where ¬ represents a logical negation in Boolean algebra,
F ^* d2 (p, q) represents a Boolean value defined by the following equation.

SAd (p, q) and SCd (p, q) are the theoretical values of the arrival times stored in the storage unit 315M of the arrival time table creation unit 315.

  The set of the information Fd2 (p, q) (the set for all values (1, 1) to (P, Q) of the coordinates (p, q)) constitutes the second direct wave domain filter information RFd2. .

  If the sound speed in water is constant, the boundary of the region defined by the above equation (10) is as shown by the line BL2 in FIG. In FIG. 16, the line BL <b> 2 passes through the midpoint of the first receiver 201 and the third receiver 301, and is perpendicular to the line Lac connecting the first receiver 201 and the third receiver 301. Indicates a surface (region boundary surface).

For direct waves, as shown in FIGS. 8A and 8C, if arrival time TC ₁ in third receiver 301 is earlier than arrival time TA ₁ in first receiver 201, As shown in FIG. 16, the sound source WS is determined to be located obliquely below the surface BL2, and accordingly, the region obliquely above the surface BL2 is determined to be an impossible region.

Concerning the direct wave, contrary to the example shown in FIGS. 8A and 8C, the arrival time TC ₁ in the third receiver 301 is later than the arrival time TA ₁ in the first receiver 201. If there is, the sound source WS is determined to be located obliquely above the surface BL2, and therefore, the region obliquely below the surface BL2 is determined to be an impossible region.

  Information indicating a region to be excluded as described above, or information indicating a region that should not be excluded (a region that may be a position of a sound source) is created by the second direct wave region filter creation unit 332. The second direct wave region filter information RFd2 is supplied to the region filter coupling unit 334.

  The third direct wave domain filter creation unit 333 is supplied via the input terminals 322 and 323, and includes information TB and TC indicating the detection timing of the incoming waves in the second and third receivers 202 and 301, and arrival. The second and third arrival time determination units receive the arrival time table data SBd (p, q) and SCd (p, q) from the time table creation unit 315 and the discrimination result DSA by the arrival wave discrimination unit 810 as inputs. Based on the direct wave arrival times TBd and TCd detected at 206 and 303, third direct wave region filter information RFd3 is created.

ここで、¬はブール代数における論理否定演算を表し、
Ｆ^＊ｄ３（ｐ，ｑ）は以下の式で定義されるブール値を表す。

ＳＢｄ（ｐ，ｑ）、ＳＣｄ（ｐ，ｑ）は到来時刻テーブル作成部３１５の記憶部３１５Ｍに記憶されている到来時刻の理論値である。 First, the third direct wave region filter creation unit 333 determines the front-to-back relationship between the direct wave arrival time TBd in the second receiver 202 and the direct wave arrival time TCd in the third receiver 301. This determination is made based on whether TBd <TCd. Then, based on the determination result, the calculation represented by the following expression (12) is performed, and information Fd3 (p, q) indicating whether each of the virtual sound source positions (p, q) is located within the possible area. q).

Where ¬ represents a logical negation in Boolean algebra,
F ^* d3 (p, q) represents a Boolean value defined by the following equation.

SBd (p, q) and SCd (p, q) are theoretical values of arrival times stored in the storage unit 315M of the arrival time table creation unit 315.

  The set of the information Fd3 (p, q) (the set for all the values (1, 1) to (P, Q) of the coordinates (p, q)) constitutes the third direct wave region filter information RFd3. .

  If the sound speed in water is constant, the boundary of the region defined by the above equation (12) is as shown by the line BL3 in FIG. In FIG. 16, a line BL3 passes through the midpoint of the second receiver 202 and the third receiver 301 and is perpendicular to the line Lbc connecting the second receiver 202 and the third receiver 301. Indicates a surface (region boundary surface).

With respect to the direct wave, as shown in FIGS. 8B and 8C, if the arrival time TC ₁ in the third receiver 301 is earlier than the arrival time TB ₁ in the second receiver 202, As shown in FIG. 16, the sound source WS is determined to be located obliquely above the surface BL3, and accordingly, the region obliquely below the surface BL3 is determined to be an impossible region.

Concerning the direct wave, contrary to the example shown in FIGS. 8B and 8C, the arrival time TC ₁ in the third receiver 301 is later than the arrival time TB ₁ in the second receiver 202. If there is, it is determined that the sound source WS is located obliquely below the surface BL3, and therefore the region obliquely above the surface BL3 is determined to be an impossible region.

  The third direct wave region filter creation unit 333 creates information indicating a region to be excluded as described above, or information indicating a region that should not be excluded (a region that may be a position of a sound source). The third direct wave region filter information RFd3 is supplied to the region filter coupling unit 334.

  The direct wave region filter coupling unit 334 includes first, second, and third direct wave region filter information RFd1, RFd2 created by the first, second, and third direct wave region filter creating units 331, 332, and 333. , And RFd3 are combined to create direct wave domain filter information RFd and output from the output terminal 338.

FIG. 16 shows not only the boundary of the region defined by the second and third direct wave region filter information RFd2 to RFd3, but also the boundary of the region defined by the first direct wave region filter information RFd1. Yes.
The first direct wave filter creation unit 331 determines that the sound source position is above the boundary surface BL1, and the second direct filter creation unit 332 determines that the sound source position is obliquely below the boundary surface BL2, When the third direct filter creation unit 333 determines that the sound source position is diagonally above the boundary surface BL3, the sound source position is above the boundary surface BL1 and the boundary by combining these determination results. The region is obliquely below the surface BL2 and obliquely above the boundary surface BL3.

The direct wave region filter coupling unit 334 performs synthesis or combination of the determination results as described above. That is, the direct wave region filter combining unit 334 combines the determination results from the first, second, and third direct wave filter creating units 331, 332, and 333, and all the conditions obtained from the respective determination results are combined. The region to be filled is output as a determination result after combining (information indicating whether or not each virtual sound source position can be a sound source position by the direct wave 311).
The combination of the determination results by the first to third direct wave region filter creation units 331 to 333 for each virtual sound source position (p, q) is expressed by the following equation (14).
Fd (p, q)
= Fd1 (p, q) ∧Fd2 (p, q) ∧Fd3 (p, q) (14)
∧ represents a logical AND operation in a Boolean algebra.
A set of determination results Fd (p, q) represented by the above formula (14) (a set of all values (1, 1) to (P, Q) of coordinates (p, q)) is a direct wave. The area filter information RFd is configured.

  The direct wave region filter information RFd generated by the direct wave region filter creation unit 311 is supplied to the region filter combination unit 314.

ここで、¬はブール代数における論理否定演算を表し、
Ｆ^＊ａ（ｐ，ｑ）は以下の式で定義されるブール値を表す。

ＳＡａ（ｐ，ｑ）、ＳＣａ（ｐ，ｑ）は到来時刻テーブル作成部３１５の記憶部３１５Ｍに記憶されている到来時刻の理論値である。 The sea surface reflected wave region filter creation unit 312 determines the front-rear relationship between the arrival time TAa of the sea surface reflected wave in the first receiver 201 and the arrival time TCa of the sea surface reflected wave in the third receiver 301. To do. This determination is made based on whether TAa <TCa. Then, based on the determination result, an operation represented by the following expression (15) is performed, and information Fa (p, q) indicating whether or not each of the virtual sound source positions (p, q) is located within the possible area. q).

Where ¬ represents a logical negation in Boolean algebra,
F ^* a (p, q) represents a Boolean value defined by the following equation.

SAa (p, q) and SCa (p, q) are theoretical values of arrival times stored in the storage unit 315M of the arrival time table creation unit 315.

  A set of the information Fa (p, q) (a set of all values (1, 1) to (P, Q) of the coordinates (p, q)) constitutes the sea surface reflected wave region filter information RFa.

  If the sound speed in water is constant, the boundary of the region defined by the above equation (15) is as shown by the line BL4 in FIG. In FIG. 17, the line BL <b> 2 passes through the midpoint between the first receiver 201 and the third receiver 301 and passes through the first receiver 201 and the third receiver 301, as in FIG. 16. A plane perpendicular to the connecting line Lac is shown, and the line BL4 is a line symmetrical to the line BL2 around the perpendicular NSA to the sea surface SA at the intersection of the line BL2 and the sea surface SA. If there is a sound source in the plane indicated by the line BL4, the reflected wave at the sea surface SA (one-time reflected wave) follows the path in the plane indicated by the line BL2.

Regarding the sea surface reflected wave, the arrival time TA ₂ in the first receiver 201 is earlier than the arrival time TC ₂ in the third receiver 301 as shown in FIGS. 8 (a) and 8 (c). For example, as shown in FIG. 17, the sound source WS is determined to be located obliquely below the surface BL4, and therefore, the region obliquely above the surface BL4 is determined to be an impossible region.

Contrary to the example shown in FIGS. 8 (a) and 8 (c), the arrival time TA _{2 at} the first receiver 201 is greater than the arrival time TC _{2 at} the third receiver 301 for the sea surface reflected wave. Later, it is determined that the sound source WS is located obliquely above the surface BL4, and accordingly, the region obliquely below the surface BL4 is determined to be an impossible region.

  The information indicating the area to be excluded in this way or the information indicating the area that should not be excluded (area that may be the position of the sound source) is generated by the sea surface reflected wave area filter creation unit 312. The wave region filter information RFa is supplied to the region filter coupling unit 314.

ここで、¬はブール代数における論理否定を表し、
Ｆ^＊ｂ（ｐ，ｑ）は以下の式で定義されるブール値を表す。

ＳＢｂ（ｐ，ｑ）、ＳＣｂ（ｐ，ｑ）は到来時刻テーブル作成部３１５の記憶部３１５Ｍに記憶されている到来時刻の理論値である。 The submarine reflected wave region filter creation unit 313 determines the front-rear relationship between the arrival time TBb of the once reflected seabed wave at the second receiver 202 and the arrival time TCb of the reflected wave once at the bottom of the third receiver 301. To do. This determination is made based on whether TBb <TCb. Then, based on the determination result, the calculation represented by the following formula (17) is performed, and information Fb (p, q) indicating whether or not each of the virtual sound source positions (p, q) is located in the possible area. q).

Where ¬ represents a logical negation in Boolean algebra,
F ^* b (p, q) represents a Boolean value defined by the following equation.

SBb (p, q) and SCb (p, q) are theoretical values of arrival times stored in the storage unit 315M of the arrival time table creation unit 315.

  A set of the above information Fb (p, q) (a set of all values (1, 1) to (P, Q) of the coordinates (p, q)) constitutes the sea surface reflected wave region filter information RFb.

  If the sound speed in water is constant, the boundary of the region defined by the above equation (17) is as shown by the line BL5 in FIG. In FIG. 18, the line BL3 passes through the midpoint between the second receiver 202 and the third receiver 301 and passes through the second receiver 202 and the third receiver 301 as in FIG. A plane perpendicular to the connecting line Lbc is shown, and the line BL5 is a line symmetrical to the line BL3 about the perpendicular NSB to the sea bottom SB at the intersection of the line BL3 and the sea bottom SB. If the sound source is in the plane indicated by the line BL5, the reflected wave (one-time reflected wave) on the seabed SB follows the path in the plane indicated by the line BL3.

Regarding the one-time reflected wave at the sea bottom, as shown in FIGS. 8B and 8C, the arrival time TC ₃ in the third receiver 301 is earlier than the arrival time TB ₃ in the second receiver 202. If there is, the sound source WS is determined to be located obliquely below the surface BL5, as shown in FIG. 18, and therefore the region obliquely above the surface BL5 is determined to be an impossible region.

Contrary to the example shown in FIGS. 8 (b) and 8 (c), the arrival time TC _{3 at} the third receiver 301 is greater than the arrival time TB _{2 at} the second receiver 202 for the sea surface reflected wave. Later, it is determined that the sound source WS is located obliquely above the surface BL5, and therefore, the region obliquely below the surface BL5 is determined to be an impossible region.

  Information indicating an area determined to be excluded in this manner, or information indicating an area that should not be excluded (an area that may be a position of a sound source) is generated by the sea surface reflected wave area filter creation unit 313. The sea surface reflected wave region filter information RFb is supplied to the region filter coupling unit 314.

  Although the description has been given on the assumption that the sea bottom is horizontal, when the sea bottom is inclined and the direction and degree of inclination vary depending on the location, the direction of the surface BL5 also changes accordingly.

ここで、Ｆ^＊ｄ１１（ｐ，ｑ）、Ｆ^＊ｄ１２（ｐ，ｑ）は以下の式で定義されるブール値を表す。

ＳＡｄ（ｐ，ｑ）、ＳＢｄ（ｐ，ｑ）は到来時刻テーブル作成部３１５の記憶部３１５Ｍに記憶されている到来時刻の理論値、Δｔは、到来時刻の検出誤差（観測値誤差）及び算出誤差（理論値誤差）を考慮に入れた時間マージンである。 Creation of the first, second and third direct wave region filter information by the first, second and third direct wave region filter creation units 331, 332 and 333, and the sea surface by the sea surface reflected wave region filter creation unit 312 The detection error (observation value error) and calculation error (theoretical value error) of the arrival time are taken into consideration in each of the creation of the reflected wave region filter information and the creation of the sea surface reflected wave region filter information by the seafloor reflected wave region filter creation unit 313. Therefore, it is desirable to narrow a region (impossible region) to be excluded because it cannot be a sound source position. In order to narrow down, for example, instead of the above formulas (8) and (9) relating to the first direct wave region filter creation unit 331, the information Fd1 is expressed by the following formulas (8B), (9B), and (9C). Find (p, q).

Here, F ^* d11 (p, q) and F ^* d12 (p, q) represent Boolean values defined by the following equations.

SAd (p, q) and SBd (p, q) are theoretical values of the arrival time stored in the storage unit 315M of the arrival time table creation unit 315, and Δt is a detection error (observation value error) and calculation of the arrival time. This is a time margin taking into account an error (theoretical value error).

  Similar modifications can be applied to the second and third direct wave region filter creating units 332 and 333, the sea surface reflected wave region filter creating unit 312 and the sea bottom reflected wave region filter creating unit 313.

  The region filter combination unit 314 combines the determination result RFd by the direct wave region filter creation unit 311, the determination result RFa by the sea surface reflected wave region filter creation unit 312, and the determination result RFb by the seabed reflected wave region filter creation unit 313. In this combination, the region filter combining unit 314 outputs information indicating a region that satisfies all the conditions according to each determination result as information (region filter information) RFL indicating a region where a sound source may exist.

An example of this combination will be described with reference to FIG. In FIG. 19, lines BL0 to BL5 have the same meaning as in FIGS.
The arrival time is as shown in FIGS. 8A to 8C, and the direct wave region filter creation unit 311 is above the surface BL1, obliquely below the surface BL2, and obliquely above the surface BL3. Is output,
The sea surface reflected wave region filter creation unit 312 outputs a determination result that is obliquely below the surface BL4,
If the determination result that the submarine reflected wave region filter creation unit 313 is obliquely below the surface BL5 is output, the region other than the region satisfying all the conditions indicated by these determination results is This is an impossible area and should be excluded.

  From another viewpoint, each virtual sound source position is defined as an impossible area defined by the direct wave area filter information, an impossible area defined by the sea surface reflected wave area filter, and an undefined area defined by the seafloor reflected wave area filter. When it belongs to one of the possible areas, it is excluded as being not a sound source position, and it is determined that it is in the possible area only when it is not.

Information F (p, q) indicating whether or not each virtual sound source position can be a sound source position as a result of the combination of such determinations for each virtual sound source position (p, q), that is, by the region filter 305. ) Is expressed by the following equation (19).
F (p, q)
= Fd (p, q) ∧Fa (p, q) ∧Fb (p, q) (19)
∧ represents a logical AND operation in a Boolean algebra.
A set of the above information F (p, q) (a set of all values (1, 1) to (P, Q) of the coordinates (p, q)) constitutes the region filter information RFL.

The region filter application unit 306 corresponds to the position in the region to be excluded based on the region filter information RFL created by the region filter creation unit 305 out of the theoretical values EU (p, q) output from the table 706. The information indicating the region narrowed down by the exclusion is supplied to the uplifting theoretical value-observed value difference calculation unit 808,
The region filter application unit 307 corresponds to the position in the region to be excluded based on the region filter information RFL created by the region filter creation unit 305 out of the theoretical values EL (p, q) output from the table 707. The information indicating the area narrowed down by the exclusion is supplied to the uplifting theoretical value-observed value difference calculation unit 809.

Hereinafter, when the multipath wave acquired by the arrival time determination units 205, 206, and 303 is composed of the sea surface once reflected wave and the seabed once reflected wave, that is, as shown in FIGS. 8 (a), 8 (b), and 8 (c). The operation of the cost function calculation unit 210 will be described in the case where TA ₃ , TB ₃ , and TC ₃ can be acquired.

In the example shown in FIGS. 8A and 8B, the first incoming wave (sound wave received at timings TA ₁ and TB ₁ ) and the second incoming wave (sound wave received at timings TA ₂ and TB ₂ ). Is received at an earlier timing by the upper receiver 201, and for the third incoming wave (acoustic waves received at TA ₃ and TB ₃ ), the lower receiver has an earlier timing. Receiving. From these facts, the first and second incoming waves arrive from above, that is, the elevation angle is upward, and the third arrival wave arrives from below, that is, the elevation angle is downward. Is determined by the elevation angle determination unit 805.
In addition, the first incoming wave (timing TA ₁ , TB ₁ , TC ₁ ) is a direct wave, and the sound wave (TA ₂ , TB ₂ , TC ₂ ) that has arrived from the first upper side after the direct wave is once on the sea surface. The arrival wave discriminating unit 810 determines that the sound wave (timing TA ₃ , TB ₃ , TC ₃ ) that is a reflected wave and arrives from the first upper side after the direct wave is a reflected wave once at the sea surface.
The sound wave determined to have the elevation angle upward is extracted by the data extraction unit 806 and supplied to the elevation angle upward theoretical value-observation value difference calculation unit 808.
The sound wave determined to be the elevation angle downward is extracted by the data extraction unit 807 and supplied to the elevation angle downward theoretical value-observed value difference calculation unit 809.

As described above, the arrival time difference in the upper receiver 201 for the multipath wave determined to be upward by the elevation angle determination unit 805 is expressed as FU _a (a = 1, 2,...)
Represented by
The arrival time difference at the lower receiver 202 for the multipath wave determined to be downward by the elevation angle determination unit 805 is _expressed as FL _b (b = 1, 2,...)
In the assumed example, a = 1 of FU _a , that is,
FU ₁
Of the arrival time difference FL _b , that is b = 1, ie
Only FL ₁ has been acquired.
That is, if each multipath propagation path is as shown in FIG. 10 (ie, as estimated), the first and second incoming waves (direct wave and first multipath wave) have an upward elevation angle. And the third incoming wave (second multipath wave) will have a downward elevation angle, so
FU ₁ = DU ₂
FL ₁ = DL ₃
It becomes.
The arrival time difference FU ₁ is supplied to the uplifting theoretical value-observed value difference calculation unit 808, and the arrival time difference FL ₁ is supplied to the uplifting theoretical value-observed value difference calculation unit 809.

The uplifting theoretical value-observation value difference calculation unit 808 is extracted by the data extraction unit 806 from the observation values of the arrival time difference supplied from the first multipath wave arrival time difference calculation unit 207 via the input terminal 803. The observation value FU ₁ of the arrival time difference for the multipath wave upward and the upward arrival time difference supplied from the elevation and arrival time difference table 706 via the input terminal 708 and further via the region filter application unit 306. The difference between the time difference theoretical values EU ₁ (p, q) ((p, q) = (1, 1) to (P, Q)) is calculated.
ΔFU ₁ (p, q) = EU ₁ (p, q) −FU ₁ (20)

In the up-and-down sum-of-squares calculation unit 812, the square of the calculation result of Expression (20) is obtained.
KU (p, q) = {ΔFU ₁ (p, q)} ² (21)
The calculation of the equation (21) corresponds to the case where I ′ = 1 in the above equation (4).

Similarly, the up-down theoretical value-observed value difference calculation unit 809 is a data extraction unit 807 among the observed arrival time difference values supplied from the second multipath wave arrival time difference calculation unit 208 via the input terminal 802. The observation value FL ₁ of the arrival time difference for the up-and-down multipath wave extracted in step ₁ and the up-down signal supplied from the up-and-down arrival time difference table 707 via the input terminal 709 and further via the region filter application unit 307. The difference between the downward arrival time difference theoretical values EL ₁ (p, q) ((p, q) = (1, 1) to (P, Q)) is calculated. That is, the calculation represented by the following formula (22) is performed.
ΔFL ₁ (p, q) = EL ₁ (p, q) −FL ₁ (22)

The up-and-down square sum calculation unit 813 obtains the square of the calculation result of Expression (22) (expressed by Expression (23) below).
KL (p, q) = {ΔFL ₁ (p, q)} ² (23)
The calculation of Expression (23) corresponds to the case where J ′ = 1 in Expression (6).

The sum calculation unit 814 performs the following calculation using the calculation results of Expression (21) and Expression (23), and outputs the calculation result from the connection terminal 815.
K (p, q) = KU (p, q) + KL (p, q) (24)
K (p, q) in equation (24) is a cost function used in the first embodiment.

  The sound source position estimation unit 211 obtains coordinates (p, q) that minimize the cost function given by the equation (24), and the obtained coordinates (p, q) represent the sound source position WS of the estimation target (target). It is estimated that.

  Both the sums of squares obtained by Expression (21) and Expression (23) become smaller as the degree of coincidence between the observed value and the theoretical value (table value) increases. Therefore, the smaller the value of the cost function of the equation (24) given by the sum of the square sums of the equations (21) and (23), the larger the degree of coincidence (the degree of inconsistency) is small. Therefore, obtaining the coordinates (p, q) that minimizes the cost function of Equation (24) means obtaining the coordinates (p, q) having the highest degree of coincidence (the smallest degree of inconsistency).

  As described above, the cost function is obtained for the virtual sound source positions narrowed down by the area filter application units 306 and 307, and the coordinates (p, q) that minimize the cost function are obtained, thereby accurately estimating the sound source position. It can be carried out.

The effects of the present invention will be described below.
First, the effect when a sufficient number of multipath waves can be obtained will be described.
The cost function of the conventional example has a small cost function difference in the region along the direction indicated by the arrow DMC in FIG. 4, and therefore, when the arrival time difference observation value has an error, the sound source position estimation error is large. However, in the first embodiment, as described above, a separate arrival time difference table is used for each of the upward and downward arrival angles of the multipath wave, and a separate degree of inconsistency (a value representing) CU CL is obtained, and the sound source position is estimated using these sums as a cost function, thereby reducing the estimation error.

  As an example, let us consider a case where the sound speed is constant and there is no change in water depth, as in the conventional example. FIG. 20 shows a cost function when there is no error in the arrival time difference observation value under the same conditions as in the conventional example. As in the conventional example, the horizontal axis represents the horizontal distance r from the receiving point Zr, and the vertical axis represents the depth (depth is greater in the downward direction) z. The lines in the figure are the cost function contours. As shown in the figure, the contour lines of the cost function are elliptical with the true value of the sound source position as the center, and the location close to the receiver changes so that the cost increases. Therefore, it can be seen that the shape of the cost function is qualitatively improved.

  Similarly to the conventional example, when an empirically known observation error is given and simulation is performed 30 times under the conditions shown in FIG. 3, the standard error of the distance is about 5% of the true distance value. Compared to the conventional example, it was about 1/3, and it was confirmed that it contributed to the reduction in the position estimation error of the sound source.

As described above, by using the configuration of FIG. 5, when the number of pulses that can be acquired is sufficiently large, the estimation error of the sound source position can be reduced even when there is an error in the arrival time difference observation value. .
However, when the number of pulses that can be acquired is small, there is a problem that the estimation error cannot be made sufficiently small.
That is, even in the case of the configuration disclosed in the prior application, when only the sea surface reflected wave and the seabed reflected wave can be acquired as multipath waves, as shown in FIG. Has a shape extending in one direction, and the cost function changes gradually in the direction of the arrow DMC in the figure, so that there is a problem that the estimation error of the sound source is large. On the other hand, as described above, the values narrowed down by the region filter application units 306 and 307 out of the theoretical values of the separate arrival time difference tables for the upward and downward arrival angles of the multipath wave are used. 22, by obtaining separate discrepancies (representing values) CU and CL for each of the upward and downward elevation angles, and estimating the sound source position using these sums as a cost function, as shown in FIG. Can be reduced.

  In the first embodiment, the difference in arrival time between the two receivers is used to determine the elevation angle. However, instead of the two receivers, for example, a vertical receiver array is used. It is also possible to use it.

Ｋｆ（ｐ，ｑ）は変換後のコスト関数である。
音源位置推定部２１１は、この変換後のコスト関数が最小となる位置（仮想音源位置のいずれか）を音源位置と推定する。 In the above example, in calculating the cost function, the cost function is calculated only for the virtual sound source positions other than the excluded area, excluding the virtual sound source positions determined not to be the positions of the sound sources. But,
Instead of doing this, a cost function may be once calculated for all virtual sound source positions, and then converted into a cost function excluding a region that cannot be the position of the sound source. For example, the cost function calculation unit 210 is configured as shown in FIG. 23, the region filter application unit 816 is inserted on the output side of the sum calculation unit 814, and the cost function K (p, q) calculated by the sum calculation unit 814 is inserted. Is converted as follows.

Kf (p, q) is a cost function after conversion.
The sound source position estimation unit 211 estimates a position (any one of the virtual sound source positions) at which the converted cost function is minimum as a sound source position.

  In the above embodiment, three receivers are used for region filtering. However, four or more receivers may be used. However, it is desirable to align the two receivers in the vertical direction so that the third and subsequent receivers do not have straight lines connecting any two receivers parallel to each other.

  In the above example, the operation in the case where only the sea surface once reflected wave and the seabed once reflected wave can be acquired as a multipath wave in addition to the direct wave has been described in detail, but three or more multipath waves could be acquired. Even in this case, the sound source position can be estimated using the filtering by region filtering in the same manner as described above.

Embodiment 2. FIG.
In the method of the first embodiment described above, there is no problem when the sound waves received by the first and second receivers 201 and 202 are composed only of signals, but noise is detected as a signal erroneously (false detection). There is a problem if a signal to be detected is not detected. In the case of false detection, there is a possibility that the difference in arrival time between receivers is calculated with an incorrect combination, and cost calculation may be performed based on the result of incorrect vertical elevation angle determination, resulting in an error in cost calculation. There is a possibility of incorrect sound source position estimation. When a signal to be detected is not detected, a pair of signals (reception timings) to be used for determining the elevation angle are not prepared, so that comparison with the arrival time difference table (calculation of the mismatch degree) can be performed appropriately. There is no problem. There is a similar problem with the creation of region filter information by the region filter creation unit 305, and it is desirable to create region filter information using only reliable arrival time data.

  The overall configuration of the sound source position estimation apparatus according to Embodiment 2 of the present invention is as shown in FIG. 5, but is different in that the configuration of the cost function calculation unit 210 is as shown in FIG.

  The cost function calculation unit 210 shown in FIG. 24 is generally the same as the cost function calculation unit 210 shown in FIG. 12, but in addition to the members shown in FIG. 24, a time difference determination unit 825 and a data extraction unit 858 are provided. Provided that data extraction units 826 and 827 are provided instead of the data extraction units 806 and 807.

  FIG. 25 shows a configuration example of the time difference determination unit 825. The illustrated time difference determination unit 825 includes an arrival time difference minimum combination calculation unit 853, an inter-receiver maximum arrival time difference calculation unit 855, and a reliability evaluation unit 854.

In the arrival time difference minimum combination calculation unit 853, the first and second receivers 201 and 202 are based on the outputs of the first and second arrival time determination units 205 and 206 supplied via the input terminals 801 and 802. Among the arrival times (timing) of the incoming waves at, a combination or pair having the smallest difference is obtained. The operation will be described with reference to FIGS. 26 (a) and 26 (b). The arrival time (timing) in the upper receiver 201 detected by the first arrival time determination unit 205 is detected by tA _a (a = 1, 2,...) And detected by the second arrival time determination unit 206. The arrival time (timing) in the lower receiver 202 is tB _b (b = 1, 2,...), And the combination having the smallest difference among the arrival times in the two receivers is a highly reliable combination. Or as a pair.

For example, for each of the time tA _a determined as a high combined to pairs reliable most things difference small among the time tB _b, time _{tB b (b = 1,2, ...} ) among the time _tA a ( Those in which the difference is not judged to be minimal for any of a = 1, 2,...) are excluded as having low reliability. Further, when it is determined that the difference between a certain time tB _b is the minimum for a plurality of times tAa, _a time tA _a having _a smaller difference from the time tB _b is selected, and the other times tA _a are low in reliability. Is excluded.

In the illustrated example, the time tB ₃ is excluded because the difference is not the minimum for any of the times tA _{1 to} tA ₆ . On the other hand, none of the time tA ₄ and time tA _5, although the difference between the time tB ₅ is minimum, since the direction of the time tA ₄ than the time tA ₅ is small difference between the time tB _5, the time tA ₄ is It remains as highly reliable (selected to be paired with time tB ₅ ) and time tA ₅ is excluded as less reliable.

As a result, time tA ₁ and time tB ₁ , time tA ₂ and time tB ₂ , time tA ₃ and time tB ₄ , time tA ₄ and time tB ₅ , time tA ₆ and time tB ₆ are mutually reliable pairs. tG ₁ , tG ₂ , tG ₃ , tG ₄ , tG ₅ are selected (thus as the same arrival wave detection timing), and time tA ₅ and time tB ₃ are due to noise (low reliability) Considered and excluded.

In the inter-receiver maximum arrival time difference calculation unit 855, the inter-receiver distances DWRa, DWRb, and DWRc included in the receiver position data WRP input from the input terminal 215, and the sound velocity profile SVP input from the input terminal 214, From the sound velocity value C in the vicinity of the receivers 201, 202 and 301,
tmab = DWRa / C (26a)
The theoretical maximum value of the time difference between the arrival timing at the first receiver 201 and the arrival timing at the second receiver 202 of the sound wave that has followed the same propagation path is calculated. Since the first and second receivers 201 and 202 are aligned in the vertical direction, the above time difference increases as the absolute value of the angle of elevation of arrival increases.

In the reliability evaluation unit 854, whether or not the time difference (observed value) dtab of the pair of arrival timings (time) determined by the combination calculation unit 853 exceeds the inter-receiver maximum arrival time difference tmab, that is,
dtab> tmab
It is determined whether or not. This determination result indicates whether or not the pair of arrival times in the first and second receivers 201 and 202 has a high probability of being the same arrival wave detection timing (high reliability), and dtab ≦ tmab If this is the case, it means that the probability is high. Therefore, it means that the data of the arrival time is suitable for use in the determination of the elevation angle, the comparison of the arrival time difference, etc. On the contrary, when dtab> tmab If it exists, it means that the probability is low and the data of the arrival time is not suitable for use in the determination of the elevation angle, the comparison of the arrival time difference, and the like.

  The reliability evaluation unit 954 further includes a pair of detection timings tA and tB of the first and second arrival waves that are determined to have high reliability for each of the detection timings tC of the arrival wave in the third receiver 301. It has a function of determining whether or not it is sufficiently close to each other, determining that only those determined to be sufficiently close are highly reliable and supplying them to the incoming wave determination 810 and the region filter creation unit 305.

  For each of the detection timings tC, whether or not the pair of detection timings tA and tB of the first and second arrival waves determined to have high reliability is close is determined by, for example, the arrival times (timing of the pair). ) Define the sum of the absolute values of the differences from each of tA and tB, or the sum of the squares of the differences from the arrival times (timing) tA and tB constituting the pair as the difference in time of the pair, This may be performed depending on whether the difference from the defined time is small (for example, compared with a predetermined value or compared with another time).

For example, in the example shown in FIGS. 26A to 26C, each of the pairs tG _g (g = 1, 2,...) Of time tA _a and time tB _b (time tA, tB constituting these pairs tGg). Are extracted as times TA and TB), the one with the smallest difference among the times tC _c (c = 1, 2,...) Is left as a highly reliable one, and the time tC _c Thus, the case where the difference is not judged to be the smallest for both the time tA _a and the time tB _b pair tG _g is excluded as having low reliability.

In the example shown in the figure, the time tC ₁ is close to the pair tG ₁ of the times tA ₁ and tB ₁ , so it is determined that the reliability is high, and is combined with the pair tG _1, and the time tC ₂ is the time tA ₂ and tB _2. since the pair tG ₂ close to it is determined to be reliable, combined with pairs tG _2, time tC ₃ is judged to be reliable it is close to the time tA ₃ and counter tG ₃ of tB _4, While combined with pairs tG _3, time tC _4, since not close with any pair, it is determined that the reliability is low. Since the time tC ₅ is close to the pair tG ₄ of the times tA ₄ and tB ₅ , it is determined that the reliability is high and is combined with the pair tG _4, and the time tC ₆ is set to the pair tG ₅ of the times tA ₆ and tB ₆ . is determined to be reliable is close, combined with pairs tG _5.

  Further, for each of the times tC determined to be highly reliable in this way, the difference between the detection times at the combined receivers 201 and 202 is the receiver 201 and the receiver 301 or the received wave. It is also possible to confirm whether the time difference is equal to or less than the theoretical maximum value of the time difference between the receiver 202 and the receiver 301 (when confirmed, the data (detection timing) is determined to be highly reliable). .

The theoretical maximum value tmac of the time difference between the arrival timing at the first receiver 201 and the arrival timing at the third receiver 301, and the arrival timing and the third reception at the second receiver 202 The theoretical maximum value tmbc of the time difference of arrival timing in the unit 301 is tmac = DWRb / C (26b), similar to the above equation (26a).
tmbc = DWRc / C (26c)
Is required.

The reliability evaluation unit 854 combines the determination result for the time tC with the determination result for the times tA and tB to determine whether the set of the times tA, tB, and tC is highly reliable. Judgment is made comprehensively, and the judgment result YVN is output.
The determination result YVN is supplied to the data extraction units 858, 806, and 807.

  In the data extraction unit 858, the determination result YVN among the outputs tA, tB, and tC of the first, second, and third arrival time determination units 205, 206, and 303 supplied via the input terminals 801, 802, and 811. Are extracted as times TA, TB, and TC, and are supplied to the arrival wave discriminating unit 810, the region filter creating unit 305, and the elevation angle determining unit 805.

  The elevation angle determination unit 805 determines whether the elevation angle is upward or downward with respect to data (TA, TB) representing the time pair supplied from the data extraction unit 858, and outputs a determination result YUL. The determination result YUL is supplied to the data extraction units 826 and 827 and the arrival wave determination unit 810.

  The data extraction unit 826 in FIG. 24 uses the determination result by the elevation angle determination unit 805 (determination result YUL whether the elevation angle is upward or downward) and the determination result by the time difference determination unit 825 (reliability of each detection timing tA, tB, tC). Of the output DU (supplied via the input terminal 803) of the first multipath wave arrival time difference calculation unit 207 at the arrival time extracted by the data extraction unit 858 according to the determination result YVN). Only the arrival time difference data value FU that corresponds to the elevation angle determination unit 805 and whose elevation angle is determined to be upward is extracted and supplied to the elevation upward theoretical value-observation value difference calculation unit 808.

  Similarly, the data extraction unit 827 determines the determination result by the elevation angle determination unit 805 (the determination result YUL whether the elevation angle is upward or downward) and the determination result by the time difference determination unit 825 (reliability of each detection timing tA, tB, tC). The arrival time extracted by the data extraction unit 858 in the output DL (supplied via the input terminal 804) of the second multipath wave arrival time difference calculation unit 208 according to the determination result YVN) , And only the arrival time difference data value FL for which the elevation angle is determined to be downward by the elevation angle determination unit 805 is extracted and supplied to the elevation downward theoretical value-observed value difference calculation unit 809.

  In the up-and-down theoretical value-observed value difference calculating unit 808 and the up-down theoretical value-observed value difference calculating unit 809, the observation values FU and FL of the arrival time differences extracted by the data extracting units 826 and 827 are used. Thus, the difference is calculated and the cost function based on the difference is calculated. The operation of other parts of the cost function calculation unit 210 is the same as that described in the first embodiment.

  In the first embodiment, the cost function that discriminates the elevation angle of the incoming sound wave can be used only in an almost ideal situation without the influence of noise, but as described in the second embodiment, the incoming sound wave (Elevation angle determination using only highly reliable data, comparison of arrival time differences (detection of mismatch)), the SN ratio is not so high, Even when there is no detection, it is possible to apply a cost function that discriminates ups and downs of incoming sound waves. For this reason, even when the SN ratio is not so high and the arrival time difference is erroneously detected or not detected, the cost function can be kept conical, and the position estimation accuracy of the sound source can be maintained.

  In addition, since only the highly reliable data is used for creating the region filter information, the highly reliable region filter information can be created.

  In the second embodiment, the time difference between the two receivers is used to measure the elevation angle. However, it is possible to determine whether the elevation angle is up or down using three or more receivers.

In the case where the sound velocity profile could not be measured, a representative example is used instead of the sound velocity value C in the vicinity of the receivers 201, 202, and 301 of the sound velocity profiles of the equations (26a), (26b), and (26c). Using a sound speed of 1500 m / s,
tmab = DWRa / 1500 (27a)
tmac = DWRb / 1500 (27b)
tmbc = DWRc / 1500 (27c)
It is also possible to obtain the theoretical maximum value of the arrival time difference using.

  201, 202, 301 Receiver, 205, 206, 303 Arrival time determination unit, 207, 208 Multipath wave arrival time difference calculation unit, 209 Arrival time difference table creation unit, 210 Cost function calculation unit, 211 Sound source position estimation unit, 305 Region filter creation unit, 306, 307 Region filter application unit, 311 Direct wave region filter creation unit, 312 Sea surface reflected wave region filter creation unit, 313 Sea bottom reflected wave region filter creation unit, 314 Region filter combination unit, 810 Arrival wave discrimination unit 816 Region filter application unit.

Claims (10)

Reception of incoming waves at the first and second receivers arranged in a line on a common vertical line and a third receiver arranged at a position shifted from the vertical line at the reception position. In the sound source estimation device that estimates the position of the sound source based on the wave time,
For each of the plurality of virtual sound source positions, the theoretical value of the arrival time difference between the direct wave of the incoming waves at the receiving position and the once reflected wave at the sea surface, assuming that a virtual sound source exists at the virtual sound source position; And an arrival time difference theoretical value storage means for storing a theoretical value of an arrival time difference between the direct wave and the reflected light at the seabed once at the receiving position;
Based on the arrival wave data observed by the first and second receivers, an observed value of the arrival time difference between the direct wave from the target sound source, the sea surface reflected wave and the seabed reflected wave is calculated. Arrival time difference observation value calculation means,
Obtaining the inconsistency between the arrival time difference between the sea surface reflected wave and the direct wave and the theoretical value of the arrival time difference for the sea surface reflected wave as the first inconsistency;
The degree of inconsistency between the arrival time difference between the seabed once reflected wave and the direct wave and the theoretical value of the arrival time difference for the seabed once reflected wave is determined as a second inconsistency, and the first inconsistency and the A cost function calculating means for obtaining a cost function based on the sum of the second mismatch degree;
Sound source position estimating means for estimating the virtual sound source position at which the cost calculated by the cost function calculating means is the minimum as the position of the target sound source,
The cost function calculating means includes
Based on the anteroposterior relationship of the arrival time observation values in the first to third receivers of each of the direct wave, the sea surface once reflected wave and the seabed once reflected wave among the incoming waves, The sound source position estimation apparatus, wherein the sound source position is estimated based on the cost function for a position other than the position determined not to be a sound source position among the virtual sound source positions. Reception of incoming waves at the first and second receivers arranged in a line on a common vertical line and a third receiver arranged at a position shifted from the vertical line at the reception position. In the sound source estimation device that estimates the position of the sound source based on the wave time,
For each of the plurality of virtual sound source positions, when the virtual sound source is assumed to exist at the virtual sound source position, the arrival angle of the first incoming wave and the second and subsequent incoming waves at the receiving position is upward. The arrival time difference theory for storing the theoretical value of the arrival time difference between each of them, and the theoretical value of the arrival time difference between each of the first arrival wave and the second and subsequent arrival waves at the receiving position where the arrival angle is downward. Value storage means;
Arrival time difference observation for calculating an observation value of the arrival time difference between the first arrival wave from the target sound source and the second and subsequent arrival waves based on the arrival wave data observed in the first and second receivers. A value calculating means;
It is determined whether the arrival elevation angle of each of the second and subsequent arrival waves is upward or downward, and among each of the second and subsequent arrival waves, each of the arrival waves having an upward elevation angle, and the first arrival wave, And the theoretical value of the arrival time difference for an incoming wave with an upward angle of elevation rising as the first inconsistency, and the incoming elevation angle of the second and subsequent incoming waves is downward A degree of inconsistency between the arrival time difference between each of the first arrival waves and the first arrival wave and the theoretical value of the arrival time difference for the arrival wave having the downward angle of elevation is defined as a second inconsistency; Cost function calculating means for obtaining a cost function based on the sum of the first mismatch degree and the second mismatch degree;
Sound source position estimating means for estimating the virtual sound source position at which the cost calculated by the cost function calculating means is the minimum as the position of the target sound source,
The cost function calculating means includes
Distinguishing the first incoming wave as a direct wave,
The first wave that arrives after the direct wave among the upward waves rising upward is determined as a sea surface reflected wave,
The first wave that arrives after the direct wave among the rising waves with the elevation angle upward is determined as a reflected wave once at the bottom of the sea,
Based on the order of arrival times at the first to third receivers of the direct wave, the once reflected wave at the sea surface, and the once reflected wave at the bottom of the sea, among the incoming waves, the virtual A sound source estimation apparatus, wherein the sound source position is estimated based on the cost function for a position other than a position determined not to be a sound source position among the sound source positions. The cost function calculating means is based on the anteroposterior relation of the arrival times at the first to third receivers of the direct wave, the sea surface reflected wave and the seabed reflected wave among the incoming waves. 3. The sound source according to claim 1, wherein the first and second inconsistencies are obtained for a position other than the position determined not to be a sound source position among the virtual sound source positions. Estimating device. The cost function calculating means is based on the anteroposterior relation of the arrival times at the first to third receivers of the direct wave, the sea surface reflected wave and the seabed reflected wave among the incoming waves. 3. The sound source according to claim 1, wherein the cost function is obtained by excluding the sum obtained for a position determined not to be a sound source position among the virtual sound source positions. Position estimation device. A straight line connecting the first receiver and the second receiver;
A straight line connecting the first receiver and the third receiver;
The first, second and third receivers are arranged so that straight lines connecting the second receiver and the third receiver are not parallel to each other. The sound source position estimation apparatus according to 1 or 2. There are four or more receivers,
Of these, the first receiver and the second receiver are aligned in the vertical direction,
A straight line connecting any two of the four or more receivers;
The sound source position estimation apparatus according to claim 1 or 2, wherein the receiver is arranged so that a straight line connecting any other two receivers is not parallel to each other. Information indicating the theoretical value of the arrival time of each of the first to third receivers of the direct wave, the sea surface reflected wave, and the seabed reflected wave from each virtual sound source position is calculated in advance, The sound source position estimation device according to claim 1 or 2, wherein the sound source position estimation device is stored as an arrival time table, and is referred to when making a determination based on a context of observation values of the arrival time. . The cost function calculating means determines whether the arrival angle of the incoming wave is upward or downward based on a difference in reception time of the same incoming wave at the first and second receivers. The sound source estimation apparatus according to claim 1, wherein: The cost function calculating means includes
Calculating the sum of the squares of the difference between the observed arrival time difference and the theoretical arrival time difference for the one-time reflected wave from the sea as the first degree of inconsistency;
2. The sound source position estimation according to claim 1, wherein a sum of squares of differences between the observed arrival time difference value and the theoretical arrival time difference value for the one-time reflected reflected seabed is calculated as the second mismatch degree. apparatus. The cost function calculating means includes
Calculating the sum of the squares of the difference between the observed arrival time difference value and the theoretical arrival time difference value for an incoming wave having an upward angle of elevation as the first mismatch degree;
The sum of the squares of the difference between the observed arrival time difference value and the theoretical arrival time difference value for an incoming wave having a downward angle of elevation is calculated as the second mismatch degree. Sound source position estimation device.

Images