Single-frequency LADAR super-resolution Doppler tomography for extended targets

Di Mo; Ning Wang; Ran Wang; Zi-Qi Song; Guang-Zuo Li; Yi-Rong Wu

doi:10.1364/OE.27.012923

1. Introduction

The nature of LADAR imaging is to invert the distribution of scattering coefficient by acquiring the spatial spectrum of the target, which is similar to the basic principle of microwave radar imaging. The acquisition of spatial spectrum mainly includes two aspects: the spatial diversity of observations [1–4], and the frequency diversity of signals [5–8]. Doppler tomography relies mainly on the spatial diversity of observations to image remote targets, which gives Doppler tomography LADAR some natural advantages. First, compared to LADAR systems that rely on frequency diversity, such as range-resolution tomography, inverse synthetic aperture LADAR (ISAL), and frequency modulated continuous wave (FMCW) LADAR, Doppler tomography LADAR avoids difficult broadband modulation of the laser over a short pulse duration [9–12]. In contrast, the Doppler tomography system is simpler and has a higher signal-to-noise ratio (SNR), making it more suitable for resource-constrained spaceborne and airborne applications [3,13,14]. Second, Doppler tomography LADAR can obtain images of spinning targets such as space debris, ballistic missile warheads, and helicopter blades [15–20]. The spin of these targets can result in a Doppler bandwidth on the order of MHz. The pulse repetition frequency (PRF) of the conventional radar system generally cannot reach this magnitude, which leads to the folding of images and the failure of target recognition. However, for Doppler tomography LADAR, only the sampling bandwidth is required to be larger than the Doppler bandwidth, which is very easy to implement. Finally, Doppler tomography is also easy to extend from two-dimensional (2-D) imaging to three-dimensional (3-D) imaging [20–22].

The concept of Doppler tomography was first proposed in the field of microwave radar [23–25]. Its imaging resolution depends on the resolution of the Doppler spectrum. Lasers and terahertz have a natural advantage in Doppler imaging compared to microwave [10,11,20,26–28]. Therefore, the research on radar systems and image reconstruction algorithms for high frequency band Doppler tomography is very promising. However, no high-quality Doppler tomograms and related imaging methods have been reported. Most research work is limited to imaging experiments and theoretical analysis for point targets. The effect of the projection model on image quality is decisive. According to the projection model, the existing Doppler tomography technologies are divided into two categories: based on complex spectra [2,13,16,21,24], and based on spectral amplitude [3,15,18–20,29,30]. The first type directly considers the Doppler spectrum as a one-dimensional (1-D) projection of the target. These methods perform coherent processing on full-angle echoes to achieve high imaging resolution. However, the side lobes of the point spread function (PSF) reach −8 dB and cannot be effectively suppressed by the window function. Due to such a serious energy leak, the dynamic range of the radar is unacceptable, and the images of extended targets have poor quality. Moreover, since the motion error of targets, and phases of scattering centers vary with the observation angle, it is practically difficult to achieve coherent processing on all angles. The other type Doppler tomography performs incoherent processing between angles. These methods can achieve good focus for point targets. However, there are more artifacts and noise in the images of extended targets, and the details inside the targets are lost. The cause of these issues is to directly consider the intensity of the Doppler spectra as 1-D projections for tomographic processing. This inaccurate processing introduces spectrum speckle into images, which has far more impact than system noise [31]. The speckle inherent in coherent imaging becomes a major factor affecting image quality, which limits the improvement in image quality.

In summary, theoretical studies and experimental validation of Doppler tomography have focused on the microwave band rather than the more promising laser and terahertz bands. Existing imaging methods do not consider the effect of speckle on the Doppler spectra and can only be applied to point targets rather than extended targets that are closer to the actual situation. These states have severely limited the practical application of Doppler tomography. To solve these problems, we propose a novel laser Doppler tomography method. The method is based on single-frequency LADAR without any form of wideband modulation of the transmitted signal. The imaging process is based on the precise relationship between the target scattering coefficient and the statistical characteristics of the Doppler spectrum. The scattering coefficient distribution of the target is inverted from the statistical characteristics of the Doppler spectrum by the maximum a posteriori estimation. This processing avoids introducing the extra artifacts and noise caused by spectral speckle into the image, which greatly improves quality of images. According to theoretical analysis, the resolution of proposed imaging method depends on the Doppler frequency resolution, which exceeds the diffraction limit of the optical aperture and is independent of the imaging distance. According to the proposed imaging method, a set of laser Doppler tomography experimental system was built. With this system, a series of 2-D images of the target with the different rotational speeds are acquired at a distance of 5 m indoors. The resolution of the fine images is 0.4 mm, which is a 200-fold increase in resolution compared to the diffraction limit. The shapes of the targets are clearly presented in the images. Even the detailed features of the targets can be distinguished. This result verifies the effectiveness of the proposed method.

2. Echo model

Doppler tomography reconstructs images by acquiring multi-angle echoes of targets. Its Imaging geometry can be simply represented as Fig. 1.

Fig. 1 Imaging geometry of laser Doppler tomography. LOS, radar line of sight.

Download Full Size | PDF

In the figure, the target spins around the z-axis and the angular velocity is ω. The radar line of sight (LOS) is perpendicular to the x-axis and the angle to the z-axis is θ. The backscatter coefficient of the target surface is defined as $σ (x, y)$ . The surface of targets is generally considered rough relative to the laser. The small undulation and non-uniformity on rough surface introduce random phase in the echo wavefront, defined as $φ (x, y)$ , which causes speckle. In order to obtain Doppler spectra at different angles, the echo signal is truncated into N sub-segments of duration $T_{p}$ . In the nth sub-segment, the internal time is defined as $t_{s} \in (- T_{p} / 2, T_{p} / 2)$ , and its relationship with real time t can be expressed as $t = t_{s} + n T_{p}$ . Due to the spin of the target, the coordinates at $t_{s} = 0$ in the nth sub-segment can be expressed as the product of the rotation matrix and the initial coordinates,

[\begin{matrix} x_{n} \\ y_{n} \end{matrix}] = [\begin{matrix} \cos n T_{p} ω & - \sin n T_{p} ω \\ \sin n T_{p} ω & \cos n T_{p} ω \end{matrix}] \cdot [\begin{matrix} x_{0} \\ y_{0} \end{matrix}] .

Then, in the nth sub-segment, the heterodyne signal of the target can be expressed as

s_{n} (t_{s}) = \iint σ (x_{n}, y_{n}) \exp [j φ (x_{n}, y_{n})] \exp [j \frac{2 π}{λ} (x_{n} \sin ω t_{s} + y_{n} \cos ω t_{s}) \sin θ] d x_{n} d y_{n},

where λ represents the laser wavelength. When

T_{p}

is small, the Doppler spectrum of the nth sub-segment can be approximated as

\begin{array}{l} S_{n} (f) = F [s_{n} (t_{s})] \\ \approx ∭ d x_{n} d y_{n} d t_{s} σ (x_{n}, y_{n}) \exp [j φ (x_{n}, y_{n})] \\ \times \exp [j \frac{2 π}{λ} (x_{n} ω t_{s} + y_{n}) \sin θ] \exp (- j 2 π f t_{s}) \\ = \int σ (\frac{λ}{ω \sin θ} f, y_{n}) \exp {j [φ (\frac{λ}{ω \sin θ} f, y_{n}) + \frac{2 π}{λ} y_{n} \sin θ]} d y_{n}, \end{array}

where

F [\cdot]

represents the Fourier transform. Similar to the 1-D projection in conventional tomography, the Doppler spectrum is also the integral of the target scattering coefficient along the y_n-axis. However, the difference is that this integral process is modulation by random phase containing

φ (x, y)

. This difference makes it impossible to obtain clear images by performing conventional tomographic processing directly on the Doppler spectra. For targets consisting of a few isolated points, the integration along the y_n-axis contains at most one non-zero value in most sub-segments. In this case, the effect of random phase can be removed by modulating the spectrum

S_{n} (f)

. Those tomography techniques based on spectral amplitude are to perform such processing to obtain clear images. Unfortunately, in reality, targets usually consist of continuous scattering surfaces. Thus, the amplitude of the Doppler spectrum cannot be directly considered as a 1-D projection in tomography.

However, it is also due to the presence of speckle (i.e., random phase $φ (x, y)$ ) that makes it possible to recover the 1-D projection required for tomography from the Doppler spectrum [32]. In Fig. 2, $F [σ (x, y)]$ represents the spatial spectrum of the scattering coefficient distribution. According to the projection-slice theorem, conventional tomography, such as X-ray transmission tomography, acquires slices of the spatial spectrum at different angles, which are Fourier transforms of 1-D projections. For laser Doppler tomography, the sampled spectrum is actually $F [σ (x, y)] \otimes F [φ (x, y)]$ due to the presence of the random phase $φ (x, y)$ , where $\otimes$ represents the convolution. Since $φ (x, y)$ is completely random, this convolution causes the target spatial spectrum to be smeared into the entire frequency domain. Therefore, even if the laser Doppler tomography detects only the annular high-frequency region in $F [σ (x, y)] \otimes F [φ (x, y)]$ , the information of the low-frequency region of $F [σ (x, y)]$ required for imaging is included. The essence of the proposed imaging process is to restore the circular target spatial spectrum through the annular original echo spectrum, as shown in Fig. 2. The advantage of this approach compared to the direct use of the original echo spectrum is that it avoids the −8 dB side lobes of the annular spectrum.

Fig. 2 Comparison of laser Doppler tomographic and conventional transmission tomographic spectra. Conventional transmission tomography samples the target spatial spectrum. The spatial spectrum sampled by laser Doppler tomography is convolved by the spectrum of the random phase.

Download Full Size | PDF

3. Imaging processing

3.1 Fast processing

The Doppler spectrum is most similar to the 1-D projection of the scattering coefficient at the same angle. Therefore, the most direct approach is to use the Doppler spectrum to recover the scattering coefficient projection at the same angle, and then use the conventional tomographic algorithm to achieve 2-D imaging. According to Eq. (3), the Doppler spectrum is actually an integral of the complex scattering coefficient composed of the scattering coefficient and the random phase. Since the target and observation conditions cannot be preset, the statistical distribution of the complex scattering coefficient cannot be modeled. However, for the extended target, the integral in Eq. (3) contains a large number of random values. According to the central limit theorem, the integral result should conform to the complex Gaussian distribution, i.e.,

S_{n} (f) \sim C N (0, \int σ {(\frac{λ}{ω \sin θ} f, y_{n})}^{2} d y_{n}) .

The integral of the square of the scattering coefficient along the y_n-axis is equal to the variance of the corresponding frequency point in the spectrum. In other words, scattering coefficients do not directly determine the amplitude of Doppler spectra, but rather variances of frequency points. Therefore, in order to obtain a 1-D projection suitable for conventional tomography algorithms, it is necessary to estimate the variance at each frequency point of the Doppler spectrum.

A convenient method is to square the amplitude of the Doppler spectrum to obtain a power spectrum which obeys a non-central chi-square distribution with a degree of freedom of 2. Then, the expectation of the non-central chi-square distribution can be expressed as

E ({| S_{n} (f) |}^{2}) = 2 \int σ {(\frac{λ}{ω \sin θ} f, y_{n})}^{2} d y_{n} .

This mean of the Doppler power spectrum of adjacent sub-segments can be considered as an estimate of the 1-D projection of the square of the scattering coefficient. Then, a 2-D distribution of the square of the scattering coefficient can be obtained by performing conventional tomographic processing on these projections. The image is the square root of this distribution. The two points that need to be explained are: (1) The 2-D distribution of the square of the scattering coefficient obtained by tomographic reconstruction can also be regarded as the distribution of the echo energy. Negative energy values are meaningless. However, the conventional tomographic processing cannot guarantee the non-negative result, complex values may appear in the square root. In the processing, the real part is taken as an image, which is equivalent to setting the negative energy values to zero before performing the square root. (2) The operation in the Eq. (5) does not need to be performed intentionally. This operation is implicit in tomographic algorithms. Taking the filtered back-projection algorithm as an example, when integrating along a sinus line, the positions of the integrated elements in the power spectrum of adjacent sub-segments are basically the same. Therefore, the sinusoidal integral itself performs the operation of averaging the power spectrum of adjacent sub-segments. The same conclusion can be obtained for the frequency domain algorithm.

The problem with fast processing is that the resolution and SNR of the image cannot be improved at the same time. The noise in the image is mainly derived from speckle in the Doppler spectra. To suppress speckle, a sufficient number of sub-segments are required to average the power spectrum. However, when the difference in viewing angle between these sub-segments is too large, the same scattering point will fall in different frequency resolution units in the power spectra. Therefore, the integral path of Eq. (3) will change beyond the resolution unit in the sub-segments with long intervals. The power spectra of these sub-segments cannot be used for averaging. The echo signal needs to be truncated into short sub-segments to ensure that there are enough spectra with similar viewing angles, which results in low image resolution.

3.2 Fine processing

In order to obtain images with high resolution and high quality, a reconstruction method based on maximum a posterior probability (MAP) estimation is proposed. Unlike fast processing, this fine processing does not attempt to obtain 1-D projections through the Doppler spectra, but directly inverts the distribution of the scattering coefficient, i.e., the image.

First, the imaging plane and the Doppler spectrum are discretized. The imaging plane $σ (x, y)$ is divided into K pixel units. These pixel units constitute a scattering coefficient vector $σ = {(σ_{i})}_{K \times 1}$ . The modulus values of the Doppler spectra of all angles are grouped into a vector, defined as $S = {(S_{j})}_{M N \times 1}$ , where M represents the number of resolution units in each spectrum.

Fine processing is based on a Bayesian framework. The magnitudes of the Doppler spectra are considered to be observation samples. The scattering coefficients $σ$ are considered as the parameters to be estimated. According to Bayes' Theorem, the posterior probability can be expressed as

P (σ | S) = \frac{P (S | σ) P (σ)}{P (S)} .

The purpose of the processing is to find the appropriate parameters

σ

so that the posterior probability is taken to the maximum. The estimates of

σ

are given by

\begin{array}{l} σ = \underset{σ \in R^{+}}{\arg \min} {- \ln P (σ | S)} \\ = \underset{σ \in R^{+}}{\arg \min} {- \ln P (S | σ) - \ln P (σ)} . \end{array}

Since the observation sample

S

is known,

P (S)

is a fixed value, so omitting it in the Eq. (7) does not affect the estimation of the parameter

σ

. The simplified expression contains two items, the likelihood function

\log P (S | σ)

and the prior probability

\log P (σ)

. Maximizing the likelihood function makes the estimate

σ

as close as possible to the observed sample

S

. The prior probability term representing the probability of the estimate itself can prevent the estimate

σ

from being over-fitting.

The likelihood function is obtained by the distribution of $S$ . According to the previous analysis, the jth element $S_{j}$ in $S$ obeys the Rayleigh distribution (its square is subject to the chi-square distribution). As shown in Eq. (4), the variance of the distribution is equal to the line integral of the square of the scattering coefficient. The integral can be approximated by a weighted summation $σ^{T} H_{j} σ$ in discrete form, where the weight matrix $H_{j}$ is a diagonal matrix, i.e.,

H_{j} = d i a g [{(h_{i, j})}_{K \times 1}] = [\begin{matrix} h_{1, j} & 0 \\ ⋱ \\ 0 & h_{K, j} \end{matrix}],

where

h_{i, j}

is the weighting factor of each resolution unit, which is zero outside the integration path. Then, the probability density function of

S_{j}

can be expressed as

f (S_{j}) = \frac{S_{j}}{σ^{T} H_{j} σ} \exp (- \frac{S_{j}^{2}}{2 σ^{T} H_{j} σ}), S_{j} > 0.

According to the probability density function, the likelihood function can be obtained as

\ln P (S | σ) = - \frac{1}{2} \sum_{j} \frac{S_{j}^{2}}{σ^{T} H_{j} σ} - \sum_{j} \ln (σ^{T} H_{j} σ) + \sum_{j} \ln S_{j},

where

S

is the determined value, so the value of the last term is fixed and can be ignored in the actual calculation without affecting the estimate of

σ

.

The a priori model of the image is approximated by the Markov Random Field (MRF) in which pixel values are only related to adjacent pixels, regardless of other pixels that are not adjacent. When only pairwise cliques are considered, the probability density function of the image can be expressed as a Gibbs distribution [33,34], i.e.,

P (σ) = \frac{1}{Z} \exp {- \sum_{{i, j} \in C} b_{i, j} ρ (\frac{σ_{i} - σ_{j}}{a})},

where

C

denotes all cliques in the image, Z is the normalization coefficient which is used to ensure that the integral of the distribution is equal to 1,

b_{i, j}

is the weight coefficient of different cliques, and a is the scaling factor which is used to adjust the intensity of regularization. For the q-Generalized Gaussian MRF (q-GGMRF), its potential function can be expressed as [34]

ρ (Δ) = \frac{{| Δ |}^{q}}{q} (\frac{{| Δ / c |}^{p - q}}{1 + {| Δ / c |}^{p - q}}), (1 \leq q \leq p \leq 2),

where c is the threshold parameter. When

Δ ≫ c

, the potential function approximates a q-norm, and when

Δ ≪ c

, the potential function approximates a p-norm. In this way, the image edges can be well preserved while smoothing the noise with small amplitudes. The probability density function of the q-GGMRF is complicated, and its partial derivatives are difficult to obtain when optimizing iterations. To simplify the calculation, an upper bound function

ρ_{s} (Δ) = K Δ^{2}

is used instead of the potential function in Eq. (12), where K is a fixed coefficient in each iteration, which is calculated by the difference

Δ_{0}

of the pixels in the current image, i.e.,

K = \frac{{| Δ_{0} |}^{p - 2}}{q c^{p - q} (1 + {| Δ_{0} / c |}^{p - q})} (p - \frac{(p - q) {| Δ_{0} / c |}^{p - q}}{1 + {| Δ_{0} / c |}^{p - q}}) .

Then, the simplified prior probability can be expressed as

\ln P (σ) = - \frac{1}{a^{2}} \sum_{{i, j} \in C} b_{i, j} K_{i, j} {(σ_{i} - σ_{j})}^{2} - \ln Z,

where the value of the last item is fixed and can be ignored in the actual calculation without affecting the estimation of

σ

.

By introducing Eqs. (10) and (14) into Eq. (7), the final optimization objective function can be obtained, which is a constrained non-convex optimization. Since the scattering coefficients are non-negative, the optimization of $σ$ has a boundary constraint. In order to remove this constraint, variable substitution $σ = \exp (η)$ can be performed, and the new optimization parameter $η$ can take any value. Moreover, although the new optimization function is non-convex (mainly in the region where the parameter is much smaller than the optimal solution), it can be approximated by the convex function in the vicinity of the optimal solution. When an appropriate initial value is given, conventional optimization algorithms, such as quasi-Newton method and expectation-maximization algorithm (EM), can solve this optimization problem quickly and stably [35,36]. An image with low resolution and high SNR can be used as the initial value, which can be obtained by fast processing. It should be noted that since the line integral only contains a small part of the entire imaging plane, most of the diagonal elements of the weight matrix $H_{j}$ are also zero. Therefore, although the dimension of $H_{j}$ is large, it is very sparse, and does not require excessive memory in the calculation.

The fast processing indirectly acquires images by inverting 1-D projections, and can only utilize information within small angles. In contrast, the fine processing directly inverts the image, and can use the information contained in the Doppler spectra of all angles. Therefore, the image obtained by fine processing has both high resolution and high SNR at the same time. The complete imaging process uses a combination of these two processing methods, the flow of which is shown in Fig. 3. With different sub-segment durations $T_{p}$ , two spectra with high angular resolution low frequency resolution, and low angular resolution high frequency resolution can be obtained. These two spectra are used for fast processing and fine processing, respectively. The low-resolution, high-SNR image obtained by the fast processing can be considered as a preliminary result of the imaging process, and it is interpolated as the initial value of the fine processing. Finally, a high-resolution, high-SNR image is obtained by fine processing.

Fig. 3 Flow chart of the complete imaging process.

Download Full Size | PDF

3.3 Image resolution

The resolution of the tomographic image depends on the resolution of the projection, i.e., the Doppler frequency resolution, which is inversely proportional to the duration $T_{p}$ of the sub-segment. According to the imaging geometry, when a rectangular window is used to truncate sub-segments, the spatial resolution $Δ r_{p}$ of the 1-D projection can be expressed as

Δ r_{p} = \frac{λ}{2 ω T_{p} \sin θ} .

The spatial spectral bandwidth of the 1-D projection is

B_{s} = 1 / Δ r_{p}

, which is equal to the diameter of the circular spatial spectrum in Fig. 2. The point spread function (PSF) can be obtained by 2-D Fourier transform and normalization of the circular spatial spectrum, which can be expressed as

P S F (x, y) = \frac{J_{1} (\frac{2 π}{λ} ω T_{p} \sqrt{x^{2} + y^{2}} \sin θ)}{\frac{π}{λ} ω T_{p} \sqrt{x^{2} + y^{2}} \sin θ},

where

J_{1} (\cdot)

represents the Bessel functions of the first kind of one order. The PSF is rotationally symmetric, and the peak sidelobe ratio is about 17.6 dB. The schematic diagram of the PSF is shown in Fig. 4. According to Eq. (16), the −3 dB width of main lobe of the PSF is approximately

1.03 Δ r_{p}

. This resolution is not limited by the optical aperture and is independent of the imaging distance.

Fig. 4 The normalized PSF of the proposed imaging process. (a) 3-D profile of the PSF, z-axis represents the normalized amplitude and x- and y-axes are normalized by $Δ r_{p}$ . (b) 2-D profile of the PSF at $y = 0$ , the amplitude of the PSF is expressed in decibels (dB).

Download Full Size | PDF

4. Results and discussion

4.1 Experiment system

To verify the effectiveness of the proposed imaging method, a single-frequency LADAR experimental system was built. The LADAR system with separated transmitter and receivers is shown in Fig. 5. The diffraction limit of the system is defined as the main lobe width of the antenna pattern which also determines the size of the imaging area. In order to obtain an suitable imaging area size, the transmit and receive lenses are used to extend −3 dB main lobe width to 16 mrad which corresponds to a resolution of 8 cm at 5 m away from the LADAR. The targets are flat patterns made of reflective material, which are driven by a motor. In the LADAR system, the output of a single-frequency laser is split into two paths. One of them is amplified by an erbium doped fiber amplifier (EDFA) and directly used as the transmitting signal. The other is the local oscillator, which is heterodyned with the echo. Quadrature detection is implemented by a 90° optical hybrid. Two balanced detectors are used to detect the optical signals of in-phase and quadrature (I and Q) channels, respectively. After sampling the electrical signal output by the detector, imaging processing is performed in the digital domain.The main parameters of the experimental system are shown in Table 1.

Fig. 5 Experimental system for demonstration of laser Doppler tomography. BS, beam splitters; EDFA, erbium doped fiber amplifier; BD, balanced detector; AD, analog-to-digital converter.

Download Full Size | PDF

Table 1. Main parameters of the experimental system

View Table

The proposed imaging method is based on the assumption that the amplitude of the echo Doppler spectrum obeys the Rayleigh distribution. This assumption is verified in an experiment. A uniform scattering plane with a diameter larger than the spot was used for the experiment. The scattering coefficients of the different parts of the plane are approximately the same. Since the spot is approximately circular, the frequency points of the same position in the Doppler spectrum of different sub-segments correspond to the integration of the same length in the circle. Therefore, it can be considered that the amplitudes of the same frequency points have the same distribution, which is convenient for obtaining data for statistical analysis. The effect of speckle on the amplitude of the Doppler spectrum is shown graphically in Fig. 6(a). The amplitude of this spectrum contains a large amount of random jitter compared to the projection sinogram in conventional tomography. These jitters are subject to the Rayleigh distribution, as shown in Fig. 6(b). The black curve in the figure is the probability density of the standard Rayleigh distribution, and the distribution of the normalized amplitude represented by the gray bars is in good agreement with the standard Rayleigh distribution.

Fig. 6 The distribution of the Doppler spectrum amplitude. (a) Speckle makes the Doppler spectrum not smooth. The vertical axis is the Doppler frequency, and the horizontal axis is the time corresponding to sub-segments. (b) At a fixed frequency (500 kHz), the Doppler spectrum amplitude obeys the Rayleigh distribution. The amplitude is normalized to its standard deviation.

Download Full Size | PDF

4.2 Comparison of fast processing with existing methods

The three differently shaped targets used for the imaging experiment are shown in Figs. 7(a)–7(c). The angular velocity of the square target shown in Fig. 7(a) is about 13.09 rad/s, and the angular velocity of targets in Figs. 7(b) and 7(c) is about 98.48 rad/s. It should be noted that the angular velocity does not need to be known in advance. It can be estimated by the periodic variation of the sub-segment Doppler spectrum. As with conventional transmission tomography, Doppler tomography enables clear imaging of the target when the target is rotated more than 180°. Therefore, for two categories of targets with different angular velocities, the time required for imaging is 250 ms and 32 ms, respectively. The echo signal needs to be intercepted as a series of sub-segments before the imaging process. For two categories of targets, the length of the sub-segments is 113 μs and 15 μs, respectively. According to Eq. (16), the resolution in both cases is 1.58 mm. The purpose of this setup is to facilitate comparison of the results of the imaging.

Fig. 7 (a)–(c) Targets for the imaging experiment. (d)–(f) Images obtained by the fast processing. (g)–(i) Images obtained from the amplitudes of spectra. (j)–(l) Images obtained from the complex spectra. For the same target, different processing methods use the same data.

Download Full Size | PDF

The results of the fast processing are shown in Figs. 7(d)–7(f). The imaging quality is the same for the slow-rotating target and the fast-rotating target. Images acquired by fast processing have few artifacts. In these images, it is easy to identify the outline of the target and distinguish the targets of different shapes. However, some existing processing methods directly perform the conventional tomographic algorithm for Doppler tomography. The Doppler spectrum amplitude is equivalent to the 1-D projection in these methods, which causes image quality degradation, as shown in the Figs. 7(g)–7(i). There are two main problems in these images. First, the amplitude of the image does not accurately reflect the intensity of the echo energy. The amplitude of the central area is reduced and the amplitude of the peripheral area is enlarged. For this set of experiments, since the energy of the wavefront is Gaussian, the amplitude of the center of the image should be larger than the periphery, but it is not reflected in Figs. 7(g)–7(i). The second problem is a lot of artifacts. These artifacts are mainly distributed in the convex envelope of the target and overlap with the image, which makes the internal edges difficult to distinguish. The images obtained by the full-angle coherent processing which is the other type of existing methods are shown in Figs. 7(j)–7(l). Obviously, the target in the image is illegible because of the −8 dB side lobes which cannot be avoided by this type of processing. The image is only relatively bright in the area where the target appears. The experimental result shows that this type of all-coherent processing cannot image extended targets.

4.3 Noise and speckle

In the experiment, the effect of noise on fast processing is analyzed. The noise mentioned here is white noise obeying Gaussian distribution such as shot noise, system thermal noise and channel noise. The result of the experiment is shown in Fig. 8. The horizontal axis is the projection SNR, i.e., the SNR of the sub-segment Doppler spectrum. The vertical axis is the image SNR, which is defined as the ratio of the average amplitude of target pixels to the average amplitude of background pixels. Obviously, when the projection SNR is attenuated from 40 dB to 10 dB, the image SNR remains essentially constant. In this interval, the dominant factor limiting the image SNR is not the SNR of the projection, but the speckle in the spectrum. When the projection SNR continues to decay, the SNR of the image begins to decrease, and the projection SNR becomes the dominant factor. The SNR of the echo is less than 10 dB for most applications of the Doppler tomography. However, since the system is narrowband, the SNR of the Doppler spectrum can usually reach 10 dB or higher. Therefore, the dominant factor affecting the image SNR in practical applications is speckle, which is consistent with the theoretical analysis in [31].

Fig. 8 The effect of noise on the image SNR. The results of the sub-segment durations of 8 μs, 15 μs, and 30 μs are obtained respectively.

Download Full Size | PDF

Speckle is an inherent phenomenon in coherent imaging. For Doppler tomography, the speckle makes the amplitude of the spectrum not equivalent to the projection of the target scattering coefficient. According to the previous analysis, the integral of the scattering coefficient is equal to the variance of the frequency point. The fast processing estimates this variance by the spectrum of multiple adjacent sub-segments. Reducing the sub-segment time length $T_{p}$ means increasing the number of spectra used for estimation, so that the estimate of the variance is more accurate. As shown in Fig. 8, under the condition that the projection SNR is the same, the decrease of $T_{p}$ leads to an improvement of the image SNR. Figure 9 shows this issue more intuitively. The $T_{p}$ corresponding to Figs. 9(a)–9(c) are 15 μs, 30 μs and 60 μs, respectively. Compared to Fig. 9(a), Figs. 9(b) and 9(c) have a higher resolution, but the quality of the image is significantly degraded. The three curves in Fig. 9(d) are the normalized power at sections s1 to 3. Compared to the other two curves, s1 has a higher SNR and the edges of the target are sharper. In the image obtained by fast processing, resolution and quality are contradictory, which is solved in the fine processing.

Fig. 9 (a)–(c) Images of the sub-segment durations of 8 μs, 15 μs, and 30 μs obtained by fast processing. (d) Comparison of cross sections s1, s2 and s3. For visualization, s1 and s2 are translated by 0.5 and 1 along the vertical axis, respectively.

Download Full Size | PDF

4.4 Results of fine processing

In our experiments, this optimization problem defined by Eq. (7) is solved by the standard Limited-memory BFGS (L-BFGS) algorithm [36]. The reconstructions are run in MATLAB using a computer with a 2.5 GHz Intel Core i7-6500U processor. As shown in Fig. 10, the computation time of the iteration is approximately linear with the size of the reconstruction. The relative change in the reconstructed image is used as the termination condition for the iteration, and the number of iterations is usually less than 100 when the threshold is set to 10⁻⁴. In comparison, fast processing does not require iteration, and its computation time is similar to conventional algorithms (such as the filtered back-projection and the polar fourier transform). Therefore, the complete processing time depends on the time of the fine processing.

Fig. 10 Time of iterative calculation under different reconstruction sizes.

Download Full Size | PDF

Through fine processing, high-resolution, high-quality images of targets can be obtained. The fine processing attempts to find the MAP estimate of the scattering coefficient distribution. Appropriate initial values can greatly speed up the iteration. This initial value can be obtained by a low-resolution, high-SNR image, which is the result of the fast processing. For slow-rotating targets with an angular velocity of 13.09 rad/s, data with a sub-segment length of 60 μs is used in fast processing. For fast-rotating targets with an angular velocity of 98.48 rad/s, the sub-segment is 8 μs. Under these conditions, images obtained by the fast processing have high SNRs of more than 23 dB. The initial values are obtained by interpolating the images onto the grid required for fine processing. It should be noted that the Doppler bandwidth of the target in the experiment only accounts for a small portion of the sampling bandwidth. In order to reduce unnecessary calculations, only the part of the spectrum containing the target is intercepted for fast processing. This interception range needs to be consistent with the fine processing. Otherwise, the image obtained by the fast processing needs to be enlarged and cropped (or reduced and zero-padded) to be used as the initial value.

The comparison of the reconstructed results of the fine processing and the conventional tomography methods is shown as Fig. 11. In the fine processing and the conventional tomography, the durations of sub-segments are 450 μs for the slow-rotating target and 60 μs for the fast-rotating target. According to Eq. (16), the theoretical resolution of all images is 0.4 mm, which is a 200-fold improvement in resolution compared to diffraction limit of 8 cm. Compared to the conventional method, the quality of these images obtained by the fine processing is satisfactory. The outline and shape of the target are clearly expressed. Even the small diamond texture with a line width of about 0.5 mm on the target can be rendered. Due to the poor SNR, this texture cannot be observed in the images obtained by the conventional method, even if the theoretical resolution is also 0.4 mm. The quality of images obtained by fine processing is greatly improved compared to fast processing or other existing imaging processing.

Fig. 11 The comparison of the reconstructed results of the fine processing and the conventional tomography methods. Images of slow-rotating target obtained by (a) the fine processing and (b) the conventional tomography methods. Images of fast-rotating target obtained by (c) the fine processing and (d) the conventional tomography methods.

Download Full Size | PDF

5. Conclusions and outlook

The proposed Doppler tomography processing method consists of two main processing flows, fast processing and fine processing. The image obtained by fast processing can be used as a preliminary imaging result. However, the resolution and quality of the image cannot be improved at the same time. By using this preliminary result as an initial value, a high-quality, high-resolution image can be obtained by fine processing. A series of laboratory experiments were carried out using the single-frequency LADAR system. Images of the slow-rotating target and the fast-rotating target are obtained. The resolution of the images is 0.4 mm, which is a 200-fold improvement in resolution compared to diffraction limit and is independent of the imaging distance. In addition, the optical system required by the proposed imaging method is consistent with a coherent LADAR, and the two complement each other: LADAR can be used for distance measurement, and the proposed method can acquire images. Therefore, this proposed method is more suitable for remote targets than conventional optical imaging methods such as cameras. The experimental results verify the effectiveness of the proposed method. The effects of noise and speckle on the image are analyzed. The inference that speckle is the dominant factor is verified. Therefore, it is precisely because the proposed method is based on the exact relationship between the statistical characteristics of the echo Doppler spectrum and the scattering coefficient that a high-quality image of the extended target can be obtained. In the experiment, clear images of targets with different rotational speeds were obtained, which indicates that the reconstruction performance of the proposed method is not affected by the target rotational speed. However, as the rotational speed increases, the duration of the sub-segment decreases accordingly, which leads to a decrease in the SNR of the Doppler spectrum. When Gaussian noise exceeds speckle as the dominant factor, the quality of the image is affected. This needs to be considered in the system design. In summary, the proposed method makes it possible for a single-frequency LADAR or a conventional coherent LADAR to acquire high-resolution, high-quality 2-D images of remote targets, including slow-rotating targets (such as space debris and warheads, etc.) and fast-rotating targets (such as blades, etc.), which is of great significance for the detection and recognition of these targets.

Funding

National Natural Science Foundation of China (NSFC) (61575198); Natural Science Foundation of Beijing Municipality (4182074).

References

1. A. Bulbul, A. Vijayakumar, and J. Rosen, “Superresolution far-field imaging by coded phase reflectors distributed only along the boundary of synthetic apertures,” Optica 5(12), 1607–1616 (2018). [CrossRef]

2. H. Sun, H. Feng, and Y. Lu, “High resolution radar tomographic imaging using single-tone CW signals,” inProceedings of IEEE Radar Conference (IEEE, 2010), pp. 975–980. [CrossRef]

3. M. S. Roulston and D. O. Muhleman, “Synthesizing radar maps of polar regions with a Doppler-only method,” Appl. Opt. 36(17), 3912–3919 (1997). [CrossRef] [PubMed]

4. J. Li, R. L. Ewing, C. Berdanier, and C. Baker, “RF tomography of metallic objects in free space: preliminary results,” Proc. SPIE 9461, 94610S (2015). [CrossRef]

5. E. W. Mitchell, M. S. Hoehler, F. R. Giorgetta, T. Hayden, G. B. Rieker, N. R. Newbury, and E. Baumann, “Coherent laser ranging for precision imaging through flames,” Optica 5(8), 988–995 (2018). [CrossRef]

6. C. L. Matson and D. E. Mosley, “Reflective tomography reconstruction of satellite features-field results,” Appl. Opt. 40(14), 2290–2296 (2001). [CrossRef] [PubMed]

7. Y. Yan, J. Sun, X. Jin, Y. Zhou, Y. Zhi, and L. Liu, “Experimental research of circular incoherently synthetic aperture imaging ladar using chirped-laser and heterodyne detection,” Chin. Opt. Lett. 10(9), 091101 (2012). [CrossRef]

8. X. Jin, J. Sun, Y. Yan, Y. Zhou, and L. Liu, “Application of phase retrieval algorithm in reflective tomography laser radar imaging,” Chin. Opt. Lett. 9(1), 12801–12804 (2011). [CrossRef]

9. R. V. Chimenti, M. P. Dierking, P. E. Powers, J. W. Haus, and E. S. Bailey, “Experimental verification of sparse frequency linearly frequency modulated ladar signals modeling,” Opt. Express 18(15), 15400–15407 (2010). [CrossRef] [PubMed]

10. M. Bashkansky, R. L. Lucke, E. Funk, L. J. Rickard, and J. Reintjes, “Two-dimensional synthetic aperture imaging in the optical domain,” Opt. Lett. 27(22), 1983–1985 (2002). [CrossRef] [PubMed]

11. S. M. Beck, J. R. Buck, W. F. Buell, R. P. Dickinson, D. A. Kozlowski, N. J. Marechal, and T. J. Wright, “Synthetic-aperture imaging laser radar: laboratory demonstration and signal processing,” Appl. Opt. 44(35), 7621–7629 (2005). [CrossRef] [PubMed]

12. Z. W. Barber, F. R. Giorgetta, P. A. Roos, I. Coddington, J. R. Dahl, R. R. Reibel, N. Greenfield, and N. R. Newbury, “Characterization of an actively linearized ultrabroadband chirped laser with a fiber-laser optical frequency comb,” Opt. Lett. 36(7), 1152–1154 (2011). [CrossRef] [PubMed]

13. J. W. McCoy, N. Magotra, and B. K. Chang, “Coherent Doppler tomography-a technique for narrow band SAR,” inProceedings of IEEE Digital Avionics Systems Conference (IEEE, 1990), pp. 200–204. [CrossRef]

14. E. Swiercz, “Doppler radar tomography of rotated object in noisy environment based on time-frequency transformation,” inProceedings of IEEE Signal Processing Symposium, (IEEE, 2015), pp. 1–6. [CrossRef]

15. T. Sato, “Shape estimation of space debris using single-range Doppler interferometry,” IEEE Trans. Geosci. Remote Sens. 37(2), 1000–1005 (1999). [CrossRef]

16. Q. Wang, M. Xing, G. Lu, and Z. Bao, “SRMF-CLEAN imaging algorithm for space debris,” Acta Electron. Sin. 55(12), 3524–3533 (2007).

17. X. Wei, S. Gong, and X. Ding, “Narrow-band tomographic radar imaging of precession cone targets,” J. Syst. Eng. Electron. 23(6), 866–874 (2012). [CrossRef]

18. G. G. Fliss, “Tomographic radar imaging of rotating structures,” Proc. SPIE 1630, 199–207 (1992). [CrossRef]

19. E. Swiercz, “Application of the reassignment of time-frequency distributions to Doppler radar tomography imaging of a rotating multi-point object,” inProceedings of IEEE International Radar Symposium, (IEEE, 2016), pp. 1–5. [CrossRef]

20. R. M. Marino, R. N. Capes, W. E. Keicher, S. R. Kulkarni, J. K. Parker, L. W. Swezey, J. R. Senning, M. F. Reiley, and E. B. Craig, “Tomographic image reconstruction from laser radar reflective projections,” Proc. SPIE 999, 248–268 (1989). [CrossRef]

21. D. J. Sego, H. Griffiths, and M. C. Wicks, “Radar tomography using Doppler-based projections,” inProceedings of IEEE Radar Conference, (IEEE, 2011), pp. 403–408.

22. D. Mo, R. Wang, N. Wang, T. Lv, K.-S. Zhang, and Y.-R. Wu, “Three-dimensional inverse synthetic aperture lidar imaging for long-range spinning targets,” Opt. Lett. 43(4), 839–842 (2018). [CrossRef] [PubMed]

23. J. L. Walker, “Range-Doppler imaging of rotating objects,” IEEE Trans. Aero. Elec. Sys. aes- 16(1), 23–52 (1980). [CrossRef]

24. D. L. Mensa, S. Halevy, and G. Wade, “Coherent Doppler tomography for microwave imaging,” Proc. IEEE 71(2), 254–261 (1983). [CrossRef]

25. D. C. Munson Jr., J. D. O’Brien, and W. Jenkins, “A tomographic formulation of spotlight-mode synthetic aperture radar,” Proc. IEEE 71(8), 917–925 (1983). [CrossRef]

26. X. C. Zhang, “Three-dimensional terahertz wave imaging,” Philos Trans A Math Phys Eng Sci 362(1815), 283–299 (2004). [CrossRef] [PubMed]

27. H. S. Lui, T. Taimre, K. Bertling, Y. L. Lim, P. Dean, S. P. Khanna, M. Lachab, A. Valavanis, D. Indjin, E. H. Linfield, A. G. Davies, and A. D. Rakić, “Terahertz inverse synthetic aperture radar imaging using self-mixing interferometry with a quantum cascade laser,” Opt. Lett. 39(9), 2629–2632 (2014). [CrossRef] [PubMed]

28. H. Guerboukha, K. Nallappan, and M. Skorobogatiy, “Exploiting k-space/frequency duality toward real-time terahertz imaging,” Optica 5(2), 109–116 (2018). [CrossRef]

29. B. Borden and M. Cheney, “Synthetic-aperture imaging from high-Doppler-resolution measurements,” Inverse Probl. 21(1), 1–11 (2005). [CrossRef]

30. E. Swiercz, “Time-frequency transform used in radar Doppler tomography,” inProceedings of IEEE International Radar Symposium, (IEEE, 2014), pp. 1–5. [CrossRef]

31. C. L. Matson, “Tomographic image quality from E-field and intensity projections,” Opt. Commun. 186(1-3), 69–82 (2000). [CrossRef]

32. D. C. Munson Jr. and J. L. C. Sanz, “Image reconstruction from frequency-offset Fourier data,” Proc. IEEE 72(6), 661–669 (1984). [CrossRef]

33. C. A. Bouman and K. Sauer, “A unified approach to statistical tomography using coordinate descent optimization,” IEEE Trans. Image Process. 5(3), 480–492 (1996). [CrossRef] [PubMed]

34. J.-B. Thibault, K. D. Sauer, C. A. Bouman, and J. Hsieh, “A three-dimensional statistical approach to improved image quality for multislice helical CT,” Med. Phys. 34(11), 4526–4544 (2007). [CrossRef] [PubMed]

35. T. K. Moon, “The expectation-maximization algorithm,” IEEE Signal Process. Mag. 13(6), 47–60 (1996). [CrossRef]

36. D. C. Liu and J. Nocedal, “On the limited memory BFGS method for large scale optimization,” Math. Program. 45(1), 503–528 (1989). [CrossRef]

Item	Value
Laser wavelength	1550 nm
Laser power	2 W
Imaging distance	5 m
Spot diameter at the target	8 cm
Equivalent aperture	0.12 mm
Spin motor speed	100~1000 r/min
Incident angle θ	20°
AD sampling rate	10 MHz

Single-frequency LADAR super-resolution Doppler tomography for extended targets

Abstract

1. Introduction

2. Echo model

3. Imaging processing

3.1 Fast processing

3.2 Fine processing

3.3 Image resolution

4. Results and discussion

4.1 Experiment system

4.2 Comparison of fast processing with existing methods

4.3 Noise and speckle

4.4 Results of fine processing

5. Conclusions and outlook

Funding

References

Cited By

Figures (11)

Tables (1)

Equations (16)

Optics Express