CN112946565B - Interferometer direction finding fuzzy error correction method, system and medium for Kalman filtering - Google Patents

Interferometer direction finding fuzzy error correction method, system and medium for Kalman filtering Download PDF

Info

Publication number
CN112946565B
CN112946565B CN202110121691.0A CN202110121691A CN112946565B CN 112946565 B CN112946565 B CN 112946565B CN 202110121691 A CN202110121691 A CN 202110121691A CN 112946565 B CN112946565 B CN 112946565B
Authority
CN
China
Prior art keywords
value
phase difference
module
ambiguity
finding
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
CN202110121691.0A
Other languages
Chinese (zh)
Other versions
CN112946565A (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.)
Shanghai Institute Of Microwave Equipment 51st Research Institute Of China Electronics Technology Group Corp
Original Assignee
Shanghai Institute Of Microwave Equipment 51st Research Institute Of China Electronics Technology Group Corp
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 Shanghai Institute Of Microwave Equipment 51st Research Institute Of China Electronics Technology Group Corp filed Critical Shanghai Institute Of Microwave Equipment 51st Research Institute Of China Electronics Technology Group Corp
Priority to CN202110121691.0A priority Critical patent/CN112946565B/en
Publication of CN112946565A publication Critical patent/CN112946565A/en
Application granted granted Critical
Publication of CN112946565B publication Critical patent/CN112946565B/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
    • G01S3/00—Direction-finders for determining the direction from which infrasonic, sonic, ultrasonic, or electromagnetic waves, or particle emission, not having a directional significance, are being received
    • G01S3/02—Direction-finders for determining the direction from which infrasonic, sonic, ultrasonic, or electromagnetic waves, or particle emission, not having a directional significance, are being received using radio waves
    • G01S3/14—Systems for determining direction or deviation from predetermined direction
    • G01S3/46—Systems for determining direction or deviation from predetermined direction using antennas spaced apart and measuring phase or time difference between signals therefrom, i.e. path-difference systems
    • G01S3/48—Systems for determining direction or deviation from predetermined direction using antennas spaced apart and measuring phase or time difference between signals therefrom, i.e. path-difference systems the waves arriving at the antennas being continuous or intermittent and the phase difference of signals derived therefrom being measured
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Data Mining & Analysis (AREA)
  • Mathematical Analysis (AREA)
  • Theoretical Computer Science (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Computational Mathematics (AREA)
  • Databases & Information Systems (AREA)
  • Algebra (AREA)
  • General Engineering & Computer Science (AREA)
  • Software Systems (AREA)
  • Probability & Statistics with Applications (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Evolutionary Biology (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Remote Sensing (AREA)
  • Operations Research (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Computing Systems (AREA)
  • Radar Systems Or Details Thereof (AREA)

Abstract

The invention provides a Kalman filtering interferometer direction finding fuzzy error correction method, which is characterized by comprising the following steps: step 1: inputting the phase of each channel of the antenna by adopting an interferometer receiver, and calculating the phase difference between every two channels; calculating the maximum value of the phase difference ambiguity of the long base line; step 2: recording input serial numbers of array element receiving data, solving the modulus of phase difference between channels, resolving the ambiguity by adopting a traversal method to obtain the non-fuzzy phase difference, and solving a direction-finding value; and step 3: adopting Kalman filtering to the input direction-finding value; and 4, step 4: if the value of the output
Figure DDA0002922254280000011
Order to
Figure DDA0002922254280000012
Otherwise, the fuzzy value is obtained according to zeta (k)
Figure DDA0002922254280000013
Figure DDA0002922254280000014
If n (k-1) ═ n (k +1) and
Figure DDA0002922254280000015
order to
Figure DDA0002922254280000016
Order to
Figure DDA0002922254280000017
Let ρ be 1; updating
Figure DDA0002922254280000018
Turning to step 3; and 5: and outputting a direction finding result value. The interferometer direction finding fuzzy error correction method, system and medium for Kalman filtering can effectively correct angle error values caused by phase measurement errors.

Description

Interferometer direction finding fuzzy error correction method, system and medium for Kalman filtering
Technical Field
The invention relates to the field of electronic technology, in particular to a Kalman filtering interferometer direction finding fuzzy error correction method.
Background
In the process of using the interferometer to carry out direction finding on the target, ambiguity resolution errors of the interferometer can be caused due to phase measurement errors, and further direction finding errors of the target are caused. To obtain stable tracking of the target requires a stable output of the interferometer direction finding.
The existing method mostly discusses an interferometer ambiguity resolution method, because the phase error of the interferometer is inevitable under the non-ideal condition, the ambiguity resolution error is further caused, and no matter which ambiguity resolution method is adopted, the ambiguity resolution error caused by the phase error is inevitable, so the existing interferometer ambiguity resolution method can not be well applied to the real environment. The interferometer direction-finding antenna usually has a plurality of array elements, and the most basic method is to set the length of the shortest base line to a value smaller than half wavelength, so that the shortest base line is not blurred, and then the longest base line is adopted for direction finding. In real-world environments, especially high-band applications, it is sometimes not practical to set the antenna distance shorter than half a wavelength, because the high-frequency wavelength is very short, such as in the order of centimeters or even millimeters, which is much smaller than the physical size of the antenna, and thus this method cannot be well applied in real-world environments. In a low-frequency environment, even if the physical size requirement is met, the method cannot overcome the ambiguity error caused by the phase measurement error. The improved interferometer ambiguity resolution method is to adopt 3 array elements to perform ambiguity resolution, set the distance between the last two array elements as the distance between the first two array elements plus a specific value, and the value is less than half wavelength, and then to perform direction finding after the ambiguity resolution by using the method.
In view of the above-mentioned related art, the inventor considers that there is a problem of ambiguity resolution error, and therefore, a technical solution is needed to improve the above technical problem.
Disclosure of Invention
Aiming at the defects in the prior art, the invention aims to provide a Kalman filtering interferometer direction finding fuzzy error correction method, system and medium.
The invention provides a Kalman filtering interferometer direction finding fuzzy error correction method, which comprises the following steps:
step 1: interferometer receiver using array element number LArray elements are arranged in sequence from 1 to L, and the phase of each channel of the input antenna is
Figure BDA0002922254260000021
i is 1 … L, and the phase difference between every two channels is calculated
Figure BDA0002922254260000022
i＝1,…,L,j＝1,…,L,i<j; the distance between the array element i and the array element j is D ij Let ρ be 0; calculating the maximum value n of the long-base-line phase difference ambiguity according to the length of the antenna and the frequency f of the signal max ：
Figure BDA0002922254260000023
Wherein the content of the first and second substances,
Figure BDA0002922254260000024
representing rounding up, so the range of long baseline phase difference ambiguity is [ -n max ，n max ]，Φ 1L And D 1L Respectively representing the phase difference and distance between array element 1 and array element L, c representing the speed of light, theta max Represents a maximum angle measurement of the interferometer;
step 2: recording the input serial number k of the array element receiving data, and comparing the phase difference phi between the channels ij (k) Modulo 2 pi to obtain
Figure BDA0002922254260000025
The value range is (-pi, pi), the blur is resolved by a traversal method to obtain the blur-free phase difference
Figure BDA0002922254260000026
Finding direction-finding values
Figure BDA0002922254260000027
Will be provided with
Figure BDA0002922254260000028
And
Figure BDA0002922254260000029
are respectively marked as
Figure BDA00029222542600000210
And
Figure BDA00029222542600000211
the corresponding ambiguity is n, and is marked as n (k); (ii) a
And step 3: for inputted direction-finding value
Figure BDA00029222542600000212
Adopting Kalman filtering, outputting a result of ζ (k), if ρ is 0, turning to the step 4, otherwise, turning to the step 5;
and 4, step 4: if the value of the output
Figure BDA00029222542600000213
Order to
Figure BDA00029222542600000214
Otherwise, the fuzzy value is obtained according to zeta (k)
Figure BDA00029222542600000215
Figure BDA00029222542600000216
Wherein<·>Denotes rounding off, if n (k-1) ═ n (k +1) and
Figure BDA00029222542600000217
order to
Figure BDA00029222542600000218
Figure BDA00029222542600000219
Order to
Figure BDA00029222542600000220
Door with alpha, beta settingLimiting value, wherein the value range is 0-10, and rho is made to be 1; updating
Figure BDA00029222542600000221
Turning to step 3;
and 5: and outputting a direction finding result value.
Preferably, the step 2 comprises the steps of:
step 2.1-calculate error E ij (n)：
Let n be [ -n max ，n max ]Traverse the value in the range according to phi 1L Derived from the degree of ambiguity n
Figure BDA00029222542600000222
Error E of ij (n) and making it take on values between (-pi, pi),
Figure BDA00029222542600000223
i＝1,…,L,j＝1,…,L,i<j，
wherein, the 2π The expression is to take the modulus of 2 pi, and the value range is (-pi, pi);
step 2.2: calculating the value of the objective function E (n):
e (n) is an objective function, let E (n) be Σ i＜j |E ij (n)| 2 ；i＝1,…,L,j＝1,…,L,i<j；
Step 2.3: taking the value of n corresponding to the minimum value of E (n):
searching the obtained objective function E (n), and searching the n value which enables the objective function E (n) to be minimum, namely the ambiguity n of the long baseline phase difference;
step 2.4: calculating the phase difference of the long base line:
calculating to obtain the long-baseline non-fuzzy phase difference according to the ambiguity n of the long-baseline phase difference
Figure BDA0002922254260000031
Figure BDA0002922254260000032
Step 2.5: estimating a single pulse direction finding value:
estimation of signal direction of arrival
Figure BDA0002922254260000033
Comprises the following steps:
Figure BDA0002922254260000034
preferably, the step 3 comprises the steps of:
step 3.1: constructing a state estimator for the kth input
Figure BDA0002922254260000035
Wherein
Figure BDA0002922254260000036
Is an estimate of sequence number k, phi (k) is
Figure BDA0002922254260000037
State estimate of (c) (. 1) T The transpose is represented by,
Figure BDA0002922254260000038
obtained by the following formula:
Figure BDA0002922254260000039
k≥2，
x (k-1) is a state quantity obtained by the k-1 th input sequence, and if k is 1, setting
Figure BDA00029222542600000310
State transition matrix
Figure BDA00029222542600000311
T represents a time length, N (k-1) is noise, and k is derived from N (k-1) to N (0, Q (k)))>0, wherein the content of the amino acid is,
Figure BDA00029222542600000312
q 2 a mean square value representing the noise acceleration;
step 3.2: predicted covariance matrix of kth input
Figure BDA00029222542600000313
Expressed as:
Figure BDA00029222542600000314
k≥2，
wherein, P (k-1) represents the state covariance at the time of k-1, when k is 1, P (0) takes the random value as the initial covariance, and the kalman gain matrix is expressed as:
Figure BDA00029222542600000315
h (k) is the Jacobian matrix of measurement errors, R is the covariance matrix of measurement errors, Jacobian matrix:
Figure BDA00029222542600000316
step 3.3: the update equation for the target state is:
Figure BDA00029222542600000317
z (k) represents the measured value,
Figure BDA00029222542600000318
Figure BDA00029222542600000319
the representation predicted value is obtained from the state predicted value at the time k-1, and is taken
Figure BDA0002922254260000041
Step 3.4: the target state covariance is updated as:
Figure BDA0002922254260000042
wherein I represents an identity matrix;
step 3.5: the filtered direction finding value ζ (k) is x (2).
The invention also provides a Kalman filtering interferometer direction finding fuzzy error correction system, which comprises the following modules:
module M1: an interferometer receiver with an array element number of L is adopted, the array element arrangement is 1 to L in sequence, and the phase of each channel of an input antenna is
Figure BDA0002922254260000043
i is 1 … L, and the phase difference between every two channels is calculated
Figure BDA0002922254260000044
i＝1,…,L,j＝1,…,L,i<j; the distance between the array element i and the array element j is D ij Let ρ be 0; calculating the maximum value n of the long-baseline phase difference ambiguity according to the length of the antenna and the frequency f of the signal max ：
Figure BDA0002922254260000045
Wherein the content of the first and second substances,
Figure BDA0002922254260000046
representing rounding up, so the range of long baseline phase difference ambiguity is [ -n max ，n max ]，Φ 1L And D 1L Respectively representing the phase difference and distance between array element 1 and array element L, c representing the speed of light, theta max Represents a maximum angle measurement of the interferometer;
module M2: recording the input serial number k of the array element receiving data, and comparing the phase difference phi between the channels ij (k) Modulo 2 pi to obtain
Figure BDA0002922254260000047
The value range is (-pi, pi), the blur is resolved by a traversal method to obtain the blur-free phase difference
Figure BDA0002922254260000048
Finding the direction-finding value
Figure BDA0002922254260000049
Will be provided with
Figure BDA00029222542600000410
And
Figure BDA00029222542600000411
are respectively marked as
Figure BDA00029222542600000412
And
Figure BDA00029222542600000413
the corresponding ambiguity is n, and is marked as n (k); (ii) a
Module M3: for inputted direction-finding value
Figure BDA00029222542600000414
Adopting Kalman filtering, outputting a result of ζ (k), if ρ is 0, turning to a module M4, otherwise, turning to a module M5;
module M4: if the value of the output
Figure BDA00029222542600000415
Order to
Figure BDA00029222542600000416
Otherwise, the fuzzy value is obtained according to zeta (k)
Figure BDA00029222542600000417
Figure BDA00029222542600000418
Wherein<·>Denotes rounding off, if n (k-1) ═ n (k +1) and
Figure BDA00029222542600000419
order to
Figure BDA00029222542600000420
Figure BDA00029222542600000421
Order to
Figure BDA00029222542600000422
Alpha and beta represent set threshold values, the value range is 0-10, and rho is set to be 1; updating
Figure BDA00029222542600000423
Transitioning to module M3;
module M5: and outputting a direction finding result value.
Preferably, the module M2 includes the following modules:
module M2.1 calculating error E ij (n)：
Let n be [ -n max ，n max ]Traverse the value in the range according to phi 1L Derived from the degree of ambiguity n
Figure BDA00029222542600000424
Error E of ij (n) and making it take on values between (-pi, pi),
Figure BDA0002922254260000051
i＝1,…,L,j＝1,…,L,i<j，
wherein, () 2π The expression is to take the module of 2 pi, and the value range is (-pi, pi);
module M2.2: calculating the value of the objective function E (n):
e (n) is an objective function, and E (n) is ∑ i＜j |E ij (n)| 2 ；i＝1,…,L,j＝1,…,L,i<j；
Module M2.3: taking the value of n corresponding to the minimum value of E (n):
searching the obtained objective function E (n), and searching the n value which enables the objective function E (n) to be minimum, namely the ambiguity n of the long baseline phase difference;
module M2.4: calculating the phase difference of the long base line:
calculating to obtain the long-baseline non-fuzzy phase difference according to the ambiguity n of the long-baseline phase difference
Figure BDA0002922254260000052
Figure BDA0002922254260000053
Module M2.5: estimating a single pulse direction finding value:
estimation of signal direction of arrival
Figure BDA0002922254260000054
Comprises the following steps:
Figure BDA0002922254260000055
preferably, the module M3 includes the following modules:
module M3.1: constructing state estimators for the kth input
Figure BDA0002922254260000056
Wherein
Figure BDA0002922254260000057
Is an estimate of sequence number k, phi (k) is
Figure BDA0002922254260000058
State estimate of (c) (. 1) T The transpose is represented by,
Figure BDA0002922254260000059
obtained by the following formula:
Figure BDA00029222542600000510
k≥2，
x (k-1) is a state quantity obtained by the k-1 th input sequence, and if k is 1, setting
Figure BDA00029222542600000511
State transition matrix
Figure BDA00029222542600000512
T represents a time length, N (k-1) is noise, and k follows from N (k-1) to N (0, Q (k)))>0, wherein the content of the compound is,
Figure BDA00029222542600000513
q 2 a mean square value representing the noise acceleration;
module M3.2: predicted covariance matrix of kth input
Figure BDA00029222542600000514
Expressed as:
Figure BDA00029222542600000515
k≥2，
wherein, P (k-1) represents the state covariance at the time of k-1, when k is 1, P (0) takes the random value as the initial covariance, and the kalman gain matrix is expressed as:
Figure BDA0002922254260000061
h (k) is the Jacobian matrix of measurement errors, R is the covariance matrix of measurement errors, Jacobian matrix:
Figure BDA0002922254260000062
module M3.3: the update equation for the target state is:
Figure BDA0002922254260000063
z (k) represents the measured value,
Figure BDA0002922254260000064
Figure BDA0002922254260000065
the representation predicted value is obtained from the state predicted value at the time k-1, and is taken
Figure BDA0002922254260000066
Module M3.4: the target state covariance is updated as:
Figure BDA0002922254260000067
wherein I represents an identity matrix;
module M3.5: the filtered direction finding value ζ (k) is x (2).
The invention also provides a computer-readable storage medium having stored thereon a computer program which, when being executed by a processor, carries out the steps of the method as set forth above.
Compared with the prior art, the invention has the following beneficial effects:
when the direction-finding angle is filtered, a direct filtering result is not adopted, but a value with correct de-blurring is kept as far as possible, an error value is eliminated, and then an output result is filtered again. Compared with the traditional ambiguity-solving direction-finding method, the method can effectively correct the angle error value caused by the phase measurement error.
Drawings
Other features, objects and advantages of the invention will become more apparent upon reading of the detailed description of non-limiting embodiments with reference to the following drawings:
FIG. 1 is a schematic flow diagram of the process of the present invention;
FIG. 2 is a schematic view of a direction-finding scene used in the method of the present invention;
FIG. 3 is a comparison graph of the direction-finding angle obtained by the method of the present invention and the direction-finding angle and the true angle obtained by other methods;
FIG. 4 is a comparison graph of the direction-finding angle error obtained by the method of the present invention and the direction-finding angle error obtained by other methods.
Detailed Description
The present invention will be described in detail with reference to specific examples. The following examples will aid those skilled in the art in further understanding the present invention, but are not intended to limit the invention in any manner. It should be noted that it would be obvious to those skilled in the art that various changes and modifications can be made without departing from the spirit of the invention. All falling within the scope of the present invention.
The technical solution of the present invention is described in further detail below with reference to the accompanying drawings.
Referring to fig. 1 and 2, the present invention uses a scenario of an array receiving antenna, which has L array elements, and the array elements receive data simultaneously, and measure the direction of a radiation source by using phase difference. The interferometer direction finding fuzzy error correction method based on Kalman filtering comprises the following steps:
step 1: an interferometer receiver with an array element number of L is adopted, the array element arrangement is 1 to L in sequence, and the phase of each channel of an input antenna is
Figure BDA0002922254260000071
i is 1 … L, and the phase difference between every two channels is calculated
Figure BDA0002922254260000072
i＝1,…,L,j＝1,…,L,i<j; the distance between the array element i and the array element j is D ij Let ρ be 0; calculating the maximum value n of the long-baseline phase difference ambiguity according to the length of the antenna and the frequency f of the signal max ：
Figure BDA0002922254260000073
Wherein the content of the first and second substances,
Figure BDA0002922254260000074
representing rounding up, so the range of long baseline phase difference ambiguity is [ -n max ，n max ]，Φ 1L And D 1L Respectively representing the phase difference and distance between array element 1 and array element L, c representing the speed of light, theta max Representing the maximum angle measurement of the interferometer.
Step 2: recording the input serial number k of the array element receiving data, and comparing the phase difference phi between the channels ij (k) Modulo 2 pi to obtain
Figure BDA0002922254260000075
The value range is (-pi, pi), the blur is resolved by a traversal method to obtain the non-blur phase difference
Figure BDA0002922254260000076
Finding the direction-finding value
Figure BDA0002922254260000077
Will be provided with
Figure BDA0002922254260000078
And
Figure BDA0002922254260000079
are respectively marked as
Figure BDA00029222542600000710
And
Figure BDA00029222542600000711
the corresponding ambiguity is n, denoted as n (k).
Step 2.1: and (3) calculating an error: let n be [ -n max ，n max ]Traverse the value in the range according to phi 1L Derived from the degree of ambiguity n
Figure BDA00029222542600000712
Error E of ij (n) and making it take on values between (-pi, pi),
Figure BDA00029222542600000713
Figure BDA00029222542600000714
i＝1,…,L,j＝1,…,L,i<j, wherein, () 2π Representing the modulo of 2 pi, and the value range is (-pi, pi).
Step 2.2: calculating the value of the objective function E (n): e (n) is a target functionLet E (n) be Σ i＜j |E ij (n)| 2 ，i＝1,…,L,j＝1,…,L,i<j。
Step 2.3: taking the value of n corresponding to the minimum value of E (n): the obtained objective function e (n) is searched for the value of n that minimizes the objective function e (n), i.e., the ambiguity n of the long baseline phase difference.
Step 2.4: calculating the phase difference of the long base line: calculating to obtain the long-baseline non-fuzzy phase difference according to the ambiguity n of the long-baseline phase difference
Figure BDA00029222542600000715
Figure BDA00029222542600000716
Step 2.5: estimating a single pulse direction finding value: estimation of signal direction of arrival
Figure BDA00029222542600000717
Comprises the following steps:
Figure BDA00029222542600000718
and step 3: for inputted direction-finding value
Figure BDA00029222542600000719
And (5) adopting Kalman filtering, outputting the result of ζ (k), and if ρ is 0, turning to the step 4, otherwise, turning to the step 5.
Step 3.1: constructing state estimators for the kth input
Figure BDA00029222542600000720
Wherein
Figure BDA00029222542600000721
Is an estimate of sequence number k, phi (k) is
Figure BDA0002922254260000081
State estimate of (c) (. 1) T The transpose is represented by,
Figure BDA0002922254260000082
is obtained by the following formula,
Figure BDA0002922254260000083
Figure BDA0002922254260000084
k is more than or equal to 2, x (k-1) is the state quantity obtained by the k-1 input sequence, if k is 1, the setting is carried out
Figure BDA0002922254260000085
Figure BDA0002922254260000086
State transition matrix
Figure BDA0002922254260000087
T represents a time length, N (k-1) is noise, and k is derived from N (k-1) to N (0, Q (k)))>0, wherein the content of the compound is,
Figure BDA0002922254260000088
q 2 represents the mean square value of the noise acceleration.
Step 3.2: predicted covariance matrix of kth input
Figure BDA0002922254260000089
Expressed as:
Figure BDA00029222542600000810
k is more than or equal to 2, wherein, P (k-1) represents the state covariance at the k-1 moment, when k is 1, P (0) takes the random value as the initial covariance, the Kalman gain matrix is expressed as,
Figure BDA00029222542600000811
h (k) is the Jacobian matrix of measurement errors, R is the covariance matrix of measurement errors, Jacobian matrix:
Figure BDA00029222542600000812
step 3.3: the update equation of the target state is：
Figure BDA00029222542600000813
z (k) represents the measured value,
Figure BDA00029222542600000814
Figure BDA00029222542600000815
the representation predicted value is obtained from the state predicted value at the time k-1, and is taken
Figure BDA00029222542600000816
Step 3.4: the target state covariance is updated as,
Figure BDA00029222542600000817
wherein I represents an identity matrix.
Step 3.5: the filtered direction finding value ζ (k) is x (2).
And 4, step 4: if the value of the output
Figure BDA00029222542600000818
Order to
Figure BDA00029222542600000819
Otherwise, the fuzzy value is obtained according to zeta (k)
Figure BDA00029222542600000820
Figure BDA00029222542600000821
Wherein<·>Denotes rounding off, if n (k-1) ═ n (k +1) and
Figure BDA00029222542600000822
order to
Figure BDA00029222542600000823
Order to
Figure BDA00029222542600000824
Figure BDA00029222542600000825
Alpha and beta represent set threshold values, the value range is 0-10, and rho is set to be 1; updating
Figure BDA00029222542600000826
Figure BDA00029222542600000827
Go to step 3.
And 5, outputting a direction finding result value zeta (k).
The effect of the present invention can be further illustrated by the following simulation results:
simulation conditions are as follows: adopt 5 array element interferometers to carry out the direction finding, linear array arranges, and overall length 2m, array element 1 is to 5 intervals in proper order: 35cm, 36.4cm, 37.8 and 90.8 cm. The maximum error of the array element phase measurement is 30 degrees, and the working frequency point is 10 GHz.
Simulation content: the interferometer is used for direction finding, and the direction finding range is from-40 degrees to 40 degrees; respectively simulating an improved half-wavelength direction finding method, a traversal ambiguity-resolving method and a traversal ambiguity-resolving and Kalman error correction direction finding method.
And (3) simulation results: FIG. 3 is a graph showing the results of direction finding obtained by the method of the present invention. The true value is a value from-40 degrees to 40 degrees, and as can be seen from the figure, the method provided by the invention can obtain higher accuracy and has the best effect compared with the other two methods. The improved half-wavelength method only adopts 3 array elements to solve ambiguity, and when phase measurement errors exist, the phase ambiguity solving errors are more and the errors are larger. While the blur is solved by adopting 5 array elements in the traversal method, compared with the improved half-wavelength method, the blur solving error is relatively small, but the actual requirement cannot be met. And the ergodic method and the Kalman error correction method have high accuracy and can meet the practical application. Fig. 4 shows the results obtained by the simulation methods and the measurement errors of the real angles, and the results correspond to fig. 3, and it can be seen from the figure that the method provided by the present patent has small errors and can better meet the practical application.
Those skilled in the art will appreciate that, in addition to implementing the system and its various devices, modules, units provided by the present invention as pure computer readable program code, the system and its various devices, modules, units provided by the present invention can be fully implemented by logically programming method steps in the form of logic gates, switches, application specific integrated circuits, programmable logic controllers, embedded microcontrollers and the like. Therefore, the system and various devices, modules and units thereof provided by the invention can be regarded as a hardware component, and the devices, modules and units included in the system for realizing various functions can also be regarded as structures in the hardware component; means, modules, units for realizing various functions can also be regarded as structures in both software modules and hardware components for realizing the methods.
The foregoing description of specific embodiments of the present invention has been presented. It is to be understood that the present invention is not limited to the specific embodiments described above, and that various changes or modifications may be made by one skilled in the art within the scope of the appended claims without departing from the spirit of the invention. The embodiments and features of the embodiments of the present application may be combined with each other arbitrarily without conflict.

Claims (5)

1. A Kalman filtering interferometer direction finding fuzzy error correction method is characterized by comprising the following steps:
step 1: an interferometer receiver with an array element number of L is adopted, the array element arrangement is 1 to L in sequence, and the phase of each channel of an input antenna is
Figure FDA0003690532100000011
Calculating the phase difference between two channels
Figure FDA0003690532100000012
Figure FDA0003690532100000013
The distance between the array element i and the array element j is D ij Let ρ be 0; according to length of antenna and signalFrequency f, calculating the maximum value n of the phase difference ambiguity of the long base line max ：
Figure FDA0003690532100000014
Wherein, the first and the second end of the pipe are connected with each other,
Figure FDA0003690532100000015
representing rounding up, so the range of long baseline phase difference ambiguity is [ -n max ，n max ]，Φ 1L And D 1L Respectively representing the phase difference and distance between array element 1 and array element L, c representing the speed of light, theta max Represents a maximum angle measurement of the interferometer;
step 2: recording input serial number k of array element receiving data, and comparing phase difference phi between channels ij (k) Modulo 2 pi to obtain
Figure FDA0003690532100000016
The value range is (-pi, pi), the blur is resolved by a traversal method to obtain the blur-free phase difference
Figure FDA0003690532100000017
Finding the direction-finding value
Figure FDA0003690532100000018
Will be provided with
Figure FDA0003690532100000019
And
Figure FDA00036905321000000110
are respectively recorded as
Figure FDA00036905321000000111
And
Figure FDA00036905321000000112
corresponding ambiguity n, notedn(k)；
And step 3: for inputted direction-finding value
Figure FDA00036905321000000113
Adopting Kalman filtering, outputting a result of ζ (k), if ρ is 0, turning to the step 4, otherwise, turning to the step 5;
and 4, step 4: if the value of the output
Figure FDA00036905321000000114
Order to
Figure FDA00036905321000000115
Otherwise, the fuzzy value is obtained according to zeta (k)
Figure FDA00036905321000000116
Figure FDA00036905321000000117
Wherein<·>Denotes rounding off, if n (k-1) ═ n (k +1) and
Figure FDA00036905321000000118
order to
Figure FDA00036905321000000119
Figure FDA00036905321000000120
Order to
Figure FDA00036905321000000121
Alpha and beta represent set threshold values, the value range is 0-10, and rho is set to be 1; updating
Figure FDA00036905321000000122
Turning to step 3;
and 5: outputting a direction finding result value;
the step 2 comprises the following steps:
step 2.1-calculate error E ij (n)：
Let n be [ -n max ，n max ]Traverse the value in the range according to phi 1L Derived from the degree of ambiguity n
Figure FDA00036905321000000123
Error E of ij (n) and is set to a value between (-pi, pi),
Figure FDA00036905321000000124
wherein, the 2π The expression is to take the module of 2 pi, and the value range is (-pi, pi);
step 2.2: calculating the value of the objective function E (n):
e (n) is an objective function, and E (n) is ∑ i＜j |E ij (n)| 2 ；i＝1,…,L,j＝1,…,L,i＜j；
Step 2.3: taking the value of n corresponding to the minimum value of E (n):
searching the obtained objective function E (n), and searching the n value which enables the objective function E (n) to be minimum, namely the ambiguity n of the long baseline phase difference;
step 2.4: calculating the phase difference of the long base line:
calculating to obtain the long-baseline non-fuzzy phase difference according to the ambiguity n of the long-baseline phase difference
Figure FDA0003690532100000021
Figure FDA0003690532100000022
Step 2.5: estimating a single pulse direction finding value:
estimation of signal direction of arrival
Figure FDA0003690532100000023
Comprises the following steps:
Figure FDA0003690532100000024
2. the kalman filter interferometer direction finding fuzzy error correction method of claim 1, wherein the step 3 comprises the steps of:
step 3.1: constructing state estimators for the kth input
Figure FDA0003690532100000025
Wherein
Figure FDA0003690532100000026
Is an estimate of sequence number k, phi (k) is
Figure FDA0003690532100000027
State estimate of (c) (. 1) T The transpose is represented by,
Figure FDA0003690532100000028
obtained by the following formula:
Figure FDA0003690532100000029
x (k-1) is a state quantity obtained by the k-1 input sequence, and if k is 1, setting
Figure FDA00036905321000000210
State transition matrix
Figure FDA00036905321000000211
T represents a time length, N (k-1) is noise, subject to N (k-1) to N (0, Q (k)), k > 0, wherein,
Figure FDA00036905321000000212
q 2 a mean square value representing the noise acceleration;
step 3.2: predicted covariance matrix of kth input
Figure FDA00036905321000000213
Expressed as:
Figure FDA00036905321000000214
wherein, P (k-1) represents the state covariance at the time of k-1, when k is 1, P (0) takes the random value as the initial covariance, and the kalman gain matrix is expressed as:
Figure FDA00036905321000000215
h (k) is the Jacobian matrix of measurement errors, R is the covariance matrix of measurement errors, Jacobian matrix:
Figure FDA0003690532100000031
step 3.3: the update equation for the target state is:
Figure FDA0003690532100000032
z (k) represents a measured value of,
Figure FDA0003690532100000033
Figure FDA0003690532100000034
the representation predicted value is obtained from the state predicted value at the time k-1, and is taken
Figure FDA0003690532100000035
Step 3.4: the target state covariance is updated as:
Figure FDA0003690532100000036
wherein I represents an identity matrix;
step 3.5: the filtered direction finding value ζ (k) is x (2).
3. The system for correcting the direction finding fuzzy error of the interferometer based on Kalman filtering is characterized by comprising the following modules:
module M1: an interferometer receiver with an array element number of L is adopted, the array element arrangement is 1 to L in sequence, and the phase of each channel of an input antenna is
Figure FDA0003690532100000037
Calculating the phase difference between two channels
Figure FDA0003690532100000038
Figure FDA0003690532100000039
The distance between the array element i and the array element j is D ij Let ρ be 0; calculating the maximum value n of the long-baseline phase difference ambiguity according to the length of the antenna and the frequency f of the signal max ：
Figure FDA00036905321000000310
Wherein the content of the first and second substances,
Figure FDA00036905321000000311
representing rounding up, so the range of long baseline phase difference ambiguity is [ -n max ，n max ]，Φ 1L And D 1L Respectively representing the phase difference and the distance between the array element 1 and the array element L, c representing the speed of light, theta max Represents a maximum angle measurement of the interferometer;
module M2: recording the input serial number k of the array element receiving data, and comparing the phase difference phi between the channels ij (k) Modulo 2 pi to obtain
Figure FDA00036905321000000312
The value range is (-pi, pi), the blur is resolved by a traversal method to obtain the blur-free phase difference
Figure FDA00036905321000000313
Finding the direction-finding value
Figure FDA00036905321000000314
Will be provided with
Figure FDA00036905321000000315
And
Figure FDA00036905321000000316
are respectively marked as
Figure FDA00036905321000000317
And
Figure FDA00036905321000000318
the corresponding ambiguity is n, and is marked as n (k);
module M3: for inputted direction-finding value
Figure FDA00036905321000000319
Adopting Kalman filtering, outputting a result of ζ (k), if ρ is 0, turning to a module M4, otherwise, turning to a module M5;
module M4: if the value of the output
Figure FDA00036905321000000320
Order to
Figure FDA00036905321000000321
Otherwise, the fuzzy value is obtained according to zeta (k)
Figure FDA00036905321000000322
Figure FDA00036905321000000323
Wherein<·>Denotes rounding off, if n (k-1) ═ n (k +1) and
Figure FDA00036905321000000324
order to
Figure FDA00036905321000000325
Figure FDA0003690532100000041
Order to
Figure FDA0003690532100000042
Alpha and beta represent set threshold values, the value range is 0-10, and rho is set to be 1; updating
Figure FDA0003690532100000043
Transitioning to module M3;
module M5: outputting a direction finding result value;
the module M2 includes the following modules:
module M2.1 calculating error E ij (n)：
Let n be [ -n max ，n max ]Traverse the value in the range according to phi 1L Derived from the degree of ambiguity n
Figure FDA0003690532100000044
Error E of ij (n) and making it take on values between (-pi, pi),
Figure FDA0003690532100000045
wherein, (.) 2π The expression is to take the module of 2 pi, and the value range is (-pi, pi);
module M2.2: calculating the value of the objective function E (n):
e (n) is an objective function, and E (n) is ∑ i＜j |E ij (n)| 2 ；i＝1,…,L,j＝1,…,L,i＜j；
Module M2.3: taking the value of n corresponding to the minimum value of E (n):
searching the obtained objective function E (n), and searching the n value which enables the objective function E (n) to be minimum, namely the ambiguity n of the long baseline phase difference;
module M2.4: calculating the phase difference of the long base line:
calculating to obtain the long-baseline non-fuzzy phase difference according to the ambiguity n of the long-baseline phase difference
Figure FDA0003690532100000046
Figure FDA0003690532100000047
Module M2.5: estimating a single pulse direction finding value:
estimation of signal direction of arrival
Figure FDA0003690532100000048
Comprises the following steps:
Figure FDA0003690532100000049
4. the kalman filter interferometer direction finding blur correction system of claim 3, wherein the module M3 comprises the following modules:
module M3.1: constructing state estimators for the kth input
Figure FDA00036905321000000410
Wherein
Figure FDA00036905321000000411
Is an estimate of sequence number k, phi (k) is
Figure FDA00036905321000000412
State estimate of (c) (. 1) T Which represents a transposition of the image,
Figure FDA00036905321000000413
obtained by the following formula:
Figure FDA00036905321000000414
x (k-1) is a state quantity obtained by the k-1 th input sequence, and if k is 1, setting
Figure FDA00036905321000000415
State transition matrix
Figure FDA00036905321000000416
T represents a time length, N (k-1) is noise, subject to N (k-1) to N (0, Q (k)), k > 0, wherein,
Figure FDA0003690532100000051
q 2 a mean square value representing the noise acceleration;
module M3.2: predicted covariance matrix of kth input
Figure FDA0003690532100000052
Expressed as:
Figure FDA0003690532100000053
wherein, P (k-1) represents the state covariance at time k-1, when k equals 1, P (0) takes the random value as the initial covariance, and the kalman gain matrix is represented as:
Figure FDA0003690532100000054
h (k) is the Jacobian matrix of measurement errors, R is the covariance matrix of measurement errors, Jacobian matrix:
Figure FDA0003690532100000055
module M3.3: the update equation for the target state is:
Figure FDA0003690532100000056
z (k) represents the measured value,
Figure FDA0003690532100000057
Figure FDA0003690532100000058
the representation predicted value is obtained from the state predicted value at the time k-1, and is taken
Figure FDA0003690532100000059
Module M3.4: the target state covariance is updated as:
Figure FDA00036905321000000510
wherein I represents an identity matrix;
module M3.5: the filtered direction finding value ζ (k) is x (2).
5. A computer-readable storage medium, in which a computer program is stored which, when being executed by a processor, carries out the steps of the method of any one of claims 1-2.
CN202110121691.0A 2021-01-28 2021-01-28 Interferometer direction finding fuzzy error correction method, system and medium for Kalman filtering Active CN112946565B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110121691.0A CN112946565B (en) 2021-01-28 2021-01-28 Interferometer direction finding fuzzy error correction method, system and medium for Kalman filtering

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110121691.0A CN112946565B (en) 2021-01-28 2021-01-28 Interferometer direction finding fuzzy error correction method, system and medium for Kalman filtering

Publications (2)

Publication Number Publication Date
CN112946565A CN112946565A (en) 2021-06-11
CN112946565B true CN112946565B (en) 2022-09-13

Family

ID=76239035

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110121691.0A Active CN112946565B (en) 2021-01-28 2021-01-28 Interferometer direction finding fuzzy error correction method, system and medium for Kalman filtering

Country Status (1)

Country Link
CN (1) CN112946565B (en)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114839587B (en) * 2022-03-25 2023-09-05 中国电子科技集团公司第二十九研究所 External correction method for interferometer system
CN115166631A (en) * 2022-08-02 2022-10-11 成都美数科技有限公司 Direction finding method and system based on interferometer direction finding receiver

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104330767A (en) * 2014-11-27 2015-02-04 中国船舶重工集团公司第七二四研究所 Interferometer direction-finding method based on phase fuzzy number search and least square fit
KR20150078539A (en) * 2013-12-31 2015-07-08 광운대학교 산학협력단 Fuzzy Prediction Error Compensation Method Control for Mobile OIS Control
CN106500588A (en) * 2016-11-18 2017-03-15 烟台职业学院 A kind of phase-interferometer inter-channel phase difference noise covariance method of estimation
CN108802670A (en) * 2018-06-05 2018-11-13 北京理工大学 A kind of phase interference angle-measuring method of robust
KR102093743B1 (en) * 2018-11-06 2020-03-26 경북대학교 산학협력단 System for lane level positioning location information of ground vehicle using sensor fusion
CN111381264A (en) * 2018-12-27 2020-07-07 北京六分科技有限公司 Long baseline ambiguity fixing method and platform in network RTK

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106842268B (en) * 2017-03-13 2020-04-14 惠州市组合科技有限公司 double-GNSS receiver carrier phase double-difference integer ambiguity floating point solution vector estimation method
CN107064980B (en) * 2017-03-24 2020-04-21 和芯星通科技(北京)有限公司 Carrier phase ambiguity fixing method and device and satellite navigation receiver

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR20150078539A (en) * 2013-12-31 2015-07-08 광운대학교 산학협력단 Fuzzy Prediction Error Compensation Method Control for Mobile OIS Control
CN104330767A (en) * 2014-11-27 2015-02-04 中国船舶重工集团公司第七二四研究所 Interferometer direction-finding method based on phase fuzzy number search and least square fit
CN106500588A (en) * 2016-11-18 2017-03-15 烟台职业学院 A kind of phase-interferometer inter-channel phase difference noise covariance method of estimation
CN108802670A (en) * 2018-06-05 2018-11-13 北京理工大学 A kind of phase interference angle-measuring method of robust
KR102093743B1 (en) * 2018-11-06 2020-03-26 경북대학교 산학협력단 System for lane level positioning location information of ground vehicle using sensor fusion
CN111381264A (en) * 2018-12-27 2020-07-07 北京六分科技有限公司 Long baseline ambiguity fixing method and platform in network RTK

Non-Patent Citations (9)

* Cited by examiner, † Cited by third party
Title
Application of neural network algorithm in fault diagnosis of mechanical intelligence;Xianzhen Xu;《Mechanical Systems and Signal Processing》;20201231;全文 *
Assessing Probability of Correct Ambiguity Resolution in the Presence of Time Correlated Errors;Kyle OKeefe;《Navigation》;20061231;全文 *
一种新的干涉测角数据处理算法;朱新国等;《北京理工大学学报》;20080615(第06期);全文 *
基于多历元递推最小二乘卡尔曼滤波方法的模糊度解算;孙红星等;《武汉大学学报(信息科学版)》;20080705(第07期);全文 *
基于数字干涉仪的机载ESM系统实现方法及测向误差分析;何晓明等;《舰船电子对抗》;20131225(第06期);全文 *
干涉仪阵列解模糊算法研究;夏韶俊等;《无线电工程》;20130805(第08期);全文 *
相差变化率的无模糊多通道检测及在单站无源定位中的应用;郁涛;《海峡科技与产业》;20160515(第05期);全文 *
相差变化率的无模糊测量;郁涛;《电子科学技术》;20140910(第02期);全文 *
运动单站干涉仪相位差直接定位方法;张敏等;《航空学报》;20130528(第09期);全文 *

Also Published As

Publication number Publication date
CN112946565A (en) 2021-06-11

Similar Documents

Publication Publication Date Title
CN112946565B (en) Interferometer direction finding fuzzy error correction method, system and medium for Kalman filtering
Forster et al. On-manifold preintegration for real-time visual--inertial odometry
CN106352876B (en) A kind of airborne distribution POS Transfer Alignments based on H ∞ and CKF mixed filterings
Wu et al. A Square Root Inverse Filter for Efficient Vision-aided Inertial Navigation on Mobile Devices.
US9098229B2 (en) Single image pose estimation of image capture devices
CN110702095B (en) Data-driven high-precision integrated navigation data fusion method
US20160305784A1 (en) Iterative kalman smoother for robust 3d localization for vision-aided inertial navigation
Öllerer et al. Robust high-dimensional precision matrix estimation
CN106802416B (en) Fast factorization back projection SAR self-focusing method
Ling et al. Modeling varying camera-imu time offset in optimization-based visual-inertial odometry
Klenske et al. Nonparametric dynamics estimation for time periodic systems
Begelfor et al. How to put probabilities on homographies
Brouri et al. Identification of series–parallel systems composed of linear and nonlinear blocks
CN112765548A (en) Covariance determination method, positioning method and device for sensor fusion positioning
Masnadi-Shirazi et al. A Step by Step Mathematical Derivation and Tutorial on Kalman Filters
CN116047459B (en) Array radar echo signal recovery method and related equipment in pulse interference environment
CN113390421B (en) Unmanned aerial vehicle positioning method and device based on Kalman filtering
CN111435168B (en) Positioning method and device
CN109840069B (en) Improved self-adaptive fast iterative convergence solution method and system
Peng et al. Ultrafast square-root filter-based VINS
Zhang et al. Maximum correntropy unbiased minimum-variance filter
CN109581284A (en) Non-market value removing method based on Interactive Multiple-Model
JP7415090B1 (en) Target state estimation device and target state estimation method
CN111093265B (en) Cooperative positioning method and device based on angle of arrival ranging
CN110988790B (en) Passive target positioning method and device

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