Spatial geometrical calibration of optoelectronic devices

. Space equipment experiences high thermal and mechanical loads that cause deformation of its structures, as well as elements of optical systems and radiation receivers, which can lead to changes in the mutual position of the optoelectronic device elements, disturbing its alignment, degrading image quality and increasing pointing errors. This paper presents a method of spatial geometric calibration that does not require additional equipment and investigates the dependence of the identification results on the viewing angle. The proposed methodology of optoelectronic instrument calibration is intended to improve efficiency, reliability and survivability of automated spacecrafts during their operation.


Introduction
Due to thermal and mechanical loads on optoelectronic devices (OEPs), there are deformations of the structure, changes in the properties of optical elements and radiation receivers, which can lead to changes in the mutual position of their optical system elements.Examples of these loads: vibrodynamic and force influences on a spacecraft in the process of its launch into orbit and during maneuvering; lack of gravity in orbit, causing redistribution of stresses and structural deformations; ionizing radiation from natural sources, affecting materials of optical parts; thermal loads, causing thermal deformations [1][2][3].
All space-based OEPs are repeatedly tested and calibrated on Earth under laboratory conditions.Testing and calibration of spacecraft orientation and navigation systems, mutual spatial orientation of its sensors and spacecraft are also carried out in laboratory conditions [4].Evaluation of the current parameters of these systems is carried out at the stage of flight design tests.For some remote sensing satellites, the estimation of current parameters of optoelectronic devices is realized at the stage of operation with the use of a limited set of test and reference objects.
However, in real conditions there is a need for systematic evaluation: for optoelectronic devices -resolution, sensitivity in different spectral ranges and at different points of the field of view, distortions of formed images, etc.; for attitude control systems -the accuracy of pointing out during maneuvers, correction of orbital parameters, correction of orientation angles, etc. [5].
Methods of spatial calibration imply special observations of a set of reference stars.The software-algorithmic complex, when used in the spacecraft control center, can provide planning of observations, collection, processing and analysis of observational data, development of recommendations for the operation of target optoelectronic equipment and elements of the spacecraft positioning control system [6][7][8].
Over the time of spacecraft operation, the mutual arrangement of its modules can change, which inevitably leads to an increase in the pointing error.If one takes the obtained frame and identifies the objects that are caught in it, one can determine the coordinates of the point of aiming with high accuracy (depending on the resolution of the camera), and then make corrections during the next observation session.It is worth noting that in order to obtain reliable results it is necessary to determine the coordinates of the photocenters of stars in the image as accurately as possible.Any distortions of the useful signal will inevitably lead to an increase in the error of determining the coordinates of the center [9][10][11].Thus, Figure 1 demonstrates the influence of matrix sensitivity inhomogeneity on the obtained result.

Identification of stars in the frame
In the proposed algorithm, the main method of star identification is to compare the angles between them, but this requires knowing the focal length of the optical system.If it is known from telemetry or instrument data, the task is simplified.In addition, it is recommended to pre-process the images to improve the quality of identification [12][13][14][15][16][17][18][19][20][21].
When the focal distance is not known, it is necessary to know which stars are present in the frame in order to identify them.To identify the stars, present in the frame, the method of similar triangles is used.First, it is necessary to find the photocenters of stars in the image (x, y), and then select the equatorial coordinates of stars from the catalog (Ra -direct ascent (α), Dec -declination (δ)).Then you need to make two lists: the first will contain all possible triangles of stars from the image (list A), the second will contain triangles from the catalog (list B).It is also necessary to translate the equatorial coordinates of stars from the catalog into their projection on the focal plane (ideal coordinates):  = ()sin (− 0 ) sin( 0 )+() cos( 0 )cos (− 0 ) ; (1)  = cos( 0 )−() sin( 0 )cos (− 0 ) sin( 0 )+() cos( 0 )cos (− 0 ) , (2) where δ0, α0 are equatorial coordinates of the center of the considered area.
After all possible triangles have been composed, it is necessary to select those of them, where one of the sides is greater than a certain threshold (Lpor, calculated as the number of pixels in the diagonal, multiplied by the coefficient k).The coefficient k is chosen from the consideration -the larger the field of view of the image, the smaller k is.So at angles of view less than 1, k = 1, and at angles ~20 0 k = 0.25.It is worth noting that if the image of stars takes from 1 to 4 pixels, then you should limit the length of the segment on the smallest side, and leave it so that it is at least several times larger than the diameters of the stars in pixels.
For the triangles in list A and B, calculate the lengths of the sides (a, b, c), with c ≥ b ≥ a.Let us use the sign of similarity of triangles by the ratio of sides: If the difference between p and q for two triangles from lists A and B satisfies the chosen error, then similar vertices are compared with each other, which are opposite to the large, middle and small sides, respectively.Each of the three similar vertices of the triangle receives a vote.After all the triangles have passed the comparison procedure, the points with the most votes are selected and they are considered to be matched in the image and in the stellar catalog.

Determining the focal distance
We introduce a coordinate system: the center (0, 0) is in the center of the corner pixel of the frame, the coordinates of the opposite (diagonal) pixel (NX -1, NY -1), where NX is the number of pixels horizontally, NY -vertically.
Calculate the coordinates of the center of the frame by the formulas: We take the identified stars and construct unit vectors from the equatorial coordinates of the stars Ra and De:   = cos  cos  ;  = cos  sin  ;  = sin . (5) We find the angular distance between pairs of vectors  ⃗  ,  ⃗  (i, j are vector numbers and i≠j) as the arc cosine of the scalar multiplication of vectors: Then the value of possible focal length of lens F (in pixels) is set and with some step.From the coordinates of the images of stars on the focal plane X and Y the vectors f optical center of the lens -focal plane are constructed:  , =   −   ; , =   −   ; , = Φ.(7) Convert vectors f to unit vectors F using normalization: where i is the number of the vector (star in the frame).The norm of the vector is equal to: Find the angular distance between pairs of vectors   ⃗⃗ ⋅   ⃗⃗ (i≠j): We calculate the inconsistency: The summation is performed in such a way that each pair of vectors enters the sum only once.The number of differences between the angles ψ and φ included in the misalignment is denoted by Np.We choose the value of focal distance F, at which the misalignment S 2 is minimal.The value ∆ = √ 2   ⁄ serves as an estimate of the error in determining the coordinates of the stars on the focal plane (in pixels).
Simplifying, we assume zero distortion, i.e. the star ray passes through the optical center of the lens and falls on the focal plane in a straight line.The star image is formed at the point of intersection of the beam and the focal plane.Angular dimensions of the frame WX, WY are calculated according to the formulas:

Recognizing star configurations, identifying these stars with the catalog, and finding potential orientation values
In this part of the work, the Gaia catalog was used, but the use of data from other catalogs is admissible [22].The following stars were selected from it and two lists were compiled: the first, up to 7.5 m in the G band (37045 stars) -the list name is G7.5; the second, up to 13.0 m in the G band (7 283 019 stars) -the list name is G13.
The catalogs contain the coordinates Ra, Deс and the mG magnitude in the G band and additional data.In the text form, the catalogs are sorted in ascending order of mG.
From the 1st catalog, the catalog of distances (D7.5) between pairs of stars, in which both stars are brighter than 7.5 m in the G-band and the angular distance between the pairs d < 10 0 is constructed.The catalog contains the following data: d is the distance between the pair in degrees; n1 is the number of the first star of the pair in the lists G7.5, G13 (numbering from 1); n2 is the number of the second star of the pair in the lists G7.5, G13; m1 is the sidereal magnitude of the first star of the pair in the list G; m2 is the sidereal magnitude of the second The input data of the stars in the frame are read.Only the data obtained directly from the frame are used: X, Y and Bri, and Bri is used only for the initial sorting of the list, and is not used in the algorithm itself.Sort the read list of stars T by decreasing Bri (the brightest stars are at the beginning of the list).
From the list T we take 3 stars with numbers i0, i1, i2 (i0<i1<i2).The angles between any pair of these stars (d01, d02, d12) must be smaller than the maximum distance value in the distance catalog.We calculate unit vectors Fi for these stars by formulas (7-9), and then find the angle between vectors F0 and F1: The resulting list will be called D0.The record structure of the list D1 includes only the numbers of stars in the lists G7.5/G13 (n0, n1).These pairs are the 1st side of possible star triangles in the sky corresponding to the images of stars with numbers i0, i1, i2 in the frame.
Find the angle between vectors F0 and F2, we denote it by d02.We select from the catalog of distances the stars for which the condition is fulfilled: where Δx is an estimate of the error in determining the coordinates of stars on the focal plane.
The resulting list will be called D1.The record structure of the list D1 (n0, n2).These pairs are the 2nd side of possible star triangles in the sky corresponding to the images of stars with numbers i0, i1, i2 in the frame.Then we have to find the entries of the lists D1 and D2 in which the numbers of the first stars in the pairs coincide.We put these entries in the list D01.The list includes the numbers of stars from both compared lists (n0, n1, n2).
These triples of stars are possible triangles of stars in the sky corresponding to the triangle of images of stars with numbers i0, i1, i2 in the frame.Check the length of its 3rd side, then find the angle between the vectors F1 and F2, and denote it as d12.We select from the catalog of distances the stars for which the condition is satisfied, taking into account the error Δx estimate: The resulting list is called D12.Compare the lists D01 and D12.We look for elements in the list D01, in which numbers 2 (n1) and 3 (n2) of stars of the triangle coincide with the corresponding numbers (n1, n2) in the list D12, and form the list D3.The structure of the list D3 coincides with that of the list D01 and includes data from both compared lists (n0, n1, n2).
These triples of stars are the possible triangles of stars in the sky corresponding to the triangle of star images with the numbers i0, i1, i2 in the frame, in which the lengths of all three sides within the errors coincide with the distances between the images of stars in the frame.
Then it is necessary to check the identifications.If there are no entries in the D3 list, the identification has failed.If there is one entry, it is most likely a correct identification, but this identification must still be checked.
Choose the entry from the list D3.For the stars with numbers i0, i1, read from the list G7.5 their equatorial coordinates Ra0, Dec0 and Ra1, Dec1.Using formulas (14-15) we calculate the unit vectors E0 and E1.Knowing the pairs of vectors E0 E1 and F0, F1, we construct a rotation matrix R, which transforms the coordinate system of the focal plane to the equatorial coordinate system and back.From the pair of vectors E0 and E1 we construct an ⃗⃗⃗⃗⃗⃗⃗ ).As a vector  0 ⃗⃗⃗⃗⃗⃗⃗ we take vector E0,  2 ⃗⃗⃗⃗⃗⃗⃗ -unit vector co-directional to vector product E0×E1, -vector product BE2×BE0.
Vectors BE are unit, mutually orthogonal and form the right triplet ( 0 ⃗⃗⃗⃗⃗⃗⃗ ,  1 ⃗⃗⃗⃗⃗⃗⃗⃗ ,  2 ⃗⃗⃗⃗⃗⃗⃗ ).From these vectors we construct the basis matrix BE (3×3, the vectors go by columns of the matrix).Similarly, we construct an orthonormal basis BF from the pair of vectors F0 and F1.These vectors are also unit, mutually orthogonal, and form the right triplet ( 0 ⃗⃗⃗⃗⃗⃗⃗ ,  1 ⃗⃗⃗⃗⃗⃗⃗⃗ ,  2 ⃗⃗⃗⃗⃗⃗⃗ ), and the basis matrix has unit determinant equal to 1 and is orthogonal.Then the rotation matrix R translating the basis BE into the basis BF is defined as  =  ⋅  T .After this, the "celestial" vector e is translated into the unit vector of the focal plane f by multiplying the matrix and vector (16).The inverse transition is represented in formula (17): =   •  .(17) Verification is done by multiplying the vector in both directions and comparing it with the original vector.We rotate all unit "celestial" vectors of stars Ei to the basis of the focal plane -we obtain vectors fi.From them we find the expected coordinates of stars XE, YE on the focal plane by the formulas: Find differences XE, YE with X and Y and memorize.Determine the position of the center of the frame in the sky.The center of the frame, namely the point (XC, YC) corresponds to a unit vector   ⃗⃗⃗ = (0, 0, 1).We convert it to the equatorial basis, and then calculate its angular coordinates: ⃗  =   •   ;  = atan 2 ( , ,  , );  = asin( , ).(19) We calculate the angular diameter of the circle W inside which the whole frame fits:  = √  2 +   2 .We select all stars from the catalog G13 which lie on the celestial sphere inside the circle with the center pointed by the unit vector EC and with the angular diameter W. We denote the list by G.
The implementation of sampling is possible in several ways (to speed up the procedure).It depends on the implementation of data structures and databases.Additionally, the stars included in the list G can be limited in magnitude by  <   .In this work, we used for wide-field lenses and for narrow-field lenses.
For stars from the list G, we do the following: we calculate the unit vectors E by formulas (5) using Ra and De values, convert E vectors to unit vectors F in the focal plane coordinate system using R matrix, and calculate the expected X, Y coordinates of stars from the list G on the focal plane by formulas (18)(19).
We exclude from the list G the stars which go beyond the frame boundaries, i.e. for which at least one of the following conditions is fulfilled: X < 0, X > NX, Y < 0, Y > NY.Then we calculate the distances on the focal plane between the stars in the frame (list T) and from the Gaia catalog (list G): where t is the number of a star from the list T, g is the number of a star from the list G. Stars from lists T and G are considered to be identified if   < ∆.The threshold ∆ can be assumed to be equal to the error of   ∆x.We can circumscribe a circle with a radius ∆ around each star in one list and count how many stars from the other list fall inside it.Three situations are possible: 1) no stars of the other list fall inside the circle -the star is not identified; 2) 1 star of the other list -unambiguous identification; 3) several stars of the other list -ambiguous (multiple) identification.Then we can calculate the following values: Taking as the main list T: T0 is the number of unidentified stars in list T; T1 is the number of uniquely identified stars in list T; T2 is the number of repeatedly identified stars in list T; T0+T1+T2 = N(T), here N(T) is the number of stars in list T; Taking list G as the main one: G0 is the number of unidentified stars in list G; G1 is the number of uniquely identified stars in list G; G2 is the number of multiply identified stars in list G; N(G) = G0 + G1 + G2 is the total number of stars in list G; GT value: the number of pairs of stars from the lists T and G for which Rtg < ∆R.This number is symmetric with respect to both lists and can exceed the number of stars in the shorter one.
The criterion for correct identification can be considered if N(T) < N(G) and T1 + T2 > 50% of N(T).Similarly, if N(G) < N(T) and G1 + G2 > 50% N(G), the identification is correct.We can also look for the largest value of T1 + T2 when N(T) < N(G) or the largest value of G1 + G2 when N(G) < N(T).
At the stage of testing it was found out that the best performance of the algorithm is at NϬ < 5. On frames with field view less than two degrees, it is possible to identify from 80% to 95% of objects, and with increasing angle of view the influence of geometric distortions becomes critically high.For example, at angles of view of ~30 degrees it is possible to identify stars that are only in the central region, which is 10%-30% of the total number of stars.

Conclusions
Since the destabilizing effects of the space environment cause the optoelectronic instruments to change their characteristics, there is a need for their periodic calibration.In the present work, the methods of calibration are considered, which do not require installation of additional equipment, to work with the frames obtained onboard the spacecraft.Algorithms for spatial sensitivity calibration of spacecraft matrix and algorithms for identification of objects in the image are presented.As a result of testing on test frames with field of view up to two degrees, the identification algorithm recognized up to 95% of objects, which is quite a satisfactory result.Also, studies have found that increasing the viewing angle has a strong impact on the quality of the identification during calibration.