Three-coordinate air radar height measurement value correction method

Abstract

The invention relates to a three-coordinate air radar height measurement value correction method, and belongs to the technical field of multi-radar data fusion. According to the method, a relative system error idea is adopted, observation data of main and auxiliary station radars on a typical air route target are selected, and a pitching relative system error resolving model derived from radar asynchronous measurement data is established through three-dimensional unified rectangular coordinate conversion, linear parameter estimation, time registration point calculation and the like; and the reliability of an error accumulation estimation result is ensured through operations such as singular value elimination. Experimental results show that the follow-up height measurement value of the auxiliary station radar is corrected according to the cumulative estimation result, the average difference of the height values of the track lines is reduced by one order of magnitude, the error correction effect is very ideal, and the requirements for air radar height estimation and fusion precision in engineering practice are met. 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 invention belongs to the technical field of multi-radar data fusion, and in particular relates to a method for correcting three-coordinate air radar altitude measurement values. Background Art

[0002] Three-coordinate air-to-air radars are generally deployed on the ground, and obtain the distance, azimuth and altitude measurement data of aerial targets through mechanical scanning in the horizontal direction (azimuth) and electrical scanning in the vertical direction (pitch angle). In actual applications, the altitude measurement values ​​of aerial targets by two radars often differ by more than a kilometer, which seriously affects the target association recognition and data fusion effects, and at the same time, puts special situation handling work at risk.

[0003] The systematic error in the pitch angle measurement is the main factor for the inaccurate height measurement of three-coordinate air radar. Generally, the systematic error refers to the inherent error of the measured value relative to a certain reference system. If the absolute height value of the target is specified as the reference system, the absolute systematic error is obtained. In the multi-radar air network observation, the absolute height value of the target is mostly unknown. What we can easily obtain is the discrete observation value of the same aerial moving target by different radars. Traditional data calibration methods all attempt to estimate the absolute systematic error in radar height measurement through such measurement data. The simulation test shows that the target height after the systematic error correction depends on the relative height of the single radar measurement value participating in the calculation, and has little correlation with the true height of the target. In actual engineering applications, the absolute systematic error estimation result of radar A will be different due to the different measurement data of radar B or radar C. If it is recognized that the height systematic error calculated by A and B data is the absolute systematic error of radar A, it may cause a contradiction between the difference between the corrected target height values ​​of A and B and the difference between the corrected height values ​​of A and C. Therefore, the idea of ​​calculating the absolute systematic error of height is difficult to grasp in engineering practice.

[0004] In most cases, the air radar network measurement system can only provide observation data of multiple radars on the same target in the same period of time, and usually a set of measurement values ​​(distance, azimuth, altitude) on a typical route (the target maintains a certain altitude and flies in a straight line) can be obtained. Based on such a data environment, we propose: ① Compared with other radars in the network, the observation of a certain radar is accurate. At this time, the height measurement value of this radar (named as the main station) can be used as a true description of the target height. Other radars (named as auxiliary stations) use this as a reference to obtain the height relative system error of the auxiliary station radar relative to the main station radar. For a regional radar network, other radars can be corrected based on the main station. In this way, while achieving consistency in the observation results of each radar, the complexity of the estimation method is simplified, which facilitates engineering implementation. ② Taking the pitch measurement error as the main component of the height system error, a pitch relative system error solution model derived from radar asynchronous measurement data is established. ③ A three-coordinate air radar pitch relative system error cumulative estimation method based on multiple pitch relative system error solution results is proposed, and the influence of abnormal error solution results is shielded with the help of historical experience data. ④ Based on the cumulative estimation results, the altitude measurement values ​​in the subsequent observations of the auxiliary station radar are corrected.

[0005] The present invention selects the observation data of the main and auxiliary station radars on typical route targets, and establishes a pitch relative system error solution model derived from radar asynchronous measurement data through three-dimensional unified rectangular coordinate conversion, linear parameter estimation and time registration point calculation, and ensures the reliability of the cumulative estimation results of the pitch relative system error through operations such as singular value removal. Experimental results show that by using this cumulative estimation result to correct the subsequent measurement of the auxiliary station radar, the degree of splitting of the main and auxiliary station height measurements is reduced by nearly an order of magnitude, and the height correction effect is very ideal. The estimation method provided by the present invention is scientific, the implementation steps of the scheme are reasonable, and the operability and practicality are very strong. Summary of the invention

[0006] 1. Technical issues to be resolved

[0007] The technical problem to be solved by the present invention is how to provide a method for correcting the altitude measurement value of a three-coordinate air-to-air radar so as to achieve the best superposition of the altitude observation results of multiple three-coordinate air-to-air radars on the same target.

[0008] (II) Technical solution

[0009] In order to solve the above technical problems, the present invention proposes a method for correcting three-coordinate air radar altitude measurement values, which comprises the following steps:

[0010] Step 1: Initialize the elevation relative system error list EP of the auxiliary station radar B;

[0011] 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, and 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 ;

[0012] 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};

[0013] 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};

[0014] Step 5: Use the linear parameter estimation model to estimate the parameters of the main station radar A altitude observation data and the auxiliary station radar B three-dimensional observation data to obtain the main station radar A altitude observation track parameters and the auxiliary station radar B three-dimensional observation track parameters;

[0015] Step 6: According to the observation time of the auxiliary station radar B, calculate the height time registration point of the main station radar A and the three-dimensional time registration point of the auxiliary station radar B; calculate the distance between the auxiliary station radar B time registration point and the radar site and the corresponding target height value, expressed as {(R′ Bi ,h′Bi ), i=1,2,…m};

[0016] Step 7: Set the relative elevation system error of the auxiliary radar B to α, use the time registration point data of the main radar A and the auxiliary radar B to construct the relative elevation system error solution equation with α as the unknown, and use the Newton method to find a real root of the nonlinear equation to solve the equation; if the equation has 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 elevation system error estimation; if the equation has a solution, set the solution of the equation to Dp;

[0017] Step 8: Perform singular value judgment on Dp, include non-singular values ​​into the elevation relative system error list EP of the auxiliary station radar B, and sort out EP;

[0018] Step 9: If the length of the elevation relative system error list EP is greater than or equal to the basic statistical number PN, the mean value β of the list EP elements 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 elevation relative system error of the auxiliary station radar B are not met at present;

[0019] Step 10: 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;

[0020] Step 11: Based on the cumulative estimation result of the pitch relative system error output in step 9, correct all subsequent height measurement values ​​of the auxiliary station radar B.

[0021] (III) Beneficial effects

[0022] The present invention proposes a method for correcting the altitude measurement value of a three-coordinate airborne radar. The present invention adopts the idea of ​​relative system error. By selecting the observation data of the main and auxiliary station radars on typical route targets, a pitch relative system error solution model derived from radar asynchronous measurement data is established through three-dimensional unified rectangular coordinate conversion, linear parameter estimation, and time registration point calculation. The reliability of the error accumulation estimation result is ensured by operations such as singular value removal. Experimental results show that by using this cumulative estimation result to correct the subsequent altitude measurement values ​​of the auxiliary station radar, the average difference of the track height value is reduced by one order of magnitude, and the error correction effect is very ideal, which meets the requirements of airborne radar altitude estimation and fusion accuracy 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 practicality are very strong. BRIEF DESCRIPTION OF THE DRAWINGS

[0023] Figure 1 The main implementation steps of the method for correcting the altitude measurement value of the air radar in the technical solution of the present invention;

[0024] Figure 2 This is a comparison diagram of the altitude measurement values ​​of the first simulated observation track line in an embodiment of the present invention;

[0025] Figure 3 The results of 100 simulation calculations of the relative system error of the radar elevation of the auxiliary station B in the embodiment of the present invention;

[0026] Figure 4 This is a schematic diagram of the XY plane projection of the first simulated observation track line of the verification experiment in an embodiment of the present invention;

[0027] Figure 5 This is a comparison diagram of the radar track height measurement values ​​of the primary and secondary stations simulated in the first verification experiment in an embodiment of the present invention;

[0028] Figure 6 This is a comparison chart of the altitude measurement values ​​of the primary and secondary station track lines and the altitude correction values ​​of the secondary station track lines in the first simulation of the verification experiment in the embodiment of the present invention;

[0029] Figure 7 Comparison diagram of the average value of the height difference between the main and auxiliary stations observed in the verification experiment in the embodiment of the present invention and the average value of the height difference between the main and auxiliary stations after correction.

[0030] Figure 8 The embodiment of the present invention verifies the correction rate of the average value of the difference between the auxiliary station track line altitude and the main station observed track line altitude before and after correction in the experiment. DETAILED DESCRIPTION

[0031] In order to make the purpose, content and advantages of the present invention more clear, the specific implementation methods of the present invention are further described in detail below in conjunction with the drawings and examples.

[0032] To achieve the above-mentioned purpose, the present invention provides a method for correcting the altitude measurement value of a three-coordinate airborne radar, which is applied to the early data preprocessing process of a multi-radar data fusion system. The method is applied to obtain the altitude system error estimation result, and to correct the subsequent altitude measurement value of the auxiliary station, which has a significant promoting effect on improving the accuracy of multi-radar target association recognition, enhancing the accuracy of target state estimation, and reducing the risk of special situation handling. The method comprises the following steps:

[0033] Step 1: Initialize the elevation relative system error list EP of the auxiliary station radar B.

[0034] 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 ,RBj ,θ Bj ,h Bj ), j = 1, 2, ... m}, n represents the number of observation data points of the main station radar A, and 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 ;

[0035] 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}.

[0036] 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}.

[0037] Step 5: Use the linear parameter estimation model to estimate the parameters of the main station radar A altitude observation data and the auxiliary station radar B three-dimensional observation data to obtain the main station radar A altitude observation track parameters and the auxiliary station radar B three-dimensional observation track parameters.

[0038] Step 6: According to the observation time of the auxiliary radar B, calculate the height time registration point of the main radar A and the three-dimensional time registration point of the auxiliary radar B. Calculate the distance between the auxiliary radar B time registration point and the radar site and the corresponding target height value, expressed as {(R′ Bi ,h′ Bi ), i=1,2,…m}.

[0039] Step 7: Set the relative elevation system error of the auxiliary radar B to α, use the time registration point data of the main radar A and the auxiliary radar B to construct the relative elevation system error solution equation with α as the unknown, and use the Newton method to find a real root of the nonlinear equation to solve the equation. If the equation has 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 elevation system error estimation. If the equation has a solution, let the solution of the equation be Dp.

[0040] Step 8: Perform singular value judgment on Dp, include non-singular values ​​into the elevation relative system error list EP of the auxiliary station radar B, and sort out EP.

[0041] Step 9: 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 list EP elements 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 elevation relative system error of the auxiliary station radar B are not met at present.

[0042] Step 10: Go to step 2, select the 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.

[0043] Step 11: Based on the cumulative estimation result of the pitch relative system error output in step 9, correct all subsequent height measurement values ​​of the auxiliary station radar B.

[0044] The step 1 comprises the following steps:

[0045] Step 1.1: Create a list EP with real number elements to store the relative system error estimates of the auxiliary radar B.

[0046] The step 2 comprises the following steps:

[0047] Step 2.1: Select the target observation data reported by the main radar and the auxiliary radar at the same time when the aerial target is on a straight track. The number of observation data of each radar is generally not less than 10 points. "Same time" means that the time difference of the first point of the main and auxiliary radar observation data is not greater than 1 radar detection cycle T (generally 10 or 20 seconds), and the time difference of the last point is also not greater than T.

[0048] 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 ,hAi ) 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.

[0049] 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}. Among them, (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 radar B. And |t SB1 -t SA1 |≤T,|t SBm -t SAn |≤T.

[0050] The step 3 comprises the following steps:

[0051] 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}. Among them:

[0052]

[0053] 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}.

[0054] in:

[0055]

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

[0057] The step 4 comprises the following steps:

[0058] 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}. Among them:

[0059]

[0060] 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}. Among them:

[0061]

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

[0063] The step 5 comprises the following steps:

[0064] Step 5.1: Use the linear parameter estimation model to estimate the parameters of the height observation data of the master radar A and obtain the height observation track parameters of the master radar (k AZ ,d AZ ), comprising the following steps:

[0065] Step 5.1.1: 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}. Assign {(x i , y i ), i = 1, 2, … n} to the structure array XY. The length of the array is n, and the array elements are structures. The structure members are (x, y). Call the linear parameter estimation function XYT_TO_kb(n, XY, k, d) to obtain the optimal straight-line track parameters (k AZ , d AZ ) = (k, d) for the Z-axis measurement of the master station A. Among them, the implementation process of the function XYT_TO_kb(n, XY, k, d) includes the following steps:

[0066] Step X.1: Function initialization. Define variables tx = 0, tx2 = 0, ty = 0, ty2 = 0, txy = 0, ii = 0.

[0067] 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.

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

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

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

[0071] Step X.6: Let:

[0072] d1 = b1 - a1 * k1, d2 = b1 - a1 * k2.

[0073] Step X.7: Let:

[0074]

[0075] 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.

[0076] 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 ends.

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

[0079] 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}. 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). Call the linear parameter estimation function XYT_TO_kb (m, XY, k, d) to obtain the best straight line trajectory parameters (k BX ,d BX )=(k,d).

[0080] 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}. 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). Call the linear parameter estimation function XYT_TO_kb (m, XY, k, d) to obtain the best straight line trajectory parameters (k BY ,d BY )=(k,d).

[0081] 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}.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). Call the linear parameter estimation function XYT_TO_kb (m, XY, k, d) to obtain the optimal straight line trajectory parameters (k BZ ,d BZ )=(k,d).

[0082] The step 6 comprises the following steps:

[0083] Step 6.1: Calculate the coordinates of the height time registration point {z′ SAi ,i=1,2,…m}. Among them: z′ SAi =k AZ *t SBi +d AZ .

[0084] 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}. Among them:

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

[0086] Step 6.3: Calculate {(x′ SBi ,y′ SBi ,z′ SBi ), i = 1, 2, ... m} and the distance from the auxiliary radar station B and the corresponding target height value {(R′ Bi ,h′ Bi ), 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 , ZSB ) Assign the three-dimensional unified rectangular coordinates of the site of the secondary station radar B represented to (XO, YO, ZO). Call the function XYZ_TO_RAh(m, XYZ, XO, YO, ZO, RAh) for converting rectangular coordinates to polar coordinates to obtain the two-dimensional polar coordinate array RAh centered on the secondary station radar B. Then, assign the members (RR, hh) of the element with subscript i - 1 in the array RAh to (R′ Bi , h′ Bi ).

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

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

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

[0090] Where, 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.

[0091] Step Z.3 Assign RAh[ii].RR to

[0092] Step Z.4 Assign RAh[ii].hh to XYZ[ii].z.

[0093] Step Z.5 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.

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

[0095] Step 7.1: Let:

[0096] A z1i = R′ Bi sinα′ Bi , A z2i = -R′ Bi cosα′ Bi .

[0097] Step 7.2: Set the relative system error of the pitch of auxiliary station B to α, and construct the relative system error solution equation with α as the unknown as follows:

[0098]

[0099] Step 7.3: Use Newton's method to solve equation (1) for a real root of a nonlinear equation. If the equation has no solution, go to step 2, select the other straight track observation data of the main station A and the auxiliary station B radar on the same air target, and start a new round of system error estimation. If the equation has a solution, let the solution of the equation be Dp. Among them, the Newton method for finding a real root of a nonlinear equation requires calculating the derivative of the function on the left side of equation (1), and its expression is:

[0100]

[0101] The step 8 comprises the following steps:

[0102] Step 8.1: Perform singular value judgment on Dp, include non-singular values ​​into the elevation relative system error list EP of auxiliary station B, and organize the error list EP. It includes the following steps:

[0103] Step 8.1.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.1.2. The basic statistical number PN represents the minimum number of times to obtain the system error estimate, which is set by the user, and generally PN≥5.

[0104] Step 8.1.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. MU is set by the user, generally 2≥MU≥3.

[0105] Step 8.1.3: If the length of the list EP is greater than the maximum number of retentions PM, recalculate the mean AveP of the elements in the list EP and delete the element in the list EP with the largest absolute difference from the mean AveP. The maximum number of retentions PM represents the maximum number of stored system error estimates, which is set by the user. Generally, PM ≥ 30, and PM > PN.

[0106] The step 11 comprises the following steps:

[0107] Step 11.1: Assume that the data of an arbitrary subsequent observation point of the auxiliary station radar B is: (t SB ,R B ,θ B ,h B ), indicating tSB The target distance R measured by the auxiliary station radar B at the moment B 、Directionθ B and altitude h B .

[0108] Step 11.2: Correction value h′ of target height measurement B The calculation formula is:

[0109]

[0110] Among them, Z SB represents the Z coordinate of the auxiliary radar B in the unified coordinate system, that is, the altitude; β is the element mean of the elevation relative system error list EP obtained in step 9, that is, the cumulative estimation result of the elevation relative system error of the auxiliary radar B.

[0111] Embodiment 1:

[0112] This embodiment specifically describes a three-coordinate air radar altitude measurement value correction method proposed by the present invention, wherein the simulation calculation process is implemented based on C# language and can be applied to the early data preprocessing process of a multi-radar networking system. The embodiment includes the following steps:

[0113] Step 1: Initialize the relative system error list EP of the auxiliary station B. That is, create a list EP with real number elements to store the relative system error estimates of the pitch angles of the auxiliary station B. The C# code is as follows: List <double>EP=newList <double>();

[0114] Step 2: Simulate the observation data of a straight line track of the same aerial target by the main radar A and the auxiliary 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}, n represents the number of observation data points of the main station radar A, and m represents the number of observation data points of the auxiliary station radar B. 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 air target are shown in Table 2; the first simulation observation data of radar A and radar B are shown in Tables 3 and 4, n = m = 25.

[0115] Table 1: Radar basic parameters

[0116] parameter Main station radar A Auxiliary station radar B x-axis coordinate 1580km 1500km y-axis coordinate 2680km 2550km Altitude h 500m 200m Distance system error 0.01km -0.2km Azimuth system error 0.002 degrees 1.8 degrees Pitch system error 0.001 degrees 0.5 degrees Random distance error 0.08km 0.1km Random error of orientation 0.08 degrees 0.1 degrees Random pitch error 0.008 degrees 0.01 degrees Random error of time scale 0.03 seconds 0.05 seconds Measuring cycle 10 seconds 10 seconds Antenna initial phase 260 degrees 50 degrees

[0117] Table 2: Track parameters

[0118] parameter Single value Start time (seconds) 50 seconds Target starting point x-axis coordinate (km) Evenly distributed within 1700±3% Target starting point y-axis coordinate (km) Uniform distribution within 2540±3% Flight altitude (remains constant, m) Evenly distributed within 1500±10% Heading (North is 0 degrees) Evenly distributed within 300±30% Target speed (km / h) Uniform distribution within 1200±5% Track points 25

[0119] Table 3: Observation data of radar A at the main station

[0120]

[0121]

[0122] Table 4: Observation data of auxiliary station radar B

[0123] Serial number <![CDATA[Time t SB (seconds)]]> <![CDATA[Distance R B (km)]]> <![CDATA[Azimuth θ B (degrees)]]> <![CDATA[Height h B (m)]]> 1 54 168.69 89.64 3.07 2 64 165.92 89.12 3.09 3 74 162.30 89.64 2.99 4 84 159.11 88.97 2.98 5 94 155.51 89.33 2.95 6 104 152.60 89.75 2.96 7 114 148.78 89.71 2.93 8 124 145.87 89.21 2.92 9 134 142.34 89.26 2.81 10 144 139.29 89.11 2.85 11 154 136.29 89.00 2.81 12 164 132.98 89.19 2.76 13 174 129.51 88.67 2.70 14 184 126.05 89.20 2.66 15 194 123.12 88.77 2.70 16 204 119.66 89.14 2.65 17 214 116.25 88.80 2.60 18 224 112.84 89.01 2.60 19 234 109.78 89.26 2.61 20 244 106.40 88.81 2.55 21 254 103.36 89.15 2.51 22 264 100.01 89.03 2.50 23 274 95.97 88.21 2.44 24 284 93.43 88.17 2.42 25 294 90.24 88.44 2.41

[0124] 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}, the calculation results of the first simulation observation data are shown in Table 5. It includes the following steps:

[0125] 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}. Among them:

[0126]

[0127] 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}.

[0128]

[0129] (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 (1580km, 2680km, 0.5km).

[0130] Table 5: Three-dimensional unified rectangular coordinates of the observation data of radar A at the master station

[0131]

[0132]

[0133] 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}, the calculation results of the first simulation observation data are shown in Table 6. It includes the following steps:

[0134] 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}. Among them:

[0135]

[0136] 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}. Among them:

[0137]

[0138] (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 (1500km, 2550km, 0.2km).

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

[0140]

[0141]

[0142] Figure 2 The height measurement comparison diagram of the first simulation track line of radar A and B drawn according to Table 5 and Table 6 is given. It can be seen from the figure that due to the pitch system error in the measurement of radar A and B, the height of the two observed track lines for the same target is split by 800 meters to 1450 meters; at the same time, the random error in the pitch angle measurement makes the single height track line appear jagged to varying degrees.

[0143] Step 5: Use the linear parameter estimation model to estimate the parameters of the main radar altitude observation data and the auxiliary radar 3D observation data to obtain the main radar altitude observation track parameters and the auxiliary radar 3D observation track parameters. It includes the following steps:

[0144] Step 5.1: Use the linear parameter estimation model to estimate the parameters of the height observation data of the master radar A and obtain the height observation track parameters of the master radar (k AZ ,d AZ ), comprising the following steps:

[0145] Step 5.1.1: 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}. 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). Call the linear parameter estimation function XYT_TO_kb (n, XY, k, d) to obtain the best observed track line parameters (k AZ ,d AZ )=(k, d). The C# implementation code of the function XYT_TO_kb(n, XY, k, d) is as follows:

[0146]

[0147]

[0148] According to step 5.1, the height observation track parameters (k AZ ,d AZ )=(0.0001,1.5948).

[0149] Step 5.2: Use the linear 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 line parameters (k BX ,d BX )、(k BY ,d BY )、(k BZ ,d BZ ), comprising the following steps:

[0150] 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}. 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). Call the straight line trajectory parameter estimation function XYT_TO_kb (m, XY, k, d) to obtain the optimal straight line trajectory parameter (k BX ,d BX )=(k,d).

[0151] 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}. 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). Call the straight line trajectory parameter estimation function XYT_TO_kb (m, XY, k, d) to obtain the best straight line trajectory parameter (k BY ,d BY )=(k,d).

[0152] 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}. 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). Call the straight line trajectory parameter estimation function XYT_TO_kb (m, XY, k, d) to obtain the optimal straight line trajectory parameter (k BZ ,d BZ )=(k,d).

[0153] According to step 5.2, the observation track parameters (k BX ,d BX )=(-0.3286,1686.5839),(k BY ,d BY )=(0.0034,2551.3949),(k BZ ,d BZ )=(-0.0028,3.2313).

[0154] Step 6: According to the observation time of the auxiliary station radar, calculate the main station radar height time registration point and the auxiliary station radar three-dimensional time registration point. Calculate the distance between the auxiliary station radar time registration point and the radar site and the corresponding target height value, expressed as {(R′ Bi ,h′ Bi ), i=1,2,…m}. It includes the following steps:

[0155] Step 6.1: Calculate the coordinates of the height time registration point {z′ SAi ,i=1,2,…m}. Among them: z′ SAi =k AZ *t SBi +d AZ .

[0156] The coordinates of the height time registration point {z′ SAi , i=1,2,…m} The calculation results are shown in Table 7.

[0157] Table 7: Coordinates of the altitude time registration points of the main station radar A

[0158] Serial number <![CDATA[Time t SB (seconds)]]> <![CDATA[Z coordinate z′ SA (kilometer)]]> 1 54 1.60 2 64 1.60 3 74 1.60 4 84 1.60 5 94 1.61 6 104 1.61 7 114 1.61 8 124 1.61 9 134 1.61 10 144 1.61 11 154 1.61 12 164 1.61 13 174 1.61 14 184 1.61 15 194 1.62 16 204 1.62 17 214 1.62 18 224 1.62 19 234 1.62 20 244 1.62 21 254 1.62 22 264 1.62 23 274 1.62 24 284 1.63 25 294 1.63

[0159] 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}. Among them:

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

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

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

[0163]

[0164] Step 6.3: Calculate {(x′ SBi ,y′ SBi ,z′ SBi ), i = 1, 2, ... m} and the distance from the auxiliary radar station B and the corresponding target height value {(R′ Bi ,h′ Bi ), 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). Call the function XYZ_TO_RAh(m, XYZ, XO, YO, ZO, RAh) to convert the rectangular coordinates to polar coordinates, and get the two-dimensional polar coordinate array RAh centered on the auxiliary radar B. Then assign the members (RR, hh) of the elements in the array RAh with the subscript i-1 to (R′ Bi ,h′ Bi ).

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

[0166]

[0167] The distance between the time registration point of the first simulated observation track line of auxiliary station radar B and the station site of auxiliary station radar B and the corresponding target height value {(R′ Bi ,h′ Bi ), i=1,2,…m} The calculation results are shown in Table 9.

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

[0169]

[0170]

[0171] Step 7: Set the relative elevation system error of auxiliary station B to α, use the time registration point data of the main station and auxiliary station to construct the relative elevation system error solution equation with α as the unknown, and use Newton's method to find a real root of the nonlinear equation to solve the equation. If the equation has no solution, go to step 2, select the other straight track observation data of the main station A and auxiliary station B radar on the same aerial target, and start a new round of altitude system error estimation. If the equation has a solution, set the solution of the equation to Dp. It includes the following steps:

[0172] Step 7.1: Order:

[0173] A z1i =R′ Bi sinα′ Bi , A z2i =-R′ Bi cosα′ Bi .

[0174] According to the first simulated observation data of auxiliary station radar B (see Table 9), A z1 , A z2 The parameter values ​​are shown in Table 10.

[0175] Table 10: Calculation parameters A of the first simulated observation data of auxiliary station radar B z value

[0176]

[0177]

[0178] Step 7.2: Set the relative system error of the pitch of auxiliary station B to α, and construct the relative system error solution equation with α as the unknown as follows:

[0179]

[0180] Step 7.3: Use Newton's method to solve equation (1) for a real root of a nonlinear equation. If the equation has no solution, go to step 2, select the other straight track observation data of the main station A and the auxiliary station B radar on the same air target, and start a new round of system error estimation. If the equation has a solution, let the solution of the equation be Dp. Among them, the Newton method for finding a real root of a nonlinear equation requires calculating the derivative of the function on the left side of equation (1), and its expression is:

[0181]

[0182] In this embodiment, the dnewtf(rp,ff,m) function is used to calculate the function value f(α) and its derivative value f′(α) on the left side of equation (1). rp is a real number that stores the value of the unknown number α; m=25, which represents the number of observation data points of auxiliary station radar B, that is, the number of time registration points; ff is a one-dimensional array with a length of 2, which returns the values ​​of f(α) and f′(α). The C# implementation code of the dnewtf() function is as follows:

[0183]

[0184] The implementation function of Newton's method for finding a real root of a nonlinear equation in this embodiment is named dnewt_H(rp,eps,kk,m). Among them, eps=0.0000001, indicating the control accuracy requirement; kk=3000, indicating the maximum number of iterations; m=25, indicating the number of observation data points of auxiliary station radar B, that is, the number of time registration points; rp is a real number, storing the initial value of the unknown number α 0.0; the function return value is the actual number of iterations, and a real number solution of the equation is stored in rp. The C# implementation code of the dnewt_H() function is as follows:

[0185]

[0186] In this embodiment, according to the first simulated observation data, the dnewt_H() function actually iterates and calculates three times, and the solution of the equation is 0.008701, which is assigned to the variable Dp, indicating that the relative system error of the pitch obtained in this round of calculation is 0.4985 degrees. Figure 3 The relative system error of the pitch of auxiliary station B is given by randomly generating track lines according to the parameters in Table 2 and calculating the simulated observation data for 100 consecutive times. The actual number of iterative calculations of the 100 dnewt_H() function does not exceed 3 times.

[0187] Step 8: Perform singular value judgment on Dp, include non-singular values ​​into the relative system error list EP of the height of auxiliary station B, and sort out EP. It includes the following steps:

[0188] Step 8.1: Perform singular value judgment on Dr, include non-singular values ​​into the distance relative system error list ER of auxiliary station B, and organize the error list ER. It includes the following steps:

[0189] Step 8.1.1: If the length of the list EP is less than the basic statistical number PN, Dp is added to the end of the list EP and the process goes to step 9. Otherwise, the process goes to step 8.1.2. The basic statistical number PN is set by the user and represents the minimum number of times to obtain the system error estimate. Generally, PN≥10. In this embodiment, PN is set to 10.

[0190] Step 8.1.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. MU is set by the user, generally 2≥MU≥3. In this embodiment, MU is set to 2.5.

[0191] Step 8.1.3: If the length of the list EP is greater than the maximum number of retention times PM, recalculate the mean AveP of the elements in the list EP, and delete the element in the list EP with the largest absolute difference from the mean AveP. The maximum number of retention times PM is set by the user, indicating the maximum number of system error estimates stored, generally PM≥30, and PM>PN. In this embodiment, PM is set to 45.

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

[0193]

[0194] Among them, the implementation code of the sample standard deviation function CalSTDev() for calculating the list is as follows:

[0195]

[0196]

[0197] Step 9: If the length of the elevation relative systematic error list EP is greater than or equal to the basic statistical number PN, the mean β of the list EP elements is calculated and output as the cumulative estimation result of the elevation relative systematic error of the auxiliary station B. Otherwise, it is considered that the conditions for estimating the elevation relative systematic error of the auxiliary station B are not met at present.

[0198] Step 10: Go to step 2, select the other straight track observation data of the main station A and the auxiliary station B radar on the same air target, and start a new round of system error estimation.

[0199] According to step 2, this embodiment simulates the observation results of 100 different straight track lines of the main station radar A and the auxiliary station radar B on the same aerial target. After processing from step 3 to step 7, the single pitch relative system error calculation result can be solved. Finally, singular value removal, error list sorting and discrimination are performed in step 8 and step 9, and the length of the pitch relative system error list EP is 45, which is greater than the specified basic statistical number PN=15. The output list mean is used as the cumulative estimation result of the pitch relative system error:

[0200] β=0.008716, i.e., 0.4994 degrees, which is close to 99.88% compared with the true value of the elevation angle system error of the auxiliary station radar B set in Table 1 (0.5 degrees).

[0201] Step 11: Based on the cumulative estimation result of the relative elevation system error output in step 9, correct all subsequent height measurement values ​​of the auxiliary station radar B. This includes the following steps:

[0202] Step 11.1: Assume that the data of an arbitrary subsequent observation point of the auxiliary station radar B is: (t SB ,R B ,θ B ,h B ), indicating t SB The target distance R measured by the auxiliary station radar B at the moment B 、Directionθ B and altitude h B .

[0203] Step 11.2: Correction value h′ of target height measurement B The calculation formula is:

[0204]

[0205] Among them, Z SB represents the Z coordinate of the radar station B in the unified coordinate system, that is, the altitude; β=0.008716 is the element mean of the elevation relative system error list EP obtained in step 9, that is, the cumulative estimation result of the elevation relative system error of auxiliary station B.

[0206] Step 12: Verification experiment. In order to further verify the effectiveness of the error estimation, we use the cumulative estimation result of the relative system error of the auxiliary station radar B obtained in step 9, β = 0.008716, and the correction formula in step 11.2 to correct the multi-point height measurement data on other track lines of the auxiliary station radar B. By comparing the change in the interval between the track line height value before and after correction and the height measurement value of the main station radar A, the effectiveness of the error estimation and correction method of the present invention is verified. The verification experiment includes the following steps:

[0207] Step 12.1: Keep the basic parameters of the main radar A and the auxiliary radar B set in Table 1 unchanged, and simulate the track observation data of the main and auxiliary radars for the same air target according to a set of composite (straight line + arc) track parameters of the air target set in Table 11.

[0208] Table 11: Composite track parameters

[0209]

[0210]

[0211] The observation data are generated by simulation according to Table 1 and Table 11. In the first simulation, 111 measurement points are obtained for each of the main and auxiliary radar stations. The projections of the main and auxiliary radar observation tracks on the XY plane are as follows: Figure 4 , the comparison of the radar track height measurement values ​​of the primary and secondary stations is as follows Figure 5 After calculation, the minimum height difference between the radar observation track lines of the main and auxiliary stations is 866.31 meters, the maximum is 2006.79 meters, and the average is 1332.84 meters.

[0212] Step 12.2: According to the cumulative estimation result of the relative system error of the auxiliary station radar B, β = 0.008716, i.e. 0.4997 degrees, the correction formula of step 11.2 is used to correct the height measurement value of the auxiliary station B track point simulated in step 12.1.

[0213] After correction, the measured values ​​of the main and auxiliary station track heights generated by the first simulation are compared with the corrected values ​​of the auxiliary station track heights as shown in the figure below: Figure 6 The corrected altitude value is very close to the main station measurement and is almost indistinguishable. After calculation, the minimum height difference between the main and auxiliary station radar observation track lines after correction is 0.02 meters, the maximum is 81.58 meters, and the average is 20.99 meters, which is 1332.84 meters higher than the average before correction, and the correction rate is 98.43%.

[0214] Step 12.3: Repeat steps 12.1 and 12.2 100 times. After calculation, the average height difference between the main and auxiliary station track lines and the average height difference between the main and auxiliary station track lines after correction are compared as follows: Figure 7 As shown, the correction rate is Figure 8 As shown. According to statistics, the minimum correction rate is 87.76%, the maximum is 98.58%, and the average is 96.64%. In 100 corrections, the correction rate was kept above 90% in 89 times, and the average difference of the altitude value of the main and auxiliary station track line was reduced by an order of magnitude, and the correction effect was very ideal. Therefore, the effectiveness of the method described in the present invention for altitude system error estimation and correction is further verified.

[0215] Without departing from the technical principle of the present invention, the present invention may also be improved and modified, and these improvements and modifications shall also be considered as the protection scope of the present invention. For example, but not limited to the following points:

[0216] (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 elements in the set can change dynamically. 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 a single linked list manner, and in C#, it can be expressed in a generic List manner.

[0217] (2) Regarding the selection of the master radar. The present invention specifies the measurement value of the most accurate air radar in the radar network as the true description of the target altitude, and names it the master. In engineering practice, target altitude measurement values ​​with higher accuracy and data rate can be obtained through air traffic control radars, flight detection equipment, ADS-B equipment, etc., which can be used as the true target altitude. Selecting measurement equipment with higher accuracy and data rate as the master will help improve the error estimation and track correction effects of the present invention, and such improvements and variations should also be considered as the scope of protection of the present invention.

[0218] (3) Regarding the problem of solving the pitch relative system error solution equation. In step 7 of the present invention, the Newton method for finding a real root of a nonlinear equation is used to solve equation (1). Simulation experiments have shown that this method is time-saving and reliable. Exploring the use of other numerical calculation methods to solve equation (1) should be regarded as an improvement and deformation of step 7.

[0219] (4) The cumulative estimation of the relative system error in the present invention is an offline, post-process. Compared with the time complexity of the program implementation, we are more concerned with the accuracy of the valuation results and the effect of online, post-track correction. Therefore, the improvement and improvement of the execution efficiency of the method or embodiment program of the present invention should also be regarded as the scope of protection of the present invention.

[0220] The present invention adopts the idea of ​​relative system error. By selecting the observation data of the main and auxiliary station radars on typical route targets, a pitch relative system error solution model derived from radar asynchronous measurement data is established through three-dimensional unified rectangular coordinate conversion, linear parameter estimation and time registration point calculation, and the reliability of the error accumulation estimation result is ensured through operations such as singular value removal. The experimental results show that the subsequent altitude measurement values ​​of the auxiliary station radar are corrected with this cumulative estimation result, and the average difference of the track height values ​​is reduced by an order of magnitude. The error correction effect is very ideal and meets the requirements for air radar altitude estimation and fusion accuracy 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 practicality are very strong.

[0221] The above is only a preferred embodiment of the present invention. It should be pointed out that for ordinary technicians in this technical field, several improvements and modifications can be made without departing from the technical principles of the present invention. These improvements and modifications should also be regarded as the scope of protection of the present invention.< / double> < / double>

Claims

1. A three-coordinate air radar altitude measurement value correction method, characterized in that: The method comprises the following steps: Step 1: Initialize the 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, and 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 linear parameter estimation model to estimate the parameters of the main station radar A altitude observation data and the auxiliary station radar B three-dimensional observation data to obtain the main station radar A altitude observation track parameters and the auxiliary station radar B three-dimensional observation track parameters; Step 6: According to the observation point time of the auxiliary station radar B, calculate the height time registration point of the main station radar A and the three-dimensional time registration point of the auxiliary station radar B; Calculate the distance between the auxiliary radar B time registration point and the radar station site and the corresponding target height value, expressed as {(R′ Bi ,h′ Bi ), i=1,2,…m}; Step 7: Set the relative elevation system error of the auxiliary radar B to α, use the time registration point data of the main radar A and the auxiliary radar B to construct the relative elevation system error solution equation with α as the unknown, and use the Newton method for finding a real root of the nonlinear equation to solve the equation; If the equation has 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 pitch system error estimation; if the equation has a solution, set the solution of the equation to Dp; Step 8: Perform singular value judgment on Dp, include non-singular values ​​into the elevation relative system error list EP of the auxiliary station radar B, and sort out EP; Step 9: If the length of the elevation relative system error list EP is greater than or equal to the basic statistical number PN, the mean value β of the list EP elements 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 relative system error estimation condition of the auxiliary station radar B is not met at present; Step 10: 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; Step 11: Based on the cumulative estimation result of the pitch relative system error output in step 9, correct all subsequent height measurement values ​​of the auxiliary station radar B.

2. The three-coordinate air radar altitude measurement value correction method as claimed in claim 1, characterized in that: The step 1 comprises the following steps: Step 1.1: Create a list EP with real number elements to store the relative system error estimates of the auxiliary radar B. The step 2 comprises the following steps: Step 2.1: Select the target observation data reported by the main radar and the auxiliary radar 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 between the first point of the main and auxiliary radar observation data is not greater than 1 radar detection cycle T, and the time difference between 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 altitude measurement value correction 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 radar A in the unified coordinate system, that is, the site coordinates of the master radar A.

4. The three-coordinate air radar altitude measurement value correction 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 radar B in the unified coordinate system, that is, the site coordinates of auxiliary radar B.

5. The three-coordinate air radar altitude measurement value correction method as claimed in claim 4, characterized in that: The step 5 comprises the following steps: Step 5.1: Use the linear parameter estimation model to estimate the parameters of the height observation data of the master radar A and obtain the height observation track parameters of the master radar (k AZ ,d AZ ), comprising the following steps: Step 5.1.1: 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} 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 linear parameter estimation function XYT_TO_kb (n, XY, k, d) to obtain the optimal straight line trajectory parameters (k AZ ,d AZ ) = (k, d); Step 5.2: Use the linear parameter estimation model to perform three-dimensional linear parameter estimation on the observation data of the auxiliary station radar B to 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} 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); call the linear parameter estimation function XYT_TO_kb (m, XY, k, d) to obtain the optimal straight line trajectory parameters (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); call the linear parameter estimation function XYT_TO_kb (m, XY, k, d) to obtain the best straight line trajectory parameters (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} 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); call the linear parameter estimation function XYT_TO_kb (m, XY, k, d) to obtain the optimal straight line trajectory parameters (k BZ ,d BZ )=(k,d).

6. The three-coordinate air radar altitude measurement value correction 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 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; 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 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; 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 three-coordinate air radar altitude measurement value correction method as described in claim 5 or 6, characterized in that: The said Step 6 includes the following steps: Step 6.1: Calculate the coordinates of the height time registration point {z′ SAi ,i=1,2,…m};where: 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: Calculate {(x′ SBi ,y′ SBi ,z′ SBi ), i = 1, 2, ... m} and the distance from the auxiliary radar station B and the corresponding target height value {(R′ Bi ,h′ Bi ), 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 two-dimensional polar coordinate array RAh centered on the auxiliary radar B; then the members (RR, hh) of the array RAh with the element subscript i-1 are assigned to (R′ Bi ,h′ Bi ); where: RAh array length is m, array element is a structure, structure member is (RR, hh), representing distance and altitude value; 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 Assign RAh[ii].RR to Step Z.4 Assign RAh[ii].hh to XYZ[ii].z; Step Z.5 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.

8. The three-coordinate air radar altitude measurement value correction method as claimed in claim 7, characterized in that: The said Step 7 includes the following steps: Step 7.1: Let: A z1i =R′ Bi sinα′ Bi ,A z2i =-R′ Bi cosα′ Bi ; Step 7.2: Set the pitch relative system error difference of the secondary station B to be α, and construct the following pitch relative system error solution equation with α as the unknown: Step 7.3: Use Newton's method for finding a real root of a non - linear equation to solve Equation (1); If the equation has no solution, go to Step 2, select the straight - line track observation data of the main station A and the secondary station B radar for the same airborne target again, and start a new round of system error estimation; if the equation has a solution, let the solution of the equation be Dp; Among them, Newton's method for finding a real root of a non - linear equation needs to calculate the derivative of the left - hand function of Equation (1), and its expression is:

9. The three-coordinate air radar altitude measurement value correction method as claimed in claim 8, characterized in that: The said Step 8 includes the following steps: Step 8.1: Perform singular value judgment on Dp, incorporate non-singular values ​​into the elevation relative system error list EP of the auxiliary station B, and organize the error list EP; including the following steps: Step 8.1.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.1.2; where the basic statistical number PN represents the minimum number of times to obtain the system error estimate, which is set by the user; Step 8.1.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; MU is set by the user; Step 8.1.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; the maximum retention number PM represents the maximum number of stored system error estimates, which is set by the user, and PM>PN.

10. The three-coordinate air radar altitude measurement value correction method as claimed in claim 8, characterized in that: The step 11 comprises the following steps: Step 11.1: Assume that the data of an arbitrary subsequent observation point of the auxiliary station radar B is: (t SB ,R B ,θ B ,h B ), indicating t SB The target distance R measured by the auxiliary station radar B at the moment B 、Directionθ B and altitude h B ; Step 11.2: Correction value h′ of target height measurement B The calculation formula is: Among them, Z SB represents the Z coordinate of the auxiliary radar B in the unified coordinate system, that is, the altitude; β is the element mean of the elevation relative system error list EP obtained in step 9, that is, the cumulative estimation result of the elevation relative system error of the auxiliary radar B.

Images