Abstract:
Concurrent determinations of errors for application of corrections to sensor data signals from an orbiting space vehicle to reduce or to eliminate the three components of boresight error. To locate more accurately the geographical position of targets detected by the sensors of an orbiting satellite, the pitch, roll, and yaw boresight errors caused by distortions or misalignments in the focal plane of the sensors are compensated for by comparing known position coordinates of an object (such as a star) to the observed position coordinates. The known position coordinates are converted to sensor coordinates by a series of intermediate coordinate conversions, essentially from the celestial frame of reference to an earth-centered frame of reference to a satellite attitude frame of reference to the satellite frame of reference to the sensor frame of reference. The corrections are applied frequently enough to take into account the change of the satellite&#39;s position relative to the sun to compensate for the thermal effects, inter alia, causing the boresight error. An optimal estimator increases the accuracy of the results.

Description:
GOVERNMENT RIGHTS IN THE INVENTION 
     The United States Government has rights in this invention pursuant to Contract No. F05603-90-C-0002 with the Air Force. 
    
    
     DOCUMENTS INCORPORATED BY REFERENCE 
     U.S. Pat. No. 5,051,751 is incorporated by reference to show other types of Cartesian transformation matrices and uses of Kalman filtering. 
     U.S. Pat. No. 4,823,134 is incorporated by reference to show the implementation of a Kalman filter. 
     BACKGROUND AND SUMMARY OF THE INVENTION 
     This invention relates generally to correction of satellite observation data. It pertains particularly to correcting errors caused by inaccuracies in the focal plane of the satellite&#39;s sensors caused by, inter alia, thermal distortion of the focal plane, employing an optimal estimator, e.g., a Kalman Filter, to estimate the boresight state from the centroid of measurement of known images. 
     To record geodetic survey information, it is necessary that the observation data from observation satellites be as accurate and as precise as possible. Observation data from satellite sensors include errors from several sources. The data is transformed through several frames of reference, some of which rotate with respect to time. Time dependent varying errors in one frame of reference are imposed by static errors in another frame of reference which cannot be corrected by static correction in another coordinate frame of reference. 
     One error source is boresight irregularities caused by inaccuracies in the focal plane, e.g., distortion due to thermal changes as the satellite&#39;s aspect relative to the sun changes. The thermal distortion causes the boresight errors to change with respect to time. The boresight error has three components, measured in angular values--pitch, roll, and yaw. Pitch is rotation about the y-axis. Roll is rotation about the x-axis. Yaw is rotation about the z-axis. 
     Thermal distortion can be at least partially compensated for using thermal control. This approach is expensive and complicated but does not correct for all the temperature-caused errors. 
     The prior art discloses several methods for maintaining a satellite in its correct attitude and position. In this description, it is assumed that the satellite is correctly oriented according to attitude and location. The following prior art disclosures are directed to corrections in satellite attitude or to errors caused by incorrect attitude or position. 
     U.S. Pat. No. 5,051,751 shows the use of Kalman filtering for estimating the position and velocity of a tracked object. A line-of-sight to a moving object (target) is measured and transformed to a derived measurement Cartesian frame of reference having one axis pointing at an estimated target position. A Kalman filter, initially set by position and velocity error states, is updated with transformed sensor measurements to estimate the target&#39;s position and velocity. 
     U.S. Pat. No. 5,184,790 is directed to correction of satellite attitude errors caused by deviations in the inclination of the satellite orbit from its normal orbit plane. Although the error can be determined from observations of known points, e.g., a star, the corrections for boresight errors are not addressed or disclosed. 
     U.S. Pat. No. 4,823,134 focuses on problems of accurately directing a shipboard-mounted antenna toward a satellite, for example. Errors are caused by motion of the ship, anomalies in the antenna structure, and distortion in the deck shape caused by temperature changes. It discloses the use of a Kalman filter to correct static and time-varying errors to compensate for errors in several of the coordinate frames of reference which the data traverses. Actual point direction is determined by comparison to a reference which is not a known point but a corrected line-of-sight relative to the ship&#39;s frame of reference. 
     U.S. Pat. No. 5,107,434 is directed to controlling the attitude of a satellite. A polar star sensor with a line of sight parallel to the pitch axis of a geosynchronous orbiting satellite is employed for a two-axis location relative to the pitch axis. An earth sensor is also employed to find the satellite pitch about the pitch axis. The result from processing the polar star data and the data from the earth sensor is used to control torquers to maintain the proper satellite attitude. 
     U.S. Pat. No. 4,896,032 employs a two-wavelength beam pointer system compensator. Boresight alignment drift is measured and then used for beam pointing compensation. It is necessary to generate a beam and requires distance measuring. 
     U.S. Pat. No. 4,884,771 is directed to the calibration of attitude control gyroscopes of a satellite. It solves the problem of determining gyroscopic drift after the gyroscopes have drifted an unknown amount. In such a case, apogee insertion (firing of thrusters during the transfer phase of a geosynchronous earth satellite) will include unacceptable directional errors. A single two-axis sensor observes a radiating reference object, e.g., the sun or a star. Calibration is performed in two stages. In the first stage, a first reference attitude command signal is transmitted to the satellite and the sensor output signals are recorded at the beginning and at the end of the first stage time period. The second stage is similar but with a different attitude command signal being given. The time integrals of the gyroscope output signals produced during the first and second calibration stages (time intervals) can be used to find the gyroscopic drift as described in the patent. 
     U.S. Pat. No. 5,177,686 discloses a method of locating stars to be observed for determining the location of a spacecraft where it cannot be predicted which stars will appear in the visual field of the sensors. It selects patterns of stars expected to be in the visual field and matches them to stored templates to determine which stars are being observed in the field. 
     U.S. Pat. No. 4,944,587 is directed to detecting attitude information of a spacecraft by comparing star image signals such as brightness and position with stored information for generating an attitude signal substantially corresponding to the orientation of the sensors. 
     U.S. Pat. No. 5,259,577 relates to roll and yaw control of a satellite employing kinetic moments. It senses the attitude of the satellite using a terrestrial sensor and a star sensor for locating Polaris (the North Star). 
     The prior art systems, while useful for the purposes for which they were designed, do not suggest boresight error corrections. The use of optimal estimators such as Kalman Filters is mainly applied to correcting satellite attitude errors. 
     In the invention being disclosed, corrections for boresight errors are determined by ascertaining the position, i.e., on which specific sensing element of a sensor array its image should appear, of a known spatial point, e.g., a star or known earth location.The actual location of the spatial point, i.e., the specific sensor on which its image actually appears, is compared to the predicted location to provide data for error estimation and corrections to be applied to the observations from the satellite. The errors are processed through an optimizing filter, typically a Kalman Filter, to improve the error estimates. 
     In accordance with the invention, the known position coordinates of an object which can be sensed are supplied in earth-centered coordinates which are converted into coordinates in a sensor frame of reference via a series of intermediate (mesne) conversions in other coordinate systems. The observed position coordinates of the known object, averaged to determine the centroid of measurement to select the data to be processed, are compared with the known position coordinates to determine concurrently the roll, pitch, and yaw errors. The errors thus determined are processed by an optimal estimator and applied as corrections to satellite mission observation data. 
     The specific utility of the invention is to make concurrent corrections to satellite survey input data to compensate for boresight errors. This is accomplished by observing a known target and transforming its coordinates to the same frame of reference as the raw input data. The difference between the observed target vector and the known target vector is applied to an optimal estimator, implemented as a Kalman Filter, to calculate the boresight state vector. The latter is then used to correct raw mission data for the boresight errors. 
    
    
     BRIEF DESCRIPTION OF THE DRAWING 
     The invention is described in detail by referring to the various figures which illustrate specific embodiments of the invention, and wherein like numerals refer to like elements. 
     FIG. 1 is a diagram of the sensor frame of reference (IRF) orientation. 
     FIG. 2 is a diagram of the satellite frame of reference (SRF) orientation. 
     FIG. 3 is a diagram of the attitude frame of reference (ARF) orientation with the earth-centered frame of reference (ECF) orientation showing their relationship to one another. 
     FIG. 4 is a block diagram of a boresight calibration apparatus. 
     FIG. 5 is a process flowchart of the boresight calibration. 
     FIG. 6 is a block diagram of a system for practicing the invention. 
     FIG. 7 is a program flowchart of the procedure for processing the input data to calculate the boresight correction state vector. 
     FIG. 8 is a program flowchart of the subroutine incorporating the Kalman Filter implementation of an optimal estimator. 
     FIG. 9 is a program flowchart for correcting the mission data for boresight errors. 
     FIG. 10 is a flowchart showing the details of making boresight corrections. 
     FIG. 11 is a graph of sample data over a twenty-four hour period showing the boresight pitch, roll, and yaw for a three-state implementation of the invention. 
     FIG. 12 is a graph of sample data over a twenty-four hour period showing the boresight pitch, roll, and yaw for a six-state implementation of the invention. 
    
    
     DESCRIPTION OF THE PREFERRED EMBODIMENT 
     Summary of Abbreviations and Symbology 
     ε--elevation angle 
     α--azimuth angle 
     γ z  --θ z  +ψ z   
     ARF--attitude frame of reference 
     ECF--earth-centered frame of reference 
     ID--source identifier 
     IR--infrared (sensors) 
     IRF--infrared (sensor) frame of reference 
     LOS--line of sight 
     NID(0,σ)--normally and independently distributed with mean of zero and standard deviation of σ, i.e., random 
     θ y ,θ x ,θ z  --boresight state variables, viz., boresight pitch, roll, and yaw 
     θ y , θ x , θ z  --boresight state velocity variables representing, respectively, ##EQU1## (See augmented boresight state matrix.) RRF--Representative Return Formation 
     SCM--satellite center of mass 
     SRF--satellite frame of reference 
     i, j, k--unit vectors along the x-, y-, z-axes, respectively, of a right-handed three-dimensional vector coordinate system 
     v--general notation denoting a vector where the components are normally expressed as a sum of scalar values times each of the unit vectors 
       !--encloses a matrix when it is not expressed in full matrix form nor identified as a matrix ##EQU2## E{m}--expected value of m, usually the mean 
     General Overview 
     The following description discloses a calibration filter for an infrared (IR) sensor array aboard a surveillance satellite. The calibration output from the filter is termed a boresight and comprises three Euler angle rotations that are applied to three-dimensional line-of-sight (LOS) vectors. 
     The boresight calibration filter calculations involve four coordinate frames of reference, each of which is based on a system of three mutually perpendicular dextral axes, x, y, and z. Unit vectors along the x-, y-, and z-axes are designated i, j, and k, respectively. 
     One coordinate frame of reference is the IR (sensor) reference frame (IRF). The IRF is referenced to the focal plane in which the sensors are arranged and upon which a surveillance telescope is focussed. 
     A second frame of reference is the satellite frame of reference (SRF). The sensor array focal plane is in the x-y plane of the SRF and the z-axis of the SRF is parallel to the line-of-sight (LOS) of the satellite&#39;s telescope, usually directed toward the earth&#39;s center. 
     A third frame of reference is the satellite attitude frame of reference (ARF) and the fourth is the earth-centered frame of reference (ECF) which is considered fixed. 
     Each of these frames of reference will now be discussed in more detail. 
     The IRF 
     The IRF is referenced to the satellite&#39;s center of mass (SCM) 103 as illustrated in FIG. 1. The T-shaped element 101 represents the sensor focal plane which contains a plurality of sensors, infrared sensors in this description. This description contemplates that the sensors provide a scan of the earth once per rotation of the satellite. The orientation of this focal plane with respect to the satellite is the subject of the calibration filter. The x-axis 105 is parallel to the center line of the focal plane 101 and the z-axis 107 intersects the bottom center point 109 of the focal plane 101. The focal plane 101 rotates about the point 109 while the IRF rotates about the z-axis. The y-axis 109 is perpendicular to the x- and z-axes. The nominal location of each of the sensing elements is specified by azimuth, angle α, and elevation, angle ε. In the IRF, the azimuth is designated α irf  and the elevation, ε irf . Each sensor in the focal plane 101 is uniquely identified by a LOS such as the vector 111 to one of the sensors 102 as shown. 
     Unit vectors to the sensors are given according to 
     
         u.sub.irf =x.sub.irf i.sub.irf +y.sub.irf j.sub.irf +z.sub.irf k.sub.irf 
    
     where 
     
         x.sub.irf =cos (α.sub.irf)sin (ε.sub.irf) ,  (1a) 
    
     
         y.sub.irf =sin (α.sub.irf)sin (ε.sub.irf) ,  (1b) 
    
     and 
     
         z.sub.irf =COS (ε.sub.irf) ,                       (1c) ##EQU3## with ##EQU4## 
    
     The SRF 
     The SRF is illustrated in FIG. 2. It is also centered at the SCM 103. The satellite 201 rotates right-handedly about the z-axis 203. The x-axis 205 is perpendicular to the LOS of the scanning telescope. The y-axis 207 is mutually perpendicular to the x- and z-axes. The azimuth and elevation angles are shown but not used in the SRF which is an intermediate (mesne) frame of reference for converting between the IRF and the ARF. The focal plane 101 is shown as it nominally relates to the z-axis 203 of the satellite. 
     The ARF and ECF 
     The ARF and ECF are illustrated in FIG. 3 to show their relationship to one another. That is, the ARF orientation is derived from the ECF. The ARF z-axis 301 points to the center of the earth 303, i.e., to the origin of the ECF. The ARF x-axis 305 is directed parallel to the earth&#39;s easterly direction. The y-axis 307 is mutually perpendicular to the x- and z-axes. A unit vector 310 to a known target, e.g., a star 309, within the field of view of the satellite having an azimuth aarf and an elevation earf is given as 
     
         u.sub.arf =x.sub.arf i.sub.arf +y.sub.arf j.sub.arf +z.sub.arf k.sub.arf 
    
     where 
     
         x.sub.arf =-sin (α.sub.arf)sin (ε.sub.arf) ,tm (2a) 
    
     
         y.sub.arf =cos (α.sub.arf)sin (ε.sub.arf) ,  (2b) 
    
     and 
     
         z.sub.arf =COS (ε.sub.arf)                         (2c) ##EQU5## with ##EQU6## and ε.sub.arf =cos.sup.-1 (z.sub.arf). 
    
     The ECF is considered the fixed frame of reference in this description and is the base reference frame for all the calibration measurement sources. The measurement sources include earth sources and celestial sources. The earth sources are known locations on or near the surface of the earth expressed in geodetic latitude, longitude, and altitude. The celestial sources are typically stars located by geocentric latitude, longitude, and time. Ephemerides are a tabulation of star locations relative to time. The ephemerides take into account the earth&#39;s rotation and orbit, proper motion, precession, nutation, and diurnal aberration. 
     The ECF x-axis 311 is oriented from the ECF origin at the center of the earth through the crossing of the Greenwich meridian at the equator, i.e., latitude 0° and longitude 0°. The z-axis 315 is oriented through the north pole, i.e., latitude 90° and longitude indeterminate. The y-axis 317 is mutually perpendicular to the x-and z-axes,i.e., latitude 0° and longitude 90° East. The ECF axes therefore rotate with the earth. Unit vectors to all ground sources are computed from geodetic latitude and longitude. 
     The filter calculations depend on the capability to express vectors in each of the frames of reference transforming between frames of reference for rotation purposes only. 
     FIG. 4 is a block diagram of an implementation of the invention. The Right Ascension (RA) and declination (dec) of one or more stars in an ephemeris 401 are converted by a CCE (conversion from celestial parameters to ECF) process. The ECF coordinates 403 are transformed to ARF coordinates 405 by a transformation process TEA, described below as a transformation matrix T erf →arf. 
     The ARF coordinates 405 are transformed to SRF coordinates 407 by a transformation process TAS, described below as a transformation matrix T arf →srf. The SRF coordinates 407 are transformed to IRF coordinates 409 by a transformation process TSI, described below as a transformation matrix T srf →irf. 
     Data returns 411 of a known object from a geodetic survey satellite are processed to generate boresight data in the IRF which are applied to one input of a difference (comparator) unit 417, the other input signals being the IRF coordinates 409 of the predicted location of the image. The output signals from the difference unit 417 are errors which are supplied as input data to a Kalman Filter 421. The Kalman Filter updates the boresight&#39;s state vector with terms θ y , θ x , and θ z  (pitch, roll, and yaw, respectively). 
     The boresight state vector and the attitude state vector, P A , R A , and Y A , are applied as input signals to a processor 415. The raw data returns representing mission data are processed by the processor 415 to supply mission data corrected for boresight errors. 
     In FIG. 5, the position coordinates of a known object to be observed, e.g., a star within the observation field of the satellite telescope, are supplied. The coordinates can be computed or can be taken from stored catalog of coordinates as represented by a process block 501. The coordinates are typically updated every 30 seconds from a list of candidates for the observation. 
     The star returns are collected, block 503, but the image of the target is not necessarily confined to a single sensor. The sensor image overlaps several of the sensor elements, producing a blurred image. That is, the star return data forms a cluster of images on the sensor array. To obtain an accurate position of the image of the known object, the returns are averaged to determine a centroid of measurement, block 505. The Representative Return Formation determines a centroid of measurement of the image, effectively cleaning up the input data and providing performance measure. The Kalman Filter optimizes the measurement by predicting the centroid of the image based on present measurements and the current state. 
     The representative returns are in terms of an augmented averaged azimuth and elevation. For greater accuracy averaging, the ε values are weighted as explained below according to intensity. The returns are used for performance measure and to collect raw data for an optimal estimator. 
     The output values (raw returns) from the process block 505 are supplied one at a time to a residual computation 507 that compares and supplies the error between the measured data and the known data of the known object prediction. The known data, supplied by the ephemeris 501 is converted from a celestial or earth coordinate frame of reference to the sensor frame of reference by mesne conversions 511. The other input information to the conversion process 511 include the satellite attitude state, i.e., the pitch, roll, and yaw terms of the satellite attitude, and the boresight state, i.e., the boresight pitch, roll, and yaw from an optimal estimator 509. 
     The optimal estimator 509 is preferably implemented as a Kalman Filter because the latter accepts multiple input data and supplies multiple output state data. The input to the optimal estimator 509 are the error data from the residual (error) computations 507. An optimal multivariable estimator, such as a Kalman Filter, basically employs iterative process for optimizing linear measurement estimates. 
     The discrete process equation for a Kalman is usually portrayed as 
     
          x(n+1)!= Φ! x(n)!+ Γ! u(n)!+ Γ.sub.1 ! w(n)! 
    
     with measurement  y(n)!= H! x(n)!+ v(n)! 
     where 
      x!=state vector, 
      u!=control vector, 
      w!=process noise NID(0,σ), 
      v!=measurement noise NID(0,σ), 
      Φ!=discrete system matrix, 
      Γ!=discrete control input matrix, 
      Γ 1  !=discrete process noise, and 
      H!=continuous output matrix. 
     The basic ideas of a Kalman Filter are that it combines previous estimates with the current measurement, and it uses the history of data as 
     
          x(n)!= Φ! x(n-1)!+ Γ! u(n-1)! 
    
     In the present application, the output from the Kalman Filter, i.e., the x state vector, is the boresight state, θ y , θ x , and θ z . The significance of the states are explained in more detail below. 
     The refined boresight state data is then used to correct mission returns, i.e., the observation data which the satellite was intended to collect. The corrected mission returns are computed from the raw data returns, the boresight state vector from the optimal estimator 509, and the converted satellite attitude vector. The corrected mission returns are converted from the sensor frame of reference to the attitude frame of reference in the correction process 515. 
     The details of the conversions (or transformations) among the different frames of reference are described next. 
     Transformation from ECF to ARF coordinates 
     The ECF to ARF transformation is accomplished by defining a set of unit vectors in the ECF that are parallel to the corresponding ARF axes. The vector R s  (with magnitude r s ) is an ECF radius vector to the SCM with coefficients a, b, and c. These are known values since they represent the position of the satellite with respect to the ECF coordinates. The value of r s  is √a 2  +b 2  +c 2  . A magnitude r xy  is |R s ,xy | and equal to √a 2  +b 2  , that is, the magnitude of the x and y ECF components of the R s  vector. The quantities r s  and r xy  are scalars. 
     The transformed components are, therefore, ##EQU7## 
     The ECF coordinates in the ARF are the vector dot products of each of the basis unit vectors. To convert a vector in the ECF to a vector in the ARF, The unit vector v erf  is found as illustrated in equation (2). The transform matrix is denoted T ecf →arf. ##EQU8## 
     To transform a vector from the ECF to the ARF, the inverse of T ecf →tarf can be used. (T ecf →arf is orthogonal so its inverse is equal to its transpose.) Therefore, T arf →ecf =T T   ecf →arf. 
     Transformation from ARF to SRF coordinates 
     Successive rotations about the pitch (ARF y), the roll (ARF x), and the yaw (ARF z) axes are used to transform from the ARF to the SRF. These rotational transforms are denoted as Ψ y , Ψ x , and Ψ z , respectively, and use the pitch angle ψ y , the roll angle ψ x , and the yaw angle ψ z  : ##EQU9## 
     The transform matrix is given by T arf →srf =Ψ z  Ψ x  Ψ y . Transformation of a vector from the ARF to the SRF is, therefore, 
     
         V.sub.srf =T.sub.arf→srf V.sub.arf =Ψ.sub.z Ψ.sub.x Ψ.sub.y V.sub.arf. 
    
     The inverse matrix transforms SRF coordinates to ARF coordinates. The-transformation matrix T srf →arf is given by Ψ T   y  Ψ T   x  Ψ T   z  since each of the Ψ matrices are orthogonal and the order of multiplication is reversed. 
     Approximations for ARF to SRF transformations can be made for ψ y  and ψ x  when the angles are small, which is usually the case. The value of ψ z  takes on angular measurements from 0 to 2π radians as the satellite rotates about the z-axis. The small angle transformation are simplified by taking a first order approximations of the sine and cosine functions (a value of 1 for the cosine function). The simplified result is ##EQU10## 
     Multiplying the matrices yields a simplified transformation matrix, ##EQU11## 
     Transformation from SRF to IRF Coordinates 
     Transformation from the satellite frame of reference to the sensor IR frame is accomplished by successive rotations about the yaw (SRF z), roll (SRF x), and pitch (SRF y) axes. These rotations are denoted by, Θ z , Θ x , and Θ y , respectively, using the yaw angle θ z , the roll angle θ x , and the pitch angle θ z , ##EQU12## The matrix product of the Θ transform matrices gives the SRF to IRF transformation matrix, T srf →irf =Θ y  Θ x  Θ z  so that a vector in the SRF is transformed to a vector in the IRF by 
     
         V.sub.irf =T.sub.srf→irf V.sub.srf =Θ.sub.y Θ.sub.x Θ.sub.z V.sub.srf. 
    
     The inverse transformation matrix is T irf →srf =Θ T   z  Θ T   x  Θ T   y  since each rotation matrix Θ is orthogonal and the order of matrix multiplication is reversed. 
     As in the ARF to SRF transformation, the θ y  and θ x  angles are usually small while θ z  takes on a constant value between 0 and 2π with small variations around the constant value. Using a first order approximation for the sine and cosine functions, the small angle transform becomes ##EQU13## Multiplying the matrices, ##EQU14## 
     Composite ARF to IRF Small Angle Transformation 
     For faster and more efficient computation, a composite transformation matrix for small angles from the ARF directly to the IRF can be derived as ##EQU15## Performing the matrix multiplication and using angle sum and difference formulae, the transformation matrix ##EQU16## where γ z  =θ z  +ψ z . 
     The application of the details to the embodiment of the invention are now described. 
     System Model 
     To accomplish the boresight calibration, the boresight pitch angle θ y , the boresight roll angle θ x , and the boresight yaw angle θ z  are determined. That is, the orientation of the IRF to the SRF is ascertained. These angles change on a diurnal basis and an annual basis, primarily due to thermal effects from the position of the sun on the satellite. Variance measurements using known diurnal boresight curves show that the variance increases according to the square of time up to several hours. Based on several measurements every ten seconds, the model chosen for the system is a double integrator for each state driven by Gaussian white noise NID(0,σ). Only the output from the second integrator is filtered. The state definition is, therefore, ##EQU17## The system model is then x(t)= w(t)! with x(t) t=0  =0. 
     The state transition matrix Φ(t) is found by solving the state equation and Φ(t)=I where I is the identity matrix, i.e., all matrix elements are zero except for the main diagonal which is all ones. 
     The second order statistic (covariance equation) matrix is P xx  (t)= Q! so the mean squared value grows according to the time. squared, i.e., ##EQU18## 
     Measurement Model 
     The measurement model is a function of the state variables plus a random component of measurement noise or disturbance which is assumed to be Gaussian white noise with zero mean and covariance  R!. That is,  z k  != h k  (x(t k ))!+ v k  !. The measurement is  z k  ! at time t k , and  h k  ! is a nonlinear function.  V k  ! is the measurement noise. 
     The measurements are provided by discrete sensor elements in the focal plane illuminated by stars and ground sources on earth. Each sensor element is identified by azimuth and elevation coordinates in the IRF. The nonlinear measurement function is obtained in the process of transforming star data from the ARF to the IRF. The measurement equations are ##EQU19## 
     
         ε.sub.irf =cos.sup.-1 z.sub.irf.                   (3c) 
    
     Since T arf →irf contains the state variables, the given measurement equations is the nonlinear measurement function  h k  (x(t k ))!. 
     Observation Matrix 
     The Kalman filter requires a linear form of  h k  (x(t k ))! to compute the Kalman gain equation and the covariance matrix update. Linearization is accomplished by taking the partial derivative of each of the measurement variables with respect to each of the state variables: ##EQU20## 
     The steps for taking the partial derivatives of α (omitting the subscripts) with respect to θ (also omitting the subscripts) are as follows: ##EQU21## by eliminating common terms. Substituting the terms in the IRF for x, y, and z as given in equations (1a-1c), ##EQU22## Eliminating common terms, using the Pythagorean identity, and restoring the subscripts, ##EQU23## Therefore, (5) ##EQU24## 
     The above equations provide expressions for the IRF azimuth and elevation changes for small changes in each of the three state variables, viz., boresight pitch θ y , boresight roll θ x , and boresight yaw θ z . They are, however, functions of the partial derivatives of the x, y, and z components in the IRF and are therefore functions of the same state variables. Expressions for these partial derivatives can.be obtained by differentiating.the following set of equations expressed in matrix form: ##EQU25## T arf →irf is derived above and the ARF x, y, and z components are not functions of the state variables. The partial derivative of the above equations with respect to the generalized variable θ is ##EQU26## 
     The following three matrices are the partial derivatives with respect to each of the state variables and the following are used with the above equation to produce nine equations: ##EQU27## 
     Using the above equations, the H matrix can be derived as follows. For each measurement, the nine partial derivatives for x irf , y irf , and z irf  with respect to the boresight state variables using the above three equations are substituted into equation (6). Those results are substituted into the six equations represented by equations (4) and (5). 
     The implementation details are now set forth to summarize the complete implementation. 
     Boresight Filter 
     The boresight filter uses a standard modified Kalman filter with a nonlinear measurement model. The filter algorithm is set forth below to facilitate a practical implementation and to illustrate the entire process of producing a boresight state vector. 
     Fixed IR sources (targets) are used to determine the boresight state vector and can be stars primarily but also stable ground earth sources with a known location. The ephemeris for star sources are stored in ECF fixed coordinates for every thirty-second interval, for example. Every thirty seconds the ECF star coordinates and ground coordinates are transformed to ARF coordinates using the above transformation equations. These require knowledge of the satellite location, i.e., the ECF radius vector, also stored in an ephemeris. Only those stars and ground sources that would be in the sensor field of view every thirty seconds need be stored. Sensor observations are made in the vicinity (±500 μrad) of the predicted ARF location of the source. The collected observations are processed by a Representative Return Formation (RRF) procedure once per scan time of the sensor, i.e., one complete rotation of the satellite. 
     The Representative Return Formation 
     The RRF eliminates noisy returns and measures filter performance. The representative return is not filtered itself. Instead, each raw return comprising the representative return is processed through the filter with measurement noise based on each return&#39;s intensity as compared to the brightest (primary) return. 
     The following steps are followed to process data from each scan: 
     1. Sort by source collection identifier (ID) for each collection ID for each star or ground source in decreasing ARF elevation order; 
     2. Beginning with the brightest intensity return for an ID, accumulate all returns until a specified gap between adjacent returns is exceeded or there is an increase in intensity; 
     3(a). If at least one adjacent return for the primary return is found, then the process is finished. 
     3(b). Otherwise, tag the primary return as unavailable and repeat step 2 using the next brightest intensity return available. 
     4. If a representative return is found, compute the representative weighted average of azimuth and elevation as follows: ##EQU28##  where n=total number of raw returns used in the RRF, and M i  =intensity of the i-th raw return; and 
     5. Apply the error (difference) data to the boresight Kalman filter in the time order received. 
     The location of each representative return α irf ,rep (the average of the returns) and ε irf ,rep (a weighted average) is compared to a predicted source (target) location to determine the error statistics in step 5 above. 
     The Kalman filter processes the error statistics (in IRF coordinates) according to the following steps: 
     1. Compute Δt k  =t k  -t k-1 , the elapsed time since the previous filter update in seconds (the initial interval is Δt=0); 
     2. Propagate the state vector from the last update to the present time as x k  =Φx k-1  which, since Φ(t)=I, leaves the state vector unchanged after propagation, i.e., the state vector is propagated from its value at time t n  to time t n+1  by multiplying it by the transition matrix; (This step is trivial but included for completeness since it is necessary when estimation is improved by including velocity states as described below where Φ(t)≠I.) 
     3. Transform the ARF location to a unit vector using equations (2a-c); 
     4. Compute the measurement residual as the difference between the measured and predicted location in the IRF as 
     
          r.sub.k != z.sub.k !- h.sub.k (x.sub.k) ! 
    
     where  h k  (x k )! is found from equations (3a-c). These equations must use the propagated boresight state vector and the current best estimate of the attitude state vector in T arf →irf and the full form should be used instead of the small angle approximation; 
     5. Compute the observation matrix H k  using equations (4), (5), (6), and (7). These equations must use the propagated boresight state vector and the current best estimate of the attitude state vector; 
     6. Propagate the covariance matrix using (all uppercase letters are matrices) 
     P k  =ΦP k-1  Φ T  +Q k-1  and since Q(t) t=0  grows with the square of time, ##EQU29## 7. Compute the Kalman gain K k  using K k  =P k  H k   T   H k  P k  H k   T  +R k  ! -1  where all upper case letters represent matrices; 
     8. Update the estimate of the boresight state using the Kalman gain and residual as x k  =x k  +K k  r k , i.e., updating is performed by applying the probability measures to correct the state using the residuals; 
     9. Update the covariance matrix using the Kalman gain and the observation matrix using P k  = I-K k  H k  !P k . 
     The state vector output from step 7 is the filtered boresight state vector and is used to correct the mission-related returns from the sensors which are transformed to the ARF for mission data correction and then to the ECF for earth location and direction. The filter optimizes the accuracy of the result. 
     The block diagram of FIG. 6 show a system embodying the invention. A data processor 607 includes a computer section 637, a data storage 627, an input section 647, and an output section 657. The processing is executed by a computer 637 controlled by a stored program 617. 
     A sensor array 601 is located in a survey satellite with a transmitter 603. The scanner data is downloaded from the transmitter 603 to a receiver 605 which couples the scanning data signals to the input section 647 and the data is stored in the data storage section 627 from where it is retrieved for processing. 
     The raw data from the survey satellite is processed for corrections as described in the preceding explanations and, via the output section 657, is supplied to a utilization device, shown in FIG. 6 as a mapper. Other utilization devices are available but are not material to the present invention and will not be discussed further. 
     A flow chart for the program 617 is shown in FIG. 7. 
     In the following description, references are made to the flowcharts depicting the sequence of operations performed by the program. The symbols used are standard flowchart symbols accepted by the American National Standards Institute and the International Standards Organization. In the explanation, an operation may be described as being performed by a particular block in the flowchart. This is to be interpreted as meaning that the operations referred to are performed by programming and executing a sequence of instructions that produces the result said to be performed by the described block. The actual instructions used depend on the particular hardware used to implement the invention. Different processors have different instruction sets but persons of ordinary skill in the art are familiar with the instruction sets with which they work and can implement the operations set forth in the blocks of the flowchart. 
     In FIG. 7, the boresight filter program that controls the data processor 607 of FIG. 6 is detailed in a flowchart. It is depicted as a subroutine that is typically called at 30-second intervals and is invoked by a hardware or software timer 701 via the enter terminal 703. 
     The raw data from the scanner or sensor array 601 of FIG. 6 is stored via the transmitter 603 and receiver 605 and the input section 647 of the data processor of FIG. 6 as shown in an input block 705 of FIG. 7. 
     Next, the coordinates of the known image, i.e., the target, candidates (those targets within the sensible area of the sensor array) are retrieved in ECF according to a process block 707. The ECF coordinates are transformed in a process block 709 to IRF coordinates as explained in detail above. 
     The raw data returns are processed in a process block 711 in the order received from the sensor array one at a time. A database gap or an increase in intensity is detected in a decision block 713. If either condition is sensed, then the primary return is tagged as being unavailable as shown in a process block 714. If, according to a decision block 717, the raw data return is the last, then the program execution returns to read in more raw data returns from the scanner at the input block 705. Otherwise, the execution returns to process the next return at process block 711. 
     If the decision block 713 finds neither a database gap nor an intensity increase, the centroid of the raw return is computed in a process block 719 by averaging as described above. A process block 721 computes the difference between the IRF coordinates of the target and centroid of measurement. These differences are processed by a subroutine 723 in the time order received to supply an optimal estimation of the differences. The optimal estimator is embodied as a Kalman Filter. The subroutine 723 is shown in detail in FIG. 8. 
     At the end, the procedure is exited via a terminal block 725. 
     The flowchart of an optimal estimator embodied as a Kalman Filter illustrates the sequence of program instructions employed to implement the Kalman Filter. It is shown as a subroutine that is entered at a terminal block 801. The methods of calculation of the various expressions are explained in detail above. In most cases, programs to make the individual process computations are commercially available. Many are readily implemented by computer programmers of ordinary skill in the art. 
     First, the time interval At from the last update of the boresight state vector is calculated by a process block 803. The state vector is propagated to the present time from its prior update in a process block 805. 
     Next, the ARF coordinates of the unit earth vector are computed in a process block 807. In a process block 809, the measurement residual is computed. The H k  matrix is computed in a process block 811 and used to propagate the previous covariance matrix according to a process block 813. 
     These computations permit the Kalman Gain to be computed by a process block 815 and to update the boresight state vector by a process block 817. The covariance matrix is updated in a process block 819. Then the subroutine is exited via a terminal block 821. 
     The program flowchart of FIG. 9 illustrates the steps for applying the boresight corrections to the mission data. It is depicted as a subroutine entered at a terminal block 901. The raw mission data is read by an input block 903 similar to the way the target data is read from the scanner (sensor array) as described for the input block 705 in FIG. 7. 
     The boresight state vector, updated by the process block 817 of the flowchart of FIG. 8, is applied to the raw mission data according to a process block 905. The corrected mission data is supplied via an output block 907 and the subroutine is exited via a terminal block 909. 
     The correction technique is explained in connection with FIG. 10. The equations with their derivation have been set forth above but will be repeated for each step in the following explanation for convenience of reference. The routine is entered via a terminal block 1001 and begins with reading in the index (id) numbers of each sensor in an input/output block 1003. In a process block 1005, the azimuth and elevation of the particular sensor is obtained using, for example, a look-up table technique. The focal plane table is accessed with the individual sensor (cell) index number and the result is a and e of the sensor. 
     A process block 1007 is supplied with a and e which are converted to IRF coordinates with the equations noted above, viz., 
     x irf  =cos α sin ε, 
     y irf  =sin α sin ε, and 
     z irf  =cos ε. 
     The IRF coordinates are then converted to corrected SRF coordinates in a process block 1009. The inverse transform, T irf →srf (the boresight state vector), is multiplied times the IRF coordinates. 
     This is ##EQU30## 
     Next, in a process block 1011, the attitude state vector, T srf →arf, is applied to compute the ARF corrected coordinates by ##EQU31## 
     For those uses of the corrected mission data that need polar coordinates, i.e., corrected α and ε, a process block 1015 converts the Cartesian coordinates to polar. As explained above, ##EQU32## 
     
         ε.sub.arf =cos.sup.-1 z.sub.arf. 
    
     The ECF coordinates of the earth frame unit vector are computed in a process block 1017 using the satellite unit vector and T arf →ecf (=T ecf →arf T ) viz., ##EQU33## 
     The ECF coordinates can be converted to geodetic latitude and longitude if necessary or desired. 
     The accuracy of the conversion from the input cell index to the earth centered frame of reference, i.e., the sensor data corrected for boresight errors, is directly related to the accuracy of the focal plane vector (look-up) table, the satellite attitude state vector, and the satellite ephemeris. The invention improves the accuracy of the computed earth frame location corresponding to a given sensor cell index by improving the accuracy of the boresight state vector estimate. 
     To improve accuracy at the expense of longer computing time, the boresight state vector can be augmented to a six-state version vice the three-state version described above. In addition to the pitch, roll, and yaw state variables, their derivatives (representing velocity in the three directions) can be included. The augmented state vector is ##EQU34## 
     The covariance matrix of the Kalman Filter is modified to comport with the six-state boresight state vector. That is, ##EQU35## 
     The transition matrix for the augmented boresight state vector is no longer a unit (or identity) matrix. Instead, it is ##EQU36## where u is replaced by y, x, or z for each of the coordinate axes. 
     Finally, the H k  matrix must be augmented for compatibility. Zeroes are simply added so that ##EQU37## 
     The improvement using the augmented matrices in the Kalman Filter is illustrated in FIG. 11 for a sample case over a twenty-four hour period. FIG. ll(A) shows the boresight pitch, roll, and yaw for the three state boresight state vector. FIG. o 11 (B) illustrates the same values for the six-state result. The graph line 1101 is the boresight yaw for the three-state implementation and 1102 is the boresight yaw for the six-state implementation. Similarly, the lines 1103 and 1105 are the boresight roll and pitch values, respectively, for the three-state implementation and the lines 1104 and 1106 are the same values, respectively, for the six-state implementation. 
     The boresight state vector values are smoother using the augmented matrices and more accurate but at the expense of more computation time. 
     FIG. 11 is a graph of sample data over a twenty-four hour period showing the boresight pitch, roll, and yaw for a three-state implementation of the invention. The boresight yaw is shown as line 1101, the boresight roll as line 1103, and the boresight pitch as line 1105. 
     FIG. 12 is a graph of the same sample data over a twenty-four hour period showing the boresight pitch, roll, and yaw for a sixstate augmented implementation of the invention. The boresight yaw is shown as line 1201, the boresight roll as line 1203, and the boresight pitch as line 1205. The smoothing due to inclusion of the velocity vector states is readily discernible. 
     While the invention has been particularly shown and described with reference to a preferred embodiment thereof, it will be understood by those skilled in the art that various changes and modifications in form and details may be made therein without departing from the spirit and scope of the invention according to the following claims.