## 3 Compensation Method  
### STFT-Based Rough Estimation  
First, we select the dominant point target from the coarsely focused image. The signal of interest (SoI) is obtained through performing inverse FFT to convert the coarsely focused image to the range-compressed phase history domain. To eliminate the linear phase in SoI, similar to PGA [29], the dominant point target is circularly shifted to the image center before performing inverse FFT. Thus, after accomplishing the above operations, we obtain the SoI, which is expressed as  
\[s(t)=\exp\left\{-j\frac{4\pi}{\lambda}\sum_{i=1}^{I}a_{i}\sin(2\pi f_{i}t+\phi _{i})\right\} \tag{11}\]  
The signal is sampled with the sampling time \(\Delta t=1/f_{PRF}\), where \(f_{PRF}\) is the pulse repetition frequency. The discretized SoI signal \(s(n)\) is represented as  
\[s(n)=s(n\cdot\Delta t),\ n\in[-N/2,N/2] \tag{12}\]  
where \(N\) represents the sample point number.  
On the one hand, the IF can be estimated by using the STFT. That is to say  
\[\hat{f}_{h}(n)=\arg\max_{f}\text{STFT}_{h}(n,f)| \tag{13}\]  
where \(\text{STFT}_{h}(n,f)\) is the STFT of the signal with window width \(h\).  
On the other hand, the IF of the signal can be acquired by taking the derivative of the phase in (11):  
\[f(t)=\frac{2}{\lambda}\sum_{i=1}^{I}2\pi f_{i}a_{i}\cos(2\pi f_{i}t+\phi_{i}) \tag{14}\]  
After discretization, the IF in (14) can be rewritten as  
\[f(n)=\sum_{i=1}^{I}[C_{0,i}\cos(2\pi f_{i}n\Delta t)+S_{0,i}\sin(2\pi f_{i}n \Delta t)] \tag{15}\]  
where  
\[\left\{\begin{array}{l}C_{0,i}=\frac{2}{\lambda}2\pi f_{i}a_{i}\cos\varphi_ {i}\\ S_{0,i}=\frac{2}{\lambda}2\pi f_{i}a_{i}\sin\varphi_{i}\end{array}\right. \tag{16}\]  
After estimating the IF \(\hat{f}_{h}(n)\) through (13), the high-frequency vibration parameters are obtained by solving the following linear equations, namely:  
\[\mathbf{\hat{f}}=\mathbf{Mx} \tag{17}\]  
where  
\[\left\{\begin{array}{l}\mathbf{\hat{f}}=\left[\hat{f}(-Q/2),\ \hat{f}(-Q/2+1),\cdots,\hat{f}(Q/2-1),\hat{f}(Q/2)\right]^{\text{T}}\\ \mathbf{x}=[C_{0,1},S_{0,1},C_{0,2},S_{0,2},\cdots C_{0,I},S_{0,I}]^{\text{ T}}\\ \mathbf{M}=[\mathbf{a}_{0,1},\mathbf{b}_{0,1},\mathbf{a}_{0,2},\mathbf{b}_{0, 2},\cdots\mathbf{a}_{0,I},\mathbf{b}_{0,I}]\\ \mathbf{a}_{0,i}=[cos(2\pi f_{i}(-Q/2)\Delta t),\cdots cos(2\pi f_{i}(Q/2) \Delta t)]^{\text{T}}\\ \mathbf{b}_{0,i}=[sin(2\pi f_{i}(-Q/2)\Delta t),\cdots sin(2\pi f_{i}(Q/2) \Delta t)]^{\text{T}}\end{array}\right. \tag{18}\]  
and \(Q+1\) is the number of azimuth instants.  
According to (14), the frequencies of the IF are the same as those of the vibration. Thus, through performing the spectral analysis on the estimated IF, the vibration frequencies are obtained. Based on the estimated vibration frequency, the LS solution of (17) is  
\[\mathbf{\hat{x}}=\left(\mathbf{M}^{\text{T}}\mathbf{M}\right)^{-1}\mathbf{M} ^{\text{T}}\mathbf{\hat{f}} \tag{19}\]Thus, the other high-frequency vibration parameters are obtained by  
\[a_{i}=\frac{\sqrt{\hat{\kappa}^{2}(2i-1)+\hat{\kappa}^{2}(2i)}}{\frac{\lambda}{2} 2\pi f_{i}} \tag{20}\]  
\[\phi_{i}=\arctan\left(\frac{\hat{\kappa}(2i)}{\hat{\kappa}(2i-1)}\right) \tag{21}\]  
Although the obtained parameters are insensitive to noise, the method's accuracy is affected by the window width. It is due to that the STFT has the limitation that it cannot satisfy the optimal frequency resolution and time resolution at the same time [30; 31]. The window width \(h\) of STFT will bring errors to the estimation results. Meanwhile, the accuracy of the IF estimation samples in the LS regression also influences the method's accuracy, especially in weak targets surrounded by clutter. Owing to these two factors, there may be a large deviation between the estimated values and the actual values, which cannot meet the requirement of THz-SAR vibration error compensation. Thus, THz-SAR images are still influenced by residual phase errors.