0% found this document useful (0 votes)
3 views11 pages

Multiexposure Imaging for Star Trackers

This research article presents a multiexposure imaging approach for intensified star trackers, which enhances the attitude update rate by recording multiple groups of star positions within a single exposure. The proposed method effectively improves the signal-to-noise ratio for dim stars and optimizes key parameters for better performance. Simulations and experiments validate the feasibility and effectiveness of this approach in overcoming limitations in pixel data transmission and processing time.

Uploaded by

mahi.zoubida27
Copyright
© All Rights Reserved
We take content rights seriously. If you suspect this is your content, claim it here.
Available Formats
Download as PDF, TXT or read online on Scribd
0% found this document useful (0 votes)
3 views11 pages

Multiexposure Imaging for Star Trackers

This research article presents a multiexposure imaging approach for intensified star trackers, which enhances the attitude update rate by recording multiple groups of star positions within a single exposure. The proposed method effectively improves the signal-to-noise ratio for dim stars and optimizes key parameters for better performance. Simulations and experiments validate the feasibility and effectiveness of this approach in overcoming limitations in pixel data transmission and processing time.

Uploaded by

mahi.zoubida27
Copyright
© All Rights Reserved
We take content rights seriously. If you suspect this is your content, claim it here.
Available Formats
Download as PDF, TXT or read online on Scribd

Research Article Vol. 55, No.

36 / December 20 2016 / Applied Optics 10187

Multiexposure imaging and parameter


optimization for intensified star trackers
WENBO YU, JIE JIANG,* AND GUANGJUN ZHANG
Key Laboratory of Precision Opto-Mechatronics Technology, Ministry of Education, School of Instrumentation Science
and Opto-electronics Engineering, Beihang University, 37 Xueyuan Rd., Haidian District, Beijing 100191, China
*Corresponding author: jiangjie@[Link]

Received 18 July 2016; revised 29 October 2016; accepted 11 November 2016; posted 14 November 2016 (Doc. ID 270696);
published 13 December 2016

Due to the introduction of the intensified image detector, the dynamic performance of the intensified star tracker
is effectively improved. However, its attitude update rate is still seriously restricted by the transmission and
processing of pixel data. In order to break through the above limitation, a multiexposure imaging approach
for intensified star trackers is proposed in this paper. One star image formed by this approach actually records
N different groups of star positions, and then N corresponding groups of attitude information can be acquired.
Compared with the existing exposure imaging approach, the proposed approach improves the attitude update
rate by N times. Furthermore, for a dim star, the proposed approach can also accumulate the energy of its N
positions and then effectively improve its signal-to-noise ratio. Subsequently, in order to obtain the optimal
performance of the proposed approach, parameter optimization is carried out. First, the motion model of
the star spot in the image plane is established, and then based on it, all the key parameters are optimized.
Simulations and experiments demonstrate the feasibility and effectiveness of the proposed approach and
parameter optimization. © 2016 Optical Society of America
OCIS codes: (120.4640) Optical instruments; (120.6085) Space instrumentation; (110.4190) Multiple imaging; (100.4145) Motion,
hyperspectral image processing.

[Link]

1. INTRODUCTION intensified image detector into the star-tracker field [15,16].


Taking stars as the observation objects, a star tracker can pro- It can significantly amplify the weak starlight and thus obtain
vide absolute attitude information with respect to the inertial high sensitivity in a very short exposure time. In this way, the
coordinate system. Among all the attitude determination devi- dynamic performance of the star tracker is effectively improved.
ces, the star tracker is by far the most accurate one [1–3]. With From the perspective of signal processing, the star-tracker
the development of space exploration technologies, such as attitude measurement is equivalent to the discrete sampling
agile satellites and space weapons, the requirements for the dy- of the continuous varying attitude processes, and the attitude
namic performance of a star tracker are rapidly increasing [4,5]. update rate is the sampling frequency. Therefore, according to
Under dynamic conditions, a star spot continuously moves and the sampling theorem, the attitude update rate of the star
forms a smeared star streak during the exposure time, resulting tracker should be improved with the increase of the dynamic
in the energy dispersion, the signal-to-noise ratio (SNR) performance of the star tracker. However, in practice, the atti-
decreasing and eventually degrading the centroiding accuracy. tude update rate is not only related to the exposure time, but
This will have adverse effects on the success rate of star iden- also to the time of pixel data transmission and processing as well
tification, tracking stabilization, and attitude precision [6–13]. as the time of star tracking and attitude estimation. Research
Yan et al. established a detailed dynamic star spot-imaging [17,18] shows that, when the parallel and pipeline structure is
model and gave the optimal exposure time on dynamic condi- adopted, compared with the sequence structure, the attitude
tions [14]. They pointed out that long exposure time will bring update rate can be improved. In this case, the attitude update
long tailing and affect the centroiding accuracy. However, short rate is determined by the most time-consuming stage in the
exposure time means that the sensitivity is low and the number pipeline structure. For the traditional star tracker, a long expo-
of detected stars is reduced, which in turn affects the improve- sure time is essential to obtain sufficient sensitivity. Therefore,
ment of dynamic performance. In order to solve the problem of exposure time is the main bottleneck of the attitude update
low sensitivity, Katake and Bruccoleri first introduced the rate of a traditional star tracker. In addition, for the existing

1559-128X/16/3610187-11 Journal © 2016 Optical Society of America


10188 Vol. 55, No. 36 / December 20 2016 / Applied Optics Research Article

intensified star tracker, because of the introduction of the image


intensifier, the exposure time is much shorter. Therefore, the
stage of transmission and processing of pixel data becomes
the most time-consuming one in the whole pipeline structure.
In particular, when the larger-array image detectors are used in
the star tracker, this stage consumes more time, and it will seri-
ously restrict the improvement of the update rate of the star
tracker.
In order to break through the limitation of the transmission
and processing time of pixel data on the attitude update rate, a
multiexposure imaging approach for intensified star trackers is
proposed in this paper. This approach utilizes the two-level ex-
posure control structure, consisting of the image detector and
the image intensifier. First, the exposure time of the image de-
tector is broadened, and then, within the broadened exposure Fig. 1. Imaging principle of the intensified image detector.
time, the image intensifier is controlled for N short-time ex-
posures. Therefore, the one star image formed in this way ac-
tually records N different groups of star positions. When they regenerate the photon signal. This photon signal is transmitted
are, respectively, performed with star tracking and attitude es- to the image detector by a fiber optic taper bundle. Finally, the
timation, we can obtain N corresponding groups of attitude image detector completes the imaging of the transmitted
information. This is equivalent to the attitude update rate of photon signal and forms a digital image that can be identified
star tracker being improved by N times. Subsequently, in order and processed. During the whole imaging process, the starlight
to obtain the optimal performance of the proposed approach, signal totally passes through six links, namely, the input
parameter optimization is carried out. First, the motion model window, the photocathode, the microchannel plate, the
of the star spot in the image plane is established, and then based phosphor screen, the fiber optic taper bundle, and the image
on it, all the key parameters are optimized. Compared with the detector. Among them, the first four constitute the image in-
existing exposure imaging approach, the approach proposed in tensifier. Therefore, the image intensifier plays the role of the
this paper greatly improves the attitude update rate without any photon-photoelectron-photon conversion and amplification.
increase of data amount of star images. Furthermore, for a dim Acceleration and multiplication of the microchannel plate is
star, the proposed approach can also accumulate the energy of the key to starlight amplification, and the amplification effect
its N positions and then effectively improve its SNR. is positively correlated with the gain voltage and the exposure
This paper is organized as follows. In Section 2, the imaging time of the image intensifier.
principle of the intensified image detector is described, and Different from the traditional image detector, the intensified
then based on it, the principle of the multiexposure imaging image detector possesses the two-level exposure control struc-
approach is presented in detail. Parameter optimization is ac- ture. On one hand, the exposure shutter Ct of the image
complished in Section 3, so as to achieve the optimal perfor- detector is completely the same as that of the traditional image
mance of the proposed approach. Simulations carried out in detector. They are both used to directly control the imaging
Section 4 and experiments in Section 5 both verify the multi- process of the whole image. When the exposure shutter
exposure imaging approach as well as parameter optimization, Ct closes, the imaging process is ended, and then the trans-
and the results show that they are feasible and effective. Finally, mission and processing of the star image pixels is started.
conclusions are drawn in Section 6. According to the above functions of Ct, it can be named
as the imaging shutter. On the other hand, the exposure shutter
I t of the image intensifier is used to control the amplification
2. MULTIEXPOSURE IMAGING APPROACH of the incident starlight, but not to directly control the imaging
The imaging principle of the intensified image detector is process of the whole image. When the exposure shutter I t is
described in this section. Then, the principle of the multiexpo- open, the starlight can be amplified and transmitted to the im-
sure imaging approach as well as the existing exposure imaging age detector. Otherwise, all the starlight is prohibited from
approaches are all presented. being transmitted to the image detector. Similarly, according
to the above functions of I t, it can be named as the gain shut-
A. Imaging Principle of the Intensified Image ter or sampling shutter. If and only if both Ct and I t are
Detector open, can stars be imaged by the intensified image detector.
Figure 1 [15,19,20] shows the imaging principle of the inten-
sified image detector. Approximately infinitely far starlight en- B. Principle of Multiexposure Imaging Approach
ters the optical lens as parallel light. Then the starlight Figure 2 shows sequence diagrams of parallel and pipeline
converges into a spot by the lens and passes through the input processing of star trackers with different exposure modes. In
window. By employing the photocathode, the photon of the Fig. 2, the pipeline structure is composed of three stages,
starlight is transformed into the photoelectron. Subsequently, namely, exposure, transmission and processing of pixels, as well
the photoelectron is accelerated and multiplied in the micro- as star tracking and attitude estimation. The time required for
channel plate and then bombards the phosphor screen to three stages are, respectively, T e , T p , and T q in seconds. The
Research Article Vol. 55, No. 36 / December 20 2016 / Applied Optics 10189

this will seriously restrict the further improvement of the


update rate F of the star tracker.
In order to break through the limitation of the transmission
and processing time T p of pixel data on the update rate F , a
multiexposure imaging approach for intensified star trackers is
proposed in this paper; its sequence diagram of parallel and
pipeline processing is shown in Fig. 2(c). Based on the charac-
teristic that time T e is much less than time T p , the core idea of
the proposed approach is that there are multiple exposures in
time T p , rather than only one. If N is the number of times of
multiexposure imaging, a single star image actually contains N
different groups of star position information at different mo-
ments. When they are, respectively, performed with star
tracking and attitude estimation, N corresponding groups of
attitude information can be acquired from only one star image.
In other words, the attitude update rate of the star tracker is
improved by N times.
In order to insert multiple exposures into time T p , the key
here is that the intensified image detector possesses the two-
level exposure control structure. In Fig. 2(b), as the exposure
time of the sampling shutter I t and that of the imaging shut-
ter Ct are exactly the same (i.e., both time T e ), each period T
of the star tracker can only have one exposure imaging. When
the exposure time of the imaging shutter Ct is broadened
from T e to T L , as shown in Fig. 2(c), due to T e ≪ T L , N
sampling shutter I t can be inserted into one exposure time
T L of the imaging shutter Ct. Each exposure time T e of the
sampling shutter It is equivalent to an exposure imaging in
the time T L . In this way, the multiexposure imaging can be
Fig. 2. Sequence diagrams of parallel and pipeline processing of star achieved within one period T of the star tracker. In particular,
trackers with different exposure modes. (a) Sequence diagram of par- when N  1, the multiexposure imaging approach in Fig. 2(c)
allel and pipeline processing of the traditional star tracker; (b) sequence is exactly the same as that in Fig. 2(b), and their effective ex-
diagram of parallel and pipeline processing of the existing intensified posure times are both the exposure time T e of the sampling
star tracker; (c) sequence diagram of parallel and pipeline processing of shutter It.
the intensified star tracker with the multiexposure imaging mode. In practice, as shown in Fig. 2(c), the time T L cannot be
Stage1, exposure; Stage2, transmission and processing of pixels; broadened without limit, so there is an upper limit value.
Stage3, star tracking and attitude estimation; T e , exposure time; When the exposure time T L meets T L ≤ T , its broadening
T p , the time of pixel data transmission and processing; T q , the time
does not change the size of the original period T of the star
of star tracking and attitude estimation; T , working period of star
tracker. In this case, the period T is still determined by the
tracker; F , attitude update rate of star tracker; T L , broadened exposure
time; T i , interval time between two adjacent short exposures. time T p and satisfies T ≥ T p . Under ideal conditions, if the
interval time between two adjacent T p is ignored, T , T p ,
and T L will meet Eq. (1), as follows:
working period and the attitude update rate of the star tracker T  T p  T L: (1)
are expressed as T in seconds and F in hertz, respectively.
Figure 2(a) shows the sequence diagram of parallel and pipeline According to Eq. (1), the above three times are no longer
processing of the traditional star tracker. In Fig. 2(a), a long distinguished and all are expressed as T L in seconds.
exposure time is essential for the traditional star tracker to ob- Compared with Fig. 2(b), the multiexposure imaging ap-
tain sufficient sensitivity. Therefore, the long exposure time T e proach in Fig. 2(c) improves the utilization efficiency of the
is the main bottleneck of the attitude update rate F . Figure 2(b) whole period T L without changing its size. If the interval time
shows the sequence diagram of parallel and pipeline processing T i (in seconds) between two adjacent I t is known, the atti-
of the existing intensified star tracker. As shown in Fig. 2(b), tude update rate F of the proposed approach is expressed as
due to the introduction of the image intensifier, the exposure 1
time T e is much shorter. In this way, the attitude update rate F F : (2)
Ti
is improved to a certain extent. Nevertheless, with the decrease
of the exposure time T e , the stage of transmission and process- As shown in Eq. (2), the update rate F is completely deter-
ing of pixel data becomes the most time-consuming one in the mined by the interval time T i . Due to T i ≪ T p , the proposed
whole pipeline structure. For larger-array image detectors, T p approach breaks through the limitation of the transmission and
will be longer. Since the working period T must meet T ≥ T p , processing time T p of pixel data on the update rate F .
10190 Vol. 55, No. 36 / December 20 2016 / Applied Optics Research Article

Fig. 3. Effect of multiexposure imaging approach. (a) Simulated star image of the proposed approach; (b) imaging process of the proposed
approach.

When the star tracker moves with the angular velocity ω, ~ optimal exposure time under dynamic conditions that decrease
the effect of the multiexposure imaging is shown in Fig. 3. with the increase of the angular velocity of the star tracker
Figure 3(a) shows the simulated star image of the proposed ap- [13,14]. Moreover, when the angular velocity increases, the lin-
proach, and Fig. 3(b) shows the imaging process of a group of star ear velocity V star of the star spot in the image plane increases
spots belonging to the same star. As is seen in Fig. 3, the selection accordingly, thus resulting in overall reduction of the second
criterion for T i is that two adjacent star spots of the same star are item. In summary, with the increase of the angular velocity,
not overlapped with each other. Otherwise, the corresponding the lower bound of the interval time T i decreases, and the
attitude information cannot be calculated from each exposure corresponding upper bound of the update rate F gradually in-
imaging, thus losing the significance of multiexposure. creases. Therefore, in theory, the update rate F of the proposed
In Fig. 3(b), the total moving distance L (in meters) of the approach is characterized by the adaptive increase with the
star spot within the period T L yields increase of the angular velocity.
L  T L V star ; (3) In practice, the working mode of multiexposure imaging
within a period T L is shown in Fig. 4, where T on and T off
where V star is the linear velocity of the star spot in the image indicate the open time and the close time (in seconds) of
plane in m · s−1 . The width L1 (in meters) of the star spot the image intensifier, respectively, and their expressions are
formed in the exposure time T e of I t is expressed as given by
L1  T e V star  2 × 3σa  Le  6σa; (4) 
T on  T e
: (8)
where σ is the Gaussian radius in pixels, which represents the T off  T i − T e
spread scale of the optical lens, a is the length (in meters) of the
In Fig. 4, the imaging shutter Ct and the sampling shutter
square pixel of the image detector, and Le is the tailing length
I t within a period T L are, respectively, expressed as
(in meters) in the exposure time T e . The interval distance L2
(in meters) between two adjacent star spots is expressed as
Ct  1; 0 ≤ t ≤ T L; (9)
L2  T i V star : (5)
As shown in Fig. 3(b), in order for the two adjacent star
spots of the same star to be completely separated, the interval
distance L2 between them should yield C(t) TL
L2 ≥ L1 : (6)
By substituting Eqs. (4) and (5) into Eq. (6), we obtain
Toff / 2 Ton(Te) Toff / 2
6σa
Ti ≥ Te  : (7) Toff
V star I(t) ...
As shown in Eq. (7), the interval time T i exists at a lower 1 2 3 N
limit value; that is to say, the update rate F cannot be improved Ti
without limitation. However, it should be noted that the first t=0 t
item of Eq. (7) [i.e., the exposure time T e of I t ] adopts the Fig. 4. Working mode of multiexposure imaging.
Research Article Vol. 55, No. 36 / December 20 2016 / Applied Optics 10191

8 T off
< 1; 2  T i · n − 1 ≤ t ≤ T on  T2off  T i · n − 1 in the image plane. In this model, the analytical expression and
I t  n  1; 2; …; N  ; the variation of the linear velocity V star of the star spot are de-
: duced in detail. Then, based on it, three parameters (i.e., T on ,
0; t  other; and 0 ≤ t ≤ T L
T off , and N ) are optimized so as to obtain the optimal perfor-
(10) mance of the proposed approach.
where 1 indicates exposure opening and 0 indicates exposure
A. Motion Model of the Star Spot in the Image Plane
closing. As shown in Fig. 4, the relation between the exposure
time T L of Ct and the parameters of I t (i.e., the open time By observing the positions of stars, a star tracker can provide
T on and the close time T off ) is given by absolute attitude information with respect to the inertial coor-
dinate system. Let W i and V i , respectively, represent the ob-
T L  N · T on  T off   N · T i ; (11) servation unit vector and the reference unit vector of the ith
where N is the number of times of multiexposure imaging. For guide star in the field of view (FOV), and then both of them
a given period T L , the update rate F of the proposed approach are expressed as
is expressed as 8 " #
N >
> −x i
F
1
 : (12) >
> W i  pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi −y
>
>
1
i
Ti TL < x 2i y 2i f 2
f
" # ; (13)
Equations (3–12) describe the principle of the multiexpo- >
> cos δi cos αi
sure imaging approach. Among them, the exposure time T e >
>
>
> V  cos δi sin αi
and the interval time T i can be derived by Eq. (8). : i
sin δi
Furthermore, the total moving distance L, the width L1 ,
and the interval distance L2 are determined by Eqs. (3–5). where (x i , y i ) are the coordinates (in meters) of the ith star spot
In summary, as shown in Fig. 4, the proposed approach can in the image plane, f is the focal length (in meters) of the star
be completely determined by four parameters, namely, the tracker, and (αi , δi ) are the right ascension and declination in
period T L , the open time T on , the close time T off , and the degrees (°) of the ith guide star in the inertial coordinate system.
number of times N of multiexposure imaging. If A indicates the attitude matrix of the star tracker, the relation
Furthermore, with the increase of the angular velocity of the between W i and V i is given by
star tracker, the linear velocity V star of the star spot in the image
plane increases accordingly. Given that, the time for a star spot W i  AV i i  1; 2; …; M ; (14)
to remain in one pixel shortens gradually, resulting in energy
obtained by each pixel reducing and the SNR decreasing. In where M is the number of stars in the FOV. Here, the QUEST
this case, for the dim star, we can accumulate the energy of algorithm can be used to solve the optimal estimation of the
its N positions gained by the proposed approach. In this attitude matrix A [21].
way, the energy of the dim star spot can be effectively enhanced, Figure 5 shows the motion of the star spot in the image
and its SNR will be increased accordingly. plane. In Fig. 5, the star tracker works on dynamic conditions,
and the angular velocity ω ~ is expressed as
3. PARAMETER OPTIMIZATION ~   ωx
ω ωy ωz T
As mentioned before, the proposed approach can be completely  ω cos φ cos β cos φ sin β sin φ T ; (15)
determined by four parameters, namely, the period T L , the
open time T on , the close time T off , and the number of times
N of multiexposure imaging. Among them, the period T L is
generally a fixed value so as to facilitate the parallel and pipeline
processing. For a given period T L , the long close time T off will
reduce the number of times N , resulting in not achieving the
optimal performance of the proposed approach. On the other
hand, although the short close time T off increases the number
of times N , it may cause two adjacent star spots of the same star
to be overlapped with each other. In this case, the correspond-
ing attitude information cannot be calculated from each expo-
sure imaging, thus losing the significance of multiexposure.
Furthermore, the open time T on also has an impact on the
number of times N ; meanwhile, it determines whether each
exposure imaging can achieve the best effect. In summary,
the other three parameters (i.e., T on , T off , and N ) directly
affect the performance of the proposed approach.
By combining Eqs. (4), (7), (8), and (11), it can be seen that
the above three parameters are all determined by the linear
velocity V star of the star spot in the image plane. Therefore,
this section first establishes the motion model of the star spot Fig. 5. Motion of the star spot in the image plane.
10192 Vol. 55, No. 36 / December 20 2016 / Applied Optics Research Article

where ω is the length in degree/second (°/second) of ω,~ φ is the as well as the motion time Δt. Furthermore, the item
angle in degrees (°) between ω ~ and the image plane O 0 X Y , and [i.e., −x it ωyt Δt  y it ωxt Δt∕f ] is usually very small.
β is the angle in degrees (°) between the projection of ω~ on the According to the parameter calculation of the star tracker in
image plane O 0 X Y and X axis. Since the motion of the star this paper, this item is less than or equal to 0.01. Ignoring
tracker is equivalent to the position change of the guide star the above small item, Eq. (20) can be rewritten as
relative to the star tracker, the position of the star spot in 
x itΔt  x it  y it ωzt Δt  f ωyt Δt
the image plane moves from position P at time t to position : (21)
y itΔt  y it − x it ωzt Δt − f ωxt Δt
P 0 at time (t  Δt), where Δt (in seconds) is sufficiently small.
Furthermore, θ indicates the angle in degrees (°) between the When Δt is sufficiently small, the instantaneous linear
! velocity of the star spot in the image plane is derived as
observation vector PO of the guide star at time t and the bore-
! 8
sight O 0 Z of the star tracker, and γ indicates the angle in < V x  lim x itΔt Δt
−x it
 y it ωzt  f ωyt
! Δt→0
y itΔt −y it ; (22)
degrees (°) between the vector O 0 P and the X axis, as shown : V y  lim Δt  −x it ωzt − f ωxt
Δt→0
in Fig. 5.
According to Eq. (14), observation unit vectors of the ith where (V x , V y ) are the components in m · s−1 of the instanta-
guide star at different times t and (t  Δt) can be, respectively, neous linear velocity along the X axis and Y axis, respectively.
derived as follows: As shown in Eq. (22), the instantaneous linear velocity is the

W it  At V it function of the focal length f , the angular velocity ω,~ and the
: (16) coordinates (x i , y i ) of the star spot in the image plane.
W itΔt  AtΔt V itΔt
Composing (V x , V y ) and then simplifying the formulation
Since the reference vectors do not change with time, namely, by using the relation of angles in Fig. 5, the instantaneous linear
V it  V itΔt , the relation between the two observation unit velocity V star is given by
vectors is expressed as qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
W itΔt  AtΔt ATt W it  R tΔt W it ; (17) V star  V 2x  V 2y
t

where ATt is the rotation matrix at time t from the star tracker fω
coordinate system to the inertial coordinate system, and R tΔt
t
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
is the transfer matrix of the star tracker coordinate system from · tan2 θ sin2 φ  cos2 φ  tan θ sin 2φ cosβ − γ:
time t to time (t  Δt). When Δt is sufficiently small, the (23)
transfer matrix R tΔt
t satisfies [22] Here, the instantaneous linear velocity V star is expressed as
R tΔt
t  I  Δξ  OjΔξj Δξ  ωt Δt; (18) the function of the focal length f , the length ω of angular
velocity, and the angles including θ, γ, φ, and β. Among
where I is the identity matrix, Δξ is the angle variation within
them, the ranges of angles are, respectively, 0° ≤ θ ≤ θmax ,
the interval time Δt, and Δξ is the cross-product matrix of
0° ≤ γ < 360°, −90° ≤ φ ≤ 90°, and 0° ≤ β < 360°.
Δξ. Ignoring the higher-order infinitesimal of Δξ, Eq. (18) can
Furthermore, as the FOV of the star tracker is 20° in this
be rewritten as
2 3 paper, the maximum value θmax of the angle θ satisfies
1 ωzt Δt −ωyt Δt θmax  FOV∕2  10°.
R tΔt
t  I  ωt Δt  4 −ωzt Δt 1 ωxt Δt 5: As shown in Eq. (23), the instantaneous linear velocity V star
ωytΔt −ωxt Δt 1 of the star spot is mainly determined by the angular velocity
(19) ~ It consists of three items, namely, the item related to the
ω.
vertical component ω sin φ of the angular velocity ω, ~ the item
By substituting Eqs. (13) and (19) into Eq. (17), and related to the horizontal component ω cos φ, and the cross
considering that f is time-invariant, the simplified two- item of the above two. It should be noted that the first two
dimensional form of the relation between two observation unit items, namely, tan2 θ sin2 φ and cos2 φ, are both even functions
vectors at different times t and (t  Δt) can be rewritten as about the angle φ. Moreover, although the last item [i.e.,
 
x itΔt 1 tan θ sin 2φ cosβ − γ] is an odd function about the angle φ,

y itΔt −x it ωyt Δt  y it ωxt Δt∕f  1 it still has the character of symmetry about φ owing to
  cosβ − γ ∈ −1; 1. In summary, the linear velocity V star is
x  y it ωzt Δt  f ωyt Δt symmetrical about the angle φ. When φ ∈ −90°; 90°, or
× it : (20)
yit − x it ωzt Δt − f ωxt Δt jφj ∈ 0°; 90°, for a given jφj, there is a minimum value
In Eq. (20), the position of the star spot at time (t  Δt) is V starmin and a maximum value V starmax of the linear velocity
determined by the position and the angular velocity ω ~ at time t, V star , and they are, respectively, derived as

8  pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
<V  min f ω · tan2 θ sin2 φ  cos2 φ − tan θ sin 2jφj
starmin
: V
0≤θ≤θmax
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi : (24)
starmax  f ω · tan2 θmax sin2 φ  cos2 φ  tan θmax sin 2jφj
Research Article Vol. 55, No. 36 / December 20 2016 / Applied Optics 10193

In order to calculate the minimum value V starmin , let


hθ; φ be expressed as
hθ; φ  tan2 θ sin2 φ − tan θ sin 2φ
 tan θ sin φtan θ sin φ − 2 cos φ; (25)
where 0 ≤ θ < 90°, and 0 ≤ φ ≤ 90°. When φ ≠ 0°, let
hθ; φ  0; we can obtain θ1  0, and θ2  arc tan
2 cotφ, and they satisfy θ1 ≤ θ2 < 90°. According to
Eq. (25), it can be easily known that when θ > θ2 ,
hθ; φ > 0, and when 0  θ1 < θ < θ2 , hθ; φ < 0.
Therefore, for a given angle φ > 0, there is a minimum value
hθ; φmin  hθm ; φ < 0, where θ1 < θm < θ2 .
Solving the partial derivative of hθ; φ relative to θ and
letting it be 0, it can be expressed as follows:
δhθ; φ 1 1
 2 tan θ · 2 · sin2 φ − 2 · sin 2φ  0: (26) Fig. 6. Variation diagram of V starmin and V starmax with the
δθ cos θ cos θ
increase of the angle jφj.
From Eq. (26), we can obtain θm  90 − φ. Therefore,
when 0° ≤ θ ≤ θmax , for a given angle φ, the minimum value
0
hθ; φmin of hθ; φ is given by
 2
0 tan θmax sin2 φ − tan θmax sin 2φθm ≥ θmax  B. Open Time T on of the Image Intensifier
hθ; φmin  : As shown in Eq. (8), the open time T on is equal to the optimal
tan2 θm sin2 φ − tan θm sin 2φθm < θmax 
exposure time T e of a single exposure, and in this way each
(27)
sampled star spot can achieve the highest centroiding accuracy,
As can be seen from Eqs. (24) and (27), if hθ; φ gets the theoretically.
minimum value hθ; φmin0 within 0° ≤ θ ≤ θmax , V star gets the Furthermore, from Eq. (4), it should be noted that the op-
minimum value V starmin accordingly, which is expressed as timal exposure time T e can be derived from the optimal tailing
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
V starmin  f ω · cos2 φ  hθ; jφjmin 0
: (28) length Le , so it is only necessary to analyze Le . Yan et al. carried
out a detailed study on the optimal tailing length Le and the
By combining Eqs. (24), (27), and (28), the variation dia- optimal exposure time T e of a single exposure star spot [14].
gram of V starmin and V starmax of the linear velocity V star with From their study, two conclusions can be obtained. First, in
the increase of the angle jφj can be obtained, which is shown in theory, the optimal tailing length exhibits considerable varia-
Fig. 6. Here, in order to simplify the expression, let the coef- tion under conditions of different stellar magnitudes and differ-
ficient f ω be equal to 1 mm/s. As shown in Fig. 6, for a given ent angular velocities. It can be achieved only for a group of
angle φ, the linear velocities of different star spots in the FOV given parameters (i.e., a given stellar magnitude and a given
lie between V starmin and V starmax . In particular, when φ  0° angular velocity). However, the actual FOV contains many dif-
(i.e., the angular velocity ω ~ is totally located within the image ferent magnitude stars whose distribution is completely ran-
plane O 0 X Y ) the item related to the vertical component and dom, and the angular velocity of the star tracker may also
the cross item are both equal to zero, and then V star is equal to be changed. Therefore, there is no determined analytical ex-
f ω. In this case, the linear velocities of star spots in the image pression of the optimal tailing length Le that can be fully appli-
plane are all the same, regardless of their positions. Moreover, cable to all the working conditions of the star tracker. Second,
when jφj  90° (i.e., the angular velocity ω ~ is perpendicular to although Le does not have a determined analytical expression,
the image plane O 0 X Y ) the item related to the horizontal com- its range is determined {i.e., Le ∈ Lmin ; Lmax } where Lmin and
ponent and the cross item are both equal to zero and then V star Lmax indicate the lower limit value and the upper limit value in
is equal to f ω tan θ. In this case, the linear velocities of differ- meters, respectively, and Lmin  4.899σa, Lmax  14.289σa.
ent star spots exhibit great difference. For a star spot in the Based on the above two points, and taking into account the
center of the FOV, θ  0° and V star  0, but for a star spot fact that the linear velocity V star lies between V starmin and
in the edge of the FOV, θ  θmax and V star  f ω tan θmax . V starmax , in this paper the global optimal exposure time T e
Equations (16–28) describe the motion model of the star is defined as the exposure time within which the star spot with
spot in the image plane. Then based on it, three parameters the median linear velocity V starmiddle generates the median tail-
(i.e., T on , T off , and N ) are optimized as follows in the next ing length Lemiddle . Therefore, the global optimal exposure
subsections. time T e (i.e., the open time T on ) is expressed as

Lemiddle Lmax  Lmin ∕2 Lmax  Lmin 


T on  T e   
V starmiddle V starmax  V starmin ∕2 V starmax  V starmin
19.188σa
 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi : (29)
0
f ω ·  tan θmax sin φ  cos φ  tan θmax sin 2jφj  cos2 φ  hθ; jφjmin
2 2 2

10194 Vol. 55, No. 36 / December 20 2016 / Applied Optics Research Article

C. Close Time T off of the Image Intensifier Table 1. Design Parameters of the Star Tracker
As mentioned before, in order for the two adjacent star spots of
Parameter Value
the same star to be completely separated, the interval time T i
between them should satisfy Eq. (7). By substituting Eq. (8) Period T L (ms) 100
into Eq. (7), the close time T off yields Pixel length aμm 5.5
Focal length f mm 31.94
Half FOV angle θ max° 10
6σa Gaussian radius σ (pixel)
T off  T i − T e ≥ : (30) 0.5
V star

As shown in Eq. (30), if and only if two adjacent star spots In summary, three parameters of the multiexposure imaging
of the star with the minimum linear velocity V starmin are approach (i.e., T on , T off , and N ) are optimized by applying
completely separated, all the star spots in the FOV are not over- Eqs. (29–34), and only in this way is the proposed approach
lapped with each other. Therefore, the close time T off should able to obtain optimal performance. Here, a total of seven
satisfy parameters, namely, the period T L , the focal length f , the angle
θmax of half FOV, the length a of the square pixel, Gaussian
6σa radius σ, the length ω, and the angle φ of the angular velocity
T off ≥ T off min 
V starmin ~ are involved in the optimization process. Among them, the
ω,
6σa length ω and the angle φ are provided by the three axis strap-
 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 0 ≤ jφj ≤ 80°; (31) down gyros or the kinematic model of the star tracker, and the
0
fω· cos2 φ  hθ; jφjmin
others are the design parameters of the star tracker.
where T off min is the lower limit value of T off . It should be
noted that, when 80° < jφj ≤ 90°, V starmin  0, which is 4. SIMULATIONS AND ANALYSIS
shown in Fig. 6. Thus, the range of angle φ is 0 ≤ jφj ≤ In order to verify the feasibility and effectiveness of the multi-
80° in Eq. (31); otherwise, there is no T off to ensure that exposure imaging approach and parameter optimization, sim-
all the star spots in the FOV are not overlapped. ulations are carried out. The results of parameter optimization,
D. Number of Times N of Multiexposure Imaging
the improvement of attitude update rate, and the enhancement
By combining Eqs. (11), (29), and (31), the number of times N effect of dim star energy are shown and analyzed in this section.
of multiexposure imaging is derived as A. Parameter Optimization and Update Rate
8h i Improvement
< TL
T on T off min ; T on  T off min < T L
To intuitively understand the proposed approach, first the op-
N  ; (32) timization results of three parameters (i.e., T on , T off , and N )
: 1; T T ≥T
on off min L are simulated by using Eqs. (29–34); the involved design
parameters of the star tracker are shown in Table 1. Then, ac-
where · denotes the integer round-down operation. As cording to Eq. (12), the update rate F of the proposed approach
shown in Eq. (32), when multiexposure imaging cannot be is proportional to the number of times N .
accomplished within the total exposure time T L , the proposed
approach degrades into the existing single-exposure imaging. 1. Results Related to the Length ω of the Angular Velocity
After the number of times N is determined by Eq. (32), Let the angle φ of the angular velocity satisfy jφj  30°, and
according to Eqs. (11) and (29), the actual close time T off then the optimization results of the open time T on and the
is given by close time T off change with the increase of length ω, which
is shown in Fig. 7. Here, according to Eq. (34), it can be de-
TL rived that the minimum effective angular velocity yield
T off  − T on : (33) ωmin  0.55°∕s. The optimization result of the number of
N
times N and the corresponding update rate F of the proposed
It should be noted that, when the star tracker is considered approach are shown in Fig. 8.
to work in approximate static conditions, namely, the As shown in Figs. 7 and 8, for a given angle jφj, when the
angular velocity is very small, the corresponding open time length ω increases, the optimization results of the open time
T on determined by Eq. (29) satisfies T on ≫ T L , directly, T on and the close time T off decrease, and the optimization re-
leading to the multiexposure imaging not being applicable. For sult of the number of times N increases accordingly.
this reason, according to Eq. (29), a minimum effective Furthermore, from Table 1 it should be noted that the update
angular velocity ωmin for the proposed approach is defined rate of the existing single-exposure imaging approach satisfies
as follows: F  1∕T L  10 Hz. Therefore, compared with the former,

19.188σa
ωmin  pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi : (34)
f TL tan θmax sin φ  cos φ  tan θmax sin 2jφj  cos2 φ  hθ; jφjmin
2 2 2 0
Research Article Vol. 55, No. 36 / December 20 2016 / Applied Optics 10195

Fig. 7. Optimization results of T on and T off change with the Fig. 10. Optimization result of N and corresponding update rate F
increase of length ω. change with the increase of angle jφj.

2. Results Related to the Angle φ of the Angular Velocity


Let the length ω of the angular velocity satisfy ω  5°∕s, and
then the optimization results of the open time T on and the
close time T off change with the increase of angle jφj, which
is shown in Fig. 9. The optimization result of the number
of times N and the corresponding update rate F of the
proposed approach are shown in Fig. 10.
As shown in Figs. 9 and 10, when 0 ≤ jφj ≤ 66.9°, N > 1
(i.e., the multiexposure imaging approach is always effective).
When 66.9° ≤ jφj ≤ 80°, N  1. In this case, multiple expo-
sures cannot be inserted into the period T L , and the proposed
approach degrades into the single-exposure imaging approach.
When 80° < jφj ≤ 90°, according to Eq. (31), there is no T off
to ensure that all the star spots in the FOV are not overlapped.
Fig. 8. Optimization result of N and corresponding update rate F
change with the increase of length ω. In this case, the proposed approach cannot be applied, and only
the existing single exposure imaging approach is applicable. In
summary, for a given length ω, when the angle jφj increases,
the optimization results of the open time T on and the close
the update rate F of the proposed approach is improved from time T off increase, and the optimization result of the number
10 to 210 Hz with the increase of length ω from 0.55°/s to of times N decreases accordingly. As a result, the increase of
20°/s. In other words, the update rate F of the proposed angle jφj is not conducive to the update rate improvement,
approach is characterized by the adaptive increase with the or rather, the effective range of the angle jφj for the proposed
increase of length ω of the angular velocity. approach is [0, 66.9°].
B. Enhancement Effect of Dim-Star Energy
To verify the enhancement effect of dim-star energy, let the
angular velocity ω~ be equal to 5°∕s; 0°∕s; 0°∕sT , and then
the length ω  5°∕s and the angle jφj  0°. According to
Eqs. (29–34), it can be derived that the optimal exposure time
of single exposure is T e  9.4 ms, and the optimization results
of three parameters of the proposed approach are, respectively,
T on  9.4 ms, T off  7.2 ms, and N  6. By using the
above information and referring to Liebe’s stellar magnitude
model [6], when the Gaussian noise is added, the gray images
and histograms of simulated star images in different conditions
can be obtained, which are shown in Fig. 11.
By comparing (d), (e), and (h) with (b), (a), and (f ) in
Fig. 11, it can be seen that for a dim star, the energy of its
N positions obtained by the proposed approach can be accu-
Fig. 9. Optimization results of T on and T off change with the mulated, and then its SNR is effectively improved. Further
increase of angle jφj. analysis shows that the SNR of the dim star obtained by the
10196 Vol. 55, No. 36 / December 20 2016 / Applied Optics Research Article

Fig. 11. Gray images and histograms of simulated star images in different conditions. (a) and (f) show the original star image obtained by the
existing single exposure imaging approach, and (c) and (g) show the original one obtained by the proposed approach. By accumulating the energy of
N positions of the same dim star gained by the proposed approach, the enhanced star image is shown in figures (e) and (h). Figures (b) and (d) are the
partial detailed views of stars in (a) and (c), as well as (e), respectively.

Fig. 12. Laboratory experiment system.

existing single-exposure imaging approach in (b), (a), and (f ) is Fig. 13. Experiment results of the star images. (a) Original star
3.05, and the improved SNR in (d), (e), and (h) is 17.79. The image of the stellar magnitude  5.0 obtained by the proposed ap-
proach; (b) original star image of the star magnitude  8.0 obtained
improved times of the SNR of the proposed approach is ap-
by the proposed approach; (c) enhanced star image of figure (b).
proximately proportional to the number of times N where
N is equal to 6.
single simulated star can enter the FOV of the intensified star
5. EXPERIMENTS AND DISCUSSION tracker from different directions. In our experiment, the
Laboratory experiments are conducted to verify the proposed angular velocity of the three-axis rotary table is set as ω ~
approach. Figure 12 shows the experiment system, including 5.000°∕s; 0.000°∕s; 0.000°∕sT ; the X axis of the three-axis
the intensified star tracker, the three-axis rotary table, and rotary table is set to rotate from −30° to 30°, so as to keep
the star simulator. The intensified star tracker is installed the angular velocity constant in the range of the FOV from
on the three-axis rotary table, which can generate quite −10° to 10°. Furthermore, the gain voltage of the image in-
accurate rotation angles and angular velocities; thereafter, the tensifier is set as 3.8 V; the multiexposure parameters are set as
Research Article Vol. 55, No. 36 / December 20 2016 / Applied Optics 10197

T L  100 ms, T on  9.4 ms, T off  7.2 ms, and N  6, REFERENCES


respectively, which is the same as the parameters in 1. C. C. Liebe, “Star trackers for attitude determination,” IEEE Trans.
Subsection B of Section 4. Aerosp. Electron. Syst. 10, 10–16 (1995).
Figure 13 shows the experiment results of the star images. 2. R. W. H. van Bezooijen, “SIRTF autonomous star tracker,” Proc. SPIE
4850, 108–121 (2003).
From (a) and (b), it can be demonstrated that the proposed
3. T. Sun, F. Xing, X. Wang, Z. You, and D. Chu, “An accuracy
approach and its parameter optimization are feasible and measurement method for star trackers based on direct astronomic
effective in the real application; besides, (c) indicates that observation,” Sci. Rep. 6, 22593 (2016).
the energy of the N positions of the same dim star can be 4. J. Shen, G. Zhang, and X. Wei, “Simulation analysis of dynamic work-
accumulated, thereby further improving the SNR of the ing performance for star trackers,” J. Opt. Soc. Am. A 27, 2638–2647
(2010).
dim star. 5. L. Ma, D. Zhan, G. Jiang, S. Fu, H. Jia, X. Wang, Z. Huang, J. Zheng,
F. Hu, W. Wu, and S. Qin, “Attitude-correlated frames approach for a
star sensor to improve attitude accuracy under highly dynamic
6. CONCLUSIONS conditions,” Appl. Opt. 54, 7559–7566 (2015).
6. C. C. Liebe, “Accuracy performance of star trackers: a tutorial,” IEEE
A multiexposure imaging approach for intensified star trackers Trans. Aerosp. Electron. Syst. 38, 587–599 (2002).
is proposed in this paper. By utilizing the two-level exposure 7. T. Sun, F. Xing, Z. You, and M. Wei, “Motion-blurred star acquisition
method of the star tracker under high dynamic conditions,” Opt.
control structure of the intensified image detector, the proposed
Express 21, 20096–20110 (2013).
approach breaks through the limitation of the transmission 8. W. Hou, H. Liu, Z. Lei, Q. Yu, X. Liu, and J. Dong, “Smeared star spot
and processing time of pixel data on the attitude update rate. location estimation using directional integral method,” Appl. Opt. 53,
Subsequently, in order to obtain the optimal performance of 2073–2086 (2014).
the proposed approach, parameter optimization is carried out. 9. Y. Liao, E. Liu, J. Zhong, and H. Zhang, “Processing centroids of
smearing star image of star sensor,” Math. Probl. Eng. 2014,
Simulations and experiments are implemented and analyzed 534698 (2014).
so as to verify the feasibility and effectiveness of the proposed 10. C. Liu, L. Hu, G. Liu, B. Yang, and A. Li, “Kinematic model for the
approach and parameter optimization. Compared with the space-variant image motion of star sensors under dynamical condi-
existing single-exposure imaging approach, the proposed ap- tions,” Opt. Eng. 54, 063104 (2015).
11. J. Zhang, Y. Hao, L. Wang, and D. Liu, “Studies on dynamic motion
proach improves the attitude update rate by N times where compensation and positioning accuracy on star tracker,” Appl. Opt.
N is the number of times of multiexposure imaging. For a 54, 8417–8424 (2015).
given angle jφj  30°, when the length ω of the angular veloc- 12. T. Sun, F. Xing, Z. You, X. Wang, and B. Li, “Smearing model and
ity increases from 0.55°/s to 20°/s, the optimization result of restoration of star image under conditions of variable angular velocity
the number of times N increases from one to 21 and the cor- and long exposure time,” Opt. Express 22, 6009–6024 (2014).
13. X. Wei, W. Tan, J. Li, and G. Zhang, “Exposure time optimization for
responding update rate F is improved from 10 to 210 Hz. highly dynamic star trackers,” Sensors 14, 4914–4931 (2014).
Therefore, the update rate F of the proposed approach is char- 14. J. Yan, J. Jiang, and G. Zhang, “Dynamic imaging model and param-
acterized by the adaptive increase with the increase of length eter optimization for a star tracker,” Opt. Express 24, 5961–5983
ω of the angular velocity. In addition, for a given length (2016).
15. A. B. Katake, “Modeling, image processing and attitude estimation of
ω  5°∕s, the optimization result of the number of times N high speed star sensors,” Ph.D. dissertation (Texas A&M University,
decreases with the increase of angle jφj, and the effective range 2006).
of the angle jφj for the proposed approach is [0, 66.9°]. Here, it 16. A. Katake and C. Bruccoleri, “StarCam SG100: a high update rate,
should be noted that, by adopting two or three star trackers high sensitivity stellar gyroscope for spacecraft,” Proc. SPIE 7536,
753608 (2010).
whose boresights are perpendicular to each other, at least
17. H. Zhong, M. Yang, and X. Lu, “Increasing update rate for star
one star tracker can be ensured to have a small angle jφj at sensor by pipelining parallel processing method,” Opt. Precis. Eng.
the same time. In this case, the adverse effect of the increasing 17, 2230–2235 (2009).
angle jφj on the number of times N will be weakened to a great 18. X. Mao, W. Liang, and X. Zheng, “A parallel computing architecture
extent. Furthermore, for a dim star, the energy of its N posi- based image processing algorithm for star sensor,” J. Astronaut. 32,
613–619 (2011).
tions obtained by the proposed approach can also be accumu- 19. Y. Jin, J. Jiang, and G. Zhang, “Three-step nonuniformity correction
lated, and then its SNR is effectively improved. In this way, the for a highly dynamic intensified charge-coupled device star sensor,”
improved times of the SNR is approximately proportional to Opt. Commun. 285, 1753–1758 (2012).
the number of times N of multiexposure imaging. 20. K. Xiong and J. Jiang, “Reducing systematic centroid errors induced
by fiber optic faceplates in intensified high-accuracy star trackers,”
Sensors 15, 12389–12409 (2015).
Funding. National Natural Science Foundation of China 21. M. D. Shuster and S. D. Oh, “Three-axis attitude determination from
(NSFC) (61222304); Specialized Research Fund for the vector observations,” J. Guid. Control Dyn. 4, 70–77 (1981).
Doctoral Program of Higher Education of China 22. M. D. Shuster, “A survey of attitude representations,” J. Astronaut. Sci.
(20121102110032). 41, 439–517 (1993).

You might also like