Three-coordinate air radar relative system error accumulation estimation method

Abstract

The invention relates to a three-coordinate air radar relative system error accumulation estimation method, and belongs to the technical field of multi-radar data fusion. According to the invention, through selecting observation data of main and auxiliary station radars for a typical air route target, a relative system error joint resolving model originating from radar asynchronous measurement data is established through three-dimensional unified rectangular coordinate conversion, three-dimensional linear track line parameter estimation, time registration point calculation and the like; and the reliability of a relative system error accumulation estimation result is ensured through operations such as singular value rejection. Experimental results show that follow-up measurement of the auxiliary station radar is corrected according to the cumulative estimation result, the split degree of the observed track is reduced by nearly one order of magnitude, and the error correction effect is very ideal. The estimation method provided by the invention is scientific, reasonable in scheme implementation steps and very high in operability and practicability.

Description

Technical Field

[0001] The present invention belongs to the technical field of multi-radar data fusion, and particularly relates to a method for cumulatively estimating the relative systematic error of a three-coordinate air surveillance radar system. Background Art

[0002] Three-coordinate air surveillance radars are generally deployed on the ground. By mechanically scanning in the horizontal direction (azimuth angle) and electronically scanning in the vertical direction (elevation angle), they obtain measurement data of the distance, azimuth, and altitude of airborne targets. The relative systematic error of a three-coordinate radar refers to the systematic error relative to a specified reference system when the radar measures the distance, azimuth angle, and elevation angle of a target. When the relative systematic errors of multiple radars are large, it will cause spatial splitting of the observation results of the same target. In severe cases, it will also prevent track association and fusion of the same target, and further cause incorrect association of tracks corresponding to different targets.

[0003] Generally, when we mention that there is a systematic error in a measurement, it is relative to a certain reference system. If the absolute position of a specified target is used as the reference system, the absolute systematic error is obtained. In multi-radar networked observations, the absolute position of the target is mostly unknown in most cases. What we can easily obtain are the discrete observation values of different radars for the same airborne moving target. Traditional data calibration methods all attempt to estimate the absolute systematic error existing in radar measurements through such measurement data. However, simulation tests show that the target track after systematic error correction depends on the relative positions of the single-radar measurement tracks participating in the calculation, and has little correlation with the true track position. In practical engineering applications, the estimated result of the absolute systematic error of radar A will vary depending on whether the measurement data of radar B or radar C is used. If it is recognized that the systematic error calculated from the data of A and B is the absolute systematic error of radar A, contradictions may occur, such as the reduction of the corrected track distance between A and B and the increase of the corrected track distance between A and C. Therefore, the idea of obtaining the absolute systematic error is difficult to grasp in engineering practice.

[0004] The relative systematic error in radar measurement is directly reflected as the spatial splitting of the observation track lines of multiple radars for the same target in a three-dimensional unified rectangular coordinate system. The error component factors causing this spatial splitting mainly include measurement (range, azimuth, elevation) errors, clock biases, site location (geographical coordinates, elevation, and true north) errors, coordinate transformation errors, time registration errors, and the inherent biases of the radar measurement caused by the site environment, etc. Estimating and eliminating the relative systematic error is an important means to reduce the spatial splitting of the target observation results. However, in engineering practice, not all systematic error components must be eliminated. Because some systematic error components have complex formation factors and are very difficult to estimate; some error components are very small, and whether to correct them has little impact on the total error. For example, the coordinate transformation errors brought about by different map projection methods and the inherent biases of the radar measurement caused by the site environment, etc., cause fixed translations in the X, Y, and Z directions for each single radar measurement point in the three-dimensional unified rectangular coordinate system except for the random errors. The formation factors of this error component are complex and it is very difficult to estimate; when using GPS or Beidou devices to accurately locate and measure the height of the radar site and adopting an effective projection model, the site location error can be ignored. Therefore, determining reasonable relative systematic error components of the three-coordinate radar and a scientific joint solution model for systematic errors is an urgent problem to be solved in engineering practice.

[0005] In most cases, the radar network measurement system can only provide the observation data of multiple radars for the same target in the same time period, and usually a set of measurement values (range, azimuth, height) on a typical flight path (the target flies in a straight line at a certain altitude) can be obtained. Based on such a data environment, we propose: ① Relative to other radars in the network, the observation of a certain radar is accurate. At this time, the measurement values of this radar (named the master station) can be used as the true description of the target position, and other radars (named slave stations) are referenced to this, so as to obtain the relative systematic error of the slave station radar relative to the master station radar. For a regional radar network, other radars can be corrected with the master station as the standard. In this way, while achieving the consistency of the observation results of each radar, the complexity of the estimation method is simplified and it is convenient for engineering implementation. ② Regarding the range, azimuth, and elevation measurement errors as the main components of the systematic error, a joint solution model for relative systematic errors derived from asynchronous radar measurement data is established. ③ A cumulative estimation method for the relative systematic error of the three-coordinate air defense radar based on the joint solution results of multiple relative systematic errors is proposed, and the influence of abnormal error estimation results is shielded by means of historical experience data. Summary of the Invention

[0006] (1) Technical Problems to be Solved

[0007] The technical problem to be solved by the present invention is how to provide a cumulative estimation method for the relative systematic error of a three-coordinate air defense radar to solve the problem of estimating the relative systematic error in radar measurement.

[0008] (2) Technical solution

[0009] To solve the above technical problems, the present invention proposes a method for cumulatively estimating the relative systematic error of a three-coordinate air radar, and the method includes the following steps:

[0010] Step 1: Initialize the distance relative systematic error list ER, the azimuth relative systematic error list EA, and the elevation relative systematic error list EP of the secondary station radar B;

[0011] Step 2: Select the observation data of a straight flight track of the same air target by the master station radar A and the secondary station radar B, {{(t SAi , R Ai , θ Ai , h Ai ), i = 1, 2,... n}} and {{(t SBj , R Bj , θ Bj , h Bj ), j = 1, 2,... m}}, where n represents the number of observation data points of the master station radar A, m represents the number of observation data points of the secondary station radar B, and (t SAi , R Ai , θ Ai , h Ai ) represents the target distance R SAi , azimuth θ Ai , and altitude h Ai measured by the master station radar A at time t Ai , and (t SBj , R Bj , θ Bj , h Bj ) represents the target distance R SBj , azimuth θ Bj , and altitude h Bj measured by the secondary station radar B at time t Bj ;

[0012] Step 3: Perform coordinate transformation on the observation data of the master station radar A to obtain a set of corresponding three-dimensional unified rectangular coordinates {{(t SAi , x SAi , y SAi , z SAi ), i = 1, 2,... n}};

[0013] Step 4: Perform coordinate transformation on the observation data of the secondary station radar B to obtain a set of corresponding three-dimensional unified rectangular coordinates {{(t SBj , x SBj , y SBj , z SBj ), j = 1, 2,... m}};

[0014] Step 5: Use the linear track line parameter estimation model in the three-dimensional unified rectangular coordinate system to perform three-dimensional linear parameter estimation on the main and auxiliary station radar observation data respectively, and obtain the main station radar A observation track line parameters and the auxiliary station radar B observation track line parameters;

[0015] Step 6: Calculate the time registration points of the main station radar A and the auxiliary station radar B according to the observation point time of the auxiliary station radar B; Convert the time registration points of the auxiliary station radar B into three-dimensional polar coordinate registration data centered on the radar station site, expressed as {(R′ Bi , θ′ Bi , h′ Bi ), i = 1, 2, … m};

[0016] Step 7: Set the distance, azimuth, and pitch relative system errors of the auxiliary station radar B as (ΔR, Δθ, Δα) respectively. Use the time registration point data of the main station radar A and the auxiliary station radar B to construct a relative system error joint solution system of equations with (ΔR, Δθ, Δα) as unknowns, and use the quasi-Newton method for finding a set of real roots of the nonlinear system of equations to solve this system of equations; If the system of equations has no solution, go to Step 2, select other linear track line observation data of the main station radar A and the auxiliary station radar B for the same airborne target, and start a new round of system error estimation; If the system of equations has a solution, assume the solution of the system of equations is (Dr, Da, Dp);

[0017] Step 8: Perform singular value judgment on Dr, Da, and Dp respectively, incorporate the non-singular values into the relative system error lists ER, EA, and EP of the auxiliary station radar B, and organize the error lists;

[0018] Step 9: If the list length is greater than or equal to the basic statistical number PN, calculate the mean value of the list elements and output it as the corresponding relative system error cumulative estimation result; Otherwise, it is considered that the corresponding relative system error estimation condition is not yet available;

[0019] Step 10: Go to Step 2, select other linear track line observation data of the main station radar A and the auxiliary station radar B for the same airborne target, and start a new round of system error estimation.

[0020] (III) Beneficial Effects

[0021] The present invention proposes a method for cumulative estimation of relative systematic error of a three-coordinate air surveillance radar. The present invention adopts the idea of relative systematic error. By selecting the observation data of the main and auxiliary station radars for typical route targets, through three-dimensional unified rectangular coordinate transformation, three-dimensional linear track line parameter estimation, time registration point calculation, etc., a joint solution model of relative systematic error derived from asynchronous radar measurement data is established, and the reliability of the cumulative estimation result of relative systematic error is ensured through operations such as singular value elimination. The experimental results show that by correcting the subsequent measurements of the auxiliary station radar with this cumulative estimation result, the degree of track splitting of the observation track is reduced by nearly one order of magnitude, and the error correction effect is very ideal, meeting the requirements for the accuracy of air surveillance radar data fusion in engineering practice. The estimation method provided by the present invention is scientific, the implementation steps of the scheme are reasonable, and the operability and practicability are very strong. Description of the Drawings

[0022] Figure 1 It is the main implementation steps of the relative systematic error cumulative estimation method of the technical solution of the present invention;

[0023] Figure 2 It is the XY-plane projection schematic diagram of the first simulated observation track line in the embodiment of the present invention;

[0024] Figure 3 It is the function flow chart of solving the linear algebraic equation system by the full-selection principal element Gaussian elimination method;

[0025] Figure 4 It is the 100-time simulated calculation result of the range relative systematic error in the embodiment of the present invention;

[0026] Figure 5 It is the 100-time simulated calculation result of the azimuth relative systematic error in the embodiment of the present invention;

[0027] Figure 6 It is the 100-time simulated calculation result of the elevation relative systematic error in the embodiment of the present invention;

[0028] Figure 7 It is the XY-plane projection schematic diagram of the first simulated observation track line in the verification experiment of the present invention;

[0029] Figure 8 It is the projection on the XY plane of the observation track lines of the main and auxiliary station radars and the auxiliary station correction track points generated in the first simulation after correction in the verification experiment of the present invention;

[0030] Figure 9 It is the comparison result of the change of the observed target height value of the auxiliary station radar generated in the first simulation before and after correction and the simulated height true value in the verification experiment of the present invention;

[0031] Figure 10Calculation results of the correction rate of the height value after the first simulation of the auxiliary station measurement correction relative to the simulated true value in the verification experiment of the present invention;

[0032] Figure 11 Average distance between the auxiliary station track line and the main station observation track line before and after the correction in the verification experiment in the embodiment of the present invention;

[0033] Figure 12 Correction rate of the average distance between the auxiliary station track line and the main station observation track line before and after the correction in the verification experiment in the embodiment of the present invention. Detailed implementation manners

[0034] To make the objectives, contents and advantages of the present invention clearer, the following further describes in detail the specific implementation manners of the present invention with reference to the drawings and embodiments.

[0035] The objective of the present invention is to propose a normalized radar relative system error cumulative estimation method applicable to a multi-radar asynchronous measurement environment, so as to achieve the best superposition of the observation results of multiple three-coordinate air defense radars on the same target in a three-dimensional unified rectangular coordinate system.

[0036] To achieve the above objective, the present invention provides a three-coordinate air defense radar relative system error cumulative estimation method, which is applied to the pre-data preprocessing process of a multi-radar data fusion system. The system error estimation result obtained by applying this method corrects subsequent measurements, and has an obvious promoting effect on improving the correct rate of multi-radar target association and recognition and enhancing the accuracy of target state estimation.

[0037] The method includes the following steps:

[0038] Step 1: Initialize the range relative system error list ER, azimuth relative system error list EA, and elevation relative system error list EP of the auxiliary station radar B.

[0039] Step 2: Select the observation data of a straight track line of the same air target by the main station radar A and the auxiliary station radar B, {{(t SAi ,R Ai ,θ Ai ,h Ai ), i = 1, 2,... n} and {{(t SBj ,R Bj ,θ Bj ,h Bj ), j = 1, 2,... m}, where n represents the number of observation data points of the main station radar A, m represents the number of observation data points of the auxiliary station radar B, and (t SAi ,R Ai ,θ Ai ,h Ai ) represents the target range R measured by the main station radar A at time t SAi ;Ai , azimuth θ Ai and altitude h Ai , (t SBj , R Bj , θ Bj , h Bj ) represents the target distance R SBj measured by the slave station radar B at time t Bj , azimuth θ Bj and altitude h Bj .

[0040] Step 3: Perform coordinate transformation on the observation data of the master station radar A to obtain a set of corresponding three-dimensional unified rectangular coordinates {(t SAi , x SAi , y SAi , z SAi ), i = 1, 2, … n}.

[0041] Step 4: Perform coordinate transformation on the observation data of the slave station radar B to obtain a set of corresponding three-dimensional unified rectangular coordinates {(t SBj , x SBj , y SBj , z SBj ), j = 1, 2, … m}.

[0042] Step 5: Use the straight-line track parameter estimation model in the three-dimensional unified rectangular coordinate system to perform three-dimensional straight-line parameter estimation on the observation data of the master and slave station radars respectively, and obtain the observation track line parameters of the master station radar A and the slave station radar B.

[0043] Step 6: Calculate the time registration points of the master station radar A and the slave station radar B according to the time of the observation points of the slave station radar B. Convert the time registration points of the slave station radar B into three-dimensional polar coordinate registration data centered on the radar station site, expressed as {(R′ Bi , θ′ Bi , h′ Bi ), i = 1, 2, … m}.

[0044] Step 7: Set the relative system errors of the distance, azimuth and pitch of the slave station radar B as (ΔR, Δθ, Δα) respectively. Use the time registration point data of the master station radar A and the slave station radar B to construct a joint solution system of equations for the relative system errors with (ΔR, Δθ, Δα) as unknowns, and use the quasi-Newton method for finding a set of real roots of the nonlinear equations to solve this system of equations. If the system of equations has no solution, go to Step 2, select the straight-line track line observation data of the master station radar A and the slave station radar B for another same airborne target, and start a new round of system error estimation. If the system of equations has a solution, assume the solution of the system of equations is (Dr, Da, Dp).

[0045] Step 8: Perform singular value judgment on Dr, Da, and Dp respectively, incorporate the non-singular values into the relative system error lists ER, EA, and EP of the slave radar B, and organize the error lists.

[0046] Step 9: If the list length is greater than or equal to the basic statistical number PN, calculate the mean of the list elements and output it as the corresponding relative system error cumulative estimation result. Otherwise, it is considered that the corresponding relative system error estimation condition is not yet available.

[0047] Step 10: Go to Step 2, select the observation data of the straight flight track of the same air target by the master radar A and the slave radar B again, and start a new round of system error estimation.

[0048] The said Step 1 includes the following steps:

[0049] Step 1.1: Create a list ER with the element type of real numbers to store the historical distance relative system error estimations of the slave radar B.

[0050] Step 1.2: Create a list EA with the element type of real numbers to store the historical azimuth relative system error estimations of the slave radar B.

[0051] Step 1.3: Create a list EP with the element type of real numbers to store the historical pitch relative system error estimations of the slave radar B.

[0052] The said Step 2 includes the following steps:

[0053] Step 2.1: Select the target observation data reported by the master and slave radars simultaneously when the air target is on a straight flight track. The number of observation data for each radar is generally not less than 10 points. "Simultaneous period" means that the time difference between the first points of the observation data of the master and slave radars is not greater than 1 radar detection period T (usually 10 or 20 seconds), and the time difference between the last points is also not greater than T.

[0054] Step 2.2: The observation data of the selected master radar A is: {(t SAi ,R Ai ,θ Ai ,h Ai ), i = 1, 2,... n}. Among them, (t SAi ,R Ai ,θ Ai ,h Ai ) represents the target distance R SAi measured by the master radar A at time t Ai , azimuth θ Ai and altitude h Ai , and n represents the number of observation data points of the master radar A.

[0055] Step 2.3: The observed data of the selected secondary radar B are: {(t SBj , R Bj , θ Bj , h Bj ), j = 1, 2, … m}. Where, (t SBj , R Bj , θ Bj , h Bj ) represents the target distance R SBj , azimuth θ Bj , and altitude h Bj measured by the secondary radar B at time t Bj , and m represents the number of observed data points of the secondary radar B. And |t SB1 - t SA1 | ≤ T, |t SBm - t SAn | ≤ T.

[0056] The said Step 3 includes the following steps:

[0057] Step 3.1: Convert the observed data {(t SAi , R Ai , θ Ai , h Ai ), i = 1, 2, … n} of the master radar A into three-dimensional rectangular coordinates {(t SAi , x Ai , y Ai , z Ai ), i = 1, 2, … n} centered on this station. Where:

[0058]

[0059] Step 3.2: Convert {(t SAi , x Ai , y Ai , z Ai ), i = 1, 2, … n} into three-dimensional unified rectangular coordinates {(t SAi , x SAi , y SAi , z SAi ), i = 1, 2, … n}.

[0060] Where:

[0061]

[0062] (X SA , Y SA , Z SA ) represents the three-dimensional rectangular coordinates of the master radar A in the unified coordinate system, that is, the station location coordinates of the master radar A.

[0063] Step 4 includes the following steps:

[0064] Step 4.1: Convert the observation data \(\{(t SBj ,R Bj ,\theta Bj ,h Bj ), j = 1, 2, \ldots m\}\) of the secondary station radar B into three-dimensional rectangular coordinates \(\{(t SBj ,x Bj ,y Bj ,z Bj ), j = 1, 2, \ldots m\}\) centered at the local station. Wherein:

[0065]

[0066] Step 4.2: Convert \(\{(t SBj ,x Bj ,y Bj ,z Bj ), j = 1, 2, \ldots m\}\) into three-dimensional unified rectangular coordinates \(\{(t SBj ,x SBj ,y SBj ,z SBj ), j = 1, 2, \ldots m\}\). Wherein:

[0067]

[0068] (X SB ,Y SB ,Z SB ) represents the three-dimensional rectangular coordinates of the secondary station radar B in the unified coordinate system, that is, the site coordinates of the secondary station radar B.

[0069] Step 5 includes the following steps:

[0070] Step 5.1: Use the linear track line parameter estimation model to perform three-dimensional linear parameter estimation on the observation data of the master station radar A, and obtain the observation track line parameters \((k AX ,d AX ), (k AY ,d AY ), (k AZ ,d AZ ) of the master station radar A, including the following steps:

[0071] Step 5.1.1: Abbreviate the X-axis observation data \(\{(t SAi ,x SAi ), i = 1, 2, \ldots n\}\) of the master station radar A as: \(\{(x i ,y i ), i = 1, 2, \ldots n\}\). Abbreviate \(\{(x i ,y i), assign i = 1, 2, … n} to the structure array XY, where the array length is n, the array elements are structures, and the structure members are (x, y). Call the linear track parameter estimation function XYT_TO_kb(n, XY, k, d) to obtain the best linear track parameter (k AX , d AX ) = (k, d). Among them, the implementation process of the function XYT_TO_kb(n, XY, k, d) includes the following steps:

[0072] Step X.1: Initialize the function, define variables tx = 0, tx2 = 0, ty = 0, ty2 = 0, txy = 0, ii = 0.

[0073] Step X.2: Accumulate the x member value of the element with subscript ii in the array XY to the variable tx; square the x member value of the element with subscript ii in the array XY and accumulate it to the variable tx2; accumulate the y member value of the element with subscript ii in the array XY to the variable ty; square the y member value of the element with subscript ii in the array XY and accumulate it to the variable ty2; multiply the x and y member values of the element with subscript ii in the array XY and accumulate it to the variable txy.

[0074] Step X.3: Let ii = ii + 1. If ii < n, go to Step X.2; otherwise, go to Step X.4.

[0075] Step X.4: Let: a1 = tx / n, a2 = tx2 / n, b1 = ty / n, b2 = ty2 / n, c0 = txy / n.

[0076] Step X.5: Let: aa = c0 - a1 * b1, bb = a2 - b2 - a1 * a1 + b1 * b1, cc = a1 * b1 - c0.

[0077] Step X.6: Let:

[0078] d1 = b1 - a1 * k1; d2 = b1 - a1 * k2.

[0079] Step X.7: Let:

[0080]

[0081] Among them, XY[0].x represents the x member value of the element with subscript 0 in the array XY, XY[0].y represents the y member value of the element with subscript 0 in the array XY, and |...| represents taking the absolute value.

[0082] Step X.8: If L1 > L2, take k = k2 and d = d2; otherwise, take k = k1 and d = d1. Output k and d as parameters and end the function run.

[0083] Step 5.1.2: Abbreviate the Y-axis observation data {(t SAi , y SAi ), i = 1, 2, … n} of the master station radar A as: {(x i , y i ), i = 1, 2, … n}. Assign {(x i , y i ) to the structure array XY, where the array length is n, the array elements are structures, and the structure members are (x, y). Call the linear track line parameter estimation function XYT_TO_kb(n, XY, k, d) to obtain the best linear track line parameters (k AY , d AY ) = (k, d) for the Y-axis measurement of the master station radar A.

[0084] Step 5.1.3: Abbreviate the Z-axis observation data {(t SAi , z SAi ), i = 1, 2, … n} of the master station radar A as: {(x i , y i ), i = 1, 2, … n}. Assign {(x i , y i ) to the structure array XY, where the array length is n, the array elements are structures, and the structure members are (x, y). Call the linear track line parameter estimation function XYT_TO_kb(n, XY, k, d) to obtain the best linear track line parameters (k AZ , d AZ ) = (k, d) for the Z-axis measurement of the master station radar A.

[0085] Step 5.2: Use the linear track line parameter estimation model to perform three-dimensional linear parameter estimation on the observation data of the slave station radar B to obtain the observation track line parameters (k BX , d BX ), (k BY , d BY ), (k BZ , d BZ ), including the following steps:

[0086] Step 5.2.1: Abbreviate the X-axis observation data {(t SBj , x SBj ), j = 1, 2, … m} of the slave station radar B as: {(x j , y j ), j = 1, 2, … m}. Assign {(xj , y j ), j = 1, 2, … m} is assigned to the structure array XY, the length of the array is m, the array elements are structures, and the structure members are (x, y). Call the linear track parameter estimation function XYT_TO_kb(m, XY, k, d), and obtain the best linear track parameters (k BX , d BX ) = (k, d).

[0087] Step 5.2.2: Abbreviate the Y-axis observation data {(t SBj , y SBj ), j = 1, 2, … m} of the secondary station radar B as: {(x j , y j ), j = 1, 2, … m}. Assign {(x j , y j ), j = 1, 2, … m} to the structure array XY, the length of the array is m, the array elements are structures, and the structure members are (x, y). Call the linear track parameter estimation function XYT_TO_kb(m, XY, k, d), and obtain the best linear track parameters (k BY , d BY ) = (k, d).

[0088] Step 5.2.3: Abbreviate the Z-axis observation data {(t SBj , z SBj ), j = 1, 2, … m} of the secondary station radar B as: {(x j , y j ), j = 1, 2, … m}. Assign {(x j , y j ), j = 1, 2, … m} to the structure array XY, the length of the array is m, the array elements are structures, and the structure members are (x, y). Call the linear track parameter estimation function XYT_TO_kb(m, XY, k, d), and obtain the best linear track parameters (k BZ , d BZ ) = (k, d).

[0089] The said step 6 includes the following steps:

[0090] Step 6.1: Calculate the three-dimensional rectangular coordinates {(x′ SAi , y′ SAi , z′ SAi ), i = 1, 2, … m} of the time registration points of the master station radar A according to the observation point time of the secondary station radar B. Among them:

[0091] x′ SAi = kAX *t SBi +d AX , y' SAi = k AY *t SBi +d AY , z' SAi = k AZ *t SBi +d AZ .

[0092] Step 6.2: Calculate the three-dimensional rectangular coordinates of the time registration points of the slave station radar B according to the observation point time of the slave station radar B, {(x' SBi , y' SBi , z' SBi ), i = 1, 2, … m}. Among them:

[0093] x' SBi = k BX *t SBi +d BX , y' SBi = k BY *t SBi +d BY , z' SBi = k BZ *t SBi +d BZ .

[0094] Step 6.3: Convert {(x' SBi , y' SBi , z' SBi ), i = 1, 2, … m} into three-dimensional polar coordinates {(R', Bi , θ', Bi , h' Bi ), i = 1, 2, … m} with the slave station radar B as the center. The method is: Assign {(x' SBi , y' SBi , z' SBi ), i = 1, 2, … m} to the structure array XYZ, the length of the array is m, the array elements are structures, and the structure members are (x, y, z); Assign the three-dimensional unified rectangular coordinates of the site of the slave station radar B represented by (X SB , Y SB , Z SB ) to (XO, YO, ZO). Call the rectangular coordinate to polar coordinate function XYZ_TO_RAh(m, XYZ, XO, YO, ZO, RAh) to obtain the three-dimensional polar coordinate array RAh with the slave station radar B as the center. Then, assign the members (RR, AA, hh) of the element with subscript i - 1 in the array RAh to (R', Bi , θ', Bi , h'Bi )。

[0095] Among them: The length of the RAh array is m, the array elements are structures, and the structure members are (RR, AA, hh), representing distance, azimuth, and altitude values. The implementation process of the function XYZ_TO_RAh(m, XYZ, XO, YO, ZO, RAh) includes the following steps:

[0096] Step Z.1: Let: ii = 0;

[0097] Step Z.2: Let: xx = XYZ[ii].x – XO, yy = XYZ[ii].y – YO, zz = XYZ[ii].z – ZO

[0098] Among them, XYZ[ii].x represents the value of the member x of the element with subscript ii in the structure array XYZ; XYZ[ii].y represents the value of the member y of the element with subscript ii in the structure array XYZ; XYZ[ii].z represents the value of the member z of the element with subscript ii in the structure array XYZ.

[0099] Step Z.3: If yy is equal to 0, go to Step Z.4; otherwise, go to Step Z.5.

[0100] Step Z.4: If xx is greater than or equal to 0, assign π / 2 to RAh[ii].AA; otherwise, assign -π / 2 to RAh[ii].AA. Go to Step Z.6. Among them, RAh[ii].AA represents the value of the member AA of the element with subscript ii in the structure array RAh.

[0101] Step Z.5: Assign arctan(xx / yy) to RAh[ii].AA. Among them, arctan() represents the arctangent function.

[0102] Step Z.6: If RAh[ii].AA is less than 0, then assign RAh[ii].AA = RAh[ii].AA + 2π.

[0103] Step Z.7: Assign RAh[ii].RR to

[0104] Step Z.8: Assign XYZ[ii].z to RAh[ii].hh.

[0105] Step Z.9: Let ii = ii + 1. If ii < m, go to Step Z.2. Otherwise, output the structure array RAh as a parameter, and the function runs to completion.

[0106] The said Step 7 includes the following steps:

[0107] Step 7.1: Let:

[0108] A x1i = R′ Bi sinθ′ Bi cosα′ Bi ,A x2i = -R′ Bi cosθ′ Bi cosα′ Bi ,A x3i = -sinθ′ Bi cosα′ Bi ,A x4i = cosθ′ Bi cosα′ Bi ,A x5i = R′ Bi sinθ′ Bi sinα′ Bi ,A x6i = -R′ Bi cosθ′ Bi sinα′ Bi ,A x7i = -sinθ′ Bi sinα′ Bi ,A x8i = cosθ′ Bi sinα′ Bi ;

[0109] A y1i = -A x2i = R′ Bi cosθ′ Bi cosα′ Bi ,A y2i = A x1i = R′ Bi sinθ′ Bi cosα′ Bi ,A y3i = -A x4i = -cosθ′ Bi cosα′ Bi ,A y4i = A x3i = -sinθ′ Bi cosα′ Bi ; A y5i = -A x6i = R′ Bi cosθ′ Bi sinα′ Bi ,A y6i = A x5i = R′ Bi sinθ′ Bi sinα′Bi , A y7i = -A x8i = -cosθ′ Bi sinα′ Bi , A y8i = A x7i = -sinθ′ Bi sinα′ Bi ;

[0110] A z1i = R′ Bi sinα′ Bi , A z2i = -sinα′ Bi , A z3i = -R′ Bi cosα′ Bi , A z4i = cosα′ Bi .

[0111] Step 7.2: Set the relative systematic errors of the distance, azimuth, and elevation of the secondary station radar B as (ΔR, Δθ, Δα) respectively, and construct the following joint solution equations of the relative systematic errors with (ΔR, Δθ, Δα) as unknowns:

[0112]

[0113] Where:

[0114]

[0115]

[0116] Step 7.3: Use the quasi - Newton method for finding a set of real roots of the non - linear equations to solve the equations (1). If the equations have no solution, go to Step 2, select the observation data of the straight - line track of the same airborne target by the master station radar A and the secondary station radar B again, and start a new round of systematic error estimation. If the equations have a solution, assume the solution of the equations is (Dr, Da, Dp).

[0117] The said Step 8 includes the following steps:

[0118] Step 8.1: Judge the singular values of Dr, include the non - singular values in the distance relative systematic error list ER of the secondary station radar B, and sort out the error list ER. It includes the following steps:

[0119] Step 8.1.1: If the length of the list ER is less than the basic statistical times PN, add Dr to the end of the list ER and go to Step 8.2. Otherwise, go to Step 8.1.2. Where the basic statistical times PN represents the minimum number of times to obtain the systematic error estimation value, which is set by the user himself, generally PN≥10.

[0120] Step 8.1.2: Calculate the sample standard deviation DevR and the mean AveR of the elements in the list ER. If the absolute value of the difference between Dr and the mean AveR is less than MU times the sample standard deviation DevR, add Dr to the end of the list ER. Otherwise, go to Step 8.2. Here, MU is set by the user himself, and generally 2 ≥ MU ≥ 3.

[0121] Step 8.1.3: If the length of the list ER is greater than the maximum retention times PM, recalculate the mean AveR of the elements in the list ER, and delete the element with the largest absolute value of the difference from the mean AveR in the list ER. Here, the maximum retention times PM represents the maximum number of times the system error estimate is stored, which is set by the user himself. Generally, PM ≥ 30 and PM > PN.

[0122] Step 8.2: Perform a singular value judgment on Da, incorporate the non-singular values into the azimuth relative system error list EA of the auxiliary station radar B, and organize the error list EA. The steps are as follows:

[0123] Step 8.2.1: If the length of the list EA is less than the basic statistical times PN, add Da to the end of the list EA and go to Step 8.3. Otherwise, go to Step 8.2.2.

[0124] Step 8.2.2: Calculate the sample standard deviation DevA and the mean AveA of the elements in the list EA. If the absolute value of the difference between Da and the mean AveA is less than MU times the sample standard deviation DevA, add Da to the end of the list EA. Otherwise, go to Step 8.3.

[0125] Step 8.2.3: If the length of the list EA is greater than the maximum retention times PM, recalculate the mean AveA of the elements in the list EA, and delete the element with the largest absolute value of the difference from the mean AveA in the list EA.

[0126] Step 8.3: Perform a singular value judgment on Dp, incorporate the non-singular values into the elevation relative system error list EP of the auxiliary station radar B, and organize the error list EP. The steps are as follows:

[0127] Step 8.3.1: If the length of the list EP is less than the basic statistical times PN, add Dp to the end of the list EP and go to Step 9. Otherwise, go to Step 8.3.2.

[0128] Step 8.3.2: Calculate the sample standard deviation DevP and the mean AveP of the elements in the list EP. If the absolute value of the difference between Dp and the mean AveP is less than MU times the sample standard deviation DevP, add Dp to the end of the list EP. Otherwise, go to Step 9.

[0129] Step 8.3.3: If the length of list EP is greater than the maximum retention times PM, recalculate the average value AveP of the elements in list EP, and delete the element in list EP with the largest absolute value of the difference from the average value AveP.

[0130] The said step 9 includes the following steps:

[0131] Step 9.1: If the length of the relative distance system error list ER is greater than or equal to the basic statistical times PN, calculate the average value of the elements in list ER and output it as the cumulative estimation result of the relative distance system error of the secondary station radar B. Otherwise, it is considered that the condition for estimating the relative distance system error of the secondary station radar B is not yet available.

[0132] Step 9.2: If the length of the relative azimuth system error list EA is greater than or equal to the basic statistical times PN, calculate the average value of the elements in list EA and output it as the cumulative estimation result of the relative azimuth system error of the secondary station radar B. Otherwise, it is considered that the condition for estimating the relative azimuth system error of the secondary station radar B is not yet available.

[0133] Step 9.3: If the length of the relative elevation system error list EP is greater than or equal to the basic statistical times PN, calculate the average value of the elements in list EP and output it as the cumulative estimation result of the relative elevation system error of the secondary station radar B. Otherwise, it is considered that the condition for estimating the relative elevation system error of the secondary station radar B is not yet available.

[0134] Embodiment 1:

[0135] This embodiment specifically describes a method for cumulative estimation of the relative system error of a three - coordinate air - to - air radar proposed by the present invention. The simulation calculation process in the said embodiment is implemented based on the C# language and can be applied to the pre - processing process of data in a multi - radar networking detection system.

[0136] The said embodiment includes the following steps:

[0137] Step 1: Initialize the relative distance system error list ER, relative azimuth system error list EA, and relative elevation system error list EP of the secondary station radar B. That is: create lists ER, EA, and EP with the element type of real numbers to store the historical relative system error estimates of distance, azimuth, and elevation of the secondary station radar B. The C# code is as follows:

[0138] List <double>ER = newList <double>(); List <double>EA = newList <double>();

[0139] List <double>EP = newList <double>();

[0140] Step 2: Simulate the observation data of a straight-line track of the same airborne target by the master station radar A and the slave station radar B, which are \(\{(t SAi , R Ai , \theta Ai , h Ai ), i = 1, 2, \ldots, n\}\) and \(\{(t SBj , R Bj , \theta Bj , h Bj ), j = 1, 2, \ldots, m\}\), where \(n\) represents the number of observation data points of the master station radar A, and \(m\) represents the number of observation data points of the slave station radar B. Among them, the positioning and error parameters of radar A and radar B are shown in Table 1; the parameters of a straight-line track of the airborne target are shown in Table 2; the first simulated observation data of radar A and radar B are shown in Table 3 and Table 4 respectively, and \(n = m = 25\).

[0141] Table 1: Basic parameters of the radar

[0142]

[0143]

[0144] Table 2: Track line parameters

[0145] Parameter Single value Start time (seconds) 50 seconds Target starting point x-axis coordinate (km) Uniformly distributed within 1700 ± 3% Target starting point y-axis coordinate (km) Uniformly distributed within 2540 ± 3% Flight altitude (constant, m) Uniformly distributed within 1500 ± 10% Course (0 degrees for true north) Uniformly distributed within 300 ± 30% Target speed (km / hour) Uniformly distributed within 1200 ± 5% Number of track points 25

[0146] Table 3: Observation data of the master station radar A

[0147] Serial number <![CDATA[Time t SA (seconds)]]> <![CDATA[Distance R A (km)]]> <![CDATA[Azimuth θ A (degrees)]]> <![CDATA[Height h A (m)]]> 1 50 152.35 143.71 1.62 2 60 150.54 144.81 1.61 3 70 149.00 145.98 1.60 4 80 146.82 146.82 1.61 5 90 144.90 148.32 1.61 6 100 143.58 149.07 1.61 7 110 142.11 150.30 1.60 8 120 140.82 151.59 1.61 9 130 139.21 152.74 1.60 10 140 137.98 153.94 1.62 11 150 136.43 155.62 1.60 12 160 134.86 156.96 1.64 13 170 133.97 157.55 1.61 14 180 132.87 159.21 1.60 15 190 131.65 160.54 1.59 16 200 130.64 161.90 1.59 17 210 129.90 163.14 1.62 18 220 128.97 164.83 1.62 19 230 128.12 165.98 1.61 20 240 127.68 166.98 1.61 21 250 126.93 169.03 1.64 22 260 126.38 170.73 1.64 23 270 125.93 171.87 1.65 24 280 125.67 173.27 1.62 25 290 125.52 175.12 1.63

[0148] Table 4: Observation data of the slave station radar B

[0149]

[0150]

[0151] Step 3: Perform coordinate transformation on the observation data of the master station radar A to obtain a set of corresponding three-dimensional unified rectangular coordinates

[0152] \{(t SAi , x SAi , y SAi , z SAi ), i = 1, 2, \ldots, n\}\), and the calculation results of the first simulated observation data are shown in Table 5. It includes the following steps:

[0153] Step 3.1: The observation data of the master station radar A \(\{(t SAi , R Ai , \theta Ai , h Ai ), i = 1, 2, … n} is converted into three-dimensional rectangular coordinates centered on this station {(t SAi , x Ai , y Ai , z Ai ), i = 1, 2, … n}. Wherein:

[0154]

[0155] Step 3.2: Convert {(t SAi , x Ai , y Ai , z Ai ), i = 1, 2, … n} into three-dimensional unified rectangular coordinates {(t SAi , x SAi , y SAi , z SAi ), i = 1, 2, … n}.

[0156]

[0157] (X SA , Y SA , Z SA ) represents the three-dimensional rectangular coordinates of radar A in the unified coordinate system. In this embodiment, the values are (1580 km, 2680 km, 0.5 km).

[0158] Table 5: Three-dimensional unified rectangular coordinates of the observation data of the main station radar A

[0159]

[0160]

[0161] Step 4: Perform coordinate transformation on the observation data of the auxiliary station radar B to obtain a set of corresponding three-dimensional unified rectangular coordinates {(t SBj , x SBj , y SBj , z SBj ), j = 1, 2, … m}. The calculation results of the first simulated observation data are shown in Table 6. The following steps are included:

[0162] Step 4.1: Convert the observation data of the auxiliary station radar B {(t SBj , R Bj , θ Bj , h Bj ), j = 1, 2, … m} into three-dimensional rectangular coordinates centered on this station {(t SBj , x Bj , y Bj , z Bj ), j = 1, 2, … m}. Wherein:

[0163]

[0164] Step 4.2: Convert \(\{(t SBj , x Bj , y Bj , z Bj ), j = 1, 2, … m\}\) into three-dimensional unified rectangular coordinates \(\{(t SBj , x SBj , y SBj , z SBj ), j = 1, 2, … m\}\). Wherein:

[0165]

[0166] (X SB , Y SB , Z SB ) represents the three-dimensional rectangular coordinates of radar B in the unified coordinate system. In this embodiment, the values are (1500 km, 2550 km, 0.2 km).

[0167] Table 6: Three-dimensional unified rectangular coordinates of the observation data of the auxiliary station radar B

[0168]

[0169]

[0170] Figure 2 Figure 6 shows the XY-plane projection schematic diagram of the first simulated observation track lines of radars A and B drawn according to Tables 5 and 6. It can be seen from the figure that due to the systematic errors in the measurements of radars A and B, the two observation track lines for the same target are split; at the same time, the random errors in the measurements make the single track line show serrations of different degrees.

[0171] Step 5: Use the straight track line parameter estimation model in the three-dimensional unified rectangular coordinate system to perform three-dimensional straight line parameter estimation on the observation data of the master and auxiliary station radars respectively, and obtain the observation track line parameters of the master station radar and the observation track line parameters of the auxiliary station radar. The steps are as follows:

[0172] Step 5.1: Use the straight track line parameter estimation model to perform three-dimensional straight line parameter estimation on the observation data of the master station radar A, and obtain the observation track line parameters of the master station radar \((k AX , d AX ), \((k AY , d AY ), \((k AZ , d AZ ). The steps are as follows:

[0173] Step 5.1.1: Convert the X-axis observation data \(\{(t SAi , x SAi ), i = 1, 2, … n} is briefly recorded as: {(x i , y i ), i = 1, 2, … n}. Assign {(x i , y i ), i = 1, 2, … n} to the structure array XY, the array length is n, the array elements are structures, and the structure members are (x, y). Call the linear track line parameter estimation function XYT_TO_kb(n, XY, k, d) to obtain the best linear track line parameters (k AX , d AX ) = (k, d) measured on the X-axis of the master station radar A. The C# implementation code of the function XYT_TO_kb(n, XY, k, d) is as follows:

[0174]

[0175]

[0176] Step 5.1.2: Abbreviate the Y-axis observation data {(t SAi , y SAi ), i = 1, 2, … n} of the master station radar A as: {(x i , y i ), i = 1, 2, … n}. Assign {(x i , y i ), i = 1, 2, … n} to the structure array XY, the array length is n, the array elements are structures, and the structure members are (x, y). Call the linear track line parameter estimation function XYT_TO_kb(n, XY, k, d) to obtain the best linear track line parameters (k AY , d AY ) = (k, d).

[0177] Step 5.1.3: Abbreviate the Z-axis observation data {(t SAi , z SAi ), i = 1, 2, … n} of the master station radar A as: {(x i , y i ), i = 1, 2, … n}. Assign {(x i , y i ), i = 1, 2, … n} to the structure array XY, the array length is n, the array elements are structures, and the structure members are (x, y). Call the linear track line parameter estimation function XYT_TO_kb(n, XY, k, d) to obtain the best linear track line parameters (k AZ , d AZ ) = (k, d).

[0178] Obtain the observation track line parameters (k AX , d AX ) = (-0.3285, 1686.3614), (k AY , d AY ) = (-0.0081, 2557.3880), (k AZ , d AZ ) = (0.0001, 1.5948) of the first simulated observation data of the master station radar A according to step 5.1.

[0179] Step 5.2: Use the straight-line track line parameter estimation model to perform three-dimensional straight-line parameter estimation on the observation data of the slave station radar B to obtain the observation track line parameters (k BX , d BX ), (k BY , d BY ), (k BZ , d BZ ), including the following steps:

[0180] Step 5.2.1: Abbreviate the X-axis observation data {(t SBj , x SBj ), j = 1, 2, … m} of the slave station radar B as: {(x j , y j ), j = 1, 2, … m}. Assign {(x j , y j ), j = 1, 2, … m} to the structure array XY, the array length is m, the array elements are structures, and the structure members are (x, y). Call the straight-line track line parameter estimation function XYT_TO_kb(m, XY, k, d) to obtain the best straight-line track line parameters (k BX , d BX ) = (k, d) of the X-axis measurement of the slave station radar B.

[0181] Step 5.2.2: Abbreviate the Y-axis observation data {(t SBj , y SBj ), j = 1, 2, … m} of the slave station radar B as: {(x j , y j ), j = 1, 2, … m}. Assign {(x j , y j ), j = 1, 2, … m} to the structure array XY, the array length is m, the array elements are structures, and the structure members are (x, y). Call the straight-line track line parameter estimation function XYT_TO_kb(m, XY, k, d) to obtain the best straight-line track line parameters (k BY , d BY ) = (k, d) of the Y-axis measurement of the slave station radar B.

[0182] Step 5.2.3: Denote the Z-axis observation data of the secondary station radar B \(\{(t SBj , z SBj ), j = 1, 2, … m\}\) simply as: \(\{(x j , y j ), j = 1, 2, … m\}\). Assign \(\{(x j , y j ), j = 1, 2, … m\}\) to the structure array XY with the array length of m and the array elements being structures, and the structure members being (x, y). Call the linear track line parameter estimation function XYT_TO_kb(m, XY, k, d) to obtain the optimal linear track line parameters (k BZ , d BZ ) = (k, d) measured on the Z-axis of the secondary station radar B.

[0183] According to Step 5.2, obtain the observation track line parameters of the first simulated observation data of the secondary station radar B: (k BX , d BX ) = (-0.3286, 1686.5839), (k BY , d BY ) = (0.0034, 2551.3949), (k BZ , d BZ ) = (-0.0028, 3.2313).

[0184] Step 6: Calculate the time registration points of the master station and secondary station radars according to the observation point time of the secondary station radar. Convert the time registration points of the secondary station radar into three-dimensional polar coordinate registration data centered on the radar site of this station, denoted as \(\{(R′ Bi , θ′ Bi , h′ Bi ), i = 1, 2, … m\}. It includes:

[0185] Step 6.1: Calculate the three-dimensional rectangular coordinates \(\{(x′ SAi , y′ SAi , z′ SAi ), i = 1, 2, … m\}\) of the time registration point of the master station radar A according to the observation point time of the secondary station radar B. Among them:

[0186] x′ SAi = k AX * t SBi + d AX , y′ SAi = k AY * t SBi + d AY , z′ SAi = k AZ * t SBi +d AZ 。

[0187] The three-dimensional rectangular coordinates of the time registration points of the first simulated observation track line of the main station radar A are {(x′ SAi , y′ SAi , z′ SAi ), i = 1, 2, … m}, and the calculation results are shown in Table 7 below.

[0188] Table 7: Three-dimensional rectangular coordinates of the time registration points of the main station radar A

[0189] Serial number <![CDATA[Time t SB (seconds)]]> <![CDATA[X coordinate x′ SA (km)]]> <![CDATA[Y coordinate y′ SA (km)]]> <![CDATA[Z coordinate z′ SA (kilometer)]]> 1 54 1668.62 2556.95 1.60 2 64 1665.34 2556.87 1.60 3 74 1662.06 2556.79 1.60 4 84 1658.77 2556.71 1.60 5 94 1655.49 2556.63 1.61 6 104 1652.20 2556.55 1.61 7 114 1648.92 2556.47 1.61 8 124 1645.63 2556.39 1.61 9 134 1642.35 2556.31 1.61 10 144 1639.06 2556.22 1.61 11 154 1635.78 2556.14 1.61 12 164 1632.49 2556.06 1.61 13 174 1629.21 2555.98 1.61 14 184 1625.93 2555.90 1.61 15 194 1622.64 2555.82 1.62 16 204 1619.36 2555.74 1.62 17 214 1616.07 2555.66 1.62 18 224 1612.79 2555.58 1.62 19 234 1609.50 2555.50 1.62 20 244 1606.22 2555.42 1.62 21 254 1602.93 2555.34 1.62 22 264 1599.65 2555.25 1.62 23 274 1596.36 2555.17 1.62 24 284 1593.08 2555.09 1.63 25 294 1589.79 2555.01 1.63

[0190] Step 6.2: Calculate the three-dimensional rectangular coordinates of the time registration points of the slave station radar B according to the observation point time of the slave station radar B, which are {(x′ SBi , y′ SBi , z′ SBi ), i = 1, 2, … m}. Among them:

[0191] x′ SBi = k BX * t SBi + d BX , y′ SBi = k BY * t SBi + d BY , z′ SBi = k BZ * t SBi + d BZ 。

[0192] The three-dimensional rectangular coordinates of the time registration points of the first simulated observation track line of the slave station radar B are {(x′ SBi , y′ SBi , z′ SBi ), i = 1, 2, … m}, and the calculation results are shown in Table 8 below.

[0193] Table 8: Three-dimensional rectangular coordinates of the time registration points of the slave station radar B

[0194] Serial number <![CDATA[Time t SB (seconds)]]> <![CDATA[X coordinate x′ SB (km)]]> <![CDATA[Y coordinate y′ SB (km)]]> <![CDATA[Z coordinate z′ SB (kilometer)]]> 1 54 1668.84 2551.58 3.08 2 64 1665.55 2551.62 3.05 3 74 1662.27 2551.65 3.02 4 84 1658.98 2551.68 2.99 5 94 1655.70 2551.72 2.97 6 104 1652.41 2551.75 2.94 7 114 1649.12 2551.79 2.91 8 124 1645.84 2551.82 2.88 9 134 1642.55 2551.86 2.85 10 144 1639.27 2551.89 2.82 11 154 1635.98 2551.93 2.80 12 164 1632.70 2551.96 2.77 13 174 1629.41 2551.99 2.74 14 184 1626.12 2552.03 2.71 15 194 1622.84 2552.06 2.68 16 204 1619.55 2552.10 2.65 17 214 1616.27 2552.13 2.63 18 224 1612.98 2552.17 2.60 19 234 1609.69 2552.20 2.57 20 244 1606.41 2552.24 2.54 21 254 1603.12 2552.27 2.51 22 264 1599.84 2552.30 2.48 23 274 1596.55 2552.34 2.46 24 284 1593.26 2552.37 2.43 25 294 1589.98 2552.41 2.40

[0195] Step 6.3: Convert {(x′ SBi , y′ SBi , z′ SBi ), i = 1, 2, … m} into three-dimensional polar coordinates centered on the slave station radar B, which are {(R′ Bi , θ′ Bi , h′ Bi ), i = 1, 2, … m}. The method is: Convert {(x′ SBi , y′ SBi , z′ SBi ), assign \(i = 1, 2, \ldots, m\) to the structure array XYZ. The length of the array is \(m\), and the array elements are structures. The structure members are \((x, y, z)\); assign the three - dimensional unified rectangular coordinates of the secondary station radar B represented by \((X SB , Y SB , Z SB ) to \((XO, YO, ZO)\). Call the rectangular - to - polar coordinate function XYZ_TO_RAh(m, XYZ, XO, YO, ZO, RAh) to obtain the three - dimensional polar coordinate array RAh centered on the secondary station radar B. Then, assign the members \((RR, AA, hh)\) of the element with index \(i - 1\) in the array RAh to \((R' Bi , \theta' Bi , h' Bi ) one by one.

[0196] Among them: The length of the RAh array is \(m\), the array elements are structures, and the structure members are \((RR, AA, hh)\), representing distance, azimuth, and altitude values. The C# implementation code of the function XYZ_TO_RAh(m, XYZ, XO, YO, ZO, RAh) is as follows:

[0197]

[0198]

[0199] The time registration points of the first - time simulated observation track line of the secondary station radar B are the three - dimensional polar coordinates \(\{(R' Bi , \theta' Bi , h' Bi ), \(i = 1, 2, \ldots, m\}\) centered on the secondary station radar B. The calculation results are shown in Table 9.

[0200] Table 9: Three - dimensional polar coordinates of the time registration points of the secondary station radar B

[0201]

[0202]

[0203] Step 7: Set the relative systematic errors of the distance, azimuth, and elevation of the secondary station radar B as \((\Delta R, \Delta \theta, \Delta \alpha)\) respectively. Use the time registration point data of the master station and the secondary station to construct a joint solution system of equations for the relative systematic errors with \((\Delta R, \Delta \theta, \Delta \alpha)\) as unknowns, and use the quasi - Newton method for finding a set of real roots of the non - linear equations to solve this system of equations. If the system of equations has no solution, go to Step 2, select the straight - track line observation data of the master station radar A and the secondary station radar B for the same airborne target again, and start a new round of systematic error estimation. If the system of equations has a solution, assume the solution of the system of equations is \((Dr, Da, Dp)\). It includes the following steps:

[0204] Step 7.1: Let:

[0205]

[0206] A x1i = R′ Bi sinθ′ Bi cosα′ Bi ,A x2i = -R′ Bi cosθ′ Bi cosα′ Bi ,A x3i = -sinθ′ Bi cosα′ Bi ,A x4i = cosθ′ Bi cosα′ Bi ,A x5i = R′ Bi sinθ′ Bi sinα′ Bi ,A x6i = -R′ Bi cosθ′ Bi sinα′ Bi ,A x7i = -sinθ′ Bi sinα′ Bi ,A x8i = cosθ′ Bi sinα′ Bi ;

[0207] A y1i = -A x2i = R′ Bi cosθ′ Bi cosα′ Bi ,A y2i = A x1i = R′ Bi sinθ′ Bi cosα′ Bi ,A y3i = -A x4i = -cosθ′ Bi cosα′ Bi ,

[0208] A y4i = A x3i = -sinθ′ Bi cosα′ Bi ; A y5i = -A x6i = R′ Bi cosθ′ Bi sinα′ Bi ,A y6i = A x5i = R' Bi sinθ' Bi sinα' Bi ,

[0209] A y7i = -A x8i = -cosθ' Bi sinα' Bi , A y8i = A x7i = -sinθ' Bi sinα' Bi ;

[0210] A z1i = R' Bi sinα' Bi , A z2i = -sinα' Bi , A z3i = -R' Bi cosα' Bi , A z4i = cosα' Bi .

[0211] The above parameter values calculated from the first simulated observation data of the secondary station radar B (see Table 9) are shown in Tables 10, 11, and 12.

[0212] Table 10: Parameter A calculated from the first simulated observation data of the secondary station radar B x Value

[0213]

[0214]

[0215] Table 11: Parameter A calculated from the first simulated observation data of the secondary station radar B y Value

[0216] Serial number <![CDATA[A y1 > <![CDATA[A y2 > <![CDATA[A y3 > <![CDATA[A y4 > <![CDATA[A y5 > <![CDATA[A y6 > <![CDATA[A y7 > <![CDATA[A y8 > 1 1.580958 168.840008 -0.009362 -0.999811 0.026952 2.878360 -0.000160 -0.017045 2 1.615422 165.554101 -0.009756 -0.999804 0.027810 2.850060 -0.000168 -0.017212 3 1.649886 162.268195 -0.010166 -0.999797 0.028691 2.821760 -0.000177 -0.017386 4 1.684349 158.982289 -0.010592 -0.999790 0.029596 2.793459 -0.000186 -0.017567 5 1.718813 155.696382 -0.011037 -0.999781 0.030526 2.765157 -0.000196 -0.017756 6 1.753276 152.410476 -0.011501 -0.999773 0.031484 2.736854 -0.000207 -0.017953 7 1.787740 149.124570 -0.011985 -0.999763 0.032471 2.708551 -0.000218 -0.018159 8 1.822204 145.838663 -0.012492 -0.999753 0.033489 2.680246 -0.000230 -0.018374 9 1.856667 142.552757 -0.013021 -0.999742 0.034540 2.651940 -0.000242 -0.018598 10 1.891131 139.266851 -0.013576 -0.999730 0.035627 2.623633 -0.000256 -0.018834 11 1.925594 135.980944 -0.014157 -0.999718 0.036752 2.595325 -0.000270 -0.019081 12 1.960058 132.695038 -0.014767 -0.999704 0.037918 2.567015 -0.000286 -0.019339 13 1.994522 129.409132 -0.015408 -0.999689 0.039128 2.538703 -0.000302 -0.019612 14 2.028985 126.123225 -0.016082 -0.999673 0.040385 2.510390 -0.000320 -0.019898 15 2.063449 122.837319 -0.016792 -0.999655 0.041694 2.482074 -0.000339 -0.020199 16 2.097912 119.551413 -0.017542 -0.999636 0.043059 2.453756 -0.000360 -0.020517 17 2.132376 116.265506 -0.018334 -0.999614 0.044484 2.425436 -0.000382 -0.020853 18 2.166840 112.979600 -0.019171 -0.999591 0.045974 2.397113 -0.000407 -0.021209 19 2.201303 109.693694 -0.020059 -0.999566 0.047536 2.368787 -0.000433 -0.021585 20 2.235767 106.407787 -0.021002 -0.999538 0.049176 2.340457 -0.000462 -0.021985 21 2.270230 103.121881 -0.022004 -0.999507 0.050901 2.312123 -0.000493 -0.022410 22 2.304694 99.835975 -0.023073 -0.999472 0.052721 2.283785 -0.000528 -0.022863 23 2.339158 96.550068 -0.024214 -0.999434 0.054643 2.255442 -0.000566 -0.023347 24 2.373621 93.264162 -0.025435 -0.999391 0.056681 2.227092 -0.000607 -0.023865 25 2.408085 89.978256 -0.026745 -0.999344 0.058845 2.198736 -0.000654 -0.024420

[0217] Table 12: Parameter A calculated from the first simulated observation data of the secondary station radar B z Value

[0218]

[0219]

[0220] Step 7.2: Set the relative systematic errors of the distance, azimuth, and elevation of the secondary station radar B to (ΔR, Δθ, Δα) respectively, and construct the following joint solution equations for the relative systematic errors with (ΔR, Δθ, Δα) as unknowns:

[0221]

[0222] In this embodiment, the function dnetnfun_3D(rap, ff, m) is used to calculate the function value f on the left side of the equation set (1). 1 (ΔR, Δθ, Δα), f 2 (ΔR, Δθ, Δα) and f 3 (ΔR, Δθ, Δα). rap is a one-dimensional array with a length of 3, storing the values of the unknowns (ΔR, Δθ, Δα); m = 25, representing the number of observation data points of the secondary station radar B, that is, the number of time registration points; ff is a one-dimensional array with a length of 3, returning the function value on the left side of the equation set (1). The C# implementation code of the dnetnfun_3D() function is as follows:

[0223]

[0224]

[0225] Step 7.3: Use the quasi-Newton method for finding a set of real roots of the nonlinear equation set to solve the equation set (1). If the equation set has no solution, go to Step 2, select the observation data of the straight flight track of the same air target by the master station radar A and the secondary station radar B again, and start a new round of system error estimation. If the equation set has a solution, assume the solution of the equation set is (Dr, Da, Dp).

[0226] Among them, the implementation function of the quasi-Newton method for finding a set of real roots of the nonlinear equation set in this embodiment is named dnetn_3D(nn, eps, tt, hhh, rap, kk, m). In this embodiment, nn = 3, indicating that there are 3 unknowns in the equation set (1); eps = 0.0000001, indicating the control accuracy requirement; tt = 0.2, which is a variable for controlling the change size of hhh; hhh = 1.0, indicating the initial value of the increment; kk = 3000, indicating the maximum number of allowed iterations; m = 25, representing the number of observation data points of the secondary station radar B, that is, the number of time registration points; rap is a one-dimensional array with a length of nn, storing the initial values (0.0, 0.0, 0.0) of the unknowns (ΔR, Δθ, Δα), the function return value is the actual number of iterations, and a set of real solutions of the equation set is stored in rap. The C# implementation code of the dnetn_3D() function is as follows:

[0227]

[0228]

[0229] Among them, agaus() is a function for solving linear algebraic equation sets by the full pivoting Gaussian elimination method, and the function implementation flowchart is as Figure 3 as shown

[0230] In this embodiment, according to the first simulated observation data, the dnetn_3D() function actually iterates and calculates 5 times, and the solution of the system of equations is (0.092578, 0.030991, 0.008693), which is assigned to the variables (Dr, Da, Dp), indicating that the relative systematic errors of distance, azimuth, and pitch obtained in this round of calculation are 0.092578 km, 1.775655 degrees, and 0.498073 degrees, respectively. Figure 4 and Figure 5 and Figure 6 show the results of the relative systematic errors of distance, azimuth, and pitch obtained by randomly generating flight routes according to the parameters in Table 2 and performing 100 consecutive simulated observation data calculations. The actual number of iterations of the dnetn_3D() function in 100 calculations does not exceed 5 times.

[0231] Step 8: Perform singular value judgment on Dr, Da, and Dp respectively, incorporate the non-singular values into the relative systematic error lists ER, EA, and EP of the auxiliary station radar B, and organize the error lists. The steps are as follows:

[0232] Step 8.1: Perform singular value judgment on Dr, incorporate the non-singular values into the distance relative systematic error list ER of the auxiliary station radar B, and organize the error list ER. The steps are as follows:

[0233] Step 8.1.1: If the length of the list ER is less than the basic statistical number PN, add Dr to the end of the list ER, and go to Step 8.2. Otherwise, go to Step 8.1.2. Among them, the basic statistical number PN is set by the user himself, indicating the minimum number of times to obtain the systematic error estimate, generally PN≥10. In this embodiment, PN is set to 15.

[0234] Step 8.1.2: Calculate the sample standard deviation DevR and the mean AveR of the elements in the list ER. If the absolute value of the difference between Dr and the mean AveR is less than MU times the sample standard deviation DevR, add Dr to the end of the list ER. Otherwise, go to Step 8.2. Among them, MU is set by the user himself, generally 2≥MU≥3. In this embodiment, MU is set to 2.5.

[0235] Step 8.1.3: If the length of the list ER is greater than the maximum retention number PM, recalculate the mean AveR of the elements in the list ER, and delete the element with the largest absolute value of the difference from the mean AveR in the list ER. Among them, the maximum retention number PM is set by the user himself, indicating the maximum number of times to store the systematic error estimate, generally PM≥30, and PM>PN. In this embodiment, PM is set to 45.

[0236] The C# implementation code for step 8.1 is as follows:

[0237]

[0238]

[0239] Step 8.2: Perform singular value judgment on Da, incorporate the non-singular values into the azimuth relative system error list EA of the secondary station radar B, and organize the error list EA. It includes the following steps:

[0240] Step 8.2.1: If the length of the list EA is less than the basic statistical times PN, add Da to the end of the list EA, and go to step 8.3. Otherwise, go to step 8.2.2.

[0241] Step 8.2.2: Calculate the sample standard deviation DevA and the mean AveA of the elements in the list EA. If the absolute value of the difference between Da and the mean AveA is less than MU times the sample standard deviation DevA, add Da to the end of the list EA. Otherwise, go to step 8.3.

[0242] Step 8.2.3: If the length of the list EA is greater than the maximum retention times PM, recalculate the mean AveA of the elements in the list EA, and delete the element with the largest absolute value of the difference from the mean AveA in the list EA.

[0243] The C# implementation code for step 8.2 is as follows:

[0244]

[0245] Step 8.3: Perform singular value judgment on Dp, incorporate the non-singular values into the elevation relative system error list EP of the secondary station radar B, and organize the error list EP. It includes the following steps:

[0246] Step 8.3.1: If the length of the list EP is less than the basic statistical times PN, add Dp to the end of the list EP, and go to step 9. Otherwise, go to step 8.3.2.

[0247] Step 8.3.2: Calculate the sample standard deviation DevP and the mean AveP of the elements in the list EP. If the absolute value of the difference between Dp and the mean AveP is less than MU times the sample standard deviation DevP, add Dp to the end of the list EP. Otherwise, go to step 9.

[0248] Step 8.3.3: If the length of the list EP is greater than the maximum retention times PM, recalculate the mean AveP of the elements in the list EP, and delete the element with the largest absolute value of the difference from the mean AveP in the list EP.

[0249] The C# implementation code for step 8.3 is as follows:

[0250]

[0251] Step 9: If the length of the list is greater than or equal to the basic statistical number PN, calculate the mean value of the list elements and output it as the corresponding cumulative estimation result of the relative system error. Otherwise, it is considered that the corresponding relative system error estimation condition is not yet available. It includes the following steps:

[0252] Step 9.1: If the length of the distance relative system error list ER is greater than or equal to the basic statistical number PN, calculate the mean value of the elements of the list ER and output it as the cumulative estimation result of the distance relative system error of the slave station radar B. Otherwise, it is considered that the distance relative system error estimation condition of the slave station radar B is not yet available.

[0253] Step 9.2: If the length of the azimuth relative system error list EA is greater than or equal to the basic statistical number PN, calculate the mean value of the elements of the list EA and output it as the cumulative estimation result of the azimuth relative system error of the slave station radar B. Otherwise, it is considered that the azimuth relative system error estimation condition of the slave station radar B is not yet available.

[0254] Step 9.3: If the length of the pitch relative system error list EP is greater than or equal to the basic statistical number PN, calculate the mean value of the elements of the list EP and output it as the cumulative estimation result of the pitch relative system error of the slave station radar B. Otherwise, it is considered that the pitch relative system error estimation condition of the slave station radar B is not yet available.

[0255] Step 10: Go to Step 2, select the observation data of the other straight track lines of the same airborne target by the master station radar A and the slave station radar B, and start a new round of system error estimation.

[0256] In this embodiment, 100 different straight track line observation results of the same airborne target by the master station radar A and the slave station radar B are simulated according to Step 2. After being processed by Steps 3 to 7, the single relative system error calculation results can be obtained. Finally, after the singular value rejection, error list sorting and discrimination are performed in Steps 8 and 9, the lengths of the distance, azimuth and pitch relative system error lists ER, EA and EP are all 45, which is greater than the specified basic statistical number PN = 15. The mean values of the output lists are used as the cumulative estimation results of the distance, azimuth and pitch relative system errors, which are (-0.1972, 0.031704, 0.008721), that is, (-0.1972 km, 1.8165 degrees, 0.4997 degrees). Compared with the system error true values (-0.2 km, 1.8 degrees, 0.5 degrees) of the slave station radar B set in Table 1, the degrees of closeness are (98.60%, 99.08%, 99.94%).

[0257] Step 11: Verification experiment. To further verify the effectiveness of the error estimation, we use the cumulative estimation results of the range, azimuth, and elevation relative systematic errors of the secondary radar B (-0.1972, 0.031704, 0.008721) obtained in Step 9 to correct the measurement data of other track lines of the secondary radar B. By comparing the changes in the distances between the track lines before and after correction and the observation track line of the primary radar A, the effectiveness of the error estimation is verified. The verification experiment includes the following steps:

[0258] Step 11.1: Keep the basic parameters of the primary radar A and the secondary radar B set in Table 1 unchanged, and simulate the generation of the track line observation data of the primary and secondary radars for the same airborne target according to a set of composite (linear + circular arc) track line parameters of the airborne target set in Table 13.

[0259] Table 13: Composite track line parameters

[0260]

[0261]

[0262]

[0263] Simulate the generation of the observation data according to Table 1 and Table 13. The XY plane projection of the observation track lines of the primary and secondary radars generated in the first simulation is as Figure 7 shown. After calculation, in the three-dimensional unified rectangular coordinate system, the average distance between the observation track lines of the primary and secondary radars is 4.9368 km.

[0264] Step 11.2: Correct the observed track points of the secondary radar B according to the cumulative estimation results of the range, azimuth, and elevation relative systematic errors of the secondary radar B (-0.1972, 0.031704, 0.008721), that is, (-0.1972 km, 1.8165 degrees, 0.4997 degrees). The XY plane projections of the observation track lines of the primary and secondary radars and the corrected track points of the secondary radar generated in the first simulation after correction are as Figure 8 shown. After calculation, in the three-dimensional unified rectangular coordinate system, the average distance between the corrected track line of the secondary radar and the observation track line of the primary radar is 0.1714 km, and the correction rate is 96.53%; the comparison results of the changes in the observed target height values of the secondary radar before and after correction and the simulated height true value in the first simulation are as Figure 9 shown. The corrected height value is very close to the true value and is almost indistinguishable; the calculation results of the correction rate of the corrected height value relative to the simulated true value are as Figure 10 shown. When the target is flying at a normal height, the height correction rate remains above 95%. When the target is taking off or landing and the flight height is relatively low (200 m - 500 m), the height correction rate is about 80%.

[0265] Step 11.3: Repeat Steps 11.1 and 11.2 for 100 times. After calculation, in the three-dimensional rectangular coordinate system, the average distances between the auxiliary station track lines before and after correction and the main station observation track line are as Figure 11 shown, and the correction rates are as Figure 12 shown. After statistics, the minimum correction rate is 84.02%, the maximum is 98.54%, and the average is 93.34%. Among the 100 corrections, the correction rate remained above 90% in 83 times, and the average distance between the main and auxiliary station track lines decreased by one order of magnitude, and the correction effect is very ideal. Thus, the effectiveness of using the method described in the present invention for systematic error estimation is further verified.

[0266] From the calculation and statistical results of the above embodiments, by using the method described in the present invention, the degree of approximation between the cumulative estimation results of the relative systematic errors of the auxiliary station radar in terms of distance, azimuth and pitch and the set true value of the systematic error is very high (the lowest value is 98.60%); using this result to correct the other track line measurement data of the auxiliary station radar, the average distance between the auxiliary station track lines before and after correction and the main station track line decreased by 93.34%, and the correction effect is very ideal, meeting the requirements for the accuracy of air radar data fusion in engineering practice.

[0267] For those of ordinary skill in the art in this technical field, without departing from the technical principle of the present invention, several improvements and deformations can be made, and these improvements and deformations should also be regarded as the protection scope of the present invention. For example, but not limited to the following points:

[0268] (1) Regarding the use of lists. The list in the present invention is a data set with a specific structure as the data type. The number of set elements can change dynamically, and the elements are connected by pointers, and elements can be easily inserted and deleted, and the capacity can be automatically changed. For example, in C and C++, the list in the present invention can be expressed in the form of a single linked list, and in C#, it can be expressed by a generic List.

[0269] (2) Regarding the selection of the main station radar. The present invention designates the measurement value of the most accurate air radar in the radar network as the true description of the target position, and names it the main station. In engineering practice, through air traffic control radars, flight inspection equipment, ADS-B equipment, etc., measurement data with higher accuracy and data rate can be obtained, and this can be used as the true description of the target position. Selecting a measurement device with higher accuracy and data rate as the main station will help improve the error estimation and track correction effects of the present invention, and such improvements and deformations should also be regarded as the protection scope of the present invention.

[0270] (3)Regarding the solution of the system of equations for the joint calculation of relative systematic errors. In step 7 of the present invention, the quasi-Newton method for finding a set of real roots of a non-linear system of equations is used to solve the system of equations (1). Simulation experiments have proved that this method is time-saving and reliable. Exploring the use of other numerical calculation methods to solve the system of equations (1) should be regarded as an improvement and variation of step 7.

[0271] (4)Regarding the values of PN, MU, and PM in step 8. In the present invention, the same values of PN, MU, and PM are used for the singular value processing and error list sorting of the relative systematic errors of distance, azimuth, and elevation. In practical engineering applications, they can be treated differently according to the different degrees of influence of the three errors on the track correction effect.

[0272] (5)In the present invention, the cumulative estimation of relative systematic errors is an offline and post-event process. Compared with the time complexity of program implementation, we pay more attention to the accuracy of the estimation results and the online and post-event track correction effect. Therefore, the improvement and enhancement of the execution efficiency of the method or embodiment program of the present invention should also be regarded as the protection scope of the present invention.

[0273] The present invention adopts the idea of relative systematic errors. By selecting the observation data of the master and slave station radars for typical route targets, through three-dimensional unified rectangular coordinate transformation, three-dimensional straight track line parameter estimation, and time registration point calculation, etc., a joint calculation model of relative systematic errors derived from asynchronous radar measurement data is established, and the reliability of the cumulative estimation result of relative systematic errors is ensured through operations such as singular value elimination. The experimental results show that using this cumulative estimation result to correct the subsequent measurements of the slave station radar, the degree of observation track splitting is reduced by nearly an order of magnitude, and the error correction effect is very ideal, meeting the requirements for the accuracy of airborne radar data fusion in engineering practice. The estimation method provided by the present invention is scientific, the implementation steps of the scheme are reasonable, and the operability and practicability are very strong.

[0274] The above are only the preferred embodiments of the present invention. It should be noted that for those of ordinary skill in the art in this technical field, without departing from the technical principle of the present invention, several improvements and variations can still be made, and these improvements and variations should also be regarded as the protection scope of the present invention.< / double> < / double> < / double> < / double> < / double> < / double>

Claims

1. A three-coordinate air radar relative system error accumulation estimation method, characterized in that: The method comprises the following steps: Step 1: Initialize the distance relative system error list ER, azimuth relative system error list EA and elevation relative system error list EP of the auxiliary station radar B; Step 2: Select the observation data of a straight line track of the main radar A and the auxiliary radar B on the same aerial target {(t SAi ,R Ai ,θ Ai ,h Ai ), i=1,2,…n} and {(t SBj ,R Bj ,θ Bj ,h Bj ), j = 1, 2, ... m}, n represents the number of observation data points of the main station radar A, m represents the number of observation data points of the auxiliary station radar B, (t SAi ,R Ai ,θ Ai ,h Ai ) indicates t SAi The target distance R measured by the master radar A at the moment Ai 、Directionθ Ai and altitude h Ai , (t SBj ,R Bj ,θ Bj ,h Bj ) indicates t SBj The target distance R measured by the auxiliary station radar B at the moment Bj 、Directionθ Bj and altitude h Bj ; Step 3: Perform coordinate transformation on the observation data of the master radar A to obtain a set of corresponding three-dimensional unified rectangular coordinates {(t SAi ,x SAi ,y SAi ,z SAi ), i=1,2,…n}; Step 4: Perform coordinate transformation on the observation data of the auxiliary radar B to obtain a set of corresponding three-dimensional unified rectangular coordinates {(t SBj ,x SBj ,y SBj ,z SBj ), j = 1, 2, ... m}; Step 5: Use the straight line track parameter estimation model to estimate the three-dimensional straight line parameters of the main and auxiliary station radar observation data in the three-dimensional unified rectangular coordinate system to obtain the observation track parameters of the main station radar A and the observation track parameters of the auxiliary station radar B; Step 6: Calculate the time alignment points of the main radar A and the auxiliary radar B according to the observation point time of the auxiliary radar B; The time registration point of the auxiliary radar B is converted into three-dimensional polar coordinate registration data centered on the radar station site, expressed as {(R B ' i ,θ B ' i ,h B ' i ), i=1,2,…m}; Step 7: Set the relative system errors of the range, azimuth and elevation of the auxiliary radar B as (ΔR, Δθ, Δα) respectively, use the time registration point data of the main radar A and the auxiliary radar B to construct a relative system error joint solution equation group with (ΔR, Δθ, Δα) as unknowns, and use the quasi-Newton method for finding a set of real roots of the nonlinear equation group to solve the equation group; If the equations have no solution, go to step 2, select the other straight track observation data of the main radar A and the auxiliary radar B on the same air target, and start a new round of system error estimation; if the equations have a solution, let the solution of the equations be (Dr, Da, Dp); Step 8: Perform singular value judgment on Dr, Da, and Dp respectively, incorporate the non-singular values ​​into the relative system error lists ER, EA, and EP of the auxiliary station radar B, and organize the error lists; Step 9: If the list length is greater than or equal to the basic statistical number PN, the mean of the list elements is calculated and output as the corresponding relative system error cumulative estimation result; Otherwise, it is considered that the corresponding relative system error estimation conditions are not yet met; Step 10: Go to step 2, select the straight track observation data of the main station radar A and the auxiliary station radar B on the same air target, and start a new round of system error estimation.

2. The three-coordinate air radar relative system error accumulation estimation method as claimed in claim 1, characterized in that: The step 1 comprises the following steps: Step 1.1: Create a list ER with real number element type to store the relative system error estimates of the auxiliary radar B. Step 1.2: Create a list EA with real number element type to store the relative system error estimates of the auxiliary station radar B. Step 1.3: Create a list EP with real number element type to store the relative system error estimates of the auxiliary station radar B. The step 2 comprises the following steps: Step 2.1: Select the target observation data reported by the primary and secondary radars at the same time when the aerial target is on a straight track; the number of observation data of each radar is not less than 10 points; "same time" means that the time difference of the first point of the primary and secondary radar observation data is not greater than 1 radar detection cycle T, and the time difference of the last point is not greater than T; Step 2.2: The observation data of the selected master radar A is: {(t SAi ,R Ai ,θ Ai ,h Ai ), i=1,2,…n}; Among them, (t SAi ,R Ai ,θ Ai ,h Ai ) indicates t SAi The target distance R measured by the master radar A at the moment Ai 、Directionθ Ai and altitude h Ai , n represents the number of observation data points of the master radar A; Step 2.3: The observation data of the selected auxiliary station radar B is: {(t SBj ,R Bj ,θ Bj ,h Bj ), j = 1, 2, ... m}; where (t SBj ,R Bj ,θ Bj ,h Bj ) indicates t SBj The target distance R measured by the auxiliary station radar B at the moment Bj 、Directionθ Bj and altitude h Bj , m represents the number of observation data points of the auxiliary station radar B; and |t SB1 -t SA1 |≤T,|t SBm -t SAn |≤T.

3. The three-coordinate air radar relative system error accumulation estimation method as claimed in claim 2, characterized in that: The step 3 comprises the following steps: Step 3.1: The observation data {(t SAi ,R Ai ,θ Ai ,h Ai ), i = 1, 2, ... n} into a three-dimensional rectangular coordinate {(t SAi ,x Ai ,y Ai ,z Ai ), i = 1, 2, ... n}; where: Step 3.2: Replace {(t SAi ,x Ai ,y Ai ,z Ai ), i = 1, 2, ... n} into three-dimensional unified rectangular coordinates {(t SAi ,x SAi ,y SAi ,z SAi ), i = 1, 2, ... n}; where: (X SA , Y SA , Z SA ) represents the three-dimensional rectangular coordinates of the master radar A in the unified coordinate system, that is, the site coordinates of the master radar A.

4. The three-coordinate air radar relative system error accumulation estimation method as claimed in claim 3, characterized in that: The step 4 comprises the following steps: Step 4.1: The observation data {(t SBj ,R Bj ,θ Bj ,h Bj ), j = 1, 2, ... m} into a three-dimensional rectangular coordinate {(t SBj ,x Bj ,y Bj ,z Bj ), j = 1, 2, ... m}; in: Step 4.2: Replace {(t SBj ,x Bj ,y Bj ,z Bj ), j = 1, 2, ... m} into three-dimensional unified rectangular coordinates {(t SBj ,x SBj ,y SBj ,z SBj ), j = 1, 2, ... m}; where: (X SB , Y SB , Z SB ) represents the three-dimensional rectangular coordinates of the auxiliary station radar B in the unified coordinate system, that is, the site coordinates of the auxiliary station radar B.

5. The three-coordinate air radar relative system error accumulation estimation method as claimed in claim 4, characterized in that: The step 5 comprises the following steps: Step 5.1: Use the straight line parameter estimation model to estimate the three-dimensional straight line parameters of the observation data of the master radar A, and obtain the observation track parameters of the master radar A (k AX ,d AX )(k AY ,d AY )(k AZ ,d AZ ), comprising the following steps: Step 5.1.1: Transform the X-axis observation data {(t SAi , x SAi ), i = 1, 2, ... n} is abbreviated as: {(x i ,y i ), i=1,2,…n}; replace {(x i ,y i ), i = 1, 2, ... n} is assigned to the structure array XY, the array length is n, the array element is a structure, and the structure member is (x, y); call the straight line trajectory parameter estimation function XYT_TO_kb (n, XY, k, d) to obtain the optimal straight line trajectory parameter (k AX ,d AX ) = (k, d); Step 5.1.2: Convert the Y-axis observation data {(t SAi ,y SAi ), i = 1, 2, ... n} is abbreviated as: {(x i ,y i ), i=1,2,…n}; replace {(x i ,y i ), i = 1, 2, ... n} are assigned to the structure array XY, the array length is n, the array element is a structure, and the structure member is (x, y); the straight line trajectory parameter estimation function XYT_TO_kb (n, XY, k, d) is called to obtain the best straight line trajectory parameter (k AY ,d AY ) = (k, d); Step 5.1.3: Transform the Z-axis observation data {(t SAi , z SAi ), i = 1, 2, ... n} is abbreviated as: {(x i ,y i ), i=1,2,…n}; replace {(x i ,y i ), i = 1, 2, ... n} are assigned to the structure array XY, the array length is n, the array element is a structure, and the structure member is (x, y); the straight line trajectory parameter estimation function XYT_TO_kb (n, XY, k, d) is called to obtain the optimal straight line trajectory parameter (k AZ ,d AZ ) = (k, d); Step 5.2: Use the straight line parameter estimation model to estimate the three-dimensional straight line parameters of the observation data of the auxiliary station radar B, and obtain the auxiliary station radar observation track parameters (k BX ,d BX )(k BY ,d BY )(k BZ ,d BZ ), comprising the following steps: Step 5.2.1: Substitute the X-axis observation data {(t SBj , x SBj ), j = 1, 2, ... m} is abbreviated as: {(x j ,y j ), j = 1, 2, ... m}; {(x j ,y j ), j = 1, 2, ... m} are assigned to the structure array XY, the array length is m, the array element is a structure, and the structure member is (x, y); the straight line trajectory parameter estimation function XYT_TO_kb (m, XY, k, d) is called to obtain the optimal straight line trajectory parameter (k BX ,d BX ) = (k, d); Step 5.2.2: Substitute the Y-axis observation data {(t SBj ,y SBj ), j = 1, 2, ... m} is abbreviated as: {(x j ,y j ), j = 1, 2, ... m}; {(x j ,y j ), j = 1, 2, ... m} is assigned to the structure array XY, the array length is m, the array element is a structure, and the structure member is (x, y); the straight line trajectory parameter estimation function XYT_TO_kb (m, XY, k, d) is called to obtain the optimal straight line trajectory parameter (k BY ,d BY ) = (k, d); Step 5.2.3: Submit the Z-axis observation data {(t SBj , z SBj ), j = 1, 2, ... m} is abbreviated as: {(x j ,y j ), j = 1, 2, ... m}; {(x j ,y j ), j = 1, 2, ... m} are assigned to the structure array XY, the array length is m, the array element is a structure, and the structure member is (x, y); the straight line trajectory parameter estimation function XYT_TO_kb (m, XY, k, d) is called to obtain the optimal straight line trajectory parameter (k BZ ,d BZ )=(k,d).

6. The three-coordinate air radar relative system error accumulation estimation method as claimed in claim 5, characterized in that: The implementation process of the function XYT_TO_kb(n, XY, k, d) includes the following steps: Step X.1: Initialize the function, define variables tx=0, tx2=0, ty=0, ty2=0, txy=0, ii=0; Step X.2: Accumulate the x member value of the element with subscript ii in the XY array into the variable tx; square the x member value of the element with subscript ii in the XY array and accumulate it into the variable tx2; accumulate the y member value of the element with subscript ii in the XY array into the variable ty; square the y member value of the element with subscript ii in the XY array and accumulate it into the variable ty2; multiply the x and y member values of the element with subscript ii in the XY array and accumulate it into the variable txy; Step X.3: Let ii = ii + 1. If ii < n, go to Step X.2; otherwise, go to Step X.4; Step X.4: Let: a1 = tx / n, a2 = tx2 / n, b1 = ty / n, b2 = ty2 / n, c0 = txy / n; Step X.5: Let: aa = c0 - a1 * b1, bb = a2 - b2 - a1 * a1 + b1 * b1, cc = a1 * b1 - c0; Step X.6: Order: d1 = b1 - a1 * k1; d2 = b1 - a1 * k2; Step X.7: Let: Among them, XY[0].x represents the x member value of the element with subscript 0 in the XY array, XY[0].y represents the y member value of the element with subscript 0 in the XY array, and |...| represents taking the absolute value; Step X.8: If L1 > L2, take k = k2, d = d2; otherwise, take k = k1, d = d1; Output k and d as parameters, and the function runs to completion.

7. The relative system error accumulation estimation method of three-coordinate air radar according to claim 5 or 6, characterized in that: The said Step 6 includes the following steps: Step 6.1: According to the observation time of the auxiliary radar B, calculate the three-dimensional rectangular coordinates {(x′ SAi ,y′ SAi ,z′ SAi ), i = 1, 2, ... m}; where: x′ SAi =k AX *t SBi +d AX ,y′ SAi =k AY *t SBi +d AY ,z′ SAi =k AZ *t SBi +d AZ ; Step 6.2: According to the observation point time of the auxiliary radar B, calculate the three-dimensional rectangular coordinates {(x′ SBi ,y′ SBi ,z′ SBi ), i = 1, 2, ... m}; where: x′ SBi =k BX *t SBi +d BX ,y′ SBi =k BY *t SBi +d BY ,z′ SBi =k BZ *t SBi +d BZ ; Step 6.3: Replace {(x′ SBi ,y′ SBi ,z′ SBi ), i = 1, 2, ... m} into a three-dimensional polar coordinate {(R′ Bi ,θ′ Bi ,h B ' i ), i=1,2,…m}; the method is: replace {(x′ SBi ,y′ SBi ,z′ SBi ), i = 1, 2, ... m} is assigned to the structure array XYZ, the array length is m, the array element is a structure, and the structure member is (x, y, z); (X SB , Y SB , Z SB ) is assigned to (XO, YO, ZO); the function XYZ_TO_RAh(m, XYZ, XO, YO, ZO, RAh) is called to convert the rectangular coordinates to polar coordinates to obtain the three-dimensional polar coordinate array RAh centered on the auxiliary radar B; then the members (RR, AA, hh) of the array RAh with the element subscript i-1 are assigned to (R′ Bi ,θ′ Bi ,h′ Bi ); where: RAh array length is m, array elements are structures, structure members are (RR, AA, hh), representing distance, azimuth and altitude values.

8. The three-coordinate air radar relative system error accumulation estimation method as claimed in claim 7, characterized in that: The implementation process of the function XYZ_TO_RAh(m, XYZ, XO, YO, ZO, RAh) includes the following steps: Step Z.1: Let: ii = 0; Step Z.2: Let: xx = XYZ[ii].x – XO, yy = XYZ[ii].y – YO, zz = XYZ[ii].z – ZO Among them, XYZ[ii].x represents the value of member x of the element with subscript ii in the structure array XYZ; XYZ[ii].y represents the value of member y of the element with subscript ii in the structure array XYZ; XYZ[ii].z represents the value of member z of the element with subscript ii in the structure array XYZ; Step Z.3: If yy is equal to 0, go to Step Z.4; otherwise, go to Step Z.5; Step Z.4: If xx is greater than or equal to 0, assign π / 2 to RAh[ii].AA; otherwise, assign -π / 2 to RAh[ii].AA; Go to Step Z.6; Among them, RAh[ii].AA represents the value of member AA of the element with subscript ii in the structure array RAh; Step Z.5: Assign arctan(xx / yy) to RAh[ii].AA; where arctan() represents the arctangent function; Step Z.6: If RAh[ii].AA is less than 0, then assign RAh[ii].AA + 2π to RAh[ii].AA; Step Z.7: Assign RAh[ii].RR to Step Z.8: Assign XYZ[ii].z to RAh[ii].hh; Step Z.9: Let ii = ii + 1; if ii < m, go to Step Z.2; otherwise, output the structure array RAh as a parameter and the function runs to completion.

9. The three-coordinate air radar relative system error accumulation estimation method as claimed in claim 8, characterized in that: The said Step 7 includes the following steps: Step 7.1: Order: A x1i =R′ Bi sinθ′ Bi cosα′ Bi ,A x2i =-R′ Bi cosθ′ Bi cosα′ Bi ,A x3i =-sinθ′ Bi cosα′ Bi ,A x4i =cosθ′ Bi cosα′ Bi , A x5i =R′ Bi sinθ′ Bi sinā′ Bi ,A x6i =-R′ Bi cosθ′ Bi sinā′ Bi ,A x7i =-sinθ′ Bi sinā′ Bi ,A x8i =cosθ′ Bi sinā′ Bi ; A y1i =-A x2i =R′ Bi cosθ′ Bi cosα′ Bi ,A y2i =A x1i =R′ Bi sinθ′ Bi cosα′ Bi ,A y3i =-A x4i =-cosθ′ Bi cosα′ Bi , A y4i =A x3i =-sinθ′ Bi cosα′ Bi ;A y5i =-A x6i =R′ Bi cosθ′ Bi sinα′ Bi ,A y6i =A x5i =R′ Bi sinθ′ Bi sinα′ Bi , A y7i =-A x8i =-cosθ′ Bi sinā′ Bi ,A y8i =A x7i =-sinθ′ Bi sinā′ Bi ; A z1i =R′ Bi sinα′ Bi ,A z2i =-sinα′ Bi ,A z3i =-R′ Bi cosα′ Bi ,A z4i =cosα′ Bi ; Step 7.2: Set the relative systematic errors of the distance, azimuth, and elevation of the slave station radar B as (ΔR, Δθ, Δα) respectively, and construct the following joint solution equations for the relative systematic errors with (ΔR, Δθ, Δα) as unknowns: Where: Step 7.3: Use the quasi-Newton method for finding a set of real roots of a non-linear equation system to solve the equation system (1); if the equation system has no solution, go to Step 2, select the straight track line observation data of the master station radar A and the slave station radar B for the same airborne target again, and start a new round of systematic error estimation; if the equation system has a solution, assume the solution of the equation system is (Dr, Da, Dp).

10. The three-coordinate air radar relative system error accumulation estimation method as claimed in claim 9, characterized in that: The said Step 8 includes the following steps: Step 8.1: Conduct singular value judgment on Dr, include the non-singular values in the distance relative systematic error list ER of the slave station radar B, and sort out the error list ER; including the following steps: Step 8.1.1: If the length of the list ER is less than the basic statistical number PN, add Dr to the end of the list ER and go to Step 8.2; otherwise, go to Step 8.1.2; where the basic statistical number PN represents the minimum number of times to obtain the systematic error estimate value, which is set by the user himself; Step 8.1.2: Calculate the sample standard deviation DevR and the mean AveR of the elements in the list ER. If the absolute value of the difference between Dr and the mean AveR is less than MU times the sample standard deviation DevR, add Dr to the end of the list ER; otherwise, go to Step 8.2; where MU is set by the user; Step 8.1.3: If the length of the list ER is greater than the maximum retention number PM, recalculate the mean AveR of the elements in the list ER, and delete the element with the largest absolute value of the difference from the mean AveR in the list ER; where the maximum retention number PM represents the maximum number of times to store the systematic error estimate value, which is set by the user himself, and PM > PN; Step 8.2: Conduct singular value judgment on Da, include the non-singular values in the azimuth relative systematic error list EA of the slave station radar B, and sort out the error list EA; including the following steps: Step 8.2.1: If the length of the list EA is less than the basic statistical number PN, add Da to the end of the list EA and go to Step 8.3; otherwise, go to Step 8.2.2; Step 8.2.2: Calculate the sample standard deviation DevA and the mean AveA of the elements in the list EA. If the absolute value of the difference between Da and the mean AveA is less than MU times the sample standard deviation DevA, add Da to the end of the list EA; otherwise, go to Step 8.3; Step 8.2.3: If the length of the list EA is greater than the maximum retention number PM, recalculate the mean AveA of the elements in the list EA, and delete the element with the largest absolute value of the difference from the mean AveA in the list EA; Step 8.3: Perform singular value judgment on Dp, incorporate non-singular values ​​into the elevation relative system error list EP of the auxiliary station radar B, and organize the error list EP; including the following steps: Step 8.3.1: If the length of the list EP is less than the basic statistical number PN, add Dp to the end of the list EP and go to step 9; otherwise, go to step 8.3.2; Step 8.3.2: Calculate the sample standard deviation DevP and mean AveP of the elements in the list EP. If the absolute value of the difference between Dp and mean AveP is less than MU times the sample standard deviation DevP, add Dp to the end of the list EP; otherwise, go to step 9; Step 8.3.3: If the length of the list EP is greater than the maximum retention number PM, recalculate the mean AveP of the elements in the list EP and delete the element in the list EP with the largest absolute value of the difference with the mean AveP.

11. The relative system error accumulation estimation method of three-coordinate air radar as claimed in claim 10, characterized in that: The step 9 comprises the following steps: Step 9.1: If the length of the distance relative systematic error list ER is greater than or equal to the basic statistical number PN, the mean of the elements in the list ER is calculated and output as the cumulative estimation result of the distance relative systematic error of the auxiliary station radar B; Otherwise, it is considered that the distance relative system error estimation condition of auxiliary station radar B is not met at present; Step 9.2: If the length of the azimuth relative systematic error list EA is greater than or equal to the basic statistical number PN, the mean of the elements in the list EA is calculated and output as the cumulative estimation result of the azimuth relative systematic error of the auxiliary station radar B; Otherwise, it is considered that the relative system error estimation condition of the auxiliary station radar B is not met at present; Step 9.3: If the length of the elevation relative system error list EP is greater than or equal to the basic statistical number PN, the mean of the elements in the list EP is calculated and output as the cumulative estimation result of the elevation relative system error of the auxiliary station radar B; Otherwise, it is considered that the conditions for estimating the relative system error of the auxiliary station radar B are not met at present.

Images