Frequency Invariant Beamforming for a Small-Sized Bi-Cone Acoustic Vector–Sensor Array

In this work, we design a small-sized bi-cone acoustic vector-sensor array (BCAVSA) and propose a frequency invariant beamforming method for the BCAVSA, inspired by the Ormia ochracea’s coupling ears and harmonic nesting. First, we design a BCAVSA using several sets of cylindrical acoustic vector-sensor arrays (AVSAs), which are used as a guide to construct the constant beamwidth beamformer. Due to the mechanical coupling system of the Ormia ochracea’s two ears, the phase and amplitude differences of acoustic signals at the bilateral tympanal membranes are magnified. To obtain a virtual BCAVSA with larger interelement distances, we then extend the coupling magnified system into the BCAVSA by deriving the expression of the coupling magnified matrix for the BCAVSA and providing the selecting method of coupled parameters for fitting the underwater signal frequency. Finally, the frequency invariant beamforming method is developed to acquire the constant beamwidth pattern in the three-dimensional plane by deriving several sets of the frequency weighted coefficients for the different cylindrical AVSAs. Simulation results show that this method achieves a narrower mainlobe width compared to the original BCAVSA. This method has lower sidelobes and a narrower mainlobe width compared to the coupling magnified bi-cone pressure sensor array.


Introduction
Interest in the measurement of the radiated noise has increased considerably in the last two decades, driven by a recognition of safeguarding the national marine safety. The early method used a single pressure sensor to measure the radiated noise [1]. However, as the development of vibration reduction and noise reduction technology, the radiated noise level has decreased significantly. The early method cannot extract the information of the radiated noise as it cannot provide the spatial gain. In order to obtain a higher gain, an array composed of many pressure sensors is usually used in the sonar system, which leads to a high cost. Acoustic vector sensors (AVSs) can simultaneously measure the acoustic pressure as well as two or three orthogonal particle velocity components at a single spatial point [2][3][4][5][6][7][8][9]. As the AVS can make use of the extra available acoustic particle velocity information, the arrays composed of AVSs have some attractive advantages compared to the pressure sensor arrays with the same configuration, such as higher spatial gain [10], stronger ability to suppress noise [11], lower direction-of-arrival (DOA) estimation error [3,12], etc. Owing to the above advantages, the measurement of the radiated noise by using the acoustic vector-sensor array (AVSA) has attracted broad attention.
The radiated noise produced by underwater targets usually consists of broadband and narrowband components. As the broadband signals carry more information of the underwater targets compared to the narrowband signals, the research concern has turned to the broadband signals. Beamforming is the most used spatial processing technique among the array signal processing techniques. However, the beamwidth of the traditional beamforming becomes narrower as the signal frequency increases. When the AVSA receives the broadband signal radiated from the volume element, except for the spindle direction of the mainlobe, the incident direction locating in the other directions will lead to the signal spectra distortion [13]. In order to remedy this shortcoming, the frequency invariant beamformers which could guarantee the width of the mainlobe holding invariant with the changing of the signal frequency were developed for maintaining the signal spectra non-distortion as long as the signal is incident to the AVSA from the directions of the mainlobe. There are several techniques for designing the constant beamwidth beamformer over the desired broadband. In [14], the theory and design methods of forming a frequency invariant far-field beampattern were proposed by using the relationship of a continuously distributed sensor array's aperture and frequency, which can be used for the one-, two-, and three-dimensional arrays. Based on the three-and four-dimensional Fourier transformation of their spatial and temporal parameters, the authors in Reference [15] proposed the design method of the frequency invariant beampattern for the two-and three-dimensional arrays that are composed of the linear sensor array. In [16], the method of selecting a reference beam is proposed for the fast Fourier-transformation-based frequency invariant beamforming, which is determined by the number of sensors and the ratio of the highest frequency and lowest frequency of the signal. However, these methods are concentrated at the linear arrays or the two-and three-dimensional array composed of linear arrays, where the disadvantage is that, when the direction of the mainlobe changes from the broadside to the endfire, the resolution becomes lower as the width of the mainlobe turns wider. Huang et al. used the Jacobi Anger expansion to approximate the beampattern of the uniform circular array and then designed the symmetric directivity pattern with the frequency invariant property, where the beamwidth of the mainlobe was invariant at any look direction in the sensor plane [17]. Although the beampattern function in [17] is independent from the signal frequency theoretically, it is subjected to the number of the sensors and the radius. Therefore, this method in [17] only exhibits a bandpass characteristic [18]. In order to obtain a frequency invariant beampattern over a large bandwidth, the authors in [19,20] used a class of uniform concentric circular arrays to design the frequency invariant beamforming and utilized the second-order cone programming to design the compensation filters, which led to a high running time of solving compensation filters due to the usage of the optimal toolbox. In [21][22][23], the harmonic nesting method was used to design the constant beamwidth pattern, where the principle was to design the frequency weighted coefficients to compensate the two subarrays with appropriate apertures in order to form a constant beamwidth over an octave. The researchers in References [24,25] extended the harmonic nesting method to linear AVSAs to obtain the frequency invariant beampattern over the octaves, which was superior to linear pressure sensor arrays in terms of spatial gain and anti-starboard ambiguity. Though the main disadvantage of the harmonic nesting method is the complex structure of the sensor arrays, the physical meaning is clear, and it can control the constant beamwidth beampattern perfectly.
It is worth noting that the above-mentioned methods focus on the planar arrays, the volumetric array composed of linear arrays, or the linear AVSAs. Currently, there are few studies on the volumetric AVSA. On the other hand, the above-mentioned methods of designing frequency invariant beamformers work well under appropriate intersensor spacing. However, the frequency band of the radiated noise from underwater targets generally falls within dozens to one thousand Hertz. This indicates it is a requirement for the sensor array with half of the wavelength interelement spacing to guarantee the desired beamwidth of the beamforming, resulting in a large aperture sensor array. In practice, large aperture sensor arrays are not usually used due to the disadvantages of high cost, difficult deployment, a high requirement for the application environment, and the ancillary facilities. Compared to large aperture sensor arrays, the small-sized sensor arrays are more suitable for applying to sonar systems. Thus, this paper aims to design a frequency invariant beamformer with a narrow beamwidth and low sidelobes for a small-sized volumetric AVSA.
The hearing system of the parasitoid fly Ormia ochracea gives us a hint to get out of this dilemma. To perpetuate the species, the female Ormia ochracea has to find a male cricket by using the cricket host's call. The Ormia ochracea features an advantage of locating these crickets accurately by using its two-ear differences in intensity and arrival time from an incident acoustic wave. This phenomenon was surprising to the researchers as the distance between the Ormia ochracea's two ears is only about 1.2mm, and the wavelength of the cricket's call is about 7cm, which is a significant mismatch to the distance of the Ormia ochracea's ears [26]. Studies in [27,28] explained that the reason about the Ormia ochracea creating this advantage was that the Ormia ochracea's left and right acoustical-sensory organs were not physically separated, but connected by a mechanical coupling, which increased the two-ear time delay and amplitude difference. Recently, the mechanical system of the Ormia ochracea's two ears has drawn sufficient attention and has been widely applied in the manufacture of microphones [29][30][31], the design of two-element antenna arrays [32,33], the DOA estimation based on a microphone array or an antenna array [34,35], etc. To our knowledge, this mechanical system of the Ormia ochracea's two ears has not been studied for applications of AVSAs.
In this paper, we provide a frequency invariant beamforming method for the small-sized bi-cone AVSA (BCAVSA), which is inspired by the coupling magnification mechanism of the Ormia ocharcea and the harmonic nesting method. First, we design the BCAVSA using several groups of cylindrical AVSAs with different radii and heights, where the cylindrical AVSA is composed of two uniform circular AVSAs (UCAVSAs), which is used to guide the construction of the constant beamwidth beamformer. Specifically, we derive the coupling magnified matrix for the cylindrical AVSAs in order to increase the distance of the two adjacent AVSs located in the same UCAVSA, as well as that of the two UCAVSAs. Therefore, a virtual BCAVSA is obtained. Then, the idea of the harmonic nesting is extended to the BCAVSA, and the frequency invariant beamforming method for the virtual BCAVSA is developed. This method can form the frequency invariant beampattern in the three-dimensional space, and the beamwidth in the azimuth plane is invariant when the elevation angle is fixed. Owing to the coupling magnified system, the beamwidth of the beampattern based on the proposed method is narrower than that of using the original BCAVSA. Additionally, compared to the bi-cone pressure sensor array, the proposed method has lower sidelobes and a narrower width of the mainlobe.
The remainder of this paper is organized as follows: In Section 2, the background is introduced, which includes the measurement model of the small-sized BCAVSA and the response of the mechanical coupling system of the Ormia ocharea. In Section 3, the proposed method is derived. In Section 4, the numerical simulations are presented. Finally, we provide conclusions in Section 5.

Measurement Model of BCAVSA
Incipiently, the wireless area proposed the bi-cone antenna array [36][37][38], which features the higher radiation power, flexible radiation pattern, etc. Later, the U.S. Navy developed a bi-cone sensor array measurement system to measure the radiated noise [39]. In the underwater, the broadband of the far-field radiated noise usually falls within dozens to a thousand of Hertz. In order to efficiently measure the low-frequency radiated noise, the interelement spacing of the bi-cone sensor array should be fitted to 0.5λ (λ is the signal wavelength). However, the size of the bi-cone sensor array is always restricted to some practical application platform such as the shallow sea. In this condition, the element spacing of the bi-cone sensor array is smaller or much smaller than the 0.5λ. On the other hand, the beamwidth of beamforming is wider with the signal frequency decreasing, which can lead to the signal spectra distortion [13]. Thus, this paper aims to propose a method for a small-sized bi-cone sensor array to form the constant beamwidth pattern to measure the low-frequency radiated noise.
Inspired by the bi-cone sensor array in [39], a small-sized BCAVSA is designed for forming the constant beamwidth pattern. It is noted that, in this paper, the "small-sized BCAVSA" means the interelement spacing is much smaller than the half of wavelength. As shown in Figure 1, the BCAVSA is composed of Q groups of cylindrical AVSAs. The qth cylindrical AVSA with the radius r q and height h q consists of the qth and −qth UCAVSAs, where h q = 2z q , z q is the position of the center of the qth UCAVSA on the z-axis, 1 ≤ q ≤ Q. The qth UCAVSA is constituted of N AVSs, where the mth AVS is located at {r q cos ψ m , r q nψ m , z q }, where ψ m is the angular position measured counterclockwise from the x-axis, 1 ≤ m ≤ N, 1 ≤ q ≤ Q. Directions of three acoustic particle velocity components of AVSs are oriented along the x-, y-, and z-axes, respectively. It is assumed that a planar wave impinges on the BCAVSA from (θ s , φ s ), where θ s ∈ [0 • , 360 • ] and φ s ∈ [0 • , 180 • ] denote the azimuth and elevation angles, respectively. Taking the origin of the coordinate system as the reference point, the steering vector of the qth UCAVSA at the signal frequency f can be expressed as (q = 1, 2, · · · , Q): where a p,q ( f , θ s , φ s ) = a pc,q ( f , θ s , φ s )e jkz q cos φ s is the steering vector related to the acoustic pressure components belonging to the qth UCAVSA, a pc,q ( f , θ s , φ s ) = [a pc,q,1 ( f , θ s , φ s ), · · · , a pc,q,M ( f , θ s , φ s )] T , a pc,q,m ( f , θ s , φ s ) = e jkr q nφ s cos(ψ m −θ s ) , k = 2π f /c is the wave number, c is the sound speed in the underwater, (·) T denotes the transpose operation, and j = √ −1 is the imaginary unit. ⊗ denotes the Kronecker product. u(θ s , φ s ) is the 4 × 1 the array manifold vector of the single AVS that is independent from the signal frequency and can be expressed as: u(θ s , φ s ) = [1, cos θ s nφ s , nθ s nφ s , cos φ s ] T .
At the signal frequency f , the steering vector of the BCAVSA can be expressed as: It can be seen from Equations (3) and (4) that the steering vector of the BCAVSA is composed of the Q steering vectors associated with Q cylindrical AVSAs, and the steering vector of each cylindrical AVSA is the Kronecker product of steering vectors those correspond to a two-element pressure sensor array along the z-axis, a uniform circular pressure sensor array, and a single AVS.
Note that in practice the closely spaced AVSs in passive sonar systems can encounter the undesired electromagnetic coupling. In this work, we assume that the calibration for the closely spaced AVSs is accomplished beforehand, viz., the performance of the passive sonar system will not be affected, and hence we ignore the effect of the undesired electromagnetic coupling in this work.

Mechanical Coupling System of Ormia ochracea's Ears
The female Ormia ochracea's hearing organ is so small that its two-ear distance is only about 1.2mm, but, magically, it can locate crickets quite accurately by using the cricket hosting's calling (wavelength is about 7cm). References [27,28] show that the Ormia ochracea's bilateral tympanal membranes are not physically separated but connected by a mechanical coupling which increases two-ears time delay and amplitude difference. This is the main reason that the Ormia ochracea can use their ears to find crickets. The mechanical coupling model of the Ormia ochracea's hearing organ is given in [28], as shown in Figure 2.
Cricket host's calling In this system, the spring α i and the dash-pot β i (i = 1 or 2) approximately represent the dynamical properties of the left or right tympanal membrane, sensory organs, and surrounding structures in the Ormia ochracea's two ears, the effective mass m 0 denotes all moving parts on the left/right unilateral of the intertympanal bridge, and the intertympanal bridge is composed of two rigid bars connected at the pivot by a coupling spring α 3 and the dash-pot β 3 . When a plane wave from the cricket host' call (signal frequency is about 5000Hz) is incident to the eardrums, the signals at the bilateral tympanal membranes x 1 (t) and x 2 (t) can be treated as the two input signals of the mechanical coupling system. Let y 1 (t) and y 2 (t) denote the moving displacements of the tympanal membranes, which are also the output of the mechanical coupling system. It is noted that this mechanical coupling system can be seen as a two-input two-output filter system. Except for the effect of x 1 (t) and x 2 (t), when the tympanal membranes move, tympanal membranes are affected by the inertia, damping, and recovery forces. According to Newton's second law of the motion, the dynamic equation of this mechanical structure can be built as: whereẏ i (t) andÿ i (t) (i = 1, 2) denote the speed and the acceleration of the dynamic displacements of bilateral tympanal membranes. The three terms in the left of Equation (5) denote the recovery, damping, and inertia forces, respectively. In order to solve Equation (5) and obtain the transfer function of the mechanical coupling system, the Fourier transform is used to Equation (5) with the assumption of zero initial values, and the following can be obtained: where X i ( f ) and Y i ( f ) are the Fourier transformation of the x i (t) and y i (t) with i = 1, 2, respectively. According to Ref. [35], denote the difference and the sum of the two acoustic pressure signals. Then, Equation (6) can be further written as: where Next, we analyse the performance of the response of the two ears. In [28], the distance d between two ears is approximately equal to 1.2mm. The incident plane wave impinges on the two ears from ϑ, where ϑ is the angle between the propagation direction and the normal of bilateral tympanal membranes. It is assumed that the parameters of each ear are the same, viz., β 1 = β 2 and α 1 = α 2 . The values of parameters β i , α i , and m 0 (i = 1, · · · , 3) can be found in [28]. Taking the centre of two ears as the reference point, the response of the two ears can be obtained: where τ = d cos(ϑ)/c c with c c ≈ 344m/s. To demonstrate the effect of the coupling system of the Ormia's ears, Figure 3 shows the biauricular amplitude difference of the coupling and uncoupling response (β 3 = 0 and α 3 = 0); meanwhile, Figure 4 shows the biauricular phase difference of the coupling and uncoupling response (β 3 = 0 and α 3 = 0). It can be seen from Figures 3 and 4 that the mechanical coupling system amplifies the amplitude and phase differences between the responses of the two ears. It is demonstrated that the mechanical coupling system can effectively create a larger distance between the two ears.

Proposed Method
In this section, we first introduce the mechanical coupling system of the Ormia's ears into the small-sized BCAVSA to obtain a virtual BCAVSA with a larger aperture. Then, based on the principle of the harmonic nesting, we propose a frequency invariant beamforming for the virtual BCAVSA to form the beampattern with the constant beamwidth on the azimuth and elevation planes.

Coupling Magnified BCAVSA
According to Equation (3), the steering vectors of the two-element pressure sensor array and the uniform circular pressure sensor array are two parts of constituting the steering vector of the qth cylindrical AVSA. Consequently, we first consider extending the magnification effect of the coupling system of the Ormia's ears to a two-element pressure sensor array. The two pressure sensors are identical, and they have the same spring α 1 = α 2 = α 0 and the same dash-pot β 1 = β 2 = β 0 . According to Equation (7), we can obtain the coupling magnified matrix for the two-element pressure sensor array, and it can be expressed as: (·) −1 denotes the inverse operation. We next extend the magnification effect of the coupling system of the Ormia ochracea's ears to the uniform circular pressure sensor array composed of N identical pressure sensors which own the same spring α 0 and the same dash-pot β 0 . Since a coupling magnification can be constructed among the two successive sensors, the N circular pressure sensor array can form N coupling magnifications. According to the analysis, the coupling magnified matrix for the circular array can be expressed as follows: T where . where According to Equations (9) and (10), the coupling magnified matrix for the qth cylindrical AVSA composed of the −qth and qth UCAVSAs can be deduced as (q = 1, t · · · , Q): where T l,q ( f ) is the coupling magnified matrix corresponding to the qth two-element pressure sensor array, and T c,q ( f ) is the coupling magnified matrix for the qth the uniform circular pressure sensor array. I 4 is the 4 × 4 identity matrix. According to Equation (12), the coupling magnified matrix for the BCAVSA can be obtained, and it can be expressed as: Consequently, we can apply the coupling magnified matrix T( f ) to the BCAVSA to obtain a virtual BCAVSA with the larger interelement spacing. The steering vector of the virtual BCAVSA can be expressed as: whereã It can be seen from Equation (14) that the performance of the coupling system relies on the coupling parameters. In order to obtain the good effect of the phase difference enhancement for the passive sonar system, it is necessary to set the parameters of the coupling system reasonably. According to [28], it can be known that the Ormia ochracea are sensitive to the sound signal with f 0 = 5kHz from the cricket host and the distance between the Ormia ochracea's two ears is approximately equal to d 0 = 1.2mm. When the coupling processing is performed on the signal with frequency f , the parameter m 0 remains unchanged and the parameters β i (i = 0, 3) and 2 times of the original parameters to obtain the phase difference enhancement of a pc,q ( f , θ s , φ s ), where d c,q = 2r q n(π/N) is the distance of the two adjacent AVSs of the qth UCAVSA, q = 1, · · · , Q. Similarly, in order to obtain the phase difference enhancement for a z,q ( f , φ s ) at the signal frequency f , the parameters β i (i = 0, 3) and ) 2 times of the original parameters, where f 0 = 5kHz and d 0 = 1.2mm are the sound signal from the cricket host and the distance between the Ormia ochracea's two ears, respectively.

Frequency Invariant Beamforming Method
We develop the harmonic nesting method to the virtual BCAVSA to achieve a desired frequency invariant beampattern in the azimuth and elevation planes. First, we take an easy example to understand the harmonic nesting method. We take two uniform linear arrays as the example. It is assumed that the high-frequency and the low-frequency uniform linear arrays have the same number of elements, and the interelement spacing of the high-frequency array is half of the low-frequency array. It is assumed that the desired broadband is equal to , the directivity functions P L ( f ) and P H ( f ) for the low-frequency and high-frequency linear arrays can be calculated. The weighted coefficients R L ( f ) and R H ( f ) can be calculated by using the directivity functions. Then, the R L ( f ) and R H ( f ) are used to multiply the directivity functions beamforming V( f ) can be obtained. The specific implementation process is shown in Figure 5.   Then, we introduce the above-mentioned harmonic nesting method into the virtual BCAVSA to form the constant beamwidth beampattern. The directivity pattern of the qth cylindrical AVSA of the virtual BCAVSA under the signal frequency f can be obtained by setting the weighted vector for the qth cylindrical AVSA equal to the normalized steering vector, and it can be written as: For the qth cylindrical AVSA of the virtual BCAVSA, the directivity function associated with the signal frequency f can be written as: , which is independent from the signal frequency. From Equation (16), it can be seen that the directivity function of the qth cylindrical AVSA can be treated as the product of three directivity functions which are respectively corresponding to the two-element pressure sensor array along the z-axis, the uniform circular sensor arrays, and a single acoustic vector sensor. The normalized directivity function of the qth cylindrical AVSA can be written as: where max{(·)} denotes the maximum value of the function (·).
The directivity function of the qth cylindrical AVSA is the function of the signal frequency f . As the signal frequency increases, the beamwidth of the directivity function of the qth cylindrical AVSA becomes narrower. When the AVSA receives the broadband signal from the volumetric element, the constant beamwidth pattern is quite important to obtain the signal energy of the radiated volumetric element exactly.
In order to obtain the constant beamwidth beampattern in the azimuth and elevation planes, the harmonic nesting method is extended into this paper to form the constant beamwidth pattern.
We take an octave [ f L , f H ] as an example, where f H = 2 f L . As shown in Figure 1, the qth cylindrical AVSA is composed of the −qth and qth UCAVSAs (q = 1, 2). According to Equation (17), the expressions of the directivity function of the high-frequency and low-frequency cylindrical AVSA can be written as: , We set the radius of the high-frequency cylindrical AVSA is equal to half that of the low-frequency cylindrical AVSA, viz., 2r 1 = r 2 . Meanwhile, the height of the high-frequency cylindrical AVSA is equal to half that of the low-frequency cylindrical AVSA, viz., 2z 1 = z 2 . In addition, the norm of the difference between the coupling magnified matrix T ∓1 ( f H ) associated with the high-frequency cylindrical AVSA and the coupling magnified matrix T ∓2 ( f L ) associated with the low-frequency cylindrical AVSA is equal to T ∓1 ( f H ) − T ∓2 ( f L ) 2 , and it is close to zero. Consequently, the directivity function of the high-frequency cylindrical AVSA associated with the f H is approximately equal to that of the low-frequency cylindrical AVSA corresponding to the f L : Since P u (θ, φ) is independent from the signal frequency, it can be deduced from Equation (19): To maintain a constant beamwidth of the high-and low-frequency cylindrical AVSAs under the bandwidth [ f L , f H ], the frequency weighted coefficients are necessary to introduce into the two cylindrical AVSAs. Linear combination of the directivity functions of the two cylindrical AVSAs using the frequency weighted coefficients can obtain the directivity function with a constant beamwidth under the bandwidth [ f L , f H ]. In order to obtain the frequency weighted coefficients for the two cylindrical AVSAs with different interelement spacing, we will consider the azimuth and elevation planes separately.
Firstly, to form the constant beampattern in the azimuth plane, the directivity function can be obtained by doing the following transformation, and it can be expressed as: where Γ z ( f , φ) = P z,1 ( f , φ)/P z,2 ( f , φ). R c,1 ( f ) and R c,2 ( f ) are the frequency weighted coefficients associated with the high-and low-frequency cylindrical AVSAs. θ, φ). For the azimuth plane, we set φ = φ s . Due to P z,1 ( f , φ s ) = 1. Equation (21) can be simplified as: We need two functions to solve the frequency weighted coefficients R c,1 ( f ) and R c,2 ( f ). Since P u (θ, φ s ) is independent from the signal frequency, we only calculate the frequency weighted coefficients using the directivity function V c ( f , θ, φ s ). According to [25], we can select the main axis and the half-power point of the normalized directivity function V c ( f , θ, φ s ) to calculate the R c,1 ( f ) and R c,2 ( f ): (1) when we set θ = θ s , the normalized value of V c ( f , θ, φ) is equal to 1, and it can be written as: viz., it means that the frequency response in the octave [ f L , f H ] is flat and the amplitude is equal to 1.
(2) when θ = θ s + δ θ (δ θ = BW c,0.5 /2, BW c,0.5 is of the beamwidth of the half-power point), the following can be obtained: Due to P c,1 ( f H , θ s , φ s ) = P c,2 ( f L , θ s , φ s ) = 1, by solving Equations (23) and (24), the frequency weighted coefficients can be obtained, and they can written as: Next, we consider the frequency invariant directivity function in the elevation plane. To obtain (21) can be converted into: In V 1 ( f , θ, φ) and V 2 ( f , θ, φ), V c ( f , θ, φ) and P u (θ, φ) are equivalent; meanwhile, P z,1 ( f , φ) and P z,2 ( f , φ) are the directivity functions associated with the high-and low-frequency scalar two-element arrays, respectively. Consequently, there are a couple of the frequency weighted coefficients to make the directivity function maintain a constant beamwidth in the elevation plane, and the directivity function can be expressed as: where R l,1 ( f ) and R l,2 ( f ) are the frequency weighted coefficients for the the high-and low-frequency scalar two-element arrays, respectively.
Similar to Equations (23) and (24), the frequency weighted coefficients R l,1 ( f ) and R l,2 ( f ) among the matrix V l ( f , φ) can be solved by: where δ φ = BW l,0.5 /2, BW l,0.5 is the beamwidth at the half power point. By solving Equation (28) R l,1 ( f ) and R l,2 ( f ) are equal to: According to Equations (21), (25), (26), (27), and (29), the frequency weighted coefficients for the high-and low-frequency cylindrical AVSAs P H ( f , θ, φ) and P L ( f , θ, φ) can be obtained: where It can be seen from Equation (30) that, by using the frequency weighted coefficients is obtained, which has constant beamwidth pattern in azimuth and elevation planes. Although the above-described theory and method are designed for the two cylindrical AVSAs to form the constant beamwidth over the octave, this constant beamwidth beamforming method can be extended to the multiple cylindrical AVSAs to construct the frequency invariant beampattern over multiple octaves. Consider Q cylindrical AVSAs, where the radii and heights of the 1st to the Qth cylindrical AVSAs are equal to r 1 , 2r 2 , · · · , 2 Q−1 r 1 and h 1 , 2h 1 , · · · , 2 Q−1 h 1 (h 1 = 2z 1 ), respectively. It is assumed that the broadband is equal to , the directivity functions P H,1 ( f , θ, φ) and P L,1 ( f , θ, φ) for the 1st and 2nd cylindrical AVSAs can be calculated by Equation (18); similarly, the frequency weighted coefficients R H,1 ( f , φ) and R L,1 ( f , φ) can be also calculated by Equation (30). Then, by summing R H,1 ( f , φ)P H,1 ( f , θ, φ) and R L,1 ( f , φ)P L,1 ( f , θ, φ), the constant beamwidth pattern over the [ f H /2 Q−1 , f H /2 Q−2 ] can be obtained. Similarly, using the qth and (q + 1)th cylindrical AVSAs, the directivity functions P H,q ( f , θ, φ) and P L,q ( f , θ, φ) and weighted frequency coefficients R H,q ( f , φ) and R L,q ( f , φ) can be calculated based on Equations (18) and (30). They can form the frequency invariant beampattern over the octave θ, φ), which has the same beamwidth to that of using the 1st and 2nd cylindrical AVSAs. Finally, the frequency invariant beampattern over [ f L , f H ] can be obtained by combining Q − 1 constant beamwidth pattern. The schematic diagram of forming a constant beamwidth over Q − 1 octave using a BCAVSA composed of Q cylindrical AVSAs is shown in Figure 6.
Finally, the array gain is discussed in this section. It is assumed that the signal is a unidirectional plane wave and the noise is isotropic. That is, the signal is perfectly coherent and the noise power in per unit solid angle is the same in all directions. At this time, the array gain reduces to the quantity called directivity index. In the case of a plane-wave signal in the isotropic noise impinging on the array and the BCAVSA steered in the direction of the signal, the array gain expression using the directivity function of the signal and noise can be written as [40]: where AG and DI denote the array gain and the directivity index, respectively. Ω denotes the unit solid angle.

Simulation Results
In this section, we provide some computational simulations to illustrate the frequency invariant property of the proposed method and the advantages of the proposed method compared to the original BCAVSA and the bi-cone pressure sensor array with the same configuration parameters. The BCAVSA is composed of several groups of cylindrical AVSAs. Each cylindrical AVSA consists of two UCAVSAs, and each UCAVSA is composed of eight AVSs. The directions of the acoustic particle velocity components of Each AVS is oriented along the x-, y-, and z-axes.

Original BCAVSA and Coupling Magnified BCAVSA
We consider a BCAVSA composed of two cylindrical AVSAs. The desired bandwidth is equal to f ∈ [500Hz, 1000Hz]. The radii of the first and second cylindrical AVSA are equal to r 1 = 0.225m = 0.15λ H and r 2 = 0.45m = 0.3λ H , respectively, λ H = c/ f H is the wavelength associated with f H . The two UCAVSAs among the first cylindrical AVSA respectively locates at z 1 = 0.075m = 0.05λ H and z −1 = −0.075m = −0.05λ H . The two UCAVSAs among the second cylindrical AVSA respectively locates at z 2 = 0.15m = 0.1λ H and z −2 = −0.15m = −0.1λ H . The sound speed in the underwater is equal to 1500m/s. The main beam is targeted at θ s = 180 • and φ s = 90 • . As the beampattern is a function of three variables f , θ, and φ, we only provide the beampatterns at some frequency bins. Figure 7 shows the normalized beampatterns of the coupling magnified and the original BCAVSAs at f = 500Hz and 700Hz. Figure 8 shows the beampatterns in the azimuth plane, which correspond to the coupling magnified and the original BCAVSAs, respectively, where the elevation angle is equal to Figure 9 shows the beampatterns in the elevation plane, which correspond to the coupling magnified and the original BCAVSAs, where the azimuth angle is θ = 180 • . It can be seen from Figure 7 that the original and coupling magnified BCAVSAs can both form the constant beamwidth patterns at 500Hz and 700Hz. In addition, the harmonic nesting method based on the coupling magnified BCAVSA has a much narrower mainlobe and lower sidelobes than that using the original BCAVSA. The frequency invariant property of the proposed constant beamwidth beamforming method can be again clearly discerned from Figures 8 and 9. In Figures 8 and 9, the frequency invariant beampattern of the coupling magnified BCAVSA has a 66 • beamwidth in the azimuth plane and a 20 • beamwidth in the elevation plane, while the frequency invariant beampattern of the original BCAVSA has a 98 • beamwidth in the azimuth plane and a 94 • beamwidth in the elevation plane. From Figures 7-9, we know that the coupling magnification processing inspired by the Ormia ochracea's ears can be applied in the BCAVSA to increase the interelement distance. As a result, compared to the original BCAVSA, the coupling magnified BCAVSA has lower sidelobes and a narrower mainlobe. According to Equation (31), if the mainlobe is narrower and sidelobes is lower, the denominator in Equation (31) is smaller, and thus the AG is higher. Consequently, the coupling magnified BCAVSA has a higher AG than the original BCAVSA.

Coupling Magnified BCAVSA and a Bi-Cone Pressure Sensor Array
Except for the main beam targeted at θ s = 130 • and φ s = 90 • , the other simulation conditions are the same as those in Section 4.1. The bi-cone pressure sensor array has the same configuration parameters as the BCAVSA. Figure 10 shows the beampatterns of the coupling magnified BCAVSA and bi-cone pressure sensor array in the azimuth plane, where the elevation angle is fixed to 90 • . Figure 11 shows the beampatterns of the coupling magnified BCAVSA and bi-cone pressure sensor array in the elevation plane, where the azimuth angle is fixed to 130 • . It can be seen from Figure 10 that the beampattern of the coupling magnified BCAVSA in the azimuth plane has the lower sidelobes and a narrower mainlobe compared to those of the coupling magnified bi-cone pressure sensor array. This phenomenon is caused by the inherent characteristic of the AVSA, viz., the AVSA can use the more acoustic information than the pressure sensor array, and thus the coupling magnified BCAVSA can provide lower sidelobes and a narrower mainlobe than those of the bi-cone pressure sensor array. In addition, it can be seen from Figure 11 that the beampattern of the coupling magnified BCAVSA in the elevation plane has lower sidelobes than that of the coupling magnified bi-cone pressure sensor array. In addition, the coupling magnified bi-cone pressure sensor array and the coupling magnified BCAVSA have the same width of the mainlobe. This is because when the elevation angle is equal to 90 degrees, the magnification effect of the coupling system plays a stronger role than the AVSA. In addition, from Figures 10a and 8a, we can see that when the elevation angle is invariant, and the beampatterns in Figures 10a and 8a have the same beamwidth in the azimuth plane. This is because the directivity function of the circular pressure sensor array is one part of the cylindrical AVSA. As the coupling BCAVSA has a narrower mainlobe and lower sidelobes compared to the coupling bi-cone pressure sensor array, the denominator in Equation (31) associated with coupling BCAVSA is smaller than that of the bi-cone pressure sensor array. Thus, according to Equation (31), the AG of the coupling BCAVSA is higher than the coupling bi-cone pressure sensor array.  Figure 12 shows the beampatterns of the coupling magnified BCAVSA and bi-cone pressure sensor array in the elevation plane, where the azimuth angle is equal to 130 • . The main beam is targeted at φ s = 80 • and θ s = 130 • . It can be seen from Figure 12 that the coupling magnified BCAVSA's beampattern has lower sidelobes and a narrower mainlobe than those of the coupling magnified bi-cone pressure sensor array. This is because the AVSA can use more acoustic information than the pressure sensor array. The advantage of the AVSA plays a stronger role than the magnification effect of the coupling system when the elevation angle is not equal to 90 • . In this condition, since the coupling BCAVSA has a narrower mainlobe and lower sidelobes compared to the coupling bi-cone pressure sensor array, the denominator in Equation (31) associated with coupling BCAVSA is smaller than that of the bi-cone pressure sensor array. As a result, the AG of the coupling BCAVSA is higher than the coupling bi-cone pressure sensor array.

Frequency Invariant Beampattern of BCAVSA over Multiple Octaves
A BCAVSA composed of three cylindrical AVSAs is considered. The desired bandwidth of the BCAVSA is equal to f ∈ [250Hz, 1000Hz]. The desired bandwidth can be divided into two subbands, the first and the second cylindrical AVSAs are used to form the constant beamwidth pattern in [500Hz, 1000Hz], and the second and third cylindrical AVSAs are used to construct the constant beamwidth pattern in [250Hz, 500Hz]. The radii of the three cylindrical AVSAs are equal to r 1 = 0.225m = 0.15λ H,2 , r 2 = 0.45m = 0.3λ H,2 , and r 3 = 0.9m = 0.3λ H,1 , respectively, λ H,1 = c/ f H,1 ( f H,1 = 500Hz), λ H,2 = c/ f H,2 ( f H,2 = 1000Hz). The two UCAVSAs of the three cylindrical AVSAs are located at z 1 = ∓0.075m = 0.05λ H,2 , z 2 = ∓0.15m = 0.1λ H,2 , and z 3 = ∓0.3m) = 0.1λ H,1 , respectively. The constant beamwidth pattern over this desired bandwidth is constructed by using the three cylindrical AVSAs. The sound speed in the underwater is equal to 1500m/s. The beam is targeted at θ s = 180 • , φ s = 90 • . Since the beampattern is a function of three variables f , θ, and φ, we only provide the beampatterns at some frequency bins. Figure 13 shows the normalized frequency invariant beampatterns of the coupling magnified BCAVSA, which are associated with the frequency 400Hz and 900Hz. In addition, Figure 14 shows the frequency invariant beampattern over the frequency [250Hz, 1000Hz] in the azimuth plane, where the elevation angle is fixed to 90 • . Figure 15 shows the frequency invariant beampattern over the frequency [250Hz, 1000Hz] in the elevation plane, where the azimuth angle is fixed to 180 • .  It can be seen from Figure 13a,b that the normalized beampatterns at f = 400Hz and f = 900Hz have approximately the same beamwidth. Form Figures 14 and 15, it can be seen that, by using the three cylindrical AVSAs, the proposed method can form the frequency invariant beampattern over the two octaves [250Hz, 1000Hz] in the azimuth and elevation planes. The frequency invariant property can be verified by Figures 13-15. This simulation case illustrates that, by using multiple cylindrical AVSAs constituting the BCAVSA, the proposed method can form the frequency invariant beampattern over multiple octaves in the azimuth and elevation planes.

Frequency Invariant Beamforming in the Presence of the Noise
The coupling magnified system not only amplifies the phase and amplitude differences of the incident signal but also the ambient noise [35]. In order to reduce the influence of the coupling magnified system to the noise, in this paper, we first construct the covariance matrix R vv,q of the qth cylindrical AVSA (q = 1, · · · , Q). The covariance matrix R vv,q is decomposed into signal subspace E s and noise subspace N n . We then use the signal subspace E s and the eigenvalue corresponding to the signal power to reconstruct a new covariance matrix R new = E s D s E H s , from which the influence of noise has theoretically been removed. Subsequently, we introduce the coupling magnified matrix to R new to obtain the covariance matrix of the virtual cylindrical AVSA with larger aperture. Finally, we use the frequency weighted coefficients to obtain the constant beamwidth pattern. Figures 16 and 17 show the frequency invariant beampattern, where signal-to-noise (SNR) is equal to 5dB and −5dB, respectively. The noise received by the BCAVSA is the white Gaussian noise. The other simulation conditions are as same as those in Section 4.1. It can be seen from Figure 16 that, when the SNR is equal to 5dB, the proposed method can maintain the perfect constant beamwidth beampatterns on the azimuth and elevation planes. Form Figure 17, it can be seen that, when the SNR decreases, the proposed frequency invariant beamforming can form the perfect constant beamwidth pattern in the azimuth plane. In addition, the proposed method can keep a good constant beamwidth beampattern within the mainlobe on the elevation plane, although the width outside the mainlobe of the beampattern at different frequencies are different. From Figures 16 and 17, it illustrates that the proposed frequency invariant beamforming method is not very sensitive to the noise.  Figure 18 shows AG curves of the coupling magnified BCAVSA, the coupling magnified bi-cone pressure sensor array, and the original BCAVSA, where SNR is set to be −5dB to 10dB with the interval of 1dB. The AG curves are the average of 50 times of independent computer simulations. It can be seen from Figure 18 that the coupling magnified BCAVSA has higher AG compared to the coupling magnified bi-cone pressure array and original BCAVSA. This is because, based on advantages of the coupling magnified system and AVSs, the coupling magnified BCAVSA provides the lower sidelobes and a narrower mainlobe than the coupling magnified bi-cone pressure array and original BCAVSA, and thus, according to Equation (31), the coupling magnified BCAVSA has higher AG than the other methods.

Frequency Invariant Beamforming in the Presence of the Steering Vector Error
This section considers the robustness of the proposed frequency invariant beamforming method in the presence of the steering vector error. In this simulation, the steering vector error for each cylindrical AVSA is defined as: where χ is the coefficient to modify the mismatch error. * denotes the Hadamard product.ǎ c anď a l denote the N × 1 and 2 × 1 Gaussian random process with the mean of zero and variance of 1, respectively. ϑ is the uniformly distributed random process in [0, 2π]. ϕ is the uniformly distributed random process in [0, π]. Figures 19 and 20 show the frequency invariant beampattern in the presence of the steering vector error, where χ is equal to 0.001 and 0.01, respectively. It can be seen from Figure 19 that, when the amplitude of the steering vector error is equal to 0.001, the proposed frequency invariant beamforming can hold the perfect constant beamwidth on the azimuth and elevation planes. When the amplitude of the steering vector error is equal to 0.01, it can be seen from Figure 20 that the proposed frequency invariant beamforming can hold a good constant beamwidth on the azimuth plane. However, the proposed method cannot hold the constant beamwidth on the elevation plane, and it can be also seen that, for some frequency bins, the mainlobe is not targeted at the 90 • . The main reason is that, for the elevation plane, the cylindrical AVSA is equivalent to a two-element array, and it is difficult to achieve robustness using too less number of elements. It illustrates that the proposed method can maintain the perfect constant beamwidth pattern on the azimuth plane in the presence of the steering vector error. In addition, the proposed frequency invariant beampattern on the elevation plane is sensitive to the steering vector error.

Conclusions
In this paper, we present an approach of designing a small-sized BCAVSA and a theory and method to form a characteristic of the nearly constant beamwidth in the directivity pattern. By appropriately adjusting the coupling parameters, the coupling magnified system of the Ormia ochracea's two ears is extended into the AVSA, and the coupling magnified matrix for the cylindrical AVSAs is derived in order to amplify the phase difference of the two adjacent AVSs. The cylindrical AVSAs with different radii, heights, and the same number of AVSs can be viewed as units which constitute the virtual BCAVSA. By extending the principle of the harmonic nesting to the several groups of cylindrical AVSAs, the constant beamwidth beamforming method is developed. The simulation results demonstrate the following aspects: (1) Compared to the original BCAVSA, the coupling magnified BCAVSA can provide a directivity pattern with a narrower beamwidth. This is because the coupling magnified system can convert the original BCAVSA into a virtual BCAVSA with larger interelement spacing.
(2) Since the directivity function of the circular array is one part of the cylindrical AVSA, the frequency invariant beampattern has the same beamwidth over the 360 degrees in the azimuth plane when the elevation is fixed.
(3) Since AVS can use more acoustic information than pressure sensors, the coupling magnified BCAVSA has lower sidelobes compared to the coupling magnified bi-cone pressure sensor array in both of the elevation and azimuth planes. In addition, the beampattern in azimuth plane of the coupling magnified BCAVSA has a narrower mainlobe compared to that of the coupling magnified bi-cone pressure sensor array. In the elevation plane, the coupling magnified BCAVSA has a narrower mainlobe compared to that of the coupling magnified bi-cone pressure sensor array when the mainpattern is not targeted at the elevation angle φ = 90 • since the advantage of the AVSA plays a stronger role than the effect of a coupling magnified system; when the main beam is targeted at the elevation angle φ = 90 • , the coupling magnified BCAVSA and bi-cone pressure sensor array have the same beamwidth, which is because the magnification effect of the coupling system plays a stronger role than the AVSA.
(4) By using multiple cylindrical AVSAs to constitute the BCAVSA, the proposed method can form the frequency invariant beampattern on the azimuth and elevation planes over multiple octaves.
(5) As the proposed method is based on advantages of the coupling magnified system and AVSs, the coupling magnified BCAVSA has higher AG than the other methods. (6) In order to reduce the influence of the coupling magnified system to the noise, in this paper, we use the signal subspace and the eigenvalue corresponding the signal subspace to construct a new matrix before applying the coupling magnified system. It is found that the proposed beamforming method can still obtain a perfect constant beamwidth pattern when the SNR is low.
(7) For the small steering vector error, the proposed method can hold a good constant beamwidth pattern. However, when the steering vector error is relatively large, the proposed method cannot maintain a good constant beamwidth pattern on the elevation plane, although it can keep a good constant beamwidth pattern on the azimuth plane.
Finally, in our future work, we will discuss the influence of the electromagnetic coupling for the small sensor arrays and study the frequency invariant beamforming under the condition of the electromagnetic coupling. In addition, the method to improve the robustness of this method will also be studied in the future.

Conflicts of Interest:
The authors declare no conflict of interest.