Effective Automatic Radar Bias Estimation Scheme Using Targets of Opportunity

In multi-radar air surveillance scenario bias in individual radar measurements is an important prerequisite for effective fusion of radar data. Conventional bias computation approaches depend on GPS fitted air sorties for computation of bias in radar measurements. The method proposed in this paper computes radar bias automatically using GPS data available from Automatic Dependent Surveillance Broadcast (ADS-B) measurements from targets of opportunity. The proposed method computes the bias in radar measurements without the need of costly dedicated calibration sorties or multiple radars and thereby reducing the cost and manual effort. This paper provides details about the mathematical frame works of recursive weighted least square based bias estimation algorithm and the associated data conversion formats in different coordinate systems. The proposed method preprocess the correlated radar and AD-B measurements and selects measurement pairs from linear segments for time alignment. This paper provides a detailed scheme of automatic bias estimation process with different parameters involved in controlling the accuracy and periodicity of the bias to be computed. The scheme developed in this paper also includes a bias validation module and a unique averaging method to ensure accuracy and smooth variations across the measurement batches considered for computation. The simulation results obtained in this paper suggest the utility of the proposed approach for practical bias computation applications in a cost effective way. The method developed in this paper is based on insight obtained from analysis of recorded data from field deployed radars.


I. INTRODUCTION
In the modern air surveillance scenario, large surveillance regions are covered with multiple radar systems and the radars employed may have common surveillance regions. The multiple target tracks generated by these radars are send to an information fusion center to create the collective air surveillance picture. The collective picture of the surveillance region is presented as global tracks (estimated kinematic parameters of targets) generated, at a fusion center. The global tracks are generated at fusion center using measurements obtained from different radars (centralized tracking) or by track-to-track fusion of local tracks (distributed tracking). The radar measurements used at the fusion center can have two types of errors, the first one is the random component referred as 'noise', uncorrelated in time and the second one, a systematic component referred as 'bias'. The errors due to uncorrelated noise are minimized by using appropriate The associate editor coordinating the review of this manuscript and approving it for publication was Yilun Shang. filtering mechanism, whereas the errors due to bias can be minimized only by estimating the bias under different measurement modalities. The bias in radar measurements may include scaling and offset biases in range, azimuth and elevation measurements [1]. The quality of the global tracks could deteriorate and could result in track discontinuities if the radar biases are not accounted [2], [3]. Hence for successful multi-radar integration the data from reporting radars are transformed to a common reference frame free of bias or registration error. Because of the huge impact of biases in the final output of the global tracks, bias modeling, estimation and compensation are essential steps in multi-radar target tracking systems. Much research has been carried out in the past for solving the bias estimation problem. The approaches developed for bias estimation can be classified as the techniques applicable in case only the radar data is available or in the case where radar data is available along with reference data.
The most common approach among techniques applicable in case only the networked radar data is available is to use VOLUME 9, 2021 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ an Augmented state Kalman Filter (ASKF) that stack target states and bias parameters into a single vector [4], [5]. The drawback of this approach is high computational load due to increased dimension of stacked vector. In [6], an exact solution for ASKF is given by decoupling the equations for implementation feasibility. In [7] a two-sensor bias estimation problem is discussed using local sensor tracks. A dimensionless likelihood ratio is used for the local sensor track-to-track association and a maximum likelihood estimator is used for bias estimation. An exact solution for bias estimation of two sensors using pseudo measurements is presented [8]. For pseudo measurements, computation Kalman gain from local tracks was assumed to be available which is not usually sent to fusion center in practical systems. In [9] an extension to [8] is proposed, in which a practically feasible algorithm based on tracklet computation [10]- [12] is shown. In the tracklet computation based approach, Kalman gain of local trackers are reconstructed and a sequential update for the fusion of tracks is presented after correcting the local tracks with last updated bias estimate. However these approaches are presented on the assumption of synchronous sensors, but the data reported by sensors are not usually synchronous due to different sampling times. In [13] and its extension in [14], an exact solution for the bias estimation of asynchronous sensors is provided using pseudo measurements. To account for asynchronous sampling times, a proper time slot is defined such that there exist a linear combination of measurements (from the same target) obtained from different sensors in the time slot which is independent of the target state. Then a Least Square (LS) estimator across time slots is defined to compute the bias of asynchronous sensors.
Among the techniques of radar bias computation based on reference data the most common approach is to use the target latitude, longitude and altitude obtained from Global Positioning Sensor (GPS) sensor. In this paper GPS data received using ADS-B (Automatic Dependent Surveillance Broadcast) receiver is used as reference data for computing bias in radar measurement. Legacy systems use methodologies in which registration errors are corrected using a pair of radars or multiple radars which could result in relative registration. ADS-B is a GPS dependent broadcast system in which aircraft with ADS-B transponder report their position to ground stations capable of decoding it. As GPS positions are usually fairly accurate, ADS-B opens the possibility of significant accuracy improvements in the determination of radar biases. The paired GPS data of target and radar are used for computing the bias after coordinate and time alignment. But to get the GPS data of target a dedicated transponder or offline data collection from calibration air sorties are required.
The advent of Automatic Surveillance Broadcast (ADS-B) over Mode-S data link [15], [16] and its reliability also steered its use as reference data for radar bias estimation. The ADS-B transponder provides the target latitude, longitude and altitude in a periodic way and these data collected from targets of opportunity provides an efficient alternative to GPS data collected from dedicated calibration sorties. In [17], radar data is aligned against ADS-B data using genetic algorithm that iteratively computes the bias in radar measurements such that the distance between correlated radar-ADS-B report is minimized. The paper [18] has presented a weighted LS based multi-sensor registration algorithm, taking sensors pairwise, using ADS-B data. In [19], a Probability Hypothesis Density(PHD) based solution for bias estimation of radar without prior data association using ADS-B as reference is proposed. A bias estimator for radar that evaluates offset parameters in range, azimuth measurements and time bias using ADS-B as reference is presented in [20]. But these algorithms require manual intervention in collection of paired data and computation of bias in a fixed periodic interval.
This paper presents an automated method to compute the bias in range, azimuth and elevation measurements obtained from a 3D surveillance radar using ADS-B as a reference. Quality of ADS-B is ensured by indicator Navigation Accuracy Category for Position (NACp) [21] and values NACp > 8 were used such that positional accuracy of ADS-B is limited to less than 30m. The ADS-B positional data of a target is asynchronous with respect to its radar measurements. The ADS-B data and radar measurements are synchronized using time projection method. To avoid positional errors due to time projection, correlated ADS-B data and radar measurements having time stamp difference of less than 1 sec are chosen for bias computation. The proposed method consists of a data prepossessing module that extracts nearly constant velocity portions of target trajectory, the maneuvering section of target trajectories are omitted to avoid any positional errors occur during time projections. A weighted LS based bias estimator and a statistical validation module for the estimated bias parameters are used in this paper. In the proposed approach, the bias is estimated after interpolating radar measurements to the sampling time of ADS-B measurements. The complete method is repeated after every user defined period without manual intervention. This paper also presents a new measurement model for bias estimation and also define a moving average filter to smoothen the bias estimates over time. The proposed method is equally applicable to primary radars as well as Secondary Surveillance Radar(SSR). The automatic bias module is tested and validated with field recorded radar data for primary and secondary radars and also with simulated data. The field radar data results were not presented in the paper due to restrictions on publishing Defense radar data. The simulations are carried out as close as real data to capture the features and advantages of the method developed in this paper.
The rest of the paper is organized as follows: a brief description on the new radar measurement model, and the required coordinate transformations is presented in Section II. Approach of Bias estimation algorithm is detailed in Section III. The computed bias undergoes validation and the methodology of bias validation is detailed in Section IV. The proposed automatic bias update process and periodic computation is detailed in Section V. Simulation results are shown in Section VI. The results obtained were concluded in Section VII.

II. MEASUREMENT MODELS AND COORDINATE SYSTEMS A. RADAR MEASUREMENT MODEL
The Radar measurement model with radar biases and the noise associated with the radar measurements can be modeled as, where R mG , θ mG and ϕ mG are the true range, azimuth and elevation of the target, K is the scale bias in the range (this parameter will capture the change in the bias parameter with range, same can be extended to scale bias in azimuth and elevation and present proposed paper is limited to scale bias in range), R s , θ s and ϕ s are the systematic errors (bias) in range, azimuth and elevation, n rp , n θ p and n ϕp are the measurement noise (random errors) in range, azimuth and elevation. [R mP , θ mP , ϕ mP ] T is the measurement vector in polar coordinates available from the radar.

B. GPS MEASUREMENT MODEL
GPS uses geodetic coordinates consisting of latitude φ mG , longitude λ mG , and height over ellipsoid h mG . To estimate the biases in radar measurements using GPS measurements the spherical coordinates centered in radar position is required.
To obtain the radar centered measurements the earth centered earth fixed (ECEF) coordinates of both the radar and target is computed. The modulus of vector difference between radar and target position in ECEF denoted as || r || provides the GPS measured range, Let (x , y , z ) be the ECEF coordinates of the target GPS measurement and (x , y , z ) be the ECEF coordinates of the radar location, then the vector difference (dx, dy, dz) between the target and the radar location is given by Then the modulo of vector difference between the target and the radar location is the range of the GPS measurement with respect to the radar location and is given by The GPS measured azimuth is computed as the angle between local north axis and the projection over the local radar horizontal plane. The GPS measured azimuth can be obtained as [22], where u x and u y are the local north and local west unitary vectors, obtained from geodetic radar position (φ R , λ R , h r ) in an offline manner as [22], The elevation angle in local coordinates with respect to the radar is given by where The ADS-B data follows geodetic coordinate system and the (3)-(9) are used to convert geodetic coordinate system of targets to radar centered coordinates. In this paper we are considering the GPS biases are negligible and the GPS data is the true data of the targets.

C. GEODETIC COORDINATE SYSTEM
The geodetic coordinate system characterizes a coordinate point near the earths surface in terms of longitude, latitude, and height (or altitude), which are respectively denoted by λ mG , φ mG , and h mG . The longitude measures the rotational angle (ranging from −180 • to 180 • ) between the Prime Meridian and the measured point. The latitude measures the angle (ranging from −90 • to 90 • ) between the equatorial plane and the normal of the reference ellipsoid that passes through the measured point. The height (or altitude) is the local vertical distance between the measured point and the reference ellipsoid. Coordinate vectors expressed in terms of the geodetic frame are denoted with a subscript g, i.e., the position vector in the geodetic coordinate system is denoted by

D. EARTH CENTERED EARTH FIXED (ECEF) COORDINATE SYSTEM
ECEF (Earth-Centered, Earth-Fixed), also known as ECR (Earth Centered Rotational), is a Cartesian coordinate system. It represents positions as an X, Y, and Z coordinate. The point (0, 0, 0) is defined as the center of mass of the Earth, hence VOLUME 9, 2021 the name Earth-Centered. Its axes are aligned with the International Reference Pole (IRP) and International Reference Meridian (IRM) that are fixed with respect to the surface of the Earth, hence the name Earth-Fixed. The following formula are used in converting geodetic to ECEF coordinate system.
Important parameters associated with the geodetic frame include the semi-major axis R Ea , the flattening factor f , the semi-minor axis R Eb , the first eccentricity e, the meridian radius of curvature M E , and the prime vertical radius of curvature N E .

III. BIAS ESTIMATION ALGORITHM
In this section the method of estimating radar bias parameters are detailed In the proposed method of bias estimation as shown in (1) the bias is an additive value added to the radar measurement to compensate the inaccuracies in the radar. The bias estimation algorithms are based on the differences in measurements generated by sensors (primary or secondary and GPS) [23]. From (1) it is evident that the parameters to estimate bias are in 3 dimensions (range, azimuth and elevation) and one pair of measurement will not be sufficient to solve (1) having 4 unknown parameters. To estimate radar bias, standard statistical estimation approaches could be used, such as least squares (LSE [22], [24]), best linear unbiased estimators ( [25]) and Kalman filters (KF [7], [25]). The differences in range, azimuth and elevation between the radar measurements and GPS for different targets are considered to evaluate the bias. As explained in (4)-(9) the bias is computed using the time aligned measurements of radar and GPS. But, in practice the measurements received from two sensor need not be aligned in time, hence a preprocessing module is employed to make the measurements time aligned. To project to common time the preprocessing module need to use the motion model of the target, but for maneuvering target to get the appropriate motion model is difficult task. Hence to reduce the errors in time projection the preprocessing module developed in this paper uses only linear target motion segments.
The preprocessing step is carried out for each track, it includes selecting the measurement pairs (GPS & radar measurements) for the target moving linearly (almost constant velocity). The preprocessing will detect and extract the almost constant velocity segments of the targets and neglect the maneuvering segment to avoid model mismatching errors during time alignment of the measurement pairs. GPS and radar measurements are sent to a correlator algorithm where the targets from radar are associated with the respective GPS data. The correlated measurements (radar and GPS) are used as input to the bias estimation process as shown in FIGURE1. The following steps ensures the quality of input data for bias computation: a. NACp 8 values for quality ADS-B data b. Difference between the time of detection of correlated ADS-B and radar measurements less than 1 sec.
c. Pre-processing module for selecting only nearly constant velocity segments of target trajectory for reducing the time alignment errors.
There are 4 blocks in Fig. 1. The correlated targets are the inputs to the block 1. They are converted from geodetic to ECEF coordinate system. Further the ECEF coordinates will be converted to spherical coordinates in block 2. In block 3 the radar measurements will be time aligned with the GPS measurements. The time aligned radar measurements and GPS measurements will be fed to block 4. The output of block 4 is the radar biases in 3 dimensions range, azimuth and elevation. The detail explanation of the blocks in FIGURE1 is given below.

A. BLOCK-1
The input to the block 1 is GPS measurements and Radar location in geodetic coordinate system. The geodetic coordinates of GPS data and Radar location in latitude, longitude and altitude are converted to Earth centric and Earth fixed (ECEF) coordinate system denoted as (x , y , z ) and (x , y , z ). Given a geodetic point say (10) the ECEF coordinates are calculated using (11) to (16) and is given by where Pe(x e , y e , z e ) is a point in ECEF coordinate system. After the conversion of GPS data and the radar location from geodetic coordinate system to ECEF coordinate system fed to Block 2 for ECEF to spherical conversions.

B. BLOCK-2
The output of Block1 is the ECEF coordinates of GPS data (x , y , z ) and Radar location (x , y , z ). Using these parameters the range, azimuth and elevation of the GPS measurements will be calculated with respect to Radar position. The range, azimuth and elevation are calculated from these two points using (5) to (9). The outputs of block are the spherical coordinates of the GPS targets represented as (R mG , θ mG , ϕ mG ). The spherical coordinates of target location obtained from GPS are used as input to correlation module along with radar measurements in spherical coordinates.

C. BLOCK-3
The measurements obtained from the Radar are time projected so as to time-align with the GPS data in x, y and z based on the kinematics of the targets. The time aligned Cartesian coordinates will be converted to spherical coordinates as range, azimuth and elevation angles.
The input for this Bias estimation block is range, azimuth and elevation of radar measurements and GPS data. The output of the block is range bias, azimuth bias and elevation bias of the respective radar. The bias estimation is carried out using recursive weighted least squares algorithm. The time aligned correlated data samples of Radar measurements and GPS data are the input for the Recursive weighted least square algorithm. The initial bias estimate and its covariance is computed using finite number of correlated radar and GPS measurements by batch processing. The batch processing will provide the initial estimate of the bias and initial error covariance. Below are the mathematical equations involved for using batch least squares estimation.
is a stacked vector of measurements (obtained as difference of radar measurements & ADS-B data using (1)) -of dimension kn z * 1, is a stacked measurement matrix -of dimension kn z * n x , is a measurement noise block diagonal matrix of dimension kn z * kn z and σ 2 R , σ 2 θ , σ 2 ϕ are the measurement error covariance in range, azimuth and elevation respectively. The batch least square estimate is given bŷ Thex(k) is used as the initial estimate (obtained from the computed bias estimate X (T − 1) from previous time instant) of the recursive weighted least squares. The covariance matrix is given by P(k) is used as the initial estimate of the error covariance for recursive weighted least squares.
Recursive Weighted Least Square Estimation [26]: Weighted Least Square Estimation (WLSE) is the variant of Least Squares adopted to compute the Bias correction solutions. The Recursive WLSE is computationally simple when compared to batch WLSE, as the latter requires matrix inversions of higher order. The mathematical steps are described below. Let Z(k + 1), H(k + 1),x(k + 1), denotes the measurement vector (obtained as difference of radar measurements and ADS-B data using (1)), the measurement matrix, and estimated bias corrections at instant.
The measurement matrix from (1) is given as The estimated bias corrections from (1) are given aŝ The Recursion for Covariance Update: where S(k + 1) denotes residual covariance and is given by and P(k +1), P(k) and R denotes the posterior estimated error covariance, prior estimated error covariance and measurement error covariance respectively. and the parameter update gain as W(k + 1) is given by With (23), (24), (22) the recursion for covariance can be rewritten more compactly as in (25) P(k + 1) = P(k) − W(k + 1)S(k + 1)W(k + 1) T (29) The Recursion for the Estimate The updated estimate of the bias is given bŷ The new (updated) estimatex(k + 1) is therefore equal to the previous estimate plus a correction term. This correction term consists of the gain W(k + 1) multiplying the residual difference between the observation Z(k +1) and the predicted value of this observation from the previous k measurements.

IV. BIAS VALIDATION MODULE
The quality of the estimated bias correction solutions can be measured by the Mean Square Error (MSE) or Goodness of fit. The goodness of fit can be approximated as a chi-square statistics. The chi-square statistics is a measure of data samples following the measurement model. The goodness of fit of the estimated parameters can be evaluated as given in the (31) where z(k),x(k) and R are the measurements, estimated vector, Measurement noise covariance respectively of a Least Squares problem. If the measurement noise is assumed to be Gaussian distributed, then J is Chi-square distributed with Nn z − n x degrees of freedom, where N is the number of measurements, n z and n x are the dimensions of measurements and estimated vector respectively [26]. The estimated vector is valid [22] or follows the bias measurement model if where c is obtained from Chi-square table such that the probability of a Chi-square random variable with Nn z − n x degrees of freedom exceeding it is α (usually 1% or 5%) If the goodness of fit is acceptable two other tests are carried out.
A. TEST FOR REASONABLENESS [27] Let σ R s , σ θ s and σ ϕ s denote estimated standard deviation in Range bias, azimuth bias, elevation bias respectively, σ R , σ θ and σ ϕ denotes the standard deviations of measurement noise. if Then discard the solution for σ R s . Where N is the number of measurements obtained. Max S igm is the threshold obtained statistically. Similar tests are done on θ s as well as ϕ s .

B. TEST FOR SIGNIFICANCE OF PARAMETERS [27]
Parameter significance determines whether the estimated parameter really contribute to the Bias measurement model.
the (34) implies a two sided probability region. The threshold c is obtained from the Gaussian distribution table.
If the (30) does not hold then the estimate of the parameter is statistically insignificant [2], and set R s = 0 (no detection of bias).
Similar tests are done on K , θ s , and ϕ s . Moving average of the bias computed: The bias computation is an automatic process and for every user defined time period T, the bias is calculated and aforementioned process will be carried out. The validated bias and the previous time period validated bias are weighted average to send out the final bias estimate for the current update. The final weighted average update of the bias is summation of bias estimate computed at previous time instant T − 1 and current time instant T is given bŷ where K (T ) is given by where WA is given by where n is the dimension of the measurement vector, N is the number of measurement samples obtained and J is the residual of the goodness fit. With one dimension of measurement the residual J should ideally tend to N and with the increase of dimension of measurements n, J will tend to multiplication of n * N . Thex(T − 1) is the final bias (periodically reported) reported by the Automatic bias computation tool for the previous time instant andx(T ) is the final bias reported for the present time instant.

V. AUTOMATIC BIAS UPDATE PROCESS
The present bias estimate is reported by weighted average of the previous and present bias estimates given by (35) for every time interval T the validation module will be invoked.
The computed bias will undergo bias validation and moving average with previous T − 1 and present time T instants. The final bias estimate will be released or updated to the user. The block diagram of the automatic bias computation is shown in Fig. 2. The Automatic bias computation tool has user defined tunable parameters. These Tunable Parameters: can be classified as, Parameters used for pre processing, Parameters used for bias computation, Parameters used for bias validation. The details of these parameters are provided in Appendix A along with bias computation module GUI in Fig. 15.
In case of networked radar scenarios the bias computation module will have to receive paired measurements from multiple radars. The bias of multiple radars can be estimated using the same technique but by running parallel instincts of the algorithm. The high level block diagram shown in Fig. 3 details the extension of automatic bias computation module to multiple radars. The data from multiple radars are stored based on Radar ID. The bias computation module will compute bias for each radar simultaneously and independently based on the user defined timers 1 and 2.
As shown in Fig. 2 the Automatic bias computation module consists of mainly 4 blocks, the Data Storage, Bias-preprocessing module, Bias Estimation module, and Bias Validation module. The functionality of each of these modules are as follows: A. DATA STORAGE The Data Storage module receives the correlated Radar and ADS-B data over a defined communication protocol. The received data is stored along with the Radar measurement and the ADS-B measurement based on the track number allotted to the target by the Multi sensor tracking module or single sensor tracking module RDP. The data storage continues for every different target with the unique track number, till the target goes out of the coverage of radar.

B. BIAS-PREPROCESSING MODULE
User defines time intervals for bias preprocessing at what periodicity the module should operate. When the timer-2 VOLUME 9, 2021 enables the tracks collected in the data base are passed to bias preprocessing module. The tracks which are not having updates for time interval M * T (say M = 5 and T = time interval of the radar) are considered as dead tracks and sent for processing. The value M is configurable and denoted as timer-2 in Fig. 2. Bias Pre-Processing selects the straight line sections of track, using the velocity and heading of the target provided by the track data. The track segments with nearly constant velocity and nearly steady course are filtered out from the tracks under process. The measurement pairs (ADS-B and Radar measurements) are collected and stored in the bias data base for further processing. The processed tracks are deleted from the target data base.
The functional flow chart of Bias preprocessing is shown in Fig. 4. The track selected is checked for number of scans more say N and in simulation have used 30 scans of data, if available then the track is checked for ADS-B more than 75% correlation. This check is to avoid wrong ADS-B correlations provided by input track data. The track with more than 75% of ADS-B ID will be selected for further processing. The track will be checked for across scans for constant velocity and constant heading segments using Velocity_Limit and heading_Limit parameters. If segments are available then the measurement pairs are stored in data base.

C. BIAS ESTIMATION MODULE
The bias estimation module calculates the bias whenever user defined timer-1 triggers. The module is explained in the Section-III as Bias estimation algorithm.

D. BIAS VALIDATION MODULE
The computed bias from the bias estimation module is passed for Bias validation module. If the computed bias satisfies the bias validation then the bias parameters are triggered to the user. The functionality of the module is explained in Section-IV.

VI. SIMULATION RESULTS
The simulations are carried out with simulated data and the performance is also validated successfully with real recordings. The objective of this section is to show the usefulness of the concepts developed in this paper using the simulated scenarios. The simulations are carried out as close as real data to show the various features of the novel concepts brought out in this paper such as (i). Using the prior bias with the existing set using (34) for improved estimation, (ii). Automatic computation of bias on set time periods with the tunable parameters mentioned in V (iii). With and without validation checks brought out in (33) and (34).
The error is computed between radar measurements with bias compensated and ADSB data and compared with the error between radar measurements without bias compensation and ADSB data. The error between radar measurements (compensated with bias parameters) and ADSB data in range, azimuth and elevation are denoted as RB , θB , ϕB and the error between the radar measurements (uncompensated with bias parameters) and ADSB data in range, azimuth and elevation are denoted as RP , θ P and ϕP .
where R mB , θ mB , ϕ mB are the bias compensated radar measurements and are given by The Error graphs have been plotted between range error with bias compensation ( RB ) and range error without bias compensation ( RP ), azimuth error with bias compensation ( θ B ) and without bias compensation ( θ P ) and elevation error with bias compensation ( ϕB ) and without compensation ( ϕP ) for targets from simulated. In the simulations carried out 2D scenario is considered.

A. SIMULATION SETUP
The block diagram of simulation set up is shown in Fig. 6. The multi radar multi target simulator generates the true target trajectories as ADS-B data and are sent on to the Ethernet. Here in simulation the ADS-B data are generated as true target trajectories with almost negligible error and interchangeable used in the further sections. The ADS-B data is simulated in geodetic coordinates consisting of latitude φ mG , longitude λ mG , and height over ellipsoid h mG . The geodetic coordinates are converted into local radar coordinate system consisting of Range R mG , azimuth θ mG and elevation ϕ mG using (5), (6) and (8). The true targets trajectories are received at Multi radar bias measurement generation module where in measurement noise and bias parameters are added to generate bias measurement and is passed on to the Ethernet. The bias measurement and ADS-B data are correlated at Multi radar and ADS-B data correlator. The correlated radar measurements and ADS-B data are received over Ethernet at Automatic bias computation tool. The Automatic bias computation module generates bias results (range bias, azimuth bias and elevation bias) for every T (say 1 hr) time period automatically. The simulation scenario details are provided in Table-1. For simplicity purpose only one 2D radar is considered for simulation. Time aligned ADS-B data and radar measurements are considered for simulation and in real data the time projection module shown in Fig. 1 would do the time alignment..

B. SIMULATION SCENARIO
The simulated scenario consists of 8 targets moving with constant velocity for one hour 15 minutes and taking a maneuver for 15 minutes and again going for a constant velocity segment for another one hour 30 minutes. The polar plot of the scenario is shown on Fig. 6. All the targets start at the same time and end at the same time. The maneuver segment is considered for simulation to have different motion models for target trajectory.
The advantage of large data availability at all azimuths and ranges to improve bias estimation is shown by considering   only 2 targets (out of 8 eight targets) as depicted in Fig. 7. The two targets were simulated at 30 deg and 60 deg respectively. The computed bias parameters considering only 2 targets are shown in Table-2 and considering 8 targets is shown  in Table-3. Table-3 and Table-2 depicts the bias parameters computed automatically with weighted average of the prior bias and present bias obtained by (35) for every one hour and the declares the validity of the computed result based on (32), (33) and (34). The first column in Table-3 and  Table-2 provides the Validity of the computed bias, the sec-  ond column provides the information about the time of bias computation, the third column, fourth column and the fifth column provides the information about bias computed parameters like Range scale bias, Range offset bias and Azimuth offset bias respectively. By comparing Table-3 and Table-2, it can be concluded that the large data availability at all azimuths and ranges will improve the bias estimation parameters.
The error graphs are plotted comparing the true bias and estimated bias, targets positions compared with true position, compensated with bias parameters and uncompensated data. The results show that the bias compensated measurements are overlapping with ADS-B data (true data). Also error between bias compensated measurements and ADS-B is very less compared to the error between uncompensated measurements and ADS-B. The results provided below depict the same.

C. SIMULATED RESULTS
The results obtained for every T time is shown in the Table-3. Here in the simulated case T is considered as 1 hour. The simulation is started at 14:08:45 and the prior bias is considered as zero. The simulated scenario is of 3 hours and the results are obtained at 15:08:46, 16:08:47 3 17:08:48 consecutively. The results obtained were almost matching with the true bias parameters. The results shown in Table-3 depict that after every T computation time the result is getting better and converging to true value. The improvement in bias after the weighted average with prior and present bias is shown in Table-   bias is 98.837859 and after the weighted average the result obtained is 100.175278 which is close enough to the actual bias result compared to the computed result. The result ascertains that use of the prior bias with the existing set using (34) will improve the bias estimation. Similar pattern is observed for azimuth bias and range scale bias. The First column in Table-4 shows the time at which the bias is obtained and WA refers to Weighted averaged result and PRES refers to present obtained bias result without weighted average with previous obtained bias result. At 15:08:46 the prior bias is zero, hence the WA and PRES bias parameters are the same.
The bias computed results obtained with only 2 targets is shown in Table-2. The results shows that range bias parameters range bias and range scale bias are converging to the true value but the azimuth bias is having a variation from the true value. The experiment was conducted to analyze the importance of diversified data required for computing the bias parameters. The 8 targets were distributed at all angles (0 to 360 degrees) and the bias result obtained is better compared to the 2 target scenario where the targets are distributed only in the quadrant 1 as shown in Fig. 7 and results could be seen in Table-3 and Table-2.
The bias pre-processing module selects the constant velocity segments (removes some part of the data in the trajectory which are not having constant velocity segments) as per the limits Velocity Limit and Heading Limit. To show the analysis of bias pre-processing module, the current simulation scenario is extended further for one more hour with accelerating segment. Since the aforementioned scenario of 3 hours was having almost constant velocity segment   the bias pre-processing module would not remove much samples. The accelerating segment was fully removed by the bias pre-processing module and overlap of actual trajectory Vs pre-processed trajectory is shown on Fig. 8.
The validity checks play a vital role in validating the computed bias result. If the computed bias fails any one of the valid test as described in section V , the result will be declared as invalid. For the current simulated scenario the result obtained is valid at for all instants except the accelerating segment. The samples for the accelerating segment  were removed by bias pre-processing module the number of samples collected during accelerating segment were failing to meet the Minimum Number of Input due to which, the bias result obtained for accelerating segment would be Invalid and is shown in Table-5. The invalidity due to failed criteria described in section V was obtained mostly in real time data because of model and measurement noise mismatches. To show that measurement noise mismatch would compute invalid bias, erroneous measurement noise values were used while computing bias. The measurement noise considered for simulation is 30m range and 0.037 deg in azimuth refer Table-1, while computing bias the measurement noise considered is 10m in range and 0.01 deg in azimuth in measurement error covariance R. The result obtained was invalid bias due to failure condition of (32), obtained residual J was 25596.358493 and Chi-square threshold c was 7636.804245. For invalid cases  the current bias obtained will be discarded and the prior bias only will be maintained for further processes.
Error graphs for different targets are given below for trajectories of target 1 and target 2 are shown in Fig. 9, Fig. 10, Fig. 11 and Fig. 12. It is observed that the overall error with bias compensation in range ( RB ) < ( RP ) is less than compared to error in range without bias compensation for all the targets and same is depicted in Fig. 9 and Fig. 11 for target 1 and target 2. Similar pattern is observed for azimuth ( θ B ) < ( θ P ) and is shown in Fig. 10 and Fig. 12.
In Fig. 13 and Fig. 14 it is shown that bias compensated measurements are aligned with ADS-B data. The results of other 6 targets were also have the similar pattern as target 1 and 2. The trajectories of ADS-B data and bias compensated radar measurements are overlapping with each other and uncompensated radar measurements are deviated due to bias. Remaining six targets also have shown similar behavior. The overlapping trajectories could help in validating the output of bias during real time (since there would be no true bias) testing and validation.

VII. CONCLUSION
The radar bias estimation and correction is one of the essential requirements for fusing the tracks received from multiple radar in a multi-radar air surveillance scenario. The computation of radar bias in a conventional mode using calibrated air sorties is costly and time-consuming task. The method developed in this paper details about the innovative way of computing the radar bias using GPS measurements received through ADS-B receiver from targets of opportunity. The methods developed in this paper details the mathematical framework for forming measurement model for 3D radar. The bias estimation algorithm developed in this paper based on recursive weighted least square (RWLS) has the advantage of computing the bias recursively in a computationally efficient way. The batch processing method described in this paper gives a way for initializing the bias estimate in the absence of prior information about the bias. In the case of available prior value this paper has provided a method to combine the prior value with the currently computed bias value.
The major contribution of this paper is the automatic bias computation scheme along with mathematical framework to carry out the bias computation algorithms in an automatic way. The developed method avoids the manual intervention in collecting the data and triggering the computation process and computes the bias effectively in a pre-sated periodicity. The simulation results carried out with different target trajectories shows that the computed bias converges to the true value as more data is received and attains a steady state. The simulated results also brings out the practical utility of the developed module by avoiding manual intervention in selection of the radar -GPS measurement pairs with the help of bias preprocessing module. The validation module developed in this paper finally checks the validity of the computed bias according to Chi Square test. The developed module based on the algorithms described is tested in field and the results are matching with the traditional calibrated air sortie-based exercise.

APPENDIX A GUI OF AUTOMATIC BIAS COMPUTATION
The GUI of Automatic bias computation module is shown in Fig. 15. The GUI shown in Fig. 15 has input details and output details of the bias module. The GUI consists of 5 types of major criteria entry Sections. The developed GUI is capable to select different radars for computing bias. The top push buttons are for selecting radar names from the available 'Radar List' pre entered and saved in the data base of this application. The next two rows are for entering 'Key Parameters' required for selecting the paired radar ADS-B measurements for a given radar. Next 5 rows of data entry field in two columns are for entering 'Primary and Secondary Radar' information. The next two rows are for displaying the 'Bias Results'. The column in right is the 'Display Area' to show the messages to user about the status of the application execution. The detailed list of parameters and their significance are detailed as follows: 1) Radar List • Radar Name: The list of radars considered for computing bias will be available. User can select any radar from the list for checking bias result.
• Show Result: User after selecting the 'Radar Name' can click 'Show Result' for viewing the bias result computed.
• Add New Radar: User can manually add the details of the new radar for which user wishes to compute the bias result. (In general the radars will be added automatically) • Clear Screen: User can remove the data of particular radar from data base. 2) Key Parameters: • Radar ID: Radar ID is the identity of the selected radar. In general the radar ID is the identification of the radar. a) Parameters used for Preprocessing: • Timer-2: Periodicity to trigger the bias preprocessing module for selecting the measurement pairs (ADS-B and radar measurements) (Say 60 (units in minutes)).
• Track_Out_of_Coverage: If the last update of the track is more than the parameter (Typical value is 5*time interval of the radar) then the track will be considered for preprocessing module.
• Heading Limit: The heading Limit used for selecting the constant heading segments say 5 deg/sec.
• Velocity Limit: The velocity Limit used for selecting the constant velocity segments say 20 meter/sec. b) Parameters used for Bias Computation: • Timer-1: Periodicity to trigger bias computation module (Say 12 (units in hours) denoted as 'T').
• Number_of_Measurement_Pairs: If the measurement pairs available in bias data base are greater than the parameter, then the bias estimation Computes the bias if following conditions (Timers) are satisfied along with Number of measurements (Typical value is 800). c) Parameters used for Bias Validation: • Maximum Sigma: Used in Bias Validation module. The typical value for Maximum Sigma is 3.
• Confidence Probability: The probability of confidence for significance test is taken as 90%. The threshold is computed form the Chi-square table. These parameters are useful for tuning the automatic bias computation module for various requirements. The parameters used for tuning bias preprocessing module help in selecting the most ideal pairs of radar and ADS-B measurements so that input error in bias computation is controlled. The periodicity of bias computation is set according to the urgency of obtaining the bias value with a tradeoff accuracy. The tightening of parameters in bias validation module will ensure accuracy in computed bias but with a delay in time in case of absence of sufficient paired measurements.

3) Primary & Secondary radar Info:
• Dimension: Specifies the dimension of the radar. • Timer-1: The timer for computing the bias in hours. • Latitude: Latitude information of the radar in degrees.
• Longitude: Longitude information of the radar in degrees.
• Altitude: Altitude information of the radar in feet.

4) Bias result info for Primary & Secondary radar:
• KFactor: Result of Range scale parameter of the radar. VOLUME 9, 2021 • Range: Result of Range bias of the radar in meters.
• Azimuth: Result of Azimuth bias of the radar in degrees.
• Elevation: Result of elevation bias of the radar in degrees.
• Time of the Result: Time at bias result got published.

5) Display Area:
• Bias Result Message: The preliminary details of the bias result is shown to user for checking the bias result.
• Clear Message Are: User can click and clear the message area of bias result message.
• Plot Data: User can plot the results in graphical form by overlapping the ADS-B data, radar measurements without bias compensation and also radar measurements with bias compensation.