Modeling and Analysis of an Echo Laser Pulse Waveform for the Orientation Determination of Space Debris

: Orientation information of space debris is required to improve the orbital prediction accuracy for mitigation or elimination of a signiﬁcant threat to not only human space activities but also operational satellites. Obtaining orientation information is currently achievable by applying photometry, adaptive optics (AO) and satellite laser ranging (SLR) technologies. In this study, a new method is proposed based on an echo laser pulse waveform (ELPW) for the orientation determination of space debris; its feasibility was also investigated by numerical simulations. Unlike the photometry and AO technologies available just under the sun-illumination condition and the SLR technology applicable only for cooperative targets, the ELPW is achievable by using a high power laser regardless of the above measurement constraints. A mathematical model is derived to generate the ELPW, and the beam broadening and spreading due to the atmospheric turbulence is taken into account. The Gaussian decomposition based on a genetic algorithm was employed to the ELPWs in order to analyze the orientation features. It is demonstrated from the numerical simulations that the ELPWs have distinctive shapes characterizing the orientation of space debris and therefore our approach was capable of providing orientation information.


Introduction
Space debris has become an important global problem because the exponential growth of its population poses a significant threat to not only human space activities but also operational satellites. Improving the orbital prediction accuracy of space debris is required to mitigate or eliminate the threat through avoidance of catastrophic collisions, which is attainable by using various measurements from the optical, radar and laser tracking systems. Even though these tracking systems provide accurate measurements, the dynamic model has to take into consideration the space debris characteristics such as size, shape and orientation information in order to improve the orbital prediction accuracy. Unlike the gravitational perturbations, non-gravitational perturbations such as the atmospheric drag and solar radiation pressure are dominated by these characteristics [1]. Thus the orientation determination including rotation rates and axis of space debris became a topic of interest in the space situational awareness communities.
There are several methods to achieve orientation information of space debris in terms of the ground-based observations: photometry, adaptive optics (AO) and satellite laser ranging (SLR). Photometry is literally the measurement of light intensity or brightness and provides light curves which are commonly used in the astronomical communities to extract physical information of asteroids. Recently light curves have been recognized as a promising tool to determine shapes, spin axis orientation and spin rate of space debris by analyzing amplitude variation as well as periodic difference, thus are widely used by many researchers [2][3][4][5][6][7]. The AO is a technology to improve the image quality nearly to diffraction limited resolution by compensating the incoming wavefronts distorted by the atmospheric turbulence. So the high-resolution optical images obtained from the AO system provide a means of characterizing size, shape and orientation of space debris and its identification [8][9][10].
As opposed to the passive technologies (i.e., photometry and AO) which receive sunlight reflected by space debris and so work only at night (in particular, twilight for space debris in low earth orbit), the SLR technology actively fires laser pulses into space objects and then measures the distance between the ground station and space objects equipped with corner cube reflectors (CCRs). The spinning array of the CCRs causes a millimeter-scale distance variation from which information of spin axis orientation and spin rate is achievable [11] for not only geodetic satellites with CCRs on the outer surface [12][13][14] but also defunct satellites equipped with CCRs [15][16][17]. Moreover, the SLR technology combined with photometry has been applied to space debris without CCRs, called uncooperative targets [18][19][20].
Debris laser ranging (DLR) is similar to SLR but applicable to uncooperative targets by employing a higher power laser rather than the SLR (normally 2 or 3 orders of power). It generates ranging data even under the low signal-to-noise ratio resulting from a small laser cross-section of uncooperative targets. For DLR, a Geiger-mode avalanche photodiode, called a photon-counting detector, is applied to respond to the presence of return signals and then records accurate timing of the first arrived photon [21][22][23][24]. In spite of benefits from the photon-counting detector with higher detection sensitivity and lower dark noise, it does not provide an echo laser pulse waveform (ELPW) which records the scattering signals as a function of time. The ELPWs are achievable for uncooperative targets by employing a linear-mode avalanche photodiode (LmAPD) and a much higher power laser on the DLR system because the LmAPD output current is linearly proportional to photocurrent and thus can provide intensity information as well as ranging measurements. An ELPW allows not only range information but also 3-D shape in the direction of laser beam propagation by analyzing the features of the scattering signals, without an optical scanning system [25].
In this study, a new method based on the ELPW analysis was proposed for orientation determination of space debris and its feasibility was also investigated by numerical simulations of two space targets (i.e., defunct satellite and rocket body). Unlike the passive and SLR technologies, the proposed method is free of the measurement constraints such as sun-illumination and CCRs requirement. Approximating that the reflected surface of a space target consists of small grids, a mathematical model was derived to generate the ELPW for the fully reflected surface, which takes into consideration temporal beam broadening and spatial beam spreading due to the atmospheric turbulence. The simulation took into account detector noise and amplifier noise by assuming that the noises are independent and of Gaussian distribution in the optical power domain. A genetic algorithm was also introduced to determine the Gaussian parameters of Gaussian decomposition and to analyze the relationship between ELPW features and orientation of the space target in terms of Gaussian parameters. It was shown from the numerical simulations that the proposed method was a feasible alternative to achieve orientation information because the ELPW had a distinctive shape characterizing the orientation of space debris.

Echo Laser Pulse Waveform
The intensity of a laser beam is not uniform in temporal and spatial domains along the direction of laser beam propagation, which requires that the laser pulse profile should be modeled in time as well as in space. The laser beam propagation can be approximated by assuming that the laser beam has a Gaussian shape in terms of temporal and spatial domains. So the irradiance distribution and normalized waveform of the propagated laser pulse are expressed as [26,27] where d is the radial distance from the axis of propagation (i.e., d 2 = x 2 + y 2 ), z is the axial distance from the beam's narrowest point called waist, P T is the transmitted power, W e f f (z) is the effective beam size referred to as the Gaussian beam radius, and τ 0 is the transmitted pulse width which is related to the full width at half maximum (FWHM) as follows: FWHM = 2 √ 2 ln 2τ 0 [28]. The reflected or intersected surface of space debris by the laser beam is assumed to consist of small grids in order to generate a total ELPW. Some of the laser power received by the ith grid located at (x i , y i ) is absorbed while the rest is reradiated back into space, in proportion to the effective receiving area depending on the incident angle θ i . Therefore the scattered power by the ith grid is described by where A i is the ith grid area, and ρ is the reflectivity on the assumption that all grids have the same reflectivity in this study. Large grid areas can save computation time but might provide ELPW features less accurate from the true orientation of space target. So all grids should have optimal areas so that the ELPW simulation closely resembles a practical situation with minimal computation time, smaller than 50 cm 2 in this study. If the incoming radiation to the receiver is scattered uniformly into a cone of solid angle Ω [29], the laser power coming into the receiver (i.e., received power) is given by where A R is the receiving area, and R i is the range to the ith grid. The solid angle can be replaced with π by considering the grid to be a diffuse scattering target, or a Lambertian target. Taking into consideration the normalized waveform given by Equation (2), the received power in terms of the time domain can be written by where c is the light speed in a vacuum, and τ R is the received pulse width, which is broadened by the atmospheric turbulence. The echo waveform model was also given by Hao et al. [30] for the case of specular targets such as mirrors, which increases significantly the received power compared to diffuse targets. Space debris and most satellites are considered as diffuse targets since they are non-cooperative targets without CCRs and the echo signals are diffused from their surface. In the practical DLR system, the transmitted laser beam suffers from various attenuations by the time it comes back to a detector, e.g., atmospheric attenuation, and optical transmittance. Taking into account various attenuations in Equation (5), the mathematical equation of ELPW is obtained by summing the received power corresponding to a grid as follows: where N is the total number of grids illuminated by the transmitted laser beam, T A and T C are the one-way transmittance by the atmospheric attenuation and cirrus clouds, respectively, η T is the transmitting optical efficiency, η R is the receiving optical efficiency, and η F is the optical filter transmittance.
In the near-ultraviolet to visible spectral band (0.3~0.7 µ m), the atmospheric attenuation is dominated by aerosol and molecular scattering as well as ozone absorption. In the near infrared band beyond 0.7 µm, attenuation is subjected to absorption features of atmospheric constituents including water vapor, oxygen and carbon dioxide. Typically, cirrus clouds are overhead about 50% of the time at most locations. The one-way transmittances by the atmospheric attenuation and cirrus clouds are expressed as [31] where σ(λ, V, 0) and V are the atmospheric attenuation coefficient and the visibility at sea level, respectively, λ is the laser beam wavelength, h SH is the scale height of 1.2 km, θ Z is the zenith angle of space target, h s is the station height, and C T is the cirrus cloud thickness.

Laser Beam Spreading and Broadening
When the laser beam propagates through the atmosphere, the atmospheric refractive index is randomly fluctuated due to the atmospheric turbulence, resulting in the intensity fluctuations known as scintillation, phase distortion, beam spot wandering, beam spreading and broadening. In particular, the spatial spreading and temporal broadening of laser beams are limiting factors in most applications such as remote sensing and free space optical communications [32]. The effective beam size of a Gaussian beam in a horizontal path is approximated in a far-field region as follows [33]: where: T = 4.35k 7/6 (ΛH) 5/6 sec 11/6 (θ z ) where T is the additional parameter corresponding to the beam spreading by the atmospheric turbulence, k is the optical wavenumber, W 0 is the beam size at z = 0 (called beam waist), θ D is the far-field divergence half-angle, Λ is the Fresnel ratio defined as Λ = 2z/kW 2 (z), H is the space debris altitude from the station, and C 2 n (h) is the refractive-index structure parameter which represents the intensity of atmospheric turbulence. Many models have been proposed to predict the behavior of the refractive-index structure parameter. The Hufnagel-Valley model is widely used in satellite optical communications and optical remote sensing because it allows easy variation of day-time and night-time profile by varying various site parameters such as wind speed, isoplanatic angle and altitude [34]. The model is given by where v is the mean square value of the wind speed in m/s, and A 0 is the nominal value of the refractive-index structure parameter at the ground in m −2/3 . The temporal broadening of a laser pulse through the atmosphere is caused physically by two mechanisms: scattering process and beam spot wandering; these are too important to be neglected for high bandwidth signals such as narrow pulses in the range from femtosecond to picosecond [35]. The received pulse width broadened by the atmospheric turbulence is given by [32] where: where β is the inner scale factor, 0 < β < 1, and l 0 (h) and L 0 (h) represent the inner and outer scale sizes as a function of altitude, respectively.

Noise Equivalent Power
The received signal is subject to amplifier and detector noises including background noise, shot noise and dark current noise, all of which should be considered in the ELPW simulation. The noise equivalent power (NEP) is a common metric that quantifies the sensitivity of a photodetector or the power generated by a noise source, which means the minimum optical power required to achieve a unity signal-to-noise ratio in a 1 Hz bandwidth. To take into account noise in the ELPW simulation, a Gaussian random variable with a zero mean and a standard deviation equal to NEP is generated and then added to the received power because the noises are independent and identically distributed Gaussian random variables [36]. The NEP resulting from the detector noise is given by where q is the electronic charge, B is the electrical bandwidth, I ds and I db are the surface and bulk dark currents, respectively, I b and I s are the currents due to the background noise and received signal, respectively, G is the detector gain, F A is the excess noise factor, is the detector responsivity defined as = qη Q G/hv using the quantum efficiency of the detector η Q , h is Plank's constant and the laser beam frequency is v. The excess noise factor is expressed as F A = k e f f G + 1 − k e f f (2 − 1/G) in terms of the detector gain and the ionization coefficient ratio k e f f . Two currents owing to the received signal and background noise in Equation (18) can be obtained from [37] I s (t) = P r (t), Remote Sens. 2020, 12, 1659 6 of 17 where I λ is the spectral radiance of the sky, Ω FoV is the receiver field of view in steradians, and ∆λ F is the bandwidth of the optical bandpass filter.
The NEP corresponding to the amplifier noise is given by [35] NEP where κ is Boltzmann's constant, T is the temperature in Kelvin, F N is the noise factor, and R L is the load resistance. The NEPs for the detector and amplifier are combined into the total NEP which is interpreted as the standard deviation for the Gaussian distribution of the additive noises in the optical power domain:

Genetic Algorithm for Gaussian Decomposition
The returned laser pulses are distorted by the cross-section characteristics (i.e., the scattering properties of targets) even though the transmitted pulses are Gaussian functions. However, it is assumed that the returned laser pulse from a small grid follows a Gaussian distribution in the time domain and then the ELPW summing the returned pulses is also represented by a series of Gaussian functions [29]. According to the analysis work of returned waveform [38], more than 98% of the observed waveforms with a topographic lidar could be correctly fitted with the sum of Gaussian functions. Therefore, the ELPW can be fitted with a mixture of Gaussian functions as follows [39]: where n is the number of Gaussian functions, and Gaussian parameters of X i = A i , µ i , σ i are the amplitude, position and half-width of the ith Gaussian function, respectively. The Gaussian parameters are determined by the Gaussian decomposition so that Equation (23) fits the ELPW given by Equation (7) with a prescribed accuracy ε using some number of Gaussian functions: where t k is the discrete time, and N T is the total number of discrete time intervals. There are several methods for the Gaussian decomposition, e.g., non-linear least square approach, maximum likelihood estimate, and stochastic approach [38]. The genetic algorithm (GA) is employed for the Gaussian decomposition in this study, which is a parameter optimization technique inspired by natural evolution using the concepts of mating, mutation, crossover, and selection [40,41]. The GA is briefly addressed in this section because it is only used as a tool to find Gaussian parameters. The GA works by building a population of chromosomes which is a set of potential solutions to Equation (24). The fitness ζ(X, t) of each chromosome in the population is evaluated and then a new population is created by repeating the following steps until the fitness of a new chromosome is satisfied: selection, crossover, and mutation. A probabilistic selection is performed based on the chromosome's fitness such that the better chromosomes have an increased chance of being selected to reproduce into the next generation. The selected chromosomes, generally called parents, are to be mixed and recombined for the production of offspring in the next generation. In order to facilitate the GA evolution cycle, there are two fundamental operators: crossover and mutation. Based on the crossover and mutation probabilities, crossover takes two chromosomes and produces two new offsprings while mutation alters one chromosome to produce a single new offspring. A more complete explanation of the GAs can be found in the book by Davies [40]. The Gaussian decomposition based Remote Sens. 2020, 12, 1659 7 of 17 on the GA was implemented using the Matlab toolbox [42] and the flowchart of the GA process is shown in Figure 1.

Simulation and Analysis of Laser Pulse Waveform
Two types of space debris were simulated to investigate the relationship between ELPW and orientation: defunct satellite (dimensions: 4.2 × 1.3 × 1.3 m) and rocket body (dimensions: 24.9 × 10.0 × 10.0 m). As shown in Figure 2, representing 3-dimensional models available from the NASA resources (https://nasa3d.arc.nasa.gov), the defunct satellite consists of a bus module, a solar panel, and small payloads, while the rocket body consists of five rocket engines, forward and aft skirts, tank and thrust structures, and an aft interstage. To generate the ELPW similar to the real measurement environment, the surface is made up of 113,241 grids for the defunct satellite, and 3,232,611 grids for the rocket body, respectively.

Simulation and Analysis of Laser Pulse Waveform
Two types of space debris were simulated to investigate the relationship between ELPW and orientation: defunct satellite (dimensions: 4.2 × 1.3 × 1.3 m) and rocket body (dimensions: 24.9 × 10.0 × 10.0 m). As shown in Figure 2, representing 3-dimensional models available from the NASA resources (https://nasa3d.arc.nasa.gov), the defunct satellite consists of a bus module, a solar panel, and small payloads, while the rocket body consists of five rocket engines, forward and aft skirts, tank and thrust structures, and an aft interstage. To generate the ELPW similar to the real measurement environment, the surface is made up of 113,241 grids for the defunct satellite, and 3,232,611 grids for the rocket body, respectively.

Simulation and Analysis of Laser Pulse Waveform
Two types of space debris were simulated to investigate the relationship between ELPW and orientation: defunct satellite (dimensions: 4.2 × 1.3 × 1.3 m) and rocket body (dimensions: 24.9 × 10.0 × 10.0 m). As shown in Figure 2, representing 3-dimensional models available from the NASA resources (https://nasa3d.arc.nasa.gov), the defunct satellite consists of a bus module, a solar panel, and small payloads, while the rocket body consists of five rocket engines, forward and aft skirts, tank and thrust structures, and an aft interstage. To generate the ELPW similar to the real measurement environment, the surface is made up of 113,241 grids for the defunct satellite, and 3,232,611 grids for the rocket body, respectively. We performed a numerical simulation by considering the slant range R O of 500 km from the ground station to the origin of body frame, a zenith angle θ z of 50 degrees, the ground station height h s of 1000 m, and the grid reflectivity ρ of 0.175. Taking the clear sky (V = 15 k m) and 1,064 nm wavelength of the laser beam into account, the atmospheric attenuation coefficient is σ(λ, V, 0) = 0.16, and the cirrus cloud thickness C T of 0.5 km is considered.
The following parameters related to the NEP are also used for the numerical simulation: the electrical bandwidth B = 100MHz, the surface and bulk dark currents I ds = 20nA and I db = 50nA, respectively, the detector gain G = 30, the excess noise factor F A = 0.25, the detector quantum efficiency η Q = 0.38, the spectral radiance of the sky I λ = 4.0 × 10 −8 watts/cm 2 · ster · µm, the receiver field of view Ω FoV = 25arcsec, the bandwidth of the optical bandpass filter ∆λ F = 1nm, the amplifier temperature T = 19 • C, the noise factor F N = 1.5, and the load resistance R L = 1, 000Ω. In addition, the optical parameters related to the transmitter and receiver are summarized in Table 1. Considering that the laser beam is transmitted from the ground station at sea level to the space debris located at the slant range of 500 km, Figure 3 shows the beam spreading and broadening in terms of the atmospheric turbulence condition, zenith angle, and inner scale factor. The beam spreading due to the atmospheric turbulence has a considerable influence on the received power because the beam size can increase more than 2 times under the strong turbulence condition (A 0 = 1.5 × 10 −13 m −2/3 , ν = 30m/s) and high zenith angle (θ Z = 70deg). The beam broadening in time is little affected by the inner scale factor even though it suffers from the turbulence condition as well as the zenith angle. It is shown that the beam broadening effect decreases as the transmitting pulse width increases and that the effect also increases as the turbulence strength and zenith angle increase. In contrast to the nanosecond regime of the transmitted pulse width, the beam broadening is too significant in the femtosecond regime to be neglected for the ELPW analysis. However, the beam broadening is negligible in this study because the transmitted pulse width is in the order of nanoseconds. Regarding the beam spreading, the most commonly used values (A 0 = 1.7 × 10 −14 m −2/3 , ν = 21m/s) are employed for the ELPW simulation.
Because of the small beam divergence, the DLR system controls the laser beam to the direction of a pre-determined position in space based on the predicted ephemeris. However, the target of space debris may deviate a little from the axis of beam propagation even though the receiver detects the echo signal from the target located within the coverage of beam divergence. The pointing error illustrated in Figure 4a can be reduced if an acquisition camera is employed to detect space debris illuminated by sunlight, which has an effect on the ELPW intensity as shown in Equation (1). Assuming that the space target is translated from the axis of beam propagation to the y-axis with pointing error φ, Equation (6) can be rewritten as Remote Sens. 2020, 12, 1659 9 of 17 where x b,i , y b,i is the coordinate of the ith grid on the body frame, and R O is the distance to the origin of the body frame. Figure 4b shows the ELPWs of defunct satellite in terms of the pointing error. Compared to the ELPWs for three pointing errors (φ = 0, 3, 5arcseconds), the ELPW intensity decreases with the same shapes as the pointing error increases. In addition, two pointing errors (φ = 5, −5arcseconds) cause some ELPW shape differences because the reflected surface of the defunct satellite is not symmetric with respect to the z-axis. But the difference of two shapes can be mitigated if a flattened Gaussian beam in which the energy distribution is nearly identical on the cross-section of the laser beam [43] is used. Thus the pointing error was not considered in the ELPW simulation.
Considering that the laser beam is transmitted from the ground station at sea level to the space debris located at the slant range of 500 km, Figure 3 shows the beam spreading and broadening in terms of the atmospheric turbulence condition, zenith angle, and inner scale factor. The beam spreading due to the atmospheric turbulence has a considerable influence on the received power because the beam size can increase more than 2 times under the strong turbulence condition (

). The beam broadening in time is
little affected by the inner scale factor even though it suffers from the turbulence condition as well as the zenith angle. It is shown that the beam broadening effect decreases as the transmitting pulse width increases and that the effect also increases as the turbulence strength and zenith angle increase. In contrast to the nanosecond regime of the transmitted pulse width, the beam broadening is too significant in the femtosecond regime to be neglected for the ELPW analysis. However, the beam broadening is negligible in this study because the transmitted pulse width is in the order of nanoseconds. Regarding the beam spreading, the most commonly used values (  Because of the small beam divergence, the DLR system controls the laser beam to the direction of a pre-determined position in space based on the predicted ephemeris. However, the target of space debris may deviate a little from the axis of beam propagation even though the receiver detects the echo signal from the target located within the coverage of beam divergence. The pointing error illustrated in Figure 4a can be reduced if an acquisition camera is employed to detect space debris illuminated by sunlight, which has an effect on the ELPW intensity as shown in Equation (1). Assuming that the space target is translated from the axis of beam propagation to the y-axis with pointing error  , Equation (6) can be rewritten as is the coordinate of the ith grid on the body frame, and O R is the distance to the origin of the body frame. Figure 4b shows the ELPWs of defunct satellite in terms of the pointing error. Compared to the ELPWs for three pointing errors ( arcseconds 5 , 3 , 0   ), the ELPW intensity decreases with the same shapes as the pointing error increases. In addition, two pointing errors ( arcseconds 5 , 5    ) cause some ELPW shape differences because the reflected surface of the defunct satellite is not symmetric with respect to the z-axis. But the difference of two shapes can be mitigated if a flattened Gaussian beam in which the energy distribution is nearly identical on the cross-section of the laser beam [43] is used. Thus the pointing error was not considered in the ELPW simulation. The analytic expression related to the pulse broadening for laser altimeters [44] is introduced for the ELPW analysis: where h  is the rms width of the receiver impulse response,   is the surface roughness, and  is the nadir angle of the laser beam.
S is the surface slope that consists of two components, namely a surface slope parallel to the nadir direction ( || S ) and a surface slope normal to nadir direction (  S ). 5 due to the solar panel slope parallel to the beam direction ( || S ), even though the cross-section area of the solar panel is smaller than that in Figure 5. The second curve in the ELPW was decomposed into 2 or 3 Gaussian functions, which can be also decomposed into many Gaussian functions by adding small Gaussian functions such as the third Gaussian function in Figure 9b. The complex reflected surface requires many Gaussian functions for the Gaussian decomposition and then makes the explicit analysis difficult in relation to the orientation of space debris. The analytic expression related to the pulse broadening for laser altimeters [44] is introduced for the ELPW analysis: where σ h is the rms width of the receiver impulse response, ∆ξ is the surface roughness, and ϕ is the nadir angle of the laser beam. S is the surface slope that consists of two components, namely a surface slope parallel to the nadir direction (S || ) and a surface slope normal to nadir direction (S ⊥ ). According to Equation (26), the received pulse width is also broadened by the nadir angle and the target surface characteristics: surface roughness and slope of the ground target, and beam curvature. The nadir angle is negligible in this debris study because the reflected surface of the space debris is so small compared to the footprint of the laser altimeter that the reflected surface area is the same regardless of the nadir angle and thus makes no impact on the pulse broadening. In place of the nadir angle, the reflected surface area of space debris has an influence on the pulse broadening. Figures 5-8 show the echo waveforms as well as the Gaussian decomposition results for four orientations of the defunct satellite. The starting point of the range gate is determined to be arbitrary but the same for the simulations of four orientations (i.e., t = 0 in the x-axis) in order to more easily compare and analyze the ELPWs. Unlike the Gaussian decomposition in Figure 8, the ELPWs in Figures 5-7 were decomposed into two Gaussian functions because the reflected surface is characterized by two main elements: the bus module and the solar panel. The first Gaussian function is related to the bus module and the second one corresponds to the solar panel because the range to the bus module is shorter than that to the solar panel (i.e., µ 1 < µ 2 ). The half-width of the first Gaussian function (σ 1 ) in Figure 6 is larger than that in Figure 5 because the cross-section area of the reflected bus module is larger and consists of two surfaces of the bus module which cause larger roughness. The half-width of the second Gaussian function (σ 2 ), corresponding to the solar panel in Figure 6, is also larger than that in Figure 5 even though the cross-section area of the reflected solar panel is smaller than that in Figure 5. The larger value of σ 2 results from the solar panel slope being normal to the beam direction (i.e., S ⊥ ) due to the 45 • rotation about the x-axis. In the case of the defunct satellite orientation with -90 • rotation about the x-axis as shown in Figure 7, the amplitude of the second Gaussian function is very small because the solar panel surface is perfectly parallel to the beam direction and thus the cross-section area is small. As shown in Figure 8, the second Gaussian function cannot fit the second curve in the ELPW well, which implies that the reflected surface is more complex than shown in Figures 5-7 and the second curve is not separated into the bus module and solar panel elements. Figure 9 shows two Gaussian decomposition results based on 3 and 4 Gaussian functions for the ELPW given in Figure 8, both of which fit the ELPW well. According to the orientation shown in Figure 8b, it is obvious that the first Gaussian function corresponds to the bus module and that the others result from the combination of the bus module and solar panel elements. The half-width of the Gaussian functions (σ 2 and σ 3 in Figure 9a, σ 2 and σ 4 in Figure 9b) is larger than σ 2 of Figure 5 due to the solar panel slope parallel to the beam direction (S || ), even though the cross-section area of the solar panel is smaller than that in Figure 5. The second curve in the ELPW was decomposed into 2 or 3 Gaussian functions, which can be also decomposed into many Gaussian functions by adding small Gaussian functions such as the third Gaussian function in Figure 9b. The complex reflected surface requires many Gaussian functions for the Gaussian decomposition and then makes the explicit analysis difficult in relation to the orientation of space debris.            The Gaussian decomposition results including the ELPWs are also shown in Figures 10-13 for four kinds of orientations of the rocket body. Unlike the ELPWs of the defunct satellite, the size of the rocket body is so large that the noise power is very weak compared to the received power. In addition, the reflected surface is too complex compared to the defunct satellite to fit the ELPWs well even by using several Gaussian functions. Therefore, the total number of Gaussian functions ( n ) is selected as the smallest number that satisfies the criterion of Equation (24). As shown in Figure 10, the first two Gaussian functions corresponding to the first peak represent the forward and aft skirts as well as the tank structure because their shapes are cylindrical, and the other Gaussian functions correspond to five rocket engines, the thrust structure and the aft interstage. The ELPWs in Figures  11-12 are so complex that each Gaussian function cannot be exactly in one-to-one correspondence with the components of the rocket body. Compared to the narrow distribution of ELPW in Figure 10, the wide distribution of ELPWs in Figures 11-12 results from the rocket body slope parallel to the beam direction ( || S ). However, the ELPW has the distinctive shape subjected to the orientation of the rocket body, which enables orientation information to be achieved. Regarding the Gaussian decomposition in Figure 13, the first two Gaussian functions represent five rocket engines while the next three Gaussian functions correspond to the thrust structure and the sixth small Gaussian function comes from all five rocket engines. The Gaussian decomposition results including the ELPWs are also shown in Figures 10-13 for four kinds of orientations of the rocket body. Unlike the ELPWs of the defunct satellite, the size of the rocket body is so large that the noise power is very weak compared to the received power. In addition, the reflected surface is too complex compared to the defunct satellite to fit the ELPWs well even by using several Gaussian functions. Therefore, the total number of Gaussian functions (n) is selected as the smallest number that satisfies the criterion of Equation (24). As shown in Figure 10, the first two Gaussian functions corresponding to the first peak represent the forward and aft skirts as well as the tank structure because their shapes are cylindrical, and the other Gaussian functions correspond to five rocket engines, the thrust structure and the aft interstage. The ELPWs in Figures 11 and 12 are so complex that each Gaussian function cannot be exactly in one-to-one correspondence with the components of the rocket body. Compared to the narrow distribution of ELPW in Figure 10, the wide distribution of ELPWs in Figures 11 and 12 results from the rocket body slope parallel to the beam direction (S || ). However, the ELPW has the distinctive shape subjected to the orientation of the rocket body, which enables orientation information to be achieved. Regarding the Gaussian decomposition in Figure 13, the first two Gaussian functions represent five rocket engines while the next three Gaussian functions correspond to the thrust structure and the sixth small Gaussian function comes from all five rocket engines.
It was shown that the orientation of a defunct satellite and a rocket body have a distinctive shape in the ELPW and the Gaussian functions that result from the Gaussian decomposition, which are a tool to characterize the orientation in more detail. The ELPW feature resulted from a single laser pulse that can provide sufficient information for the orientation determination of space debris. However, there is also a possibility in practical application that the ELPW features do not characterize the orientation of the space target due to its simple shape. If the ELPWs are available sequentially while the DLR is activated, the ELPW changes not only in shape but also in amplitude, which can provide more useful information for the orientation determination of space targets including even tumbling space debris. A number of small satellites, especially CubeSats which have simple and cubic shapes, have now been launched into space for scientific and commercial missions. It is difficult to determine the orientation of defunct CubeSats from the ELPW analysis even though a higher power laser is employed since the ELPW does not indicate the distinct features depending on their orientations because of a simple and cubic shape. But the orientation information of CubeSats is not relatively important in terms of the orbital prediction because the cross-sectional area variation is negligible over the entire orbit in computing the atmospheric drag and solar radiation pressure. The Gaussian decomposition results including the ELPWs are also shown in Figures 10-13 for four kinds of orientations of the rocket body. Unlike the ELPWs of the defunct satellite, the size of the rocket body is so large that the noise power is very weak compared to the received power. In addition, the reflected surface is too complex compared to the defunct satellite to fit the ELPWs well even by using several Gaussian functions. Therefore, the total number of Gaussian functions ( n ) is selected as the smallest number that satisfies the criterion of Equation (24). As shown in Figure 10, the first two Gaussian functions corresponding to the first peak represent the forward and aft skirts as well as the tank structure because their shapes are cylindrical, and the other Gaussian functions correspond to five rocket engines, the thrust structure and the aft interstage. The ELPWs in Figures  11-12 are so complex that each Gaussian function cannot be exactly in one-to-one correspondence with the components of the rocket body. Compared to the narrow distribution of ELPW in Figure 10, the wide distribution of ELPWs in Figures 11-12 results from the rocket body slope parallel to the beam direction ( || S ). However, the ELPW has the distinctive shape subjected to the orientation of the rocket body, which enables orientation information to be achieved. Regarding the Gaussian decomposition in Figure 13, the first two Gaussian functions represent five rocket engines while the next three Gaussian functions correspond to the thrust structure and the sixth small Gaussian function comes from all five rocket engines.

Conclusions
Orientation information of space debris is important to improve the orbital prediction accuracy for mitigation or elimination of a significant threat to not only human space activities but also operational satellites. The ELPW records the scattering signals from the space target surface and thus reflects the target shape in the direction of laser beam propagation. The new method based on the ELPW analysis was proposed to achieve the orientation information of space debris because it has no measurement constraints, unlike the photometry, AO and SLR technologies. This paper discussed the feasibility of using ELPW for the orientation determination of space debris by numerical simulations for a defunct satellite and a rocket body. By assuming that the reflected surface consists of small grids, the mathematical model was derived to generate the ELPW and also took into account the temporal beam broadening and spatial spreading due to the atmospheric turbulence. The beam broadening is so significant only in the femtosecond regime; it is negligible in the numerical simulation using the nanosecond pulse laser. Gaussian decomposition based on a genetic algorithm was employed to the ELPWs in order to analyze the orientation features. In particular, Gaussian functions obtained from the Gaussian decomposition were able to characterize the orientation features in more detail for a defunct satellite with a simpler configuration than a rocket body. It was demonstrated from the numerical simulations that the ELPW had a distinctive shape corresponding to the orientation of space debris, and therefore it could be considered as a feasible method to provide orientation information. Future works will focus on machine learning as a powerful tool to identify the orientation based on waveform shape recognition.