CN108845325B - Towed line array sonar subarray error mismatch estimation method - Google Patents

Towed line array sonar subarray error mismatch estimation method Download PDF

Info

Publication number
CN108845325B
CN108845325B CN201810516947.6A CN201810516947A CN108845325B CN 108845325 B CN108845325 B CN 108845325B CN 201810516947 A CN201810516947 A CN 201810516947A CN 108845325 B CN108845325 B CN 108845325B
Authority
CN
China
Prior art keywords
array
vector
matrix
error
sub
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201810516947.6A
Other languages
Chinese (zh)
Other versions
CN108845325A (en
Inventor
乔文昇
王立
雷志雄
李维科
李明兵
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Southwest Electronic Technology Institute No 10 Institute of Cetc
Original Assignee
Southwest Electronic Technology Institute No 10 Institute of Cetc
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Southwest Electronic Technology Institute No 10 Institute of Cetc filed Critical Southwest Electronic Technology Institute No 10 Institute of Cetc
Priority to CN201810516947.6A priority Critical patent/CN108845325B/en
Publication of CN108845325A publication Critical patent/CN108845325A/en
Application granted granted Critical
Publication of CN108845325B publication Critical patent/CN108845325B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S15/00—Systems using the reflection or reradiation of acoustic waves, e.g. sonar systems
    • G01S15/88—Sonar systems specially adapted for specific applications
    • G01S15/89—Sonar systems specially adapted for specific applications for mapping or imaging
    • G01S15/8906—Short-range imaging systems; Acoustic microscope systems using pulse-echo techniques
    • G01S15/8997—Short-range imaging systems; Acoustic microscope systems using pulse-echo techniques using synthetic aperture techniques

Landscapes

  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Acoustics & Sound (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • General Physics & Mathematics (AREA)
  • Measurement Of Velocity Or Position Using Acoustic Or Ultrasonic Waves (AREA)

Abstract

The invention discloses a towed line array sonar subarray error mismatch estimation method, which provides an estimation method capable of reducing positioning error, obtaining an accurate azimuth estimation value and having higher angular resolution, and is realized by the following technical scheme: in the array manifold matrix model, representing a real full array manifold matrix as a linear combination of the displacement error quantity of each sub-array and the contribution of each sub-array error quantity to the full array manifold matrix; introducing position errors among the sub-arrays into a direction-finding model, establishing a data fusion model for a full array model containing inter-sub-array displacement mismatch by applying Bayes law, solving a sub-array position error vector and a target real orientation simultaneously, and estimating sub-array displacement errors and direction of arrival simultaneously by using a Bayes algorithm according to data acquired by fusion sensor nodes to obtain a likelihood function of a multi-snapshot observation value; and obtaining an azimuth estimation value and a displacement error estimation value of the variation of the root mean square error along with the signal-to-noise ratio by using the posterior function.

Description

Towed line array sonar subarray error mismatch estimation method
Technical Field
The invention relates to a sonar which is used for embedding hydrophones on a cable to form a linear array and detecting a target in water after a towing cable is towed at a ship tail. In particular to a method for estimating the subarray displacement error and the direction of arrival of the subarray simultaneously under the condition that the subarray displacement error exists when passive aperture synthesis is carried out on a sonar towed linear array.
Background
Towed sonar is a sonar that tows a transducer array behind the tail of a vehicle platform to a target on the water side. The towed-line array sonar basic array is flexible. In the towing process, due to the influence of ship maneuvering and ocean current and self-shaking resonance, the array shape is difficult to keep stable, the array shape distortion enables the towed line array sonar to hardly reach the theoretical performance, and the problem is more serious particularly in modern array processing methods such as adaptive signal processing and space spectrum estimation. Because the conventional array processing method only accumulates the energy of the array elements, and the modern array processing method also calculates the correlation of the signals of the array elements, the characteristic value of the covariance matrix of the signals and the like, although the modern array processing method can greatly improve the precision and the resolution capability of the target, the stability of the array type is more strict. In addition, the towline array has no vertical aperture, can not distinguish water surface underwater targets, and a large number of targets generate serious interference to the performance of the towline array when in shipping intensive sea areas or formation combat, so that the phenomenon of starry sky is generated everywhere. Low-frequency active towed linear array sonar is the most effective means for detecting the quiet submarine at present. Because the towed linear array sonar array is far away from the ship and is a flexible array, the detected direction has larger error. Modern towed linear array sonar systems tend to work at low frequency, with ever-increasing detection distances, and towards higher range and detection accuracy. In such an environment, to obtain better spatial resolution, an array with a larger aperture is required, which in practice usually means more system complexity and higher equipment cost. Passive Synthetic Aperture (PSA) with an array of towed lines without changing the parameters of the towed array increases the effective aperture probability. The PSA technology utilizes the time and space correlation of signals and the motion information and the position information of an array to construct an array containing a plurality of virtual sub-arrays so as to obtain higher azimuth resolution. At present, the common moving linear array PSA methods include: yen, Carey passive synthetic aperture method, fast Fourier transform based passive synthetic aperture, and extended towed array measurement method, etc. Although the PSA with multiple virtual sub-arrays has been widely applied to the fields of sonar beam forming, direction of arrival (DOA) estimation, etc., the performance may be drastically reduced in the case of a large number of targets and the inconsistency between the motion parameters of the target ship and the real conditions. Specifically, when only a single target exists, the PSA method can work under the condition that errors occur in relative speed, a real phase correction factor is directly obtained from data, and DOA estimation can still be accurately carried out when the positions of synthesized array elements are not overlapped; however, when the relative velocity is not in error, the PSA method can work, and if the relative velocity is in error, the DOA result at this time is even worse than that of the single-array conventional beam forming, and the array expansion does not bring any benefit, but increases the positioning error. In addition to this, there are other factors that affect the accuracy of DOA estimation, such as array element spacing (usually half the operating wavelength), the angle difference between the incoming waves of the target, the coherence of the signal, etc. On the basis of PSA, accurate angle estimates of the target can be obtained, typically using high resolution DOA methods. Many high resolution methods, such as multiple Signal classification music (multiple Signal classification) multiple Signal classification methods, rotation invariant subspace methods, sparse spatial spectrum estimation, etc., can improve the spatial resolution of the array when the array parameters are precisely known. The MUSIC algorithm is a method based on matrix eigenspace decomposition, and from a geometric point of view, an observation space for signal processing can be decomposed into a signal subspace and a noise subspace, wherein the two subspaces are orthogonal. The MUSIC algorithm uses the orthogonality property between the two complementary spaces to estimate the azimuth of the spatial signal. All vectors of the noise subspace are used to construct the spectrum, and the peak positions in all spatial orientation spectra correspond to the incoming wave orientations of the signals, the basic idea is to perform eigen decomposition on the covariance matrix of any array output data, so as to obtain the signal subspace corresponding to the signal classification and the noise orthogonal to the signal components. The processing task of the MUSIC algorithm is to try to estimate the number D of the spatial signals incident on and exiting from the array, and the strength of the spatial signal sources and the direction of the incoming waves. In actual processing, the observed array output data complex vector Y results in data that is a finite number of samples in a finite time period, also known as snapshots or snapshots, but the prototype MUSIC algorithm requires that the incoming signals be incoherent. The algorithm represented by MUSIC has a disadvantage of non-ideal coherent signal processing. Among the processing schemes for coherent signal sources, spatial smoothing techniques such as Spatial Smoothing (SS) and Modified Spatial Smoothing (MSS) algorithms are more classical. However, spatial smoothing techniques are at the expense of the array effective aperture and are only applicable to equidistant Uniform Linear Arrays (ULA). When array element errors exist, an assumed clutter model is not matched with actual received data, and under the condition of model parameter mismatch, a target signal originally belonging to a signal subspace is wrongly divided into a noise subspace, so that the orthogonality of a MUSIC algorithm target function is damaged, and positioning failure is caused. Similarly, other high-resolution DOA algorithms are also sensitive to the errors of the array manifold, and the performance of the high-resolution algorithms is greatly degraded when the conditions of low signal-to-noise ratio, few fast beats and sound source coherence occur. Since the reconstructed clutter covariance matrix depends on the accuracy of the clutter model, if the assumed clutter model does not match with the actual data, the performance of the algorithm is reduced and even the algorithm fails. Therefore, in the application of the towed linear array PSA, array element errors, signal coherence, array element cross coupling, channel mismatch and the like can make an ideal model no longer be established. How to estimate the actual real array manifold vector by more effectively using the received data and perform accurate direction estimation becomes a very important topic.
Currently, a series of researches have been conducted to solve such problems both domestically and abroad. These documents correct array errors from different sides, depending on the actual error sources present in the respective application domain. For example, Weiss and Friedlander propose a method for array self-correction based on maximum likelihood estimation, which performs block-wise alternate optimization iteration on target orientation and displacement mismatch. Because the method models all array elements, when the method is applied to long-baseline, large-scale and long-time PSA, the dimensionality is increased sharply, the convergence of an objective function is very slow, and the calculation efficiency is reduced sharply. Therefore, there is a need to develop array manifold self-correction and orientation estimation methods for long baseline, large scale, long time condition towed linear array PSAs.
The towed linear array sonar is also called "towed array sonar" (may be called towed array for short). The hydrophone is embedded on a cable to form a linear array, and a towed cable is towed on a naval vessel tail to detect a target sonar in water. Mainly used for listening and detecting divingThe boat radiates noise, and remote monitoring, direction finding and identification are performed, and some of the boats can also be used for distance measurement. The device comprises a linear array matrix, a dragging cable, a winding and unwinding device, a winch, an electronic cabinet and the like. The towed linear array is composed of front, instrument, basic, back and tail segments, and has a length of tens to hundreds of meters and a variable working depth. Has large array size,Frequency of operationLow, is beneficial to line spectrum detection, can find the target in a long distance and hidden way, and the like; but has adverse effect on the maneuvering such as the turning back and backing of the towing naval vessel. The sound source, ocean channel and hydrophone array areHydroacousticsThree basic elements in the study. A sound source radiates sound signals in water, and the sound source is a source for forming a sound field; the ocean channel determines the propagation characteristics of sound waves in the ocean; the hydrophone array is used for receiving the acoustic signals and sampling the sound field distribution in water. The three parts are closely related to each other to form an inseparable unified whole. Knowing two of them, the third party can be inferred, which is the basic basis for the processing of the underwater acoustic matching field. If the hydrophone array receiving signals and ocean channel information are known, sound source information including the sound source position is to be solved, and therefore the matching field passive localization is achieved. If the hydrophone array received signal and the source information including the source position are known, the ocean channel information is to be solved, which is the Matched Field Inversion (MFI). They are all important contents of the underwater acoustic matching field processing research.
In recent years, the matching field processing technology has attracted wide attention in the aspects of underwater target detection, passive positioning, marine environment parameter inversion and the like. The Matched Field Matching Field Processing (MFP) of the underwater acoustic array signal is to calculate the amplitude and phase of the acoustic Field of a receiving matrix through an underwater acoustic Field model by utilizing marine environment parameters and acoustic propagation channel characteristics to form a copy Field vector, and match the copy Field vector with the matrix receiving data, thereby realizing the passive positioning of an underwater target and the accurate estimation of the marine environment parameters. The towed line array sonar platform based on far-near field acoustic propagation characteristics has the advantages that the plane wave propagation characteristics of far-field target signals combine a matching field positioning technology and a plane wave target direction estimation technology, and N angular frequencies exist under the assumption of far-field plane wavesAn acoustic source signal of ω is incident on a linear array having P sub-arrays. Suppose there is M inside each subarraypThe number of array elements of the whole array is
Figure GDA0003487512410000031
The position of the array element in the subarray is known accurately, and the model of the p-th subarray can be expressed as
xp(t)＝Aps(t)+ep(t), (1)
Wherein, the p sub-array manifold matrix Ap＝[ap(θ1),ap(θ2),...,ap(θN)]，θnIs the azimuth of the nth incoming wave; vector ap(θn)＝[1,exp(-jωd cos(θn)/c),...,exp(-jω(M-1)d cos(θn)/c)]TIs the p-th sub-array thetanArray manifold vectors corresponding to the directions, wherein c is sound velocity, and superscript T represents transposition; vector s (t) ═ s1(t),s2(t),...,sN(t)]TRepresenting a signal waveform vector corresponding to the t moment; vector ep(t)＝[e1(t),e2(t),...,eM(t)]TRepresenting the noise corresponding to the p-th sub-array at the time t. In general, the number N of sound sources is smaller than the number M of array elements, and incoming waves have sparseness in the space domain. Suppose a scalar rpIs the distance from the first array element of the p-th sub-array to the first array element of the first sub-array, and the array manifold vector of the full array can be expressed as
aw(θ)＝V(θ)h(θ), (2)
In the formula
Figure GDA0003487512410000041
Is an array manifold matrix within a subarray, apAnd (theta) is an array manifold vector of the p-th sub-array, h is an approximate array manifold vector between the sub-arrays, and theta is the azimuth angle of an incoming wave. Although a isp(theta) precisely known, actually measured relative position vector of subarrays
Figure GDA0003487512410000042
Generally not equal to the pre-set phase of the subarrayFor the position vector r, so the true full array manifold vector
Figure GDA0003487512410000043
With a preset array manifold vector awThere may be a deviation between (θ). An array model with displacement mismatch between sub-arrays and a target direction finding method.
Disclosure of Invention
Aiming at the defects in the prior art, the invention provides the target positioning method of the error mismatch model of the sonar towed line array subarray, which can reduce the positioning error, obtain an accurate azimuth estimation value, has higher angular resolution and can improve the robustness.
The above object of the present invention can be achieved by the following technical solutions, and a method for estimating error mismatch of towed line array sonar subarrays is characterized by comprising the following steps: in an array manifold matrix model with displacement mismatch among sub-arrays, representing a real full array manifold matrix as a linear combination of displacement error amount of each sub-array and contribution of each sub-array error amount to the full array manifold matrix; introducing the position error between the subarrays into a direction-finding model, establishing a data fusion model for a full array model containing displacement mismatch between the subarrays by applying Bayesian law, and establishing a position error vector beta of the subarrays and a real target azimuth angle alpha-1Simultaneously solving, estimating the subarray displacement error and the direction of arrival by using a Bayesian algorithm according to data acquired by fusing sensor nodes, and calculating to obtain a likelihood function of the multi-snapshot observed value; and obtaining an azimuth estimation value and a displacement error estimation value of the root mean square error along with the change of the signal-to-noise ratio by using the posterior function.
Compared with the prior art, the invention has the following beneficial effects:
the real full-array matrix under the condition of sub-array displacement error is approximated to be linear combination of the displacement error amount of each sub-array and the contribution of each sub-array error amount to the full-array manifold matrix, then a Bayes frame is used for solving, the sub-array position error vector and the target real azimuth can be simultaneously solved under the condition of low signal to noise ratio, the self-correction and the target accurate direction finding of the array can be simultaneously realized under the condition of sub-array displacement error, the position error among the sub-arrays is introduced into the direction finding model, the simultaneous estimation of the sub-array displacement error and the direction of arrival is realized through the Bayes algorithm, the error of the array manifold is subjected to insensitivity, the positioning error is reduced, and the robustness is improved.
The invention uses Bayes algorithm to calculate the variation of the root mean square error of the displacement error estimation value along with the signal-to-noise ratio as shown in FIG. 2, and the root mean square error of the displacement error estimation value can be seen to be reduced along with the increase of the signal-to-noise ratio. The root mean square error of the azimuth estimate varies with the signal-to-noise ratio as shown in figure 3 and is more or less 1 degree more stable than the root mean square error obtained using the multiple signal classification MUSIC algorithm.
The root mean square errors calculated by the Bayesian algorithm are all reduced along with the increase of the signal-to-noise ratio, are all smaller than 1 degree, and the performance is greatly improved compared with the MUSIC algorithm. The positioning results of the full-array CBF algorithm, the full-array MUSIC algorithm and the method are shown in figure 4 when the signal-to-noise ratio is 0dB, and comparison shows that the method not only can obtain an accurate azimuth estimation value, but also has higher angular resolution than other methods.
The method has obvious implementation effect in estimating the azimuth angle of the coherent sound source by using the array with the subarray displacement error. Compared with the method of directly using the full-array CBF method and the full-array MUSIC method to carry out azimuth estimation, the method has the advantages that:
(1) by introducing the subarray error into the array model, the displacement error of the subarray and the direction of arrival of the target can be estimated simultaneously;
(2) the array model containing the subarray displacement error is introduced into a Bayesian framework, so that more stable and higher-precision positioning performance can be obtained.
The invention is suitable for sensor subarray processing, radio frequency antenna array expansion, sonar passive synthetic aperture and other scenes, is mainly used for listening submarine radiation noise, and is used for carrying out remote monitoring, distance measurement, direction finding and identification array signal processing and sonar signal processing methods,
drawings
FIG. 1 is a schematic diagram of the array expansion with subarray displacement error according to the present invention.
FIG. 2 is a graphical representation of the RMS error versus SNR for the displacement error estimate of the present invention.
FIG. 3 is a graphical representation of the RMS error versus signal-to-noise ratio for the azimuth estimate of the present invention.
FIG. 4 is a graph illustrating the positioning results when the SNR is 0dB according to the present invention.
The invention is further described below with reference to the accompanying drawings.
Detailed Description
See fig. 1-4. According to the method, in an array manifold matrix model with displacement mismatch among sub-arrays, a real full array manifold matrix is represented as a linear combination of displacement error amount of each sub-array and contribution of each sub-array error amount to the full array manifold matrix; introducing the position error between the subarrays into a direction-finding model, calculating the position and the speed of a target by establishing a data fusion model and data acquired by a fusion sensor node for a full array model containing displacement mismatch between the subarrays by using a Bayesian method, and calculating a position error vector beta of the subarrays and a real azimuth angle alpha of the target-1Simultaneously solving, and calculating to obtain a likelihood function of the multi-snapshot observation values; and simultaneously estimating the subarray displacement error and the direction of arrival, and obtaining an azimuth estimation value and a displacement error estimation value of the root mean square error along with the change of the signal-to-noise ratio by using a posterior function. Comprises the following specific steps
Step 1: mismatching real inter-subarray array manifold vector
Figure GDA0003487512410000051
Performing first-order Taylor expansion approximation at a preset relative position vector r of the real relative position vector of the subarray to obtain an approximate array manifold vector between the subarrays
Figure GDA0003487512410000061
In the formula, the vector h is approximate array manifold vector between subarrays, e is Euler constant, j is imaginary unit constant, k is wave number, superscript T represents transposition, P is subarray number, and theta represents azimuth angle of incoming waveVector of motion
Figure GDA0003487512410000062
For the true relative position vector of the sub-array,
Figure GDA0003487512410000063
for the true relative position of the p-th sub-array, the vector r ═ r1,...,rP]TIs a predetermined relative position vector of the sub-array, rpFor the preset relative position of the p-th sub-array, β ═ β1,...,βP]TTo be the sub-array position error vector,
Figure GDA0003487512410000064
for the p-th sub-array displacement error, diag (β) represents a diagonal matrix with the elements of the vector β as diagonal elements.
According to the inter-subarray approximate array manifold vector obtained by first-order Taylor expansion approximation, the real full-array manifold vector a is obtainedwIs approximated to
Figure GDA0003487512410000065
Inter-subarray array manifold matrix
Figure GDA0003487512410000066
In the formula (I), the compound is shown in the specification,
Figure GDA0003487512410000067
the p-th column is v when the direction of the incoming wave is thetap(theta) the remaining columns are matrices of zero vectors, vector vpIs the p-th column, a, of the inter-subarray array manifold matrix V (theta)pAnd (theta) is an array manifold vector of the p-th sub-array.
In the full-array model, the array manifold vector of the full array is converted into the full-array manifold vector
Figure GDA0003487512410000068
Along N theta1-θNIs scanned and combined into an approximate array manifold matrix A of a full array1Then, the displacement error amount of each subarray and the error amount of each subarray are added to the full array streamThe linear combination of the contributions of the shape matrices is:
Figure GDA0003487512410000069
and the preset full-array manifold matrix A without position error between sub-arraysw＝[aw(θ1),...,aw(θN)]The p-th sub-array error projection matrix Bp＝-jk[Vp(θ1)h(θ1,r),...,Vp(θN)h(θN,r)]Where N is the number of the scanning azimuth grids, awFor a predetermined full array manifold vector, betapBpIs the first-order approximation product of the full array manifold matrix error caused by the p-th sub-array position error.
Step 2, in Bayes positioning based on displacement mismatch array model among subarrays, a Bayes algorithm is used for carrying out position error vector beta containing subarrays and target true azimuth alpha among subarrays of the full array model-1Simultaneously solving to obtain a full-array manifold matrix representing approximation
Figure GDA0003487512410000071
And an approximate full-array model full-array received signal vector x (t) ═ Φ (β) s (t) + e (t), and based on the full-array received signal vector x (t) ═ x1(t),x2(t),...,xN(t)]TSignals received by the array at the time t and a sound source vector s (t) of a waveform corresponding to the target signal at the time t1(t),s2(t),...,sN(t)]TAnd the noise vector e (t) of the full array noise at the t-th time (e)1(t),e2(t),...,eM(t)]TIn the case of multi-snapshot, the array model of the full array is rewritten into an array reception signal matrix X ═ Φ (β, θ) S + E at a plurality of snapshots, where T is the snapshot, and X ═ X (1), X (2),. once, X (T)]X (t) is the array received signal vector at time t of the array, and the matrix S ═ S (1), S (2),.., S (t)]Denotes a sound source signal matrix, s (t) is a sound source signal vector at time t, E ═ E (1), E (2),.., E (t)]Representing noise matrices, e (t) representing time tThe noise vector, θ, represents the azimuth vector of the scanned incoming wave.
Scanning azimuth vector of incoming wave
Figure GDA0003487512410000072
Scanning grid vectors for L directions, wherein L is greater than N, and signal vectors of all scanning directions at time t are
Figure GDA0003487512410000073
Matrix array
Figure GDA0003487512410000074
The method is characterized in that the method is a target signal matrix of T moments under L scanning grid conditions, and a full array model under L azimuth scanning grid conditions of an array manifold matrix with displacement mismatch among sub-arrays can be expressed as
Figure GDA0003487512410000075
In the form of matrix
Figure GDA0003487512410000076
Scanning the array manifold matrix, vector, for the full array in the presence of inter-subarray errors
Figure GDA0003487512410000077
A scanning azimuth vector representing the scanning incoming wave,
Figure GDA0003487512410000078
for the first scanning orientation,
Figure GDA0003487512410000079
scanning the manifold matrix and vector of the array for the full array without position error between sub-arrays
Figure GDA00034875124100000710
Is the true full-matrix scanning array manifold vector and matrix at the ith scanning position
Figure GDA00034875124100000711
Scanning projection matrix for p-th sub-matrix error, matrix
Figure GDA00034875124100000712
For scanning an azimuth thetalWhen the p is listed as
Figure GDA00034875124100000713
The rest columns are matrixes of zero vectors, beta is a position error vector of the submatrix,
Figure GDA00034875124100000714
is the target signal matrix for T instants under L scan grid conditions,
Figure GDA00034875124100000715
e represents the noise matrix for the signal vectors of all scan directions at time t.
The noise of each array element is independent and satisfies that the mean value is 0 and the variance is the real initial azimuth angle of the target
Figure GDA00034875124100000716
When complex Gaussian distribution is performed, the likelihood function of the multi-snapshot observed value is obtained
Figure GDA0003487512410000081
Wherein I is an L-dimensional identity matrix, det () represents a matrix determinant, exp () represents an exponentiation function, | () | electrically non-conductive2Representing the vector binorm. Noise accuracy alpha0Gamma distribution p (alpha) obeying parameters a and b0|a,b)＝Gamma(α0|a,b)
In the formula
Figure GDA0003487512410000082
Γ (a) represents a gamma function with a variable a.
Source vector
Figure GDA0003487512410000083
Obeying a complex Gaussian distribution with a probability density function of the matrix of the target signal of
Figure GDA0003487512410000084
In the formula (I), the compound is shown in the specification,
covariance matrix
Figure GDA0003487512410000085
α2The accuracy of the ith orientation. Signal precision vector parameter alpha ═ alpha1,α2,...,αL]Gamma distribution obeying parameters c and d
Figure GDA0003487512410000086
Target signal matrix
Figure GDA0003487512410000087
Has a posterior probability distribution of
Figure GDA0003487512410000088
In the formula, alpha0For the initial noise precision parameter, μ (t) represents
Figure GDA0003487512410000089
A posteriori mean vector, Σ representation
Figure GDA00034875124100000810
The posterior covariance matrix of (a).
Solving an initial noise precision parameter alpha at the t moment according to a Bayesian model updating method0The updated a posteriori mean vector μ (t) and the a posteriori covariance matrix Σ, wherein
Posterior mean vector
Figure GDA00034875124100000811
Posterior covariance matrix
Figure GDA00034875124100000812
In the formula (I), the compound is shown in the specification,
Figure GDA00034875124100000813
when the error between the subarrays is the full-array scanning array manifold matrix, H is the conjugate transpose, beta is the error vector of the subarray position,
Figure GDA00034875124100000814
for the unknown received signal vector at the time t, the matrix Λ is a prior covariance matrix, and the matrix Λ is-1Is the inverse of the matrix Λ.
Solving a signal precision parameter alpha at the t moment according to a Bayesian model updating methodiAnd a noise accuracy parameter alpha0Update noise precision parameter of
Figure GDA00034875124100000815
Updating noise accuracy parameter calculation formula
Figure GDA00034875124100000816
Updating the variance σ
Figure GDA0003487512410000091
Wherein σ is the variance, | () | calculationFFrobenius norm operator representing a vector, mean matrix H ═ μ (1), μ (2),., μ (T)]，ΣiiIs the ith diagonal element of the covariance matrix Σ, scalar γi＝1-αiΣii。
Obtaining the estimation value of the position error vector beta of the subarray according to the Bayes model iterative updating method
Figure GDA0003487512410000092
Figure GDA0003487512410000093
The first intermediate matrix T-G + Q,
second intermediate matrix
Figure GDA0003487512410000094
Third intermediate matrix
Figure GDA0003487512410000095
Fourth intermediate matrix
Figure GDA0003487512410000096
First intermediate vector
Figure GDA0003487512410000097
Second intermediate vector
Figure GDA0003487512410000098
The whole Bayesian model updating iterative process can be summarized as follows: initial noise precision parameter alpha of iteration updating method of leaf-shaped model0The signal accuracy vector alpha and the subarray position error vector beta are given initial values, the mean vector mu and the covariance matrix sigma are updated by using the updated posterior mean vector mu (t) formula (8) and the posterior covariance matrix sigma formula (9), and then the noise accuracy parameter calculation formula (10) is updated by using the formula, and the variance sigma (sigma) is updated2)newCalculating the estimate of equation (11) and the subarray position error vector beta
Figure GDA0003487512410000101
Calculation formula (12) updates initial noise precision parameter α0And repeating the above processes until convergence. Initial noise precision parameter alpha after iteration is finished0The signal precision vector alpha and the subarray position error vector beta respectively represent noise energy and characteristicsThe azimuth signal energy and the displacement error of each virtual sub-array.
The following is a specific example:
see fig. 1. The preset position of the array extension of the sub-array displacement error comprises a virtual sub-array 1, a virtual sub-array 2 and a virtual sub-array 3, the actual position deviation beta of the virtual sub-array 2 is 0.11, the actual position deviation beta of the virtual sub-array 3 is 0.2, a target number K is 2 far-field narrow-band signals are incident on an array element number M is 4-element uniform linear array, the central frequency f is 250Hz, and the incident angle theta of sound source signals is theta1And theta 260 degrees and 65 degrees respectively, and the distance between the linear array elements is 0.68 meter. The preset positions of the first array elements of the 3 virtual sub-arrays are respectively 0m, 3.4m and 6.8m away from the initial position, and the sub-array displacement error beta1,β2And beta 30 meter, 0.11m and 0.2m respectively. And acquiring a target signal matrix X by taking the snapshot number T of each position array as 200. Calculating a full array manifold matrix B caused by the position error of the p sub-array according to a formula (5)p。
2. Initial noise accuracy parameter alpha0The hyper-parameters are respectively set as a ═ b ═ 1 × 10-4C is 1, d is 0.01; initial noise precision parameter alpha of iterative process0Is set as an initial value of
Figure GDA0003487512410000102
Initial value setting of signal precision vector
Figure GDA0003487512410000103
The initial value of the sub-array position error vector β is set to β ═ 0.
3. Adding variance sigma to the signal received by each array element2Independent white Gaussian noise, defining the SNR
Figure GDA0003487512410000104
The simulated signal-to-noise ratio ranges from 0dB to 10dB, and the simulation times R at each signal-to-noise ratio are 200.
4. Updating posterior mean according to posterior mean vector mu (t) formula (8) and posterior covariance matrix sigma formula (9)Vector mu and posterior covariance matrix sigma, then updating noise precision parameter calculation formula (10) according to the formula, updating variance sigma (sigma)2)newCalculating equation (11) and an estimate of the subarray position error vector beta
Figure GDA0003487512410000105
Calculating formula (12) to update and update initial noise precision parameter alpha0Repeating the above processes until convergence to obtain the estimated value of the displacement error vector of the ith simulation
Figure GDA0003487512410000106
And orientation estimate
Figure GDA0003487512410000107
Wherein
Figure GDA0003487512410000108
The vector is formed by combining the inverses of K components with the minimum median of the signal precision vector alpha.
The root mean square error of the direction of arrival under a given signal-to-noise ratio is calculated according to a formula
Figure GDA0003487512410000111
The calculation is carried out according to the calculation,
displacement error root mean square error according to displacement error root mean square error calculation formula
Figure GDA0003487512410000112
And (4) calculating, wherein R represents the simulation times, and P represents the subarray number.
The root mean square error of the displacement error estimate, as shown in fig. 2, is calculated using a root mean square error calculation formula as a function of the signal to noise ratio, and it can be seen that the root mean square error of the displacement error estimate decreases as the signal to noise ratio increases. The root mean square error of the azimuth estimate varies with the signal-to-noise ratio as shown in figure 3 and is more or less 1 degree more stable than the root mean square error obtained using the multiple signal classification MUSIC algorithm. The root mean square errors calculated by the method are all reduced along with the increase of the signal-to-noise ratio, are all less than 1 degree, and have greatly improved performance compared with an MUSIC algorithm. The positioning results of the full-array CBF algorithm, the full-array MUSIC algorithm and the method are shown in figure 4 when the signal-to-noise ratio is 0dB, and the comparison shows that the method not only can obtain an accurate azimuth estimation value, but also has higher angular resolution than other methods.

Claims (6)

1. A towed linear array sonar subarray error mismatch estimation method is characterized by comprising the following steps: in the array manifold matrix model with the displacement mismatch among the sub-arrays, the mismatched real inter-sub-array manifold vector
Figure FDA0003585595440000011
Performing first-order Taylor expansion approximation at a preset relative position vector r of the real relative position vector of the subarray to obtain an approximate array popular vector between the subarrays
Figure FDA0003585595440000012
Figure FDA0003585595440000013
Wherein, the vector h is approximate array manifold vector between subarrays, e is Euler constant, j is imaginary unit constant, k is wave number, superscript T represents transposition, P is subarray number, theta represents azimuth angle of incoming wave, and vector
Figure FDA0003585595440000014
For the true relative position vector of the sub-array,
Figure FDA0003585595440000015
for the true relative position of the p-th sub-array, the vector r ═ r1,...,rP]TIs a predetermined relative position vector of the sub-array, rpFor the preset relative position of the p-th sub-array, β ═ β1,...,βP]TTo be the sub-array position error vector,
Figure FDA0003585595440000016
for the p-th sub-array displacement error, diag (β) represents a diagonal matrix with the elements of the vector β as diagonal elements; according to the inter-subarray approximate array manifold vector obtained by first-order Taylor expansion approximation, the real full-array manifold vector is obtained
Figure FDA0003585595440000017
Is approximately as
Figure FDA0003585595440000018
Inter-subarray array manifold matrix
Figure FDA0003585595440000019
In the formula (I), the compound is shown in the specification,
Figure FDA00035855954400000110
the p-th column is v when the direction of the incoming wave is thetap(theta) the remaining columns are matrices of zero vectors, vector vpIs the p-th column, a, of the inter-subarray array manifold matrix V (theta)p(θ) is the array manifold vector for the p-th sub-array; in the full array model, the array manifold vector of the full array is expressed
Figure FDA00035855954400000111
Along N theta1To thetaNScanning the incoming wave azimuth to form an approximate array manifold matrix A of a full array1And linearly combining the displacement error amount of each sub-array and the contribution of each sub-array error amount to the full-array manifold matrix into:
Figure FDA00035855954400000112
and preset without position error between sub-arraysFull array manifold matrix Aw＝[aw(θ1),...,aw(θN)]The p-th sub-array error projection matrix Bp＝-jk[Vp(θ1)h(θ1,r),...,Vp(θN)h(θN,r)]Where N is the number of the scanning azimuth grids, awFor a predetermined full array manifold vector, betapBpIs the first-order approximation product of the full-array manifold matrix error caused by the position error of the p-th sub-array; then introducing the position error among the subarrays into a direction-finding model, establishing a data fusion model for a full array model containing displacement mismatch among the subarrays by applying Bayes law, and in Bayesian positioning based on the displacement mismatch array model among the subarrays, using Bayes algorithm to contain a subarray position error vector beta and a target true azimuth angle alpha among the subarrays of the full array model-1Simultaneously solving to obtain a full-array manifold matrix representing approximation
Figure FDA0003585595440000021
Array model of approximate full array a full array received signal vector x (t): x (t) ═ Φ (β) s (t) + e (t), based on the full array received signal vector x (t), x (t) ═ x1(t),x2(t),...,xN(t)]TSignal received by the array at time t and sound source vector s (t) of the waveform of the target signal at time t: s (t) ═ s1(t),s2(t),...,sN(t)]TAnd the noise vector e (t) of the full array noise at time t: e (t) ═ e1(t),e2(t),...,eM(t)]TIn the case of multi-snapshot, the array model of the full array is rewritten into an array reception signal matrix X ═ Φ (β, θ) S + E at a plurality of snapshots, where T is the snapshot, and X ═ X (1), X (2),. once, X (T)]X (t) is the array received signal vector at time t of the array, and the matrix S ═ S (1), S (2),.., S (t)]Denotes a sound source signal matrix, s (t) is a sound source signal vector at time t, E ═ E (1), E (2),.., E (t)]Representing a noise matrix, e (t) representing a noise vector at the time t, and theta representing an azimuth angle vector of a scanning incoming wave; the subarray position error vector is represented as β, and the source signal matrix S is represented as S ═ S (1), S (2)...,s(T)]And E is expressed as a noise matrix, the scanning azimuth vector of the incoming wave
Figure FDA0003585595440000022
Scanning grid vectors for L directions, wherein L is greater than N, and signal vectors of all scanning directions at time t are
Figure FDA0003585595440000023
Matrix array
Figure FDA0003585595440000024
The method is characterized in that the method is a target signal matrix of T moments under L scanning grid conditions, and a full array model under L azimuth scanning grid conditions of an array manifold matrix with displacement mismatch among sub-arrays can be expressed as
Figure FDA0003585595440000025
In the form of matrix
Figure FDA0003585595440000026
Scanning the array manifold matrix, vector, for the full array in the presence of inter-subarray errors
Figure FDA0003585595440000027
A scanning azimuth vector representing the scanning incoming wave,
Figure FDA0003585595440000028
for the first scanning orientation,
Figure FDA0003585595440000029
scanning the array manifold matrix, vector for full matrix without position error between sub-matrices
Figure FDA00035855954400000210
Is a true full-matrix scanning array manifold vector and matrix at the ith scanning orientation
Figure FDA00035855954400000211
Scanning projection matrix for p-th sub-matrix error
Figure FDA00035855954400000212
For scanning an azimuth thetalWhen the p is listed as
Figure FDA00035855954400000213
The rest columns are matrixes of zero vectors, beta is a position error vector of the submatrix,
Figure FDA00035855954400000214
is the target signal matrix for T instants under L scan grid conditions,
Figure FDA00035855954400000215
e represents a noise matrix for signal vectors of all scanning directions at the moment t;
noise at each array element is independent and meets the target real initial azimuth with the mean value of 0
Figure FDA0003585595440000031
When the complex Gaussian distribution variance is obtained, the likelihood function of the multi-snapshot observed value is obtained
Figure FDA0003585595440000032
Wherein I is an L-dimensional unit matrix, det () represents a matrix determinant, exp () represents an exponentiation function, | () | circuitry2Representing the vector two norm, noise accuracy alpha0Obeying a Gamma distribution with parameters a and b, p (α)0|a,b)＝Gamma(α0|a,b)
In the formula
Figure FDA0003585595440000033
Γ (a) represents a gamma function with variable a;
source vector
Figure FDA0003585595440000034
Obeying a complex Gaussian distribution, the probability density function of the target signal matrix is
Figure FDA0003585595440000035
In the formula (I), the compound is shown in the specification,
covariance matrix
Figure FDA0003585595440000036
α2For the accuracy of the ith azimuth, the signal accuracy vector parameter α ═ α1,α2,...,αL]Gamma distribution obeying parameters c and d
Figure FDA0003585595440000037
Target signal matrix
Figure FDA0003585595440000038
Has a posterior probability distribution of
Figure FDA0003585595440000039
In the formula, alpha0For the initial noise precision parameter, μ (t) represents
Figure FDA00035855954400000310
A posteriori mean vector, Σ representation
Figure FDA00035855954400000311
A posterior covariance matrix of (a);
solving an initial noise precision parameter alpha at the t moment according to a Bayesian model updating method0And updating the a posteriori mean vector μ (t) and the a posteriori covariance matrix Σ, wherein
Posterior mean vector
Figure FDA00035855954400000312
Posterior covariance matrix
Figure FDA00035855954400000313
In the formula (I), the compound is shown in the specification,
Figure FDA00035855954400000314
when the error between sub-arrays is the full-array scanning array manifold matrix, H is the conjugate transpose, beta is the error vector of the sub-array position,
Figure FDA00035855954400000315
for the unknown received signal vector at the time t, the matrix Λ is a prior covariance matrix, and the matrix Λ is-1Is the inverse of the matrix Λ;
solving a signal precision parameter alpha at the t moment according to a Bayesian model updating methodiAnd an initial noise accuracy parameter alpha0Update noise precision parameter of
Figure FDA0003585595440000041
Updating noise accuracy parameter calculation formula
Figure FDA0003585595440000042
Updating the variance σ
Figure FDA0003585595440000043
Wherein σ is the variance, | () | non-calculationFFrobenius norm operator, mean matrix H ═ μ (1), μ (2),.., μ (T) for the vector]，ΣiiIs the ith diagonal element of the covariance matrix Σ, scalar γi＝1-αiΣii(ii) a Obtaining the estimation value of the position error vector beta of the subarray according to the Bayes model iteration updating method
Figure FDA0003585595440000044
Figure FDA0003585595440000045
The first intermediate matrix T-G + Q,
second intermediate matrix
Figure FDA0003585595440000046
Third intermediate matrix
Figure FDA0003585595440000047
Fourth intermediate matrix
Figure FDA0003585595440000048
First intermediate vector
Figure FDA0003585595440000049
Second intermediate vector
Figure FDA00035855954400000410
Estimating the subarray displacement error and the direction of arrival simultaneously by using a Bayesian algorithm according to data acquired by the fusion sensor node, and calculating to obtain a likelihood function of the multi-snapshot observation value; and obtaining an azimuth estimation value and a displacement error estimation value of the root mean square error along with the change of the signal-to-noise ratio by using the posterior function.
2. The towed line array sonar subarray error mismatch estimation method of claim 1, wherein: the preset position of the array extension of the sub-array displacement error comprises a virtual sub-array 1, a virtual sub-array 2 and a virtual sub-array 3, the actual position deviation beta of the virtual sub-array 2 is 0.11, the actual position deviation beta of the virtual sub-array 3 is 0.2, a target number K is 2 far-field narrow-band signals are incident on an array element number M is 4-element uniform linear array, the central frequency f is 250Hz, and the incident angle theta of sound source signals is theta1And theta260 degrees and 65 degrees respectively, and the distance between the linear array elements is 0.68 meter.
3. The towed line array sonar subarray error mismatch estimation method of claim 2, wherein: the preset positions of the first array elements of the 3 virtual sub-arrays are respectively 0m, 3.4m and 6.8m away from the initial position, and the sub-array displacement error beta1,β2And beta3Respectively 0m, 0.11m and 0.2m, the snapshot number T acquired by each position array is 200, and a full array manifold matrix B caused by the position error of the p-th sub-array is calculatedpAnd obtaining a target signal matrix X.
4. The towed line array sonar subarray error mismatch estimation method of claim 1, wherein: initial noise accuracy parameter alpha0The hyper-parameters are respectively set as a ═ b ═ 1 × 10-4C is 1, d is 0.01; initial noise precision parameter alpha of iterative process0Is set as an initial value of
Figure FDA0003585595440000051
Initial value setting of signal precision vector
Figure FDA0003585595440000052
The initial value of the sub-array position error vector β is set to β ═ 0.
5. The towed line array sonar subarray error mismatch estimation method of claim 1, wherein: adding variance sigma to the signal received by each array element2Independent white gaussian noise, defining the signal-to-noise ratio
Figure FDA0003585595440000053
The simulated signal-to-noise ratio ranges from 0dB to 10dB, and the simulation times R at each signal-to-noise ratio is 200.
6. The towed line array sonar subarray error mismatch estimation method of claim 1, wherein: the root mean square error of the direction of arrival under a given signal-to-noise ratio is calculated according to a formula
Figure FDA0003585595440000054
The calculation is carried out according to the calculation,
displacement error root mean square error according to displacement error root mean square error calculation formula
Figure FDA0003585595440000055
And (4) calculating, wherein R represents the simulation times, and P represents the subarray number.
CN201810516947.6A 2018-05-25 2018-05-25 Towed line array sonar subarray error mismatch estimation method Active CN108845325B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810516947.6A CN108845325B (en) 2018-05-25 2018-05-25 Towed line array sonar subarray error mismatch estimation method

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810516947.6A CN108845325B (en) 2018-05-25 2018-05-25 Towed line array sonar subarray error mismatch estimation method

Publications (2)

Publication Number Publication Date
CN108845325A CN108845325A (en) 2018-11-20
CN108845325B true CN108845325B (en) 2022-07-05

Family

ID=64213522

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810516947.6A Active CN108845325B (en) 2018-05-25 2018-05-25 Towed line array sonar subarray error mismatch estimation method

Country Status (1)

Country Link
CN (1) CN108845325B (en)

Families Citing this family (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109444896B (en) * 2018-11-21 2023-05-09 中国人民解放军91388部队 Underwater sound array positioning system and positioning method thereof
CN110008520B (en) * 2019-03-11 2022-05-17 暨南大学 Structural damage identification method based on displacement response covariance parameters and Bayesian fusion
CN111538058A (en) * 2020-04-23 2020-08-14 桂林长海发展有限责任公司 Passive positioning method, device and storage medium
CN111914641B (en) * 2020-06-30 2023-07-18 中国科学院声学研究所 Target depth identification method and system based on modal intensity matching analysis
CN112881981B (en) * 2021-01-11 2023-06-30 西北工业大学 Method for processing gain loss of wireless sensor array space under various mismatch conditions
CN113740804B (en) * 2021-08-27 2023-04-11 青岛理工大学 Hydrophone array direction finding system based on DSP and DOA estimation method thereof
CN113946955B (en) * 2021-10-14 2023-08-08 西安电子科技大学 Multi-target Bayesian direction-of-arrival estimation method based on fusion center feedback information
CN114397480B (en) * 2022-01-04 2022-10-14 湖南大学 Acoustic Doppler velocimeter error estimation method, device and system
CN114613384B (en) * 2022-03-14 2023-08-29 中国电子科技集团公司第十研究所 Deep learning-based multi-input voice signal beam forming information complementation method

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101644765A (en) * 2009-04-23 2010-02-10 中国科学院声学研究所 Amplitude and phase error correction method used for linear array of underwater acoustic transducer
CN103513238A (en) * 2012-06-15 2014-01-15 中国科学院声学研究所 A regularization least square subspace crossing target direction finding method

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101644765A (en) * 2009-04-23 2010-02-10 中国科学院声学研究所 Amplitude and phase error correction method used for linear array of underwater acoustic transducer
CN103513238A (en) * 2012-06-15 2014-01-15 中国科学院声学研究所 A regularization least square subspace crossing target direction finding method

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
Localization of low-frequency coherent sound sources with compressive beamforming-based passive synthetic aperture;Zhixiong Lei et al.;《Acoustical Society of America》;20150323;全文 *
失配状态下的双线阵波束形成研究;王华奎等;《声学学报(中文版)》;20070310(第02期);全文 *
拖线阵声纳的匹配场后置波束形成干扰抵消方法;杨坤德等;《西北工业大学学报》;20041030(第05期);全文 *
矢量水听器线列阵的被动合成孔径技术;何心怡 等;《武汉理工大学学报（交通科学与工程版）》;20031231;第27卷(第6期);全文 *

Also Published As

Publication number Publication date
CN108845325A (en) 2018-11-20

Similar Documents

Publication Publication Date Title
CN108845325B (en) Towed line array sonar subarray error mismatch estimation method
Bilik Spatial compressive sensing for direction-of-arrival estimation of multiple sources using dynamic sensor arrays
CN110197112B (en) Beam domain Root-MUSIC method based on covariance correction
RU2382379C2 (en) Method for multi-signal location of radio-frequency radiation sources at one frequency
CN109459744B (en) Steady self-adaptive beam forming method for realizing multi-interference suppression
CN107390197B (en) Radar self-adaption sum-difference beam angle measurement method based on feature space
CN109782238B (en) Combined calibration method for amplitude-phase response and array element position of sensor array element
CN111487599B (en) Bistatic MIMO radar angle estimation method and device under gain-phase error background
CN113189592B (en) Vehicle-mounted millimeter wave MIMO radar angle measurement method considering amplitude mutual coupling error
CN109541573B (en) Array element position calibration method for bending hydrophone array
CN112230215B (en) Self-adaptive monopulse direction finding method for unknown polarization parameters
CN109541526A (en) A kind of ring array direction estimation method using matrixing
CN109541572B (en) Subspace orientation estimation method based on linear environment noise model
CN111427045A (en) Underwater target backscattering imaging method based on distributed multi-input-multi-output sonar
CN115166725B (en) MUSIC algorithm-based meter wave polarization sensitive array radar height measurement method
Wang et al. Angle-polarization-range estimation using sparse polarization sensitive FDA-MIMO radar with co-prime frequency offsets
Tao et al. Oblique projection based enumeration of mixed noncoherent and coherent narrowband signals
CN114325560A (en) Super-resolution target direction finding method for beam scanning radar
CN113381793B (en) Coherent information source estimation-oriented non-grid direction-of-arrival estimation method
RU2385467C1 (en) Method for spatial polarisation-sensitive localisation of multibeam radio signals
CN113791379A (en) Orthogonal matching pursuit DOA estimation method under nested array non-Gaussian environment
Liang et al. A DOA Estimation algorithm for the vertical line array of vector hydrophone based on data fusion method
Daldorff et al. Direction of arrival estimation using single tripole radio antenna
Al Jabr et al. Modified UCA-ESPRIT for estimating DOA of coherent signals using one snapshot
Massoud et al. Angle of arrival estimation based on warped delay-and-sum (WDAS) beamforming technique

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant