Open AccessArticle

Estimation of High-Frequency Vibration Parameters for Airborne Terahertz SAR Using Chirplet Decomposition and LS Sequential Estimators

School of Electronics and Information Engineering, Beihang University, Beijing 100191, China

Department of Engineering, University of Cambridge, Cambridge CB2 1PZ, UK

School of Information and Electronics, Beijing Institute of Technology, Beijing 100081, China

Author to whom correspondence should be addressed.

Remote Sens. 2022, 14(14), 3416; https://doi.org/10.3390/rs14143416

Submission received: 27 May 2022 / Revised: 6 July 2022 / Accepted: 14 July 2022 / Published: 16 July 2022

(This article belongs to the Special Issue Advances in Synthetic Aperture Radar Remote Sensing)

Download

Browse Figures

Versions Notes

Abstract

Due to the short wavelength of the terahertz wave, airborne terahertz synthetic aperture radar (THz-SAR) suffers from echo phase errors caused by the high-frequency vibration of the platform. These errors will result in defocusing and the emergence of ghost targets, which will degrade the quality of the image. Therefore, it is necessary to compensate for phase errors in order to bring the image into focus. This paper proposes a multi-component high-frequency vibration parameter estimation method based on chirplet decomposition and least squares (LS) sequential estimators, which differs from other methods that can only be applied to simple harmonic vibrations. In particular, we first obtain the instantaneous chirp rate (ICR) of the signal by chirplet decomposition. Then, we employ the LS sequential estimators in conjunction with separable regression technique (SRT) to estimate vibration parameters. The estimated parameters are subsequently used to re-establish the ICR components for each vibration component and these parameters are further re-estimated to improve their accuracy. Based on the estimated parameters, phase compensation functions can be constructed to suppress the defocusing and ghost targets in airborne THz-SAR imaging. Simulated results on point targets and distributed imaging scenes demonstrate that the proposed method is accurate and reliable even at low signal-to-noise ratios (SNRs).

Keywords:

terahertz SAR; high-frequency vibration; chirplet decomposition; least squares sequential estimators; separable regression technique

1. Introduction

Terahertz (THz) waves are electromagnetic waves with a frequency between 100 GHz and 10 THz. Since their frequency falls between that of microwaves and infrared waves, THz waves share many characteristics with microwaves and infrared rays, including penetrability, high resolution, safety, and material identification. This has led to the widespread application of THz radar technology [1,2,3]. Especially, synthetic aperture radar (SAR) operating in the THz band has become one of the important research directions in the field of radar imaging [4,5]. Compared with conventional SAR operating in microwaves bands, THz-SAR offers high imaging resolution, good penetrability, and high frame rate imaging capability, giving it potential applications in many fields, such as video SAR, security checks, indoor object recognition, detecting micro-motion targets and stealth targets, etc. [1,2,3,4,5,6].

An important issue in THz-SAR is the compensation of the high-frequency vibration phase errors. Due to the influence of air flow and the flight characteristics of the platform, the actual airborne SAR system may deviate from the ideal flight trajectory, and the phase errors of the echo signal caused by this non-ideal motion of the platform may result in an unfocused imaging. Thus, the phase errors need to be estimated and compensated. Non-ideal motion can be divided into two parts: low-frequency components and high-frequency components. Low-frequency motion errors are usually caused by environmental changes such as air flow disturbance, while high-frequency vibration error components are specifically caused by high-speed rotation of aircraft engines or helicopter propellers. Because the physical mechanism of these two motion error components is different, it is generally believed that there is no coupling relationship between them. Several methods exist for compensating phase errors caused by the low-frequency components, such as using motion sensor measurements and autofocusing [7,8,9]; however, they are less effective for compensating errors caused by high-frequency components. Despite phase errors caused by high-frequency vibration being negligible in traditional microwave SAR imaging, these errors cannot be ignored in THz-SAR as the amplitude of high-frequency vibration is close to the wavelength of the THz-SAR signal. For the conventional motion compensation (MOCO) method based on sensor measurement data, the phase errors caused by the positioning errors of the radar antenna phase center (APC) should not be greater than π/4 to realize focusing imaging. It requires that the accuracy of the position measurement sensor must be less than λ/16, where λ is the wavelength of the radar signal. Taking the THz-SAR system with the working frequency of 100 GHz as an example, the threshold of positioning accuracy that can make the imaging results focus well is 0.188 mm. Even with very advanced filtering technology, the positioning accuracy of high-frequency vibration measurement sensors is far from this threshold. Therefore, an advanced signal processing method has to be used to estimate and compensate the phase errors caused by the high-frequency vibration of the platform in THz-SAR imaging, which is also the motivation of this paper.

There have been several methods to compensate the high-frequency vibration phase errors of THz-SAR. The most effective method is vibration isolation technology. However, the device is difficult to miniaturize, and the flexible connection between the radar and the platform may make the positioning of the radar inaccurate. Thus, the adaptive chirplet decomposition and local fractional Fourier transform (LFrFT) are applied to estimate high-frequency vibration parameters in [10,11,12,13,14,15,16]. Using advanced time–frequency representation approaches, these methods can improve the accuracy of azimuth instantaneous chirp rate (ICR) estimation and realize novel compensation of single harmonic vibration errors. In [17], a discrete sine frequency modulation transformation (DSFMT) based multi-component high-frequency vibration parameter estimation method is proposed. However, the need to search in multidimensional space imposes a significant computational burden. This computational burden is reduced in [18,19] by using intelligent optimization algorithms; nonetheless, the objective function may converge to a local extremum, which may deteriorate the estimation accuracy. In [20], the successive Doppler keystone transform (SDKT) is used to correct the range cell migration induced by vibration, and the analytical formula of the phase errors is deduced. However, the image may be unfocused when the amplitude is large.

Considering that the high-frequency vibration can be approximated as a sum of several simple harmonic vibrations with different frequencies, the echo signal is in the form of a multi-component sinusoidal phase modulation signal after eliminating the Doppler frequency modulation term caused by ideal motion and its ICR is in the form of the sum of sinusoidal models. Then, the task is transformed to estimate the ICR of the echo signal and its parameters. In [21], the ICR of the signal is estimated based on fractional Fourier transform (FrFT) combined with quasi-maximum likelihood (QML) and random sample consensus (RANSAC). Then, the parameters of high-frequency vibration can be estimated by employing the discrete Fourier transform (DFT) and a least square (LS) estimator. Nevertheless, this method assumes that the number of components of high-frequency vibration is known as a prior, and the errors may be large when the signal-to-noise ratio (SNR) is low.

In this paper, we propose a new method to estimate the THz-SAR multi-component high-frequency vibration parameters based on the chirplet decomposition and LS sequential estimators. Considering the form of the chirplet basis function, the ICR of the signal can be estimated more accurately using chirplet decomposition. It is different from the existing state of the art methods that the proposed method can estimate the parameters of high-frequency vibration under the condition of the unknown number of vibration components. Especially, we first estimate the parameters of the component with the largest amplitude using LS sequential estimators and separable regression technique (SRT). Next, we remove this estimated component in the ICR. We repeat the above two steps until there is no component with an amplitude above the set threshold in the remaining signal, then the number of vibration components can be obtained. Finally, the ICR component of each vibration component is constructed by using the estimation parameters, and these parameters are further re-estimated by using LS to obtain more accurate estimation values.

2. THz-SAR High-Frequency Vibration Error Model and Effect

The imaging geometry of the airborne THz-SAR is shown in Figure 1. The platform flies horizontally along the Y-axis at a velocity of V and a reference height of h. The point target P is in the imaging scene and its coordinate is (

x_{0}

y_{0}

0

). The slant range between P and the antenna phase center of radar is

R (t)

. Assume that

r_{0}

is the shortest slant range between the antenna phase center and P when the vibration is ignored. The slant range between antenna phase center and P is

r_{c} (t) = \sqrt{r_{0}^{2} + V^{2} {(t - t_{0})}^{2}}

(1)

where t is the slow time and

t_{0}

is the zero Doppler time.

The high-frequency vibration of the platform can be approximated as a sum of several simple harmonic vibrations with different frequencies [22], which is modeled as

r_{v} (t) = \sum_{m = 1}^{M} A_{m} \sin (2 π f_{m} t + φ_{m})

(2)

where M is the number of simple harmonic components.

A_{m}

f_{m}

, and

φ_{m}

are the amplitude, frequency, and initial phase of the m-th vibration component, respectively.

The actual antenna phase center of airborne SAR may deviate from the ideal flight trajectory because of the high-frequency vibration, and the slant range between the antenna phase center and P changes accordingly. Figure 2 is a schematic diagram of the vibration direction, where the position of the antenna phase center without vibration is E and the vibration makes it shift to E′. The angle between the vibration direction and the XEZ plane is

θ_{2}

, and the angle between the projection of the vibration direction on the XEZ plane and the X-axis is

θ_{1}

. Hence the coordinate of radar is (

Δ x (t)

V t + Δ y (t)

h + Δ z (t)

), where

{\begin{array}{l} Δ x (t) = r_{v} (t) \cos (θ_{2}) \cos (θ_{1}) \\ Δ y (t) = r_{v} (t) \sin (θ_{2}) \\ Δ z (t) = r_{v} (t) \cos (θ_{2}) \sin (θ_{1}) \end{array}

(3)

According to (3), the slant range between radar and P is

\begin{array}{l} R (t) = \sqrt{{(Δ x (t) - x_{0})}^{2} + {(V t + Δ y (t) - y_{0})}^{2} + {(h + Δ z (t))}^{2}} \\ = \sqrt{r_{c}^{2} (t) + r_{v}^{2} (t) - 2 Δ r_{v}} \end{array}

(4)

where

Δ r_{v} = x_{0} Δ x (t) - h Δ z (t) + (y_{0} - V t) Δ y (t)

(5)

For side-looking THz-SAR, because the synthetic aperture time is short, and the amplitude of high-frequency vibration is far less than the slant range, the vibration component in the flight direction can be ignored, and

R (t)

can be approximated as [22]

R (t) \approx r_{c} (t) + C r_{v} (t)

(6)

where

C = \cos (θ_{2}) \sin (\arctan (x_{0} / h) - θ_{1})

(7)

The component of high-frequency vibration along the line of sight (LOS) is the main factor affecting the signal phase according to (6). Hence, the calculation can be simplified by setting C to 1 without affecting the accuracy of the result [23].

Next, we present the phase error model for THz-SAR under high-frequency vibration and analyze the influence of the errors on THz-SAR imaging. Assume that the radar transmits a linear frequency modulated continuous wave (LFMCW) signal [24]. In one period, the radar transmission signal is expressed as

s (\tilde{τ}) = w_{r} (\tilde{τ}) \exp {j 2 π f_{0} \tilde{τ} + j π K_{r} {\tilde{τ}}^{2}}

(8)

where

\tilde{τ}

f_{0}

K_{r}

, and

w_{r}

are the fast time, center frequency, chirp rate, and envelope of range direction, respectively.

The echo signal of P received by the radar antenna is

s (\tilde{τ}, t) = w_{r} (\tilde{τ} - τ_{d}) w_{a} (t - t_{p}) \cdot \exp {j 2 π f_{0} (\tilde{τ} - τ_{d}) + j π K_{r} {(\tilde{τ} - τ_{d})}^{2}}

(9)

where

τ_{d} = 2 R (t) / c

, c is the velocity of light,

w_{a}

is the envelope of azimuth direction, and

t_{p}

is the time when the target is at the center of the azimuth beam.

After dechirp demodulation for sampling rate reduction, residual video phase (RVP) removal [7], and the time–frequency substitution

K_{r} \tilde{τ} \to f_{τ}

, the echo signal becomes

s_{0} (f_{τ}, t) = W_{r} (f_{τ}) w_{a} (t - t_{p}) \cdot \exp {- j 2 π (f_{0} + f_{τ}) τ_{d}}

(10)

where

f_{τ}

denotes the range frequency,

W_{r} (f_{τ})

is the range envelope of spectrum.

After range cell migration correction (RCMC), and the range compression through the inverse Fourier transform in range direction, we have

s_{r c} (τ, t) = p_{r} (τ - \frac{2 r_{0}}{c}) w_{a} (t - t_{p}) \cdot \exp {- j \frac{4 π}{λ} [r_{0} + \frac{V^{2} {(t - t_{p})}^{2}}{2 r_{0}}]} \cdot \exp {- j \frac{4 π}{λ} r_{v} (t)}

(11)

where

τ

is the new fast time corresponding to

f_{τ}

p_{r} (τ)

is the inverse Fourier transform of

W_{r} (f_{τ})

The phase of the echo signal is modulated by the following exponential term

s_{e r r} (t) = \exp {- j \frac{4 π}{λ} r_{v} (t)} = \exp {- j \frac{4 π}{λ} \sum_{m = 1}^{M} A_{m} \sin (2 π f_{m} t + φ_{m})}

(12)

Then, the phase errors can be expanded into a Bessel series as

s_{e r r} (t) = \sum_{n_{1} = - \infty}^{\infty} \sum_{n_{2} = - \infty}^{\infty} \dots \sum_{n_{M} = - \infty}^{\infty} \exp {- \sum_{m = 1}^{M} (j 2 n_{m} π f_{m} t + j n_{m} φ_{m})} \cdot \prod_{m = 1}^{M} J_{n_{m}} (z_{m})

(13)

where

z_{m} = 4 π A_{m} / λ

, and

J_{n_{m}} (z_{m})

is the n_mth-order Bessel function.

According to (11) and (13), we have

\begin{array}{l} s_{r c} (τ, t) = \sum_{n_{1} = - \infty}^{\infty} \sum_{n_{2} = - \infty}^{\infty} \dots \sum_{n_{M} = - \infty}^{\infty} p_{r} (τ - \frac{2 r_{0}}{c}) w_{a} (t - t_{p}) \cdot \exp {- j 2 π t \sum_{m = 1}^{M} n_{m} f_{m}} \\ \cdot \exp {- j \frac{4 π}{λ} [r_{0} + \frac{V^{2} {(t - t_{p})}^{2}}{2 r_{0}}] - j \sum_{m = 1}^{M} n_{m} φ_{m}} \cdot \prod_{m = 1}^{M} J_{n_{m}} (z_{m}) \end{array}

(14)

After the Fourier transform along the azimuth direction, (14) becomes

\begin{array}{l} S_{1} (τ, f_{t}) = \sum_{n_{1} = - \infty}^{\infty} \sum_{n_{2} = - \infty}^{\infty} \dots \sum_{n_{M} = - \infty}^{\infty} p_{r} (τ - \frac{2 r_{0}}{c}) W_{a} (f_{t} - f_{t_{p}} + \sum_{m = 1}^{M} n_{m} f_{m}) \\ \cdot \exp {- j \frac{4 π r_{0}}{λ} + j \frac{π}{K_{a}} {(f_{t} + \sum_{m = 1}^{M} n_{m} f_{m})}^{2}} \\ \cdot \exp {- j 2 π t_{p} (f_{t} + \sum_{m = 1}^{M} n_{m} f_{m}) - \sum_{m = 1}^{M} (j n_{m} φ_{m})} \cdot \prod_{m = 1}^{M} J_{n_{m}} (z_{m}) \end{array}

(15)

where

f_{t}

is the azimuth frequency,

W_{a} (f_{t})

is the azimuth envelope of spectrum,

f_{t_{p}}

is Doppler center frequency, and Doppler rate

K_{a} = - 2 V^{2} / (λ r_{0})

After the azimuth compression, we can get

\begin{array}{l} S_{a c} (τ, t) = \sum_{n_{1} = - \infty}^{\infty} \sum_{n_{2} = - \infty}^{\infty} \dots \sum_{n_{M} = - \infty}^{\infty} p_{r} (τ - \frac{2 r_{0}}{c}) p_{a} (t - t_{p} + \sum_{m = 1}^{M} \frac{n_{m} f_{m}}{K_{a}}) \\ \cdot \exp {j 2 π (t - t_{p} + \sum_{m = 1}^{M} \frac{n_{m} f_{m}}{K_{a}}) (f_{t_{p}} - \sum_{m = 1}^{M} n_{m} f_{m})} \\ \cdot \exp {- j \frac{4 π r_{0}}{λ} + j \frac{π}{K_{a}} {(\sum_{m = 1}^{M} n_{m} f_{m})}^{2}} \\ \cdot \exp {- j \sum_{m = 1}^{M} (2 π n_{m} f_{m} t_{p} + n_{m} φ_{m})} \cdot \prod_{m = 1}^{M} J_{n_{m}} (z_{m}) \end{array}

(16)

where

p_{a} (t)

is the amplitude of the azimuth impulse response and is in the form of a sinc function. It can be seen from (16) that there are paired echoes in the slow time domain when high-frequency vibration exists. The interval and intensity of paired echoes are related to the amplitude and frequency of high-frequency vibration [25]. These paired echoes make paired ghost targets appear on both sides of the real target. The distribution of these ghost targets in the slow time domain is

t = t_{p} - \sum_{m = 1}^{M} n_{m} f_{m} / K_{a}, (n_{m} = 0, \pm 1, \pm 2, \dots)

(17)

In order to eliminate the influence of paired echoes and to suppress the appearance of ghost targets, the phase errors caused by vibration must be compensated in the imaging process.

3. ICR Estimation with Chirplet Decomposition

It can be seen from (12) that the phase of the THz-SAR echo signal includes the multi-component sinusoidal terms caused by the vibration. After performing range compression and RCMC on the received signal, we can achieve a coarsely focused SAR image. Next, the image of the range bin with the strongest response is intercepted by the azimuth window. Then, the azimuth signal of the range bin with the strongest response can be obtained. The reference function is constructed according to the known azimuth velocity.

s_{r e f} (t) = \exp {j \frac{4 π}{λ} [r_{0} + \frac{V^{2} t^{2}}{2 r_{0}}]}

(18)

To eliminate the quadratic term in the phase, we multiply the azimuth signal with the above reference function. The obtained signal is in the form of a multi-component sinusoidal phase modulation signal.

s_{p h a s e} (t) = \exp [- j \frac{4 π}{λ} \sum_{m = 1}^{M} A_{m} \sin (2 π f_{m} t + φ_{m})]

(19)

The parameters of high-frequency vibration can be estimated according to this signal.

Chirplet decomposition is a method used to decompose the signal into a sum chirplet basis functions based on maximum likelihood. Each chirplet basis function corresponds to a local part of the signal, which can be used to capture the local time–frequency variation in the signal. The form of the chirplet basis function is

g (t) = {(π σ^{2})}^{- 0.25} \exp {- \frac{{(t - t_{g})}^{2}}{2 σ^{2}} + j ω (t - t_{g}) + j β {(t - t_{g})}^{2}}

(20)

where

σ

t_{g}

ω

, and

β

are the width, time center, initial frequency, and chirp rate, respectively.

The multi-component sinusoidal phase modulation signal can be decomposed into a sum of chirplet basis functions

s_{p h a s e} (t) = \sum_{k = 1}^{K} c_{k} g_{k} (t)

(21)

where

c_{k}

is the coefficient of the k-th basis function,

g_{k} (t) = \underset{g (t)}{\arg \max} {| 〈 s_{k} (t), g (t) 〉 |}^{2}

(22)

c_{k} = | 〈 s_{k} (t), g_{k} (t) 〉 |

(23)

s_{k} (t) = s_{k - 1} (t) - c_{k - 1} g_{k - 1} (t), s_{1} (t) = s_{p h a s e} (t)

(24)

It can be seen from (20) that the chirplet basis function is in the form of a linear frequency modulation signal, which is suitable for decomposing the signal whose frequency changes linearly with time. However, the frequency of the multi-component sinusoidal phase modulation signal changes nonlinearly with time. To solve this problem, the slow time can be divided into numerous short time intervals. The phase of the signal in each time interval can be approximated as a second-order polynomial, and the signal can be approximated as a linear frequency modulation signal. The division of the slow time domain can be realized by a sliding narrow Gaussian window. The width of the window is determined by

σ

, and the position of the window is determined by

t_{g}

. The phase of the signal in the i-th window

[t_{i}, t_{i} + Δ t]

can be approximately regarded as a second-order polynomial

s_{p h a s e} (t) \approx \exp {- j \frac{4 π}{λ} [ω_{i} (t - t_{i}) + β_{i} {(t - t_{i})}^{2}]}

(25)

where

ω_{i}

and

β_{i}

are the initial frequency and chirp rate of the multi-component sinusoidal phase modulation signal in this window. When

σ

is small enough,

β_{i}

can be approximately regarded as the value of ICR at the time of

t_{i}

. Therefore, the multi-component sinusoidal phase modulated signal in each window can be represented by a chirplet basis function, which has the same time center, initial frequency, and chirp rate as the signal. For the signal in the i-th window,

σ_{i}

and

t_{g, i}

of the corresponding chirplet basis function are determined by the width and position of the window, so only two parameters need to be estimated. After chirplet decomposition, the ICR of the multi-component sinusoidal phase modulation signal can be obtained by arranging the

β_{i}

in time order. To improve the estimation accuracy, adjacent windows can overlap each other.

Suppose the pulse repetition frequency (PRF) is

f_{R}

. The sampling time interval is

Δ t = 1 / f_{R}

. It can be obtained according to (19) that the discretized ICR is

\begin{array}{l} β (n) = \frac{16 π^{3}}{λ} \sum_{m = 1}^{M} f_{m}^{2} A_{m} \sin (2 π f_{m} n Δ t + φ_{m}) \\ = \sum_{m = 1}^{M} [B_{m, 1} \sin (2 π f_{m} n Δ t) + B_{m, 2} \cos (2 π f_{m} n Δ t)], n = 1, 2, \dots, N \end{array}

(26)

where N is the number of sampling points, and

{\begin{array}{l} B_{m, 1} = 16 π^{3} A_{m} f_{m}^{2} \cos (φ_{m}) / λ \\ B_{m, 2} = 16 π^{3} A_{m} f_{m}^{2} \sin (φ_{m}) / λ \end{array}

(27)

4. Parameter Estimation with LS Sequential Estimators and SRT

LS sequential estimators [26] can be used to sequentially estimate the parameters of each component of high-frequency vibration according to the discretized ICR. For the convenience of analysis, we assume that all components in (26) are arranged in descending order of amplitude, that is,

B_{1, 1}^{2} + B_{1, 2}^{2} > B_{2, 1}^{2} + B_{2, 2}^{2} > \dots > B_{M, 1}^{2} + B_{M, 2}^{2}

(28)

The parameters of the component with the largest amplitude can be estimated by LS

\begin{array}{l} ({\hat{B}}_{1, 1}, {\hat{B}}_{1, 2}, {\hat{f}}_{1}) = \underset{B_{1, 1}, B_{1, 2}, f_{1}}{\arg \min} H (B_{1, 1}, B_{1, 2}, f_{1}) \\ = \underset{B_{1, 1}, B_{1, 2}, f_{1}}{\arg \min} {[Y_{1} - X (f_{1}) α_{1}]}^{T} [Y_{1} - X (f_{1}) α_{1}] \end{array}

(29)

where

H (B_{1, 1}, B_{1, 2}, f_{1})

is the objective function and

Y_{1}

is the discretized ICR,

{\begin{array}{l} α_{1} = {[B_{1, 1}, B_{1, 2}]}^{T} \\ Y_{1} = {[β (1), β (2), \dots, β (N)]}^{T} \\ X (f_{1}) = [C (f_{1}), S (f_{1})] \\ S (f_{1}) = {[\sin (2 π f_{1} Δ t), \dots, \sin (2 π f_{1} N Δ t)]}^{T} \\ C (f_{1}) = {[\cos (2 π f_{1} Δ t), \dots, \cos (2 π f_{1} N Δ t)]}^{T} \end{array}

(30)

There are three parameters (

B_{1, 1}

B_{1, 2}

, and

f_{1}

) that need to be estimated at the same time, resulting in a huge computational burden. The estimation process can be simplified by the separable regression technique (SRT) [27].

SRT is a technique to simplify the maximum likelihood estimation problem by reducing the number of parameters in regression problems. In regression analysis, for an objective function

l (x)

with q independent variables, we define

{\begin{array}{l} x = [x^{(1)}, x^{(2)}] \\ x^{(1)} = [x_{1}, x_{2}, \dots, x_{p}] \\ x^{(2)} = [x_{p + 1}, x_{p + 2}, \dots, x_{q}] \end{array}

(31)

Assume that

{\hat{x}}^{(2)} (x^{(1)})

is the maximum likelihood value of

x^{(2)}

when

x^{(1)}

is fixed.

{\hat{x}}^{(2)}

can be derived from

\frac{\partial l (x)}{\partial x_{i}} = 0, i = p + 1, p + 2, \dots, q

(32)

Then, the objective function can be replaced by

L (x^{(1)}) = l [x^{(1)}, {\hat{x}}^{(2)} (x^{(1)})]

(33)

In this way, the number of parameters in regression analysis is reduced to p by using the new objective function.

The objective function in (29) can be replaced by the new objective function with

f_{1}

as the only independent variable using SRT. The LS estimate of

α_{1}

when

f_{1}

is fixed is

{\hat{α}}_{1} = {[X^{T} (f_{1}) X (f_{1})]}^{- 1} X^{T} (f_{1}) Y_{1}

(34)

According to Equation (34), the new objective function can be obtained as

H_{s i m} (f_{1}) = [Y_{1}^{T} (I - P_{X (f_{1})}) Y_{1}]

(35)

where I is the identity matrix,

P_{X (f_{1})} = X (f_{1}) {[X^{T} (f_{1}) X (f_{1})]}^{- 1} X^{T} (f_{1})

(36)

The number of parameters which need to be estimated in the LS process can be reduced from three to one by substituting (35) into (29), that is, Equation (29) can be simplified to

{\hat{f}}_{1} = \underset{f_{1}}{\arg \min} [Y_{1}^{T} (I - P_{X (f_{1})}) Y_{1}]

(37)

The value of

{\hat{α}}_{1}

can be obtained by substituting the estimated value of

{\hat{f}}_{1}

into (34). Then, we can obtain the estimated value of other vibration parameters according to (27), and

{\begin{array}{l} {\hat{A}}_{1} = \frac{λ}{16 π^{3} {\hat{f}}_{1}^{2}} \sqrt{{\hat{B}}_{1, 1}^{2} + {\hat{B}}_{1, 2}^{2}} \\ {\hat{φ}}_{1} = \tan^{- 1} (- \frac{{\hat{B}}_{1, 1}}{{\hat{B}}_{1, 2}}) + φ_{a d j, 1} \end{array}

(38)

where

φ_{a d j, 1} = {\begin{matrix} π, {\hat{B}}_{1, 2} > 0 \\ 0, {\hat{B}}_{1, 2} < 0 \end{matrix}

(39)

After estimating vibration parameters of one simple harmonic, we construct the corresponding ICR component and remove this component from

Y_{1}

Y_{2} = Y_{1} - X ({\hat{f}}_{1}) {\hat{α}}_{1}

(40)

Other parameters can be estimated by repeating the above processes. For cases where the number of vibration components is unknown, a threshold can be set for the amplitude. Starting from the second iteration, if the estimated amplitude is lower than the threshold, the corresponding vibration component is ignored and the iteration is ended. The constraint of the phase errors is

| Δ ϕ (t) | \leq \frac{π}{4}

(41)

According to (12), the phase errors of the echo signal caused by the m-th vibration component are

Δ ϕ_{m} (t) = \frac{4 π}{λ} A_{m} \sin (2 π f_{m} t + φ_{m})

(42)

Thus, the threshold for ending the iteration is set to

A_{m} \leq \frac{λ}{16}

(43)

The estimated high-frequency vibration after the iteration is

{\hat{r}}_{v} (t) = \sum_{m = 1}^{\hat{M}} {\hat{A}}_{m} \sin (2 π {\hat{f}}_{m} t + {\hat{φ}}_{m})

(44)

In the m-th iteration process of estimation, since

Y_{m}

contains the ICR components corresponding to several simple harmonic components, the parameter estimation of the m-th component is affected by the remaining components, causing errors in the estimation results. To reduce the errors, the ICR components corresponding to the remaining components in

Y_{m}

can be removed according to the estimated parameters

{Y^{'}}_{m} = Y_{m} - \sum_{i = m + 1}^{\hat{M}} X ({\hat{f}}_{i}) {\hat{α}}_{i}, m < \hat{M}

(45)

The accuracy of the estimation result can be enhanced by repeating the estimation step using

{Y^{'}}_{m}

instead of

Y_{m}

. The search space for frequency can be reduced according to the coarse estimation results to lessen the computational burden in the re-estimation process.

5. The Processing Flow of the Proposed Parameters Estimation Method

The processing flow of the proposed high-frequency vibration parameter estimation method based on chirplet decomposition and LS sequential estimators is summarized as follows:

Step 1: Perform range compression and RCMC on the echo signal. Intercept the azimuth signal of the range bin where the strongest response is located. Then, the multi-component sinusoidal phase modulation signal is obtained by eliminating the quadratic term.

Step 2: Estimate the ICR of the signal with chirplet decomposition.

Step 3: Perform LS to estimate the frequency of the component with the largest amplitude in the current ICR, and further estimate the amplitude and initial phase corresponding to the frequency.

Step 4: Remove the estimated components in ICR.

Step 5: Calculate the threshold and compare it with the amplitude of the currently estimated vibration component. If the latter is larger, return to step 3. Otherwise, proceed to step 6.

Step 6: Calculate the

{Y^{'}}_{m}

corresponding to each vibration component according to the estimated parameters and perform LS to re-estimate the parameters.

Figure 3 is the flow chart of proposed method. After estimating parameters of all components of high-frequency vibration, the phase error compensation function for the echo signal can be constructed as follows

s_{c o m} (t) = \exp [- j \frac{4 π}{λ} \sum_{m = 1}^{\hat{M}} {\hat{A}}_{m} \sin (2 π {\hat{f}}_{m} t + {\hat{φ}}_{m})]

(46)

After compensating the phase errors of high-frequency vibration, the focusing imaging can be completed by conventional imaging processing algorithms such as the range migration algorithm (RMA) [7].

6. Simulation Results and Discussion

To verify the effectiveness of the high-frequency vibration parameter estimation method proposed in this paper, simulation experiments were carried out on the multi-component high-frequency vibration error estimation and compensation of airborne THz-SAR, and the result of the proposed method (Chirplet-LSSE-SRT) is compared with the parameter estimation method based on FrFT combined with QML and SANSAC (FrFT-QML-SANSAC) in [21]. Firstly, the imaging simulation of a point target was carried out to verify the influence of high-frequency vibration on imaging. Then, the echo signals of point targets and distributed imaging scenes with isolated dominant scatter points were simulated, respectively. Then, the vibration parameter estimation, phase error compensation, and imaging processing were performed. The settings of system parameters in the simulation experiment are shown in Table 1. The radar transmits LFMCW signal and the farthest slant range was set to 850 m, so the bandwidth of the intermediate frequency (IF) signal after dechirp demodulation was about 188.9 MHz. The fast time length of extracted data for imaging processing was chosen as 22 μs. This means that the number of range samples is 7040.

6.1. Influence of High-Frequency Vibration on Imaging

Assume the scatter point at the center of the scene as the target. The echo signals are simulated with the platform in three different vibration states. The azimuth duration of the echo signal is 0.370 s, and there is vibration throughout. While the synthetic aperture time is 0.185 s, the azimuth resolution of 0.1 m can be obtained. The imaging results performing the RMA imaging algorithm are shown in Figure 4. Figure 4a shows the imaging results without vibration; Figure 4b shows the imaging results of the platform with simple harmonic vibration whose expression is

r_{v} (t) = A_{1} \sin (2 π f_{1} t + φ_{1})

, where

A_{1}

= 1.5 mm,

f_{1}

= 18.3 Hz,

φ_{1} = 5 π / 6

; Figure 4c shows the imaging results with dual-components high-frequency vibration whose expression is

r_{v} (t) = A_{1} \sin (2 π f_{1} t + φ_{1}) + A_{2} \sin (2 π f_{2} t + φ_{2})

, where

A_{2}

= 1.0 mm,

f_{2}

= 35 Hz,

φ_{2} = 5 π / 6

. It can be seen from Figure 4 that ghost targets occur in the imaging results of the point target when high-frequency vibration exists. The ghost targets caused by simple harmonic vibration are equally spaced, while the distribution of ghost targets caused by multi-component high-frequency vibration is more complex.

6.2. Simulation Results of Point Targets

The simulation experiment was carried out with

3 \times 3

lattice targets located in the center of the scene, and the interval distance between points in azimuth direction and ground range direction is 10 m. The high-frequency vibration errors of the platform are

r_{v} (t) = A_{1} \sin (2 π f_{1} t + φ_{1}) + A_{2} \sin (2 π f_{2} t + φ_{2})

, where

A_{1}

= 1.5 mm,

f_{1}

= 18.3 Hz,

φ_{1} = 5 π / 6

A_{2}

= 1.0 mm,

f_{2}

= 35 Hz,

φ_{2} = 5 π / 6

. We simulated the echo signal of point targets according to system parameters in Table 1 with the Gaussian white noise and signal-to-noise ratio (SNR) being 5 dB. The multi-component sinusoidal phase modulation signal was extracted, and then the ICR curves were estimated using the chirplet decomposition method and FrFT, respectively, which are shown in Figure 5.

LS sequential estimators were performed to sequentially estimate the parameters of high-frequency vibration. The objective function curves obtained in the first three iterations are shown in Figure 6a–c, where the abscissas of the points at the lowest positions of the curves are the vibration frequencies. The amplitudes and initial phases calculated according to (35) and (38) are listed in Table 2. Since the amplitude estimation value of the third component is lower than the threshold, the iteration was terminated after the third iteration and the parameters of the first two estimated vibration components are regarded as the coarse estimate of the high-frequency vibration parameters of the platform. It can be seen that this method can effectively obtain the number of vibration components. The parameters of the first component obtained by re-estimating are shown in Figure 6d and Table 2, respectively. It can be seen that LS sequential estimators have high accuracy.

We compared our method with the existing approach FrFT-QML-RANSAC [21], and the estimation results of the two methods are shown in Table 3. It can be seen that the Chirplet-LSSE-SRT method is more accurate than the FrFT-QML-RANSAC method. The phase error curve based on the estimated parameters in Table 3 is shown in Figure 7. It can be seen that the phase error curve obtained using the proposed method is closer to the true value, while the phase error curve obtained using the FrFT-QML-RANSAC method has large errors near the extreme points. The residual phase errors after phase error compensation are shown in Figure 8. The phase error corresponding to the red dotted lines is ±π/4, which represents the phase error threshold for focusing the image. It can be seen that some parts of the residual phase error curve of the FrFT-QML-RANSAC method exceed the threshold, while the residual phase error curve of our method is wholly less than the threshold, meeting the requirements for focusing imaging.

The imaging results without phase error compensation are shown in Figure 9a. It shows that the high-frequency vibration of the platform causes serious defocus and ghost targets in the imaging results. The imaging results of the two methods after compensating the phase errors are shown in Figure 9b,c. It can be seen that the imaging results obtained using the proposed method have better focusing performance. To more accurately evaluate the imaging quality of point targets, we analyzed the range and azimuth responses of all point targets in the imaging results. For this purpose, we numbered the point targets in Figure 9 as P1, P2, …, P9, in the order from left to right and from top to bottom. The measured quality metrics include peak sidelobe ratio (PSLR), integrated sidelobe ratio (ISLR), and impulse response width (IRW). All three metrics of the range responses of point targets obtained by different methods are almost identical. The measured values of PSLR, ISLR, and IRW in range direction are −13.31 dB, −8.53 dB, and 0.20 m, respectively. This is consistent with the expected result. After different vibration error compensation methods, the residual phase errors mainly affect the imaging quality in the azimuth direction but have no effect on the imaging quality in the range direction. The azimuth responses of two selected points P4 and P9 are shown in Figure 10. The imaging quality metrics of the azimuth response for all point targets are listed in Table 4. It can be seen that the PSLRs and ISLRs of the point target imaging results of the proposed method are lower than that of the FrFT-QML-RANSAC method, and they are closer to that of the perfect compensation. The IRWs obtained by the two methods are the same, which are all close to the results of perfect compensation.

To verify the effectiveness of the proposed method under different SNRs, we performed 100 Monte Carlo runs to obtain the estimation results. The root mean square error (RMSE) curves of amplitudes, frequencies, and initial phases are shown in Figure 11. Note that the RMSE curves are not smooth as outliers occasionally occur in results when noise exists, which has a huge impact on the value of the RMSE. The RMSE curve shows that the accuracy of parameter estimation decreases with the decrease of SNR. It can be seen from the RMSE curves that the estimation accuracy of the proposed Chirplet-LSSE-SRT method is better than that of the FrFT-QML-RANSAC method when the SNR is between 5 and 10 dB.

6.3. Simulation Results of Distributed Imaging Scene

A high-resolution well-focused real SAR image with an isolated dominant scatter point is used to generate echo data of distributed scenes. The simulation method is the same as that in Refs. [22,23]. In the echo simulation, the vibration errors with the same parameters as that in Section 6.2 are added. The vibration parameters are, respectively, estimated by the two methods. The obtained ICR curves are shown in Figure 12, and the parameter estimation results are listed in Table 5. The phase error curves obtained according to the parameter estimation and the residual phase errors are shown in Figure 13 and Figure 14, respectively. It can be seen that the estimation results of high-frequency vibration parameters using the proposed Chirplet-LSSE-SRT method are more accurate than those of the FrFT-QML-RANSAC method when there is an isolated dominant point in the distributed imaging scene.

Figure 15a,b show the imaging results when the high-frequency vibration does not exist and when the phase errors of high-frequency vibration are not compensated, respectively. It can be seen that the high-frequency vibration of the platform causes serious defocus and ghost targets in the imaging results. The imaging results obtained after phase error compensation using the estimates of the two methods are shown in Figure 15c,d. It can be seen that the focusing quality of Figure 15c is improved compared with Figure 15b, but the defocus phenomenon and ghost targets are not completely suppressed. In comparison, Figure 15d has no obvious defocus phenomenon or ghost target and the focusing quality of it is close to Figure 15a.

To further quantitatively evaluate the image quality of Figure 15, we calculated the image entropy and image contrast of each image. Image entropy is an effective metric used to evaluate the focusing quality of images [28]. Images with smaller image entropy usually have better focusing quality. For an image

G (p, q)

with size

P_{g} \times Q_{g}

, the calculation formula of image entropy is

E_{g} = \ln S_{g} - \frac{1}{S_{g}} \sum_{p = 1}^{P_{g}} \sum_{q = 1}^{Q_{g}} {| G (p, q) |}^{2} \ln ({| G (p, q) |}^{2})

(47)

where

S_{g} = \sum_{p = 1}^{P_{g}} \sum_{q = 1}^{Q_{g}} {| G (p, q) |}^{2}

(48)

Image contrast is another common metric used to evaluate the focusing quality of images [29]. Images with larger contrast values usually have better focusing quality. The calculation formula of image contrast is

C_{g} = \frac{s t d ({| G (p, q) |}^{2})}{m e a n ({| G (p, q) |}^{2})}

(49)

where

s t d (\cdot)

refers to the standard deviation operator, and

m e a n (\cdot)

refers to the average operator. According to Equations (47) and (49), the image entropy and image contrast calculated for each image in Figure 15 are listed in Table 6. It can be seen that compared with Figure 15c, the image entropy of Figure 15d is smaller, the image contrast is larger, and the two metric values are both closer to that of Figure 15a. It indicates that the imaging quality of Figure 15d is better than that of Figure 15c.

7. Conclusions

High-frequency vibration of the platform leads to defocus and ghost targets in the imaging results of the airborne THz-SAR. In the absence of precise measurement devices for vibration compensation, the conventional MOCO methods cannot meet the requirements of high-frequency vibration error compensation. In order to improve the focusing quality of the image, it is necessary to estimate and compensate the phase errors caused by vibration. Therefore, a novel high-frequency vibration parameter estimation method based on chirplet decomposition and LS sequential estimator is proposed in this paper. The advantage of the proposed method is that it can estimate the parameters of multi-component high-frequency vibration under the condition of the unknown number of vibration components. In addition, the proposed method can obtain sufficiently accurate phase error estimation that meets the focusing requirements under low SNR. Simulation results of multi-component high-frequency vibration of point targets and distributed imaging scenes verify the effectiveness of this method. It should be noted that the performance of the proposed method will be limited, if the actual vibration of the platform cannot be approximated due to the superposition of several simple harmonic vibrations of different frequencies, such as the broadening of the vibration spectrum.

Author Contributions

Conceptualization, Z.H. and J.S.; methodology, Z.H.; software, Z.H.; validation, J.S., Q.L. and T.S.; formal analysis, Z.H.; investigation, J.S.; resources, J.S.; data curation, Z.H.; writing—original draft preparation, Z.H.; writing—review and editing, J.S., Q.L. and T.S.; visualization, T.S.; supervision, J.S.; project administration, J.S. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the National Natural Science Foundation of China, Grant No. 62131001 and 62171029.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

References

Knipper, R.; Brahm, A.; Heinz, E.; May, T.; Notni, G.; Meyer, H.-G.; Tunnermann, A.; Popp, J. THz Absorption in Fabric and Its Impact on Body Scanning for Security Application. IEEE Trans. Terahertz Sci. Technol. 2015, 5, 999–1004. [Google Scholar] [CrossRef]
Wang, Q.; Zhou, H.; Wang, Y.; Xia, R.; Liu, Q.; Zhao, B.; Wang, Y. Super-Resolution Imaging Using Very Deep Convolutional Network in Terahertz NDT Field. In Proceedings of the International Charter: Space and Major Disasters (ICSMD), Xi’an, China, 15–17 October 2020; pp. 438–441. [Google Scholar]
Ding, X.; Xu, J.; Chen, Q.; Fu, C.; Ni, L. Retrieval of ice cloud properties with a dual frequency optimal estimation algorithm for terahertz and millimeter wave cloud radar. In Proceedings of the IET International Radar Conference, Online Conference, 4–6 November 2020; pp. 523–527. [Google Scholar]
Zuo, F.; Min, R.; Pi, Y.; Li, J.; Hu, R. Improved Method of Video Synthetic Aperture Radar Imaging Algorithm. IEEE Geosci. Remote Sens. Lett. 2019, 16, 897–901. [Google Scholar]
Kim, S.; Fan, R.; Dominski, F. ViSAR: A 235 GHz radar for airborne applications. In Proceedings of the 2018 IEEE Radar Conference, Oklahoma City, OK, USA, 23–27 April 2018; pp. 1549–1554. [Google Scholar]
Batra, A.; Sheikh, F.; Khaliel, M.; Wiemeler, M.; Göhringer, D.; Kaiser, T. Object Recognition in High-Resolution Indoor THz SAR Mapped Environment. Sensors 2022, 22, 3762. [Google Scholar] [PubMed]
Carrara, W.G.; Goodman, R.S.; Majewski, R.M. Spotlight Synthetic Aperture Radar: Signal Processing Algorithms; Artech House: Norwood, MA, USA, 1995; pp. 29, 67–70. [Google Scholar]
Chen, Z.; Zhang, Z.; Zhou, Y.; Wang, P.; Qiu, J. A Novel Motion Compensation Scheme for Airborne Very High Resolution SAR. Remote Sens. 2021, 13, 2729. [Google Scholar]
Ren, Y.; Tang, S.; Dong, Q.; Sun, G.; Guo, P.; Jiang, C.; Han, J.; Zhang, L. An Improved Spatially Variant MOCO Approach Based on an MDA for High-Resolution UAV SAR Imaging with Large Measurement Errors. Remote Sens. 2022, 14, 2670. [Google Scholar]
Wang, Y.; Wang, Z.; Zhao, B.; Xu, L. Compensation for high-frequency vibration of platform in SAR imaging based on adaptive chirplet decomposition. IEEE Geosci. Remote Sens. Lett. 2016, 13, 792–795. [Google Scholar]
Li, Y.; Ding, L.; Zheng, Q.; Zhu, Y.; Sheng, J. A novel high-frequency vibration error estimation and compensation algorithm for THz-SAR imaging based on local FrFT. Sensors 2020, 20, 2669. [Google Scholar]
Wang, Z.; Wang, Y.; Xu, L. Parameter Estimation of Hybrid Linear Frequency Modulation-Sinusoidal Frequency Modulation Signal. IEEE Signal Processing Lett. 2017, 24, 1238–1241. [Google Scholar]
Wang, Z.; Wang, Y.; Xu, L. Method for High-resolution SAR Imaging Based on Inverse Radon Transform. J. Eng. 2019, 20, 6951–6955. [Google Scholar]
Wang, Z.; Wang, Y.; Xu, L. Time-Frequency Ridge-based Parameter Estimation for Sinusoidal Frequency Modulation Signals. In Proceedings of the International Conference in Communications, Signal Processing, and Systems (CSPS), Harbin, China, 14–17 July 2017; pp. 1375–1380. [Google Scholar]
Xia, H.; Chen, Q.; Li, Y.; Fu, C.; Wang, H. A high frequency vibration compensation approach in terahertz SAR based on wavelet multi-resolution analysis. In Proceedings of the Inter-American Social Security Conference (CISS), Shanghai, China, 10–12 October 2018; pp. 1–5. [Google Scholar]
Xia, W.; Huang, L. Target vibration estimation in SAR based on phaseanalysis method. EURASIP J. Adv. Signal Processing 2016, 94, 1–11. [Google Scholar]
Sun, Z.; Chen, J.; Cao, X.; Guo, L. Parameter estimation of multicomponent SFM signals based on discrete sinusoidal frequency modulation transform. Syst. Eng. Electron. 2012, 34, 1973–1979. [Google Scholar]
Wang, Y.; Wang, Z.; Zhao, B.; Xu, L. Enhancement of azimuth focus performance in high-resolution SAR imaging based on the compensation for sensors platform vibration. IEEE Sens. J. 2016, 16, 6333–6345. [Google Scholar] [CrossRef]
Wang, Y.; Wang, Z.; Zhao, B.; Xu, L. Parameters Estimation of Sinusoidal Frequency Modulation Signal with Application in Synthetic Aperture Radar Imaging. J. Appl. Remote Sens. 2016, 10, 020502. [Google Scholar]
Zhang, Y.; Sun, J.; Lei, P.; Hong, W. SAR-based paired echo focusing and suppression of vibrating targets. IEEE Trans. Geosci. Remote Sens. 2014, 52, 7593–7605. [Google Scholar]
Li, Y.; Wu, Q.; Wu, J.; Li, P.; Zheng, Q.; Ding, L. Estimation of High-Frequency Vibration Parameters for Terahertz SAR Imaging Based on FrFT With Combination of QML and RANSAC. IEEE Access 2021, 9, 5485–5496. [Google Scholar]
Zhang, Y.; Sun, J.; Lei, P.; Wang, H. High-frequency vibration compensation of helicopter-borne THz-SAR. IEEE Trans. Aerosp. Electron. Syst. 2016, 52, 1460–1466. [Google Scholar]
Sun, J.; Hao, Z.; Li, Q.; Li, D. Vibration Compensation of Airborne Terahertz SAR Based on Along Track Interferometry. IEEE Geosci. Remote Sens. Lett. 2022, 19, 1–5. [Google Scholar] [CrossRef]
Charvat, G.L. Small and Short-Range Radar Systems; CRC Press: Boca Raton, FL, USA, 2014; pp. 29, 69–76. [Google Scholar]
He, Q.; Zhang, Q.; Luo, Y.; Sun, L. Sinusoidal frequency modulation Fourier-Bessel series for multicomponent SFM signal estimation and separation. Math. Probl. Eng. 2017, 2017, 1–14. [Google Scholar] [CrossRef] [Green Version]
Prasad, A.; Kundu, D.; Mitra, A. Sequential estimation of the sum of sinusoidal model parameters. J. Stat. Plan. Inference 2008, 138, 1297–1313. [Google Scholar] [CrossRef] [Green Version]
Richards, F.S.G. A Method of Maximum-Likelihood Estimation. J. R. Stat. Soc. Ser. B Methodol. 1961, 23, 469–475. [Google Scholar] [CrossRef]
Wang, J.; Li, Y.; Song, M.; Huang, P.; Xing, M. Noise Robust High-Speed Motion Compensation for ISAR Imaging Based on Parametric Minimum Entropy Optimization. Remote Sens. 2022, 14, 2178. [Google Scholar] [CrossRef]
Jia, X.; Song, H.; He, W. A Novel Method for Refocusing Moving Ships in SAR Images via ISAR Technique. Remote Sens. 2021, 13, 2738. [Google Scholar] [CrossRef]

Figure 1. Imaging geometry of THz-SAR.

Figure 2. Vibration direction diagram.

Figure 3. Flow chart of the proposed method.

Figure 4. Single point target imaging results. (a) Without vibration; (b) simple harmonic vibration; (c) multi-component sinusoidal vibration.

Figure 5. ICR estimation.

Figure 6. Objective function curve. (a) The first iteration; (b) the second iteration; (c) the third iteration; (d) re-estimate of the first component.

Figure 7. Estimation of phase errors.

Figure 8. Residual phase errors.

Figure 9. Imaging results. (a) Without compensation; (b) FrFT-QML-RANSAC; (c) Chirplet-LSSE-SRT. P4 and P9 are two selected points for further azimuth responses analyze.

Figure 10. Azimuth response of point targets. (a) Point target P4 of Figure 9c; (b) point target P9 of Figure 9c.

Figure 11. RMSE in different SNRs. (a) Amplitude

A_{1}

; (b) amplitude

A_{2}

; (c) frequency

f_{1}

; (d) frequency

f_{2}

; (e) initial phase

φ_{1}

; (f) initial phase

φ_{2}

Figure 11. RMSE in different SNRs. (a) Amplitude

A_{1}

; (b) amplitude

A_{2}

; (c) frequency

f_{1}

; (d) frequency

f_{2}

; (e) initial phase

φ_{1}

; (f) initial phase

φ_{2}

Figure 12. ICR estimation of the distributed imaging scene.

Figure 13. Estimate of phase errors of the distributed imaging scene.

Figure 14. Residual phase errors of the distributed imaging scene.

Figure 15. Imaging results of the distributed imaging scene. (a) Without vibration; (b) without compensation; (c) FrFT-QML-RANSAC; (d) Chirplet-LSSE-SRT.

Table 1. System parameters of the simulation.

Parameter	Value
Center frequency	216 GHz
Azimuth beam width	2.6°
Azimuth resolution	0.1 m
Range resolution	0.2 m
Slant range of scene center	800 m
Platform height	200 m
Platform velocity	30 m/s
Pulse width	30 μs
Signal bandwidth	1.0 GHz
Sample frequency	320 MHz
Pulse repetition frequency	6000 Hz

Table 2. Estimate of high-frequency vibration parameters obtained by LS sequential estimators.

Parameter	Component 1	Component 2	Component 3	Re-Estimate of Component 1
$f / Hz$	34.900	18.300	49.640	35.000
$A / mm$	0.939	1.508	0.027	1.068
$φ / rad$	2.595	2.604	6.070	2.599

Table 3. Estimate of high-frequency vibration parameters.

Parameter	True Value	Chirplet-LSSE-SRT	FrFT-QML-RANSAC
$f_{1} / Hz$	18.300	18.300	18.310
$f_{2} / Hz$	35.000	35.000	35.020
$A_{1} / mm$	1.500	1.508	1.536
$A_{2} / mm$	1.000	1.068	0.957
$φ_{1} / rad$	2.618	2.604	2.601
$φ_{2} / rad$	2.618	2.599	2.585

Table 4. Measured imaging quality metrics of azimuth response.

Target No.	Perfect Compensation			FrFT-QML-RANSAC			Chirplet-LSSE-SRT
Target No.	PSLR/dB	ISLR/dB	IRW/m	PSLR/dB	ISLR/dB	IRW/m	PSLR/dB	ISLR/dB	IRW/m
P1	−13.23	−8.67	0.10	−10.88	−6.62	0.11	−13.20	−8.10	0.11
P2	−13.27	−8.49	0.10	−10.92	−6.11	0.10	−13.23	−8.02	0.10
P3	−13.29	−8.60	0.10	−11.08	−6.35	0.11	−13.28	−8.20	0.11
P4	−13.30	−8.57	0.10	−10.81	−6.27	0.10	−13.27	−8.12	0.10
P5	−13.31	−8.58	0.10	−11.14	−6.26	0.10	−13.30	−8.18	0.10
P6	−13.31	−8.56	0.10	−11.05	−6.20	0.10	−13.27	−8.04	0.10
P7	−13.28	−8.68	0.10	−11.00	−6.69	0.11	−13.21	−8.31	0.11
P8	−13.29	−8.52	0.10	−11.17	−6.23	0.10	−13.24	−8.12	0.10
P9	−13.28	−8.51	0.10	−10.91	−6.30	0.11	−13.19	−8.16	0.11

Table 5. Estimate of high-frequency vibration parameters of the distributed imaging scene.

Parameter	True Value	Chirplet-LSSE-SRT	FrFT-QML-RANSAC
$f_{1} / Hz$	18.300	18.300	18.230
$f_{2} / Hz$	35.000	35.020	34.980
$A_{1} / mm$	1.500	1.506	1.503
$A_{2} / mm$	1.000	0.998	0.961
$φ_{1} / rad$	2.618	2.616	2.583
$φ_{2} / rad$	2.618	2.591	2.571

Table 6. Image entropy and image contrast of the imaging results.

	Image Entropy	Image Contrast
Figure 15a	9.2485	9.1317
Figure 15b	9.7678	2.8002
Figure 15c	9.6374	4.9559
Figure 15d	9.3657	8.7684

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Hao, Z.; Sun, J.; Li, Q.; Shan, T. Estimation of High-Frequency Vibration Parameters for Airborne Terahertz SAR Using Chirplet Decomposition and LS Sequential Estimators. Remote Sens. 2022, 14, 3416. https://doi.org/10.3390/rs14143416

AMA Style

Hao Z, Sun J, Li Q, Shan T. Estimation of High-Frequency Vibration Parameters for Airborne Terahertz SAR Using Chirplet Decomposition and LS Sequential Estimators. Remote Sensing. 2022; 14(14):3416. https://doi.org/10.3390/rs14143416

Chicago/Turabian Style

Hao, Zhaoxin, Jinping Sun, Qing Li, and Tao Shan. 2022. "Estimation of High-Frequency Vibration Parameters for Airborne Terahertz SAR Using Chirplet Decomposition and LS Sequential Estimators" Remote Sensing 14, no. 14: 3416. https://doi.org/10.3390/rs14143416

APA Style

Hao, Z., Sun, J., Li, Q., & Shan, T. (2022). Estimation of High-Frequency Vibration Parameters for Airborne Terahertz SAR Using Chirplet Decomposition and LS Sequential Estimators. Remote Sensing, 14(14), 3416. https://doi.org/10.3390/rs14143416

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu

Estimation of High-Frequency Vibration Parameters for Airborne Terahertz SAR Using Chirplet Decomposition and LS Sequential Estimators

Abstract

1. Introduction

2. THz-SAR High-Frequency Vibration Error Model and Effect

3. ICR Estimation with Chirplet Decomposition

4. Parameter Estimation with LS Sequential Estimators and SRT

5. The Processing Flow of the Proposed Parameters Estimation Method

6. Simulation Results and Discussion

6.1. Influence of High-Frequency Vibration on Imaging

6.2. Simulation Results of Point Targets

6.3. Simulation Results of Distributed Imaging Scene

7. Conclusions

Author Contributions

Funding

Data Availability Statement

Conflicts of Interest

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI