# Exploiting Chirp Rate Estimation Methods to Improve Image Formation Quality of Synthetic Aperture Radar

Authors: N. Rabiee; H. Azad; N. Parhizgar
Journal: International Journal of Aerospace Engineering
 (2021)

---

## Abstract

A common assumption in SAR image formation and processing algorithms is that the chirp rates of the transmitted and received radar signals are exactly the same. Dechirp processing is also done based on this common assumption. In real scenarios, the chirp rate of the received signal is different from that of the transmitted signal due to several reasons. In case the difference between the chirp rates of the transmitted and received signals is obvious, the demodulation and compression of the received pulse are not carried out precisely and defocusing the targets and the output images of the SAR processor results. In the present paper, a new technique is proposed to improve the image formation quality of SAR by exploiting chirp rate estimation methods. Based on the proposed technique, the chirp rate of the received signal is estimated, and then, dechirp is carried out by using a time-reversed complex conjugate filter constructed based on the estimated chirp rate. In this stage, the existing chirp rate estimation algorithms can be used. The quality of the output image is assessed using PSLR as a quantitative criterion and the average number of point target extension pixels along the azimuth direction. Simulation results indicated that the smaller the average number of point target extension pixels along with azimuth and the higher the PSLR average is, the better the output image quality would be. Therefore, output images obtained from the proposed method by exploiting chirp rate estimation algorithms would have a better quality with a higher PSLR average (14.1 and 13.6) and also the lower average number of point target extension pixels along the azimuth directions (2.1 and 4.9) than the common method with PSLR average (8.3) and an average number of point target extension pixels along the azimuth direction (7.1).


---

## Body

## 1. Introduction

Considering the principle restrictions of the optical imaging systems with respect to the climate, the acceptable resolution of the image, the limitation of taking images during the daytime, and other challenges formed the idea of using the radar systems. Synthetic aperture radars (SARs) are extensively used for high-resolution space borne and airborne applications for monitoring the weather, climate change, and earth resource mapping. Phased array imaging radars employ long antennas to generate a fan beam that illuminates the ground below. Track resolution of these radars is determined by the beam width, while the across resolution is determined by the pulse length. Antenna dimensions and mass of such radars need to be confined, especially in airborne and space-borne systems. This can severely limit the antenna aperture of the radar and, hence, degrade its resolution. This limitation is nowadays circumvented using signal processing techniques in synthetic aperture radars (SARs). The SAR system is mounted on a moving platform such as an aircraft or a satellite to achieve high-resolution images. In fact, SAR radar uses an antenna with small real dimensions, which distributes waves along with the platform movement trajectory and then receives the return echoes from the respective area. After receiving the return echo from the ground surface, it starts a special process to take images from the area. Thus, the main characteristic of the hardware operation of SAR is its movement over the area and transmitting and receiving the waves, so that it can process the signals after receiving enough echoes. Thus, the small antenna of the radar operates similarly to a large one. The synthetic aperture radar is considered a high-resolution imaging radar. The SAR imaging resolution is acceptable and significant with respect to both range and azimuth [1–10]. Chirp signal is one of the most useful signals in the field of wireless communication; the frequency of which varies with time [11–16]. When in the specific state that the chirp signal frequency is linear with respect to time, it is called linear frequency modulation (LFM) [17]. Due to their widespread role in the frequency domain, LFM signals are of considerable use in radar, Sonar, ISAR, communications, ultrasound, geodesy, and many other areas. The imaging of moving targets using SAR is a specific application. The LFM is determined through two of its main parameters, i.e., central frequency and chirp rate. Many studies have been done for detecting and estimating the LFM parameter. Many algorithms have been proposed for estimating LFM signal parameters. One of its special and significant applications regarding the estimation of the chirp signal parameters in the radar processing algorithm is the synthetic aperture radar (SAR). For instance, when imaging the moving targets, finding the position and accurate focus on the moving target requires exact information regarding the chirp signal phase parameters, which is determined using two parameters of signal phase (chirp rate and central frequency) of the target’s movement [18–22]. Due to the SAR platform motion, the return echo is a chirp signal with an unknown chirp rate. In order to process and form the final image of SAR, the information about the return chirp rate is inevitable and necessary. The precise estimation of the chirp rate is of great importance for precisely forming the receiver’s matched filter and, finally, improving the focus and resolution of the final image. In airborne systems, especially drones, the platform motion compensation is carried out. However, the precision of the motion compensation is limited to the precision of navigation systems on the platform (GPS and INS). Most of the time, it results in a residual error. A common assumption in the proposed algorithms in the references for the SAR image formation and processing is that the chirp rates of the radar’s transmitted and received signals are precisely the same. Carrying out the dechirp processing is based on this common assumption. While in the real scenarios, due to different reasons such as vibrations of the platform or lack of precise calibration of the systems, the chip rate of the received signal will be different from that of the transmitted signals. Therefore, in case the difference between the chirp rates of the transmitted and received signals is obvious, the demodulation and compression of the received pulse are not carried out precisely; thus, the targets and output images of the SAR processor are defocused. The aim of this paper is to enhance the image formation quality of SAR by exploiting chirp rate estimation methods. As shown, by using desirable chirp rate estimators, it is possible to estimate the chirp rate more accurately and, as a result, to achieve better focus of the SAR image. This work describes a novel technique to improve the image formation quality of SAR. The chirp rate of the received signal is estimated, and then, dechirp is carried out by using time-reversed complex conjugate filter constructed based on the estimated chirp rate. In this stage, the existing chirp rate estimation algorithms can be employed. Afterwards, the output image quality is assessed using the PSLR as the quantitative criterion and lower average number of point target extension pixels along the azimuth direction which is compared to the common method. The paper is organized as follows. Section 2 discusses the process of SAR, and in Section 3, two different algorithms for estimating the chirp rate are described. The algorithm is expressed in Section 4, and the proposed technique is expressed in Section 5. The results of the simulation are provided in Section 6. The conclusion is explained in Section 7.


## 2. SAR Processing

The image resolution is considered an important parameter among the applications of imaging. The resolution in the azimuth directionρa is equal to the radar aperture Ls, i.e., multiplication of the beam width β and the distance to a target R. The beam width is determined by the ratio of the wavelength λ and the antenna length D. Therefore, the resolution in the azimuth direction ρa increases with shorter wavelength and bigger antenna size [5]:
(1)ρa=Ls=λDR,where λ is the wavelength, D is the radar antenna length in flight direction, and R is the distance of the target away from radar. Geometric expression of the resolution in the azimuth direction and the beam width in R distance are displayed in Figure 1.Figure 1
Geometric expression of the resolution in the azimuth direction [5].The transmitted signal by the SAR radar is always a LFM pulse or chirp, which is stated as follows:(2)St=recttτp×ej2πf0t+1/2Kt2.In the above relation,rect signifies a rectangular window function; the value of which on the span of τp equals one and equals zero on any other position. f0 parameter indicates the radar carrier frequency which is modulated by the frequency rate of K (LFM modulation). The return echo signal from the target equals the following:
(3)Set=rectt−ττp×ej2πf0t−τ×ej2πKt−τ2,where τ is the delay of the transmit and return signal of the radar and depends on the distance of the radar from the target.
(4)τ=2rc≈2Rc1+x22R2,where x is the location of the radar in the movement trajectory. In fact, the radar transmits a train of the LFM modulated pulse with the pulse repetition interval PRI=1/PRF which PRF is the pulse repetition frequency. Therefore, the echo signal is stated as a two-dimensional signal as follows:
(5)Sex,t=rectt−2R/cτp×ej2πf0t−4πR/λ×e−j2π/λRx2×ejπKt−2R/c2.The received echo train demodulation that is displayed bySex,t will be as follows:
(6)Sex,t:R=recttτp×e−j4πR/λ×ejπKt−2R/c2×e−j2π/λRx2.Equation (6) signifies the response (echo) of SAR radar from the point target with R parameter as the shortest distance between the radar and the target [5]. The following is the impulse response of SAR:
(7)hx,t;R=e−j4πR/λ×τpe−jπKt2sinπt/ρrπt/ρr×Lsej2π/λRx2sinπx/ρaπx/ρa,where rt⟶ is the distance vector and is defined as follows:
(8)rt⟶=r0⟶+V⟶t+12A⟶t2,where V⟶ and A⟶ are the vectors of velocity and acceleration of the platform, respectively. The following is the distance vector domain:
(9)rt⟶≈r0⟶+r0⟶×V⟶r0⟶t+r0⟶×A⟶+V⟶22r0⟶t2.The Doppler frequency created by the relative movement between the radar and the target equals the following:(10)fdt=12π×dφtdt=−2λ×r0⟶×V⟶r0⟶−2λ×r0⟶×A⟶+V⟶22r0⟶t=fdc+fdrt.Equation (10) indicates that the Doppler signal is, in fact, a frequency modulated pulse, in which fdc is the Doppler centroid frequency and fdr is the Doppler frequency rate (ramp of frequency variation with respect to time). In accordance with the theory of linear systems, SAR is a linear system comprising two match filters of range and azimuth. The SAR input is the received (echo) data by the radar, and its output is the processed image. RDA has a better quality than other algorithms among the SAR processing algorithms. Using RCMC is one of the distinguished properties of this algorithm. In this algorithm, the received energy from the point targets with the equal range, which are placed in the separate azimuth positions, is transferred to similar positions in the azimuth frequency domain. The steps of RDA implementation are displayed in Figure 2 as blocks.Figure 2
The block diagram for implementing RDA [1].


## 3. Different Chirp Rate Estimation Algorithms

Two different algorithms employed in the present paper for chirp rate estimation, suggested in References [23, 24], are explained in this section.


### 3.1. Chirp Rate and Instantaneous Frequency Estimators for Frequency Modulation Signals

In Reference [23], the new chirp rate and instantaneous frequency estimators for frequency modulation signals were introduced, which were used for designing new recursive versions of vertically synchrosqueezed short-time Fourier transform (STFT) using the method explained before [25]. As discussed by Fourer et al. [23], it is assumed that the locally analysed signal is approximated with a Gaussian modulated linear chirp.xt=Axtejϕxt, with Axt=Axe−t−tx2/2Tx2, ϕxt=φx+ωxt+αxt2/2, andj2=−1.According to the definition, derivation ofxt can be stated as follows:
(11)dxdtt=ddtlnAxt+jdϕxdttxt,(12)dxdtt=qxt+pxxt.Withqx=−1/Tx2+jαx and px=tx/Tx2+jωx. The term qxt+px=−t−tx/Tx2+jαxt+ωx is often called instantaneous complex frequency [26, 27], and its imaginary part is the instantaneous frequency of the signal. Differentiable analysis window ht is considered to define the short-time Fourier transform (STFT) for xt with Fxht,ω. In this case, the result will be as follows:
(13)Fxht,ω=∫Rxuht−u∗e−jωudu=e−jωt∫Rxt−uhu∗ejωudu.The partial derivative of this term with respect to time can be written as follows:(14)∂Fxh∂tt,ω=∫Rxudhdtt−u∗e−jωudu=−jωFxht,ω+e−jωt∫Rdxdtt−uhu∗ejωudu.Ifqxt−u+pxxt−u is placed, instead of dx/dtt−u, then,
(15)FxDht,ω=−qxFxτht,ω+qxt+px−jωFxht,ω.If we differentiate again with respect tot, then, the result will be as follows:
(16)FxD2ht,ω=−qxFxτDht,ω+qxt+px−jωFxDht,ω.And in general, forn≥1 is as follows:
(17)∂nFxh∂tnt,w=FxDnht,ω=−qxFxτDn−1ht,ω+qxt+px−jωFxDn−1ht,ω.Assuming that amplitude is constant, frompx=jωx, qx=jαx, and Equation (15), we can conclude that
(18)FxDht,ω=−jαxFxτht,ω+jαxt+ωx−ωFxht,ω.By multiplying byFxht,ω∗ and considering the real part, we can obtain the chirp rate estimator on the basis of the first derivative of STFT (in case of nonzero denominator).
(19)α^xK1t,ω=ReFxDht,ωFxht,ω∗ImFxτht,ωFxht,ω∗.This equation that is extracted from the analytical results obtained in the special problem of the Gaussian window in [28, 29] provides a chirp rate estimator for each analytical window of h (unbiased in case T2x⟶∞). According to Equation (17), all types of chirp rate estimators can be derived on the same principle:
(20)α^xKnt,ω=ReFxDnht,ωFxDn−1ht,ω∗ImFxτDn−1ht,ωFxDn−1ht,ω∗.Regarding the estimators using time derivatives, if the amplitude is not assumed to be constant as before, Equations (15) (16) can be considered as one system of two linear equations with two variables of qxt+px−jω and qx. The response of this system (in case of a nonzero denominator) obtains the qx estimator and imaginary part of a new chirp rate estimator:
(21)q^xt2t,ω=FxD2ht,ωFxht,ω−FxDht,ω2FxDht,ωFxτht,ω−FxτDht,ωFxht,ω,α^xt2t,ω=Imq^xt2t,ω.In general, if a similar process is operated on Equations (15) and (17), all types of chirp rate estimators can be achieved:
(22)q^xtnt,ω=FxDnht,ωFxht,ω−FxDn−1ht,ωFxDht,ω2FxDn−1ht,ωFxτht,ω−FxτDn−1ht,ωFxht,ω,α^xtnt,ω=Imq^xtnt,ω.Considering that∂Fxh/∂ω=jFxτht,ω−tFxht,ω, then, by deriving Equation (15) with respect to ω, the results will be as follows:
(23)FxτDht,ω+Fxht,ω=−qxFxτ2ht,ω+qxt+px−jωFxτht,ω.In general, differentiating derivation of Equation (15) n−1 times (for n≥2) with respect to ω, then,
(24)Fxτn−1Dht,ω+n−1Fxτn−2ht,ω=−qxFxτnht,ω+qxt+px−jωFxτn−1ht,ω.If Equations (15) and (23) are considered as one system of linear equations, one qx estimator and one chirp rate estimator will be obtained:
(25)q^xω2=FxτDht,ωFxht,ω+Fxht,ω2−Fxτht,ωFxDht,ω2Fxτht,ω2−Fxτ2ht,ωFxht,ω,α^xω2=Imq^xω2t,ω.Generally, a whole class of chirp rate estimators can be obtained from Equations (15) and (24):
(26)q^xωn=Fxτn−1Dh+n−1Fxτn−2hFxh−Fxτn−1hFxDhFxτn−1hFxτh−FxτnhFxh,α^xωn=Imq^xωnt,ω.The classical spectroscopy assignment and synchrosqueezing use the time and frequency assignment operators in accordance with the definition [25, 30–33]:
(27)t~t,ω=t−Fxτht,ωFxht,ωwitht^t,ω=Ret~t,ω,(28)ω^t,ω=Imω~t,ωwithω~t,ω=jω+FxDht,ωFxht,ω.By using Equations (15), (27), and (28), the results will be as follows:
(29)ω^t,ω=αxt^t,ω+ωx+1Tx2ImFxτht,ωFxht,ω.When the frequency amplitude is constant, if we multiplyFxτht,ω∗ by Equation (18) and consider its real part, then,
(30)ϕ̇^xK1t,ω=ω−ReFxDht,ωFxτht,ω∗ImFxht,ωFxτht,ω∗.In comparison to Equation (28), ϕ̇^xK1t,ω is an instantaneous frequency estimator at time of t, while ω^t,ω is an instantaneous frequency estimator at time t^t,ω,ϕ̇^xt^t,ω,ω. With the same principle, all the instantaneous frequency estimators can be obtained from Equation (15):
(31)ϕ̇^xKnt,ω=ω−ReFxDnht,ωFxτDn−1ht,ω∗ImFxDn−1ht,ωFxτDn−1ht,ω∗.By combining Equations (15) and (16), the qxt+px−jω estimator is achieved; the imaginary part of which results in an instantaneous frequency estimator. It can be generalized to any order from Equations (15) and (24) (in case of the nonzero denominator).
(32)ϕ̇^xt2t,ω=ω+ImFxD2hFxτh−FxτDhFxDhFxDhFxτh−FxτDhFxh,ϕ̇^xtnt,ω=ω+ImFxDnhFxτh−FxτDn−1hFxDhFxDn−1hFxτh−FxτDn−1hFxh.By combining Equations (15) and (24), another class of instantaneous frequency estimator is obtained as follows:
(33)ϕ̇^xωnt,ω=ω+ImFxτn−1Dh+n−1Fxτn−2hFxτh−FxτnhFxDhFxτn−1hFxτh−FxτnhFxh.The results stated that regarding the unbiased instantaneous frequency estimators can be employed for estimating the instantaneous frequency of a signal component, which is located in the vicinity of a TF pointt,ω. Furthermore, it can be used for extracting an improved synchrosqueezing process, called vertical synchrosqueezing. This process operates based on a general signal reconstruction formula [25]:
(34)xt−t0=1ht0∗∫RFxht,ωejωt−t0dω2π.For eacht0 where ht0≠0, the vertically synchrosqueezed STFT can be defined as follows [34, 35].
(35)VSFxht,ω=∫RFxht,ω′ejω′t−t0δω−ϕ̇^xt,ω′dω′.In whichϕ̇^xt,ω is one of the proposed instantaneous frequency estimators. This expression is totally different from the classical synchrosqueezed STFT. Since in that expression ϕ̇^xt,ω′ is used instead of ω^t,ω′, its squared modulus results in a precise TFR which can be inverted as follows:
(36)x^t−t0=1ht0∗∫RVSFxht,ωdω2π.In which the intervals of integration can be limited to the vicinity of a signal ridge to carry out mode extraction [36]. In order to assess their sensitivity in relation to the model inefficiencies, the log amplitude and phase are assumed to be third-order polynomial [23]:
(37)lnAxt=lnAx−t−tx22Tx2−ΔAt−tx36Tx2,ϕxt=φx+ωxt+αxt22+Δϕt36.Thus, Equations (12) and (15) will be transformed intodx/dtt=rxt2+qxt+pxxt, with rx=−ΔA/2Tx2+jΔϕ/2, qx=−1−ΔAtx/Tx2+jαx,px=−2tx−ΔAtx2/2Tx2+jωx, and
(38)FxDht,ω=rxFxτ2ht,ω−2rxt+qxFxτht,ω+rxt2+qxt+px−jωFxht,ω.Differentiating Equation (38) with respect to t and ω results as follows:
(39)FxD2h=rxFxτ2Dh−2rxt+qxFxτDh+rxt2+qxt+px−jωFxDh,(40)FxτDh+Fxh=rxFxτ3h−2rxt+qxFxτ2h+rxt2+qxt+px−jωFxτh.Equations (38), (39), and (40) can be considered as a set of three linear equations, with three variablesrx,2rxt+qx, and rxt2+qxt+px−jω. In case of a nonzero denominator, the imaginary parts of the response of this system will result in robust estimators from angular jerk Δϕ for the instantaneous chirp rate and instantaneous frequency.


### 3.2. Efficient Algorithm for Estimating the Parameters of a Chirp Signal

In [24], an efficient algorithm was suggested for estimating different parameters of a chirp signal in the presence of the stationary noise. The main advantage of this algorithm is that by starting from the reasonable initial values for the frequency and frequency rate, only in four stages, estimations will be generated that are asymptotically equal to LSE. As discussed by Lahiri et al. [24], the chirp signal in the presence of additive noise is modelled as follows:
(41)yn=A0cosα0n+β0n2+B0sinα0n+β0n2+Xn,n=1,⋯,N,where A0 and B0 are the nonzero values with restriction A02+B02≤M for constant M. α0 and β0 signify the frequency and frequency rate, respectively, and their values range between zero and π. Xn is a stationary noise sequence, which is as follows:
(42)Xn=∑j=−∞∞ajεn−j,∑j=−∞∞aj≺∞,where εn is a sequence of independent and identically distributed (i.i.d) random variables with zero mean and variance σ2. LESs of unknown parameters of Equation (41) can be obtained by minimizing SΘ with respect to variable Θ=A,B,α,β.
(43)SΘ=∑n=1Nyn−Acosαn+βn2−Bsinαn+βn22=Y−Wα,βABTY−Wα,βAB,where Y=y1,⋯,yNT is the N×1 data vector and Wθ is the matrix N×2, as demonstrated as follows:
(44)Wθ=cosα+βsinα+βcos2α+4βsin2α+4β⋮⋮cosNα+N2βsinNα+N2β.If the values ofα and β are known, LSEs of A0 and B0 can be obtained as A^θ and B^θ, in which θ=α,β and
(45)A∧θ,B∧θT=WTθWθ−1WTθY.Thus, LSEs of α0 and β0 can be obtained by minimizing the value of Qα,β with respect to α and β, in which
(46)Qα,β=SA^θ,B^θ,α,β=YTWθWTθWθ−1WTθY.After determining the LSEs ofα0 and β0 called α^ andβ^, the value of the LSEs of A0 and B0 can be simply achieved as A^α^,β^ andB^α^,β^, respectively. Kundu and Nandi [37] obtained the properties of LSEs to be as follows. LSEs of the unknown parameters of Equation (41) with the compatible corresponding parameters and asymptotic distribution are as follows:
(47)NA^−A01/2,NB^−B01/2,N3/2α^−α0,N5/2β^−β0⟶dN40,4cσ2∑,where c=∑j=−∞∞aj2 and ⟶d symbol indicate the distribution convergence, N40,4cσ2∑ indicates a 4-variable normal distribution with zero mean vector, the dispersion matrix is 4cσ2∑, and
(48)∑=1A02+B0212A02+9B02−4A0B018B0−15B0−4A0B0129A02+B02−18A015A018B0−18A096−90−15B015A0−9090.Ifα~ is an α0 estimator, such that α~−α0=opN−1−δ1 for some δ1≻0, and β~ is a β0 estimator, such that β~−β0=opN−2−δ2 for some δ2≻0. Hence, an improved estimator can be obtained as follows:
(49)α~~=α~+48N2ImPNαQN,with
(50)PNα=∑n=1Nynn−N2e−iα~n+β~n2,QN=∑n=1Nyne−iα~n+β~n2.Accordingly, an improved estimator ofβ0 can be achieved as follows:
(51)β~~=β~+45N4ImPNβQN,with
(52)PNβ=∑n=1Nynn2−N23e−iα~n+β~n2,and QN is defined in the same as Equation (50). The following theorems provide a justification for improved estimators.Theorem 1.
Ifα~−α0=opN−1−δ1 for δ1≻0, then,
(a)
α~~−α0=opN−1−2δ1 if δ1≤1/4(b)
N3/2α~~−α0⟶dN0,σ12 if δ1>1/4
In whichσ12=384cσ12/A02+B02 is the asymptotic variance of the LSE of α0.Theorem 2.
Ifβ~−β0=opN−2−δ2 forδ2≻0, then,
(a)
β~~−β0=opN−2−2δ2 if δ2≤1/4(b)
N5/2β~~−β0⟶dN0,σ22 if δ2≻1/4
In whichσ22=360cσ12/A02+B02 is the asymptotic variance of the LSE of β0.Theorem 3.
Ifα~−α0=opN−1−δ1 and β~−β0=opN−2−δ2, then, for δ1≻1/4 and δ2≻1/4, there is CovN3/2α~~−α0,N5/2β~~−β0⟶360cσ12/A02+B02.The efficient estimators can be obtained using the above method and by starting from the initial guessesα~ and β~ by the convergence rate α~−α0=opN−1 and β~−β0=opN−2, respectively. In accordance with Rice and Rosenblatt in Reference [38], it is concluded that if it can be demonstrated that the Nα~−α0⟶0 and N2β~−β0⟶0 are established for α~ and β~, and Qα,β minimizer, then, when searching over the grids πj/N,πk/N2,j=1,⋯,N and k=1,⋯,N2, there will be such an initial guess. If Θ0=A0,B0,α0,β0 is the value of real parameter and D is a diagonal matrix, such that D=1/N,1/N,1/NN,1/N2N. Note that (see [37])
(53)−S′Θ0D=Θ~−Θ0D−1DS″Θ¯D,where Θ¯ defines a point on line joining Θ~ and Θ0.
(54)DS″Θ¯D⟶limN⟶∞DS″Θ0D.Now, Equation (53) gives the following:
(55)Θ~−Θ0ND−1=−1NS′Θ0DDS″Θ¯D−1.Besides,−S′Θ0D/N⟶0 manifests that Θ~−Θ0ND−1⟶0, i.e., Nα~−α0⟶0 and Θ~−Θ0ND−1⟶0.The algorithm can be described as follows. The estimation obtained forα0 and β0 at the iteration of i-th is α~i and β~i, respectively.


## 3.1. Chirp Rate and Instantaneous Frequency Estimators for Frequency Modulation Signals

In Reference [23], the new chirp rate and instantaneous frequency estimators for frequency modulation signals were introduced, which were used for designing new recursive versions of vertically synchrosqueezed short-time Fourier transform (STFT) using the method explained before [25]. As discussed by Fourer et al. [23], it is assumed that the locally analysed signal is approximated with a Gaussian modulated linear chirp.xt=Axtejϕxt, with Axt=Axe−t−tx2/2Tx2, ϕxt=φx+ωxt+αxt2/2, andj2=−1.According to the definition, derivation ofxt can be stated as follows:
(11)dxdtt=ddtlnAxt+jdϕxdttxt,(12)dxdtt=qxt+pxxt.Withqx=−1/Tx2+jαx and px=tx/Tx2+jωx. The term qxt+px=−t−tx/Tx2+jαxt+ωx is often called instantaneous complex frequency [26, 27], and its imaginary part is the instantaneous frequency of the signal. Differentiable analysis window ht is considered to define the short-time Fourier transform (STFT) for xt with Fxht,ω. In this case, the result will be as follows:
(13)Fxht,ω=∫Rxuht−u∗e−jωudu=e−jωt∫Rxt−uhu∗ejωudu.The partial derivative of this term with respect to time can be written as follows:(14)∂Fxh∂tt,ω=∫Rxudhdtt−u∗e−jωudu=−jωFxht,ω+e−jωt∫Rdxdtt−uhu∗ejωudu.Ifqxt−u+pxxt−u is placed, instead of dx/dtt−u, then,
(15)FxDht,ω=−qxFxτht,ω+qxt+px−jωFxht,ω.If we differentiate again with respect tot, then, the result will be as follows:
(16)FxD2ht,ω=−qxFxτDht,ω+qxt+px−jωFxDht,ω.And in general, forn≥1 is as follows:
(17)∂nFxh∂tnt,w=FxDnht,ω=−qxFxτDn−1ht,ω+qxt+px−jωFxDn−1ht,ω.Assuming that amplitude is constant, frompx=jωx, qx=jαx, and Equation (15), we can conclude that
(18)FxDht,ω=−jαxFxτht,ω+jαxt+ωx−ωFxht,ω.By multiplying byFxht,ω∗ and considering the real part, we can obtain the chirp rate estimator on the basis of the first derivative of STFT (in case of nonzero denominator).
(19)α^xK1t,ω=ReFxDht,ωFxht,ω∗ImFxτht,ωFxht,ω∗.This equation that is extracted from the analytical results obtained in the special problem of the Gaussian window in [28, 29] provides a chirp rate estimator for each analytical window of h (unbiased in case T2x⟶∞). According to Equation (17), all types of chirp rate estimators can be derived on the same principle:
(20)α^xKnt,ω=ReFxDnht,ωFxDn−1ht,ω∗ImFxτDn−1ht,ωFxDn−1ht,ω∗.Regarding the estimators using time derivatives, if the amplitude is not assumed to be constant as before, Equations (15) (16) can be considered as one system of two linear equations with two variables of qxt+px−jω and qx. The response of this system (in case of a nonzero denominator) obtains the qx estimator and imaginary part of a new chirp rate estimator:
(21)q^xt2t,ω=FxD2ht,ωFxht,ω−FxDht,ω2FxDht,ωFxτht,ω−FxτDht,ωFxht,ω,α^xt2t,ω=Imq^xt2t,ω.In general, if a similar process is operated on Equations (15) and (17), all types of chirp rate estimators can be achieved:
(22)q^xtnt,ω=FxDnht,ωFxht,ω−FxDn−1ht,ωFxDht,ω2FxDn−1ht,ωFxτht,ω−FxτDn−1ht,ωFxht,ω,α^xtnt,ω=Imq^xtnt,ω.Considering that∂Fxh/∂ω=jFxτht,ω−tFxht,ω, then, by deriving Equation (15) with respect to ω, the results will be as follows:
(23)FxτDht,ω+Fxht,ω=−qxFxτ2ht,ω+qxt+px−jωFxτht,ω.In general, differentiating derivation of Equation (15) n−1 times (for n≥2) with respect to ω, then,
(24)Fxτn−1Dht,ω+n−1Fxτn−2ht,ω=−qxFxτnht,ω+qxt+px−jωFxτn−1ht,ω.If Equations (15) and (23) are considered as one system of linear equations, one qx estimator and one chirp rate estimator will be obtained:
(25)q^xω2=FxτDht,ωFxht,ω+Fxht,ω2−Fxτht,ωFxDht,ω2Fxτht,ω2−Fxτ2ht,ωFxht,ω,α^xω2=Imq^xω2t,ω.Generally, a whole class of chirp rate estimators can be obtained from Equations (15) and (24):
(26)q^xωn=Fxτn−1Dh+n−1Fxτn−2hFxh−Fxτn−1hFxDhFxτn−1hFxτh−FxτnhFxh,α^xωn=Imq^xωnt,ω.The classical spectroscopy assignment and synchrosqueezing use the time and frequency assignment operators in accordance with the definition [25, 30–33]:
(27)t~t,ω=t−Fxτht,ωFxht,ωwitht^t,ω=Ret~t,ω,(28)ω^t,ω=Imω~t,ωwithω~t,ω=jω+FxDht,ωFxht,ω.By using Equations (15), (27), and (28), the results will be as follows:
(29)ω^t,ω=αxt^t,ω+ωx+1Tx2ImFxτht,ωFxht,ω.When the frequency amplitude is constant, if we multiplyFxτht,ω∗ by Equation (18) and consider its real part, then,
(30)ϕ̇^xK1t,ω=ω−ReFxDht,ωFxτht,ω∗ImFxht,ωFxτht,ω∗.In comparison to Equation (28), ϕ̇^xK1t,ω is an instantaneous frequency estimator at time of t, while ω^t,ω is an instantaneous frequency estimator at time t^t,ω,ϕ̇^xt^t,ω,ω. With the same principle, all the instantaneous frequency estimators can be obtained from Equation (15):
(31)ϕ̇^xKnt,ω=ω−ReFxDnht,ωFxτDn−1ht,ω∗ImFxDn−1ht,ωFxτDn−1ht,ω∗.By combining Equations (15) and (16), the qxt+px−jω estimator is achieved; the imaginary part of which results in an instantaneous frequency estimator. It can be generalized to any order from Equations (15) and (24) (in case of the nonzero denominator).
(32)ϕ̇^xt2t,ω=ω+ImFxD2hFxτh−FxτDhFxDhFxDhFxτh−FxτDhFxh,ϕ̇^xtnt,ω=ω+ImFxDnhFxτh−FxτDn−1hFxDhFxDn−1hFxτh−FxτDn−1hFxh.By combining Equations (15) and (24), another class of instantaneous frequency estimator is obtained as follows:
(33)ϕ̇^xωnt,ω=ω+ImFxτn−1Dh+n−1Fxτn−2hFxτh−FxτnhFxDhFxτn−1hFxτh−FxτnhFxh.The results stated that regarding the unbiased instantaneous frequency estimators can be employed for estimating the instantaneous frequency of a signal component, which is located in the vicinity of a TF pointt,ω. Furthermore, it can be used for extracting an improved synchrosqueezing process, called vertical synchrosqueezing. This process operates based on a general signal reconstruction formula [25]:
(34)xt−t0=1ht0∗∫RFxht,ωejωt−t0dω2π.For eacht0 where ht0≠0, the vertically synchrosqueezed STFT can be defined as follows [34, 35].
(35)VSFxht,ω=∫RFxht,ω′ejω′t−t0δω−ϕ̇^xt,ω′dω′.In whichϕ̇^xt,ω is one of the proposed instantaneous frequency estimators. This expression is totally different from the classical synchrosqueezed STFT. Since in that expression ϕ̇^xt,ω′ is used instead of ω^t,ω′, its squared modulus results in a precise TFR which can be inverted as follows:
(36)x^t−t0=1ht0∗∫RVSFxht,ωdω2π.In which the intervals of integration can be limited to the vicinity of a signal ridge to carry out mode extraction [36]. In order to assess their sensitivity in relation to the model inefficiencies, the log amplitude and phase are assumed to be third-order polynomial [23]:
(37)lnAxt=lnAx−t−tx22Tx2−ΔAt−tx36Tx2,ϕxt=φx+ωxt+αxt22+Δϕt36.Thus, Equations (12) and (15) will be transformed intodx/dtt=rxt2+qxt+pxxt, with rx=−ΔA/2Tx2+jΔϕ/2, qx=−1−ΔAtx/Tx2+jαx,px=−2tx−ΔAtx2/2Tx2+jωx, and
(38)FxDht,ω=rxFxτ2ht,ω−2rxt+qxFxτht,ω+rxt2+qxt+px−jωFxht,ω.Differentiating Equation (38) with respect to t and ω results as follows:
(39)FxD2h=rxFxτ2Dh−2rxt+qxFxτDh+rxt2+qxt+px−jωFxDh,(40)FxτDh+Fxh=rxFxτ3h−2rxt+qxFxτ2h+rxt2+qxt+px−jωFxτh.Equations (38), (39), and (40) can be considered as a set of three linear equations, with three variablesrx,2rxt+qx, and rxt2+qxt+px−jω. In case of a nonzero denominator, the imaginary parts of the response of this system will result in robust estimators from angular jerk Δϕ for the instantaneous chirp rate and instantaneous frequency.


## 3.2. Efficient Algorithm for Estimating the Parameters of a Chirp Signal

In [24], an efficient algorithm was suggested for estimating different parameters of a chirp signal in the presence of the stationary noise. The main advantage of this algorithm is that by starting from the reasonable initial values for the frequency and frequency rate, only in four stages, estimations will be generated that are asymptotically equal to LSE. As discussed by Lahiri et al. [24], the chirp signal in the presence of additive noise is modelled as follows:
(41)yn=A0cosα0n+β0n2+B0sinα0n+β0n2+Xn,n=1,⋯,N,where A0 and B0 are the nonzero values with restriction A02+B02≤M for constant M. α0 and β0 signify the frequency and frequency rate, respectively, and their values range between zero and π. Xn is a stationary noise sequence, which is as follows:
(42)Xn=∑j=−∞∞ajεn−j,∑j=−∞∞aj≺∞,where εn is a sequence of independent and identically distributed (i.i.d) random variables with zero mean and variance σ2. LESs of unknown parameters of Equation (41) can be obtained by minimizing SΘ with respect to variable Θ=A,B,α,β.
(43)SΘ=∑n=1Nyn−Acosαn+βn2−Bsinαn+βn22=Y−Wα,βABTY−Wα,βAB,where Y=y1,⋯,yNT is the N×1 data vector and Wθ is the matrix N×2, as demonstrated as follows:
(44)Wθ=cosα+βsinα+βcos2α+4βsin2α+4β⋮⋮cosNα+N2βsinNα+N2β.If the values ofα and β are known, LSEs of A0 and B0 can be obtained as A^θ and B^θ, in which θ=α,β and
(45)A∧θ,B∧θT=WTθWθ−1WTθY.Thus, LSEs of α0 and β0 can be obtained by minimizing the value of Qα,β with respect to α and β, in which
(46)Qα,β=SA^θ,B^θ,α,β=YTWθWTθWθ−1WTθY.After determining the LSEs ofα0 and β0 called α^ andβ^, the value of the LSEs of A0 and B0 can be simply achieved as A^α^,β^ andB^α^,β^, respectively. Kundu and Nandi [37] obtained the properties of LSEs to be as follows. LSEs of the unknown parameters of Equation (41) with the compatible corresponding parameters and asymptotic distribution are as follows:
(47)NA^−A01/2,NB^−B01/2,N3/2α^−α0,N5/2β^−β0⟶dN40,4cσ2∑,where c=∑j=−∞∞aj2 and ⟶d symbol indicate the distribution convergence, N40,4cσ2∑ indicates a 4-variable normal distribution with zero mean vector, the dispersion matrix is 4cσ2∑, and
(48)∑=1A02+B0212A02+9B02−4A0B018B0−15B0−4A0B0129A02+B02−18A015A018B0−18A096−90−15B015A0−9090.Ifα~ is an α0 estimator, such that α~−α0=opN−1−δ1 for some δ1≻0, and β~ is a β0 estimator, such that β~−β0=opN−2−δ2 for some δ2≻0. Hence, an improved estimator can be obtained as follows:
(49)α~~=α~+48N2ImPNαQN,with
(50)PNα=∑n=1Nynn−N2e−iα~n+β~n2,QN=∑n=1Nyne−iα~n+β~n2.Accordingly, an improved estimator ofβ0 can be achieved as follows:
(51)β~~=β~+45N4ImPNβQN,with
(52)PNβ=∑n=1Nynn2−N23e−iα~n+β~n2,and QN is defined in the same as Equation (50). The following theorems provide a justification for improved estimators.Theorem 1.
Ifα~−α0=opN−1−δ1 for δ1≻0, then,
(a)
α~~−α0=opN−1−2δ1 if δ1≤1/4(b)
N3/2α~~−α0⟶dN0,σ12 if δ1>1/4
In whichσ12=384cσ12/A02+B02 is the asymptotic variance of the LSE of α0.Theorem 2.
Ifβ~−β0=opN−2−δ2 forδ2≻0, then,
(a)
β~~−β0=opN−2−2δ2 if δ2≤1/4(b)
N5/2β~~−β0⟶dN0,σ22 if δ2≻1/4
In whichσ22=360cσ12/A02+B02 is the asymptotic variance of the LSE of β0.Theorem 3.
Ifα~−α0=opN−1−δ1 and β~−β0=opN−2−δ2, then, for δ1≻1/4 and δ2≻1/4, there is CovN3/2α~~−α0,N5/2β~~−β0⟶360cσ12/A02+B02.The efficient estimators can be obtained using the above method and by starting from the initial guessesα~ and β~ by the convergence rate α~−α0=opN−1 and β~−β0=opN−2, respectively. In accordance with Rice and Rosenblatt in Reference [38], it is concluded that if it can be demonstrated that the Nα~−α0⟶0 and N2β~−β0⟶0 are established for α~ and β~, and Qα,β minimizer, then, when searching over the grids πj/N,πk/N2,j=1,⋯,N and k=1,⋯,N2, there will be such an initial guess. If Θ0=A0,B0,α0,β0 is the value of real parameter and D is a diagonal matrix, such that D=1/N,1/N,1/NN,1/N2N. Note that (see [37])
(53)−S′Θ0D=Θ~−Θ0D−1DS″Θ¯D,where Θ¯ defines a point on line joining Θ~ and Θ0.
(54)DS″Θ¯D⟶limN⟶∞DS″Θ0D.Now, Equation (53) gives the following:
(55)Θ~−Θ0ND−1=−1NS′Θ0DDS″Θ¯D−1.Besides,−S′Θ0D/N⟶0 manifests that Θ~−Θ0ND−1⟶0, i.e., Nα~−α0⟶0 and Θ~−Θ0ND−1⟶0.The algorithm can be described as follows. The estimation obtained forα0 and β0 at the iteration of i-th is α~i and β~i, respectively.


## 4. Algorithm

(1)
First step: SelectN1=N8/9. Consequently, α~0−α0=OPN−1=OPN1−1−1/8 and β~0−β0=OPN−2=OPN1−2−1/4.Perform Equation (49) and Equation (51). Thus, after the first iteration, α~1−α0=OPN1−1−1/4=OPN−10/9 and β~1−β0=OPN1−2−1/2=OPN−20/9.
(2)
Second step: SelectN2=N80/81. Consequently, α~1−α0=OPN2−1−1/8 and β~1−β0=OPN2−2−1/4.Perform Equation (49) and Equation (51). Thus, after the second iteration, α~2−α0=OPN2−1−1/4=OPN−100/81 and β~2−β0=OPN2−2−1/2=OPN−200/81.
(3)
Third step: SelectN3=N. Consequently, α~2−α0=OPN3−1−19/81 and β~2−β0=OPN3−2−38/81.Perform Equation (49) and Equation (51). Thus, after the third iteration, there is α~3−α0=OPN−1−38/81 and β~3−β0=OPN−2−76/81.
(4)
Fourth step: SelectN4=N and perform Equation (49) and Equation (51). Now, you can obtain the required convergence rate, i.e., α~4−α0=OPN−3/2 and β~4−β0=OPN−5/2.In fact, this iteration-based algorithm is for finding efficient estimators of frequency and frequency rate.


## 5. Proposed Technique

Due to its high time-bandwidth product, the LFM signal can simultaneously obtain acceptable resolution and the signal-to-clutter-plus-noise ratio (SCNR). Therefore, it is the most applicable signal in different types of SAR radars. In the SAR radar’s receiver, the time-reversed complex conjugate version from the transmitted pulse is multiplied by the received signal. Thus, it will demodulate and compress it. This processing is called dechirp on receive. By doing so, in the ideal state, the LFM modulation is demodulated and merely one sinusoidal signal is remained and its frequency corresponds to the delay of transmission and return of radar pulse. On the basis of this delay, the range of the target or the illuminated area can be estimated and employed in other processing. A common assumption in the proposed algorithms in the references for the SAR image formation and processing is that the chirp rates of the radar’s transmitted and received signals are precisely the same. Carrying out the dechirp processing is based on this common assumption. While in the real scenarios, due to different reasons such as vibrations of the platform (in the airborne systems, especially drones) or lack of precise calibration of the systems, the chip rate of the received signal will be different from the transmitted signals. The amount of this difference depends on the intensity of the above-said factors [39]. Therefore, in case the difference between the chirp rates of the transmitted and received signals is obvious, the demodulation and compression of the received pulse are not carried out precisely and defocusing the targets and the output images of the SAR processor is resulted (this problem is demonstrated in Simulation Result). For instance, in airborne systems, especially drones, platform motion compensation is carried out. However, the precision of the motion compensation is limited to the precision of navigation systems on the platform (GPS and INS); most of the time, it results in a residual error. In the present paper, at first, the received chirp rate is estimated, and then, dechirp is carried out by using time-reversed complex conjugate filter constructed based on the estimated chirp rate. In this stage, various existing chirp rate estimation algorithms can be employed and their operation quality can be assessed. Here, we used two types of chirp rate estimation algorithms proposed in [23, 24]. Finally, the quality of the output images of the proposed technique and the common technique was compared and assessed (using the chirp rate of the transmitted signal for dechirp) with the average PSLR as a quantitative criterion and the average number of point target extension pixels along with azimuth. As in Figure 3, it is processed to form the radar output image. So, it is only needed to add the chirp rate estimation to the block diagram of RDA and use an additional processing of chirp rate estimation. Simulation results of the aforementioned method will be presented in the next section.Figure 3
The block diagram of range-Doppler algorithm (RDA) containing the block proposed (the proposed plan).


## 6. Simulation Result

In the simulation section, SAR in the strip map imaging mode is transmitting pulses and receiving the return echoes from the area surface. The respective geometry for SAR system motion is shown in Figure4. In this simulation, 9-point targets are considered as seen in Figure 4. Besides, RDA [1] is employed for forming the image. It is worth noting that, in the stage of generating the received signal by the radar, the 3D vibrations will be applied to the platform motion. Then, in the processor, the unacceptable platform motions will be eliminated using MOCO methods provided in [40]. In Table 1, the values of the radar’s parameters are illustrated. In Figure 5(a), the simulation result of SAR image formation from the total area using the algorithm of [23] is displayed. For instance, if we zoom on any of the three targets in the middle along with azimuth, a magnified image of a part of the area is demonstrated in Figure 5, in which Figure 5(b) is the middle target, Figure 5(c) is the left target, and Figure 5(d) is the right target. Figure 6(a) indicates the simulation result of the SAR image formation by employing the algorithm in [24] from the total area. For instance, if we zoom on any of the three targets in the middle along with azimuth, in Figure 6, we can see the magnified image of a part of the area, in which Figure 6(b) is the middle target, Figure 6(c) is the left target, and Figure 6(d) is the right target.Figure 4
Respective geometry for motion of SAR system.Table 1
Parameters of the respective systems.
ParameterValueDescriptionh1e3 mHeight of platformθinc50 degreesIncident anglec3e8 m/sPropagation speedPRF500 HzPulse repetition frequencyfo9.6e9 HzCarrier frequencyλc/fomWavelengthτp3e−6 secChirp pulse durationBW180e6 HzBaseband bandwidthLe0.1 mAntenna length in elevationθaz5 degreesAntenna azimuth beam widthvp50Platform velocity (m/s)np1301Number of processed pulsesSNR10Signal-to-noise ratioFigure 5
(a) Image of the total area. Magnified image of a part of the area. (b) Middle target. (c) Left target. (d) Right target.
(a)(b)(c)(d)Figure 6
(a) Image of the total area. Magnified image of a part of the area. (b) Middle target. (c) Left target. (d) Right target.
(a)(b)(c)(d)In Figure7(a), the result of the SAR image formation simulation of the total area is displayed in the state of using the chirp rate of the transmitted signal in order to perform dechirp in the receiver (the common method without chirp estimation). For instance, if we zoom on any of the three targets in the middle along with azimuth, we can see the magnified image of a part of the area, in which Figure 7(b) is the middle target, Figure 7(c) is the left target, and Figure 7(d) is the right target.Figure 7
(a) Image of the total area. Magnified image of a part of the area. (b) Middle target. (c) Left target. (d) Right target.
(a)(b)(c)(d)In Table2, the average of PSLR obtained for 9 targets in the final images produced by using two different algorithms described in Section 3 of the paper is provided and compared to the average of PSLR using the common method (using the chirp rate of the transmitted signal in order to perform dechirp). Besides in Table 3, the average number of point target extension pixels along with azimuth in the images is generated using two different algorithms explained in Section 3 of the present paper and is provided for 9 targets and compared to the average number of point target extension pixels of the common method. The simulation results indicate that, with the lower average of point target extension pixels along the azimuth direction and the higher average of PSLR, the quality of the output image will be better.Table 2
Average of PSLR obtained for 9 targets in the final images produced by using different algorithms.
AlgorithmAverage PSLR (dB)Algorithm in [23]14.1Algorithm in [24]13.6Without chirp estimation(common method)8.3Table 3
Average number of point target extension pixels along with azimuth for 9 targets in the final images produced by using different algorithms.
AlgorithmAzimuth target extension (m)Algorithm in [23]2.1Algorithm in [24]4.9Without chirp estimation (common method)7.1According to Tables2 and 3, the proposed technique is shown to enhance the image formation quality of SAR by exploiting chirp rate estimation method with higher PSLR average (14.1 and 13.6) and lower average number of point target extension pixels along the azimuth direction (2.1 and 4.9) than the common method with PSLR average (8.3) and average number of point target extension pixels along the azimuth direction (7.1).


## 7. Conclusions

This paper proposed a novel technique for enhancing the image formation quality of SAR by exploiting chirp rate estimation method. Taking into account that carrying out the dechirp processing is on the basis of this common assumption in the algorithms of processing and formation of SAR image, the chirp rates of the radar’s transmitted and received signals are exactly the same. While in the real scenarios, due to different reasons, the chip rate of the received signal will be different from the transmitted signal. In case the difference between the chirp rates of the transmitted and received signals is obvious, defocusing the targets and the output images of the SAR processor is resulted. Based on the proposed technique, the received signal chirp rate was estimated. Afterwards, by using the time-reversed complex conjugate filter constructed based on the estimated chirp rate, the dechirp was carried out. In this process, different chip rate estimation algorithms could be used and, as a result, to achieve better focus of the SAR image. The simulation results indicated that the quality of the output image produced by the proposed technique by exploiting chirp rate estimation algorithms with higher PSLR average and lower average number of point target extension pixels along the azimuth direction was better than that of the common method (using the transmitted signal chip rate in order for dechirp).



---

*Source: 6636286-2021-09-28.xml*