Spatial Feature-Based ISAR Image Registration for Space Targets

Abstract

Image registration is essential for applications requiring the joint processing of inverse synthetic aperture radar (ISAR) images, such as interferometric ISAR, image enhancement, and image fusion. Traditional image registration methods, developed for optical images, often perform poorly with ISAR images due to their differing imaging mechanisms. This paper introduces a novel spatial feature-based ISAR image registration method. The method encodes spatial information by utilizing the distances and angles between dominant scatterers to construct translation and rotation-invariant feature descriptors. These feature descriptors are then used for scatterer matching, while the coordinate transformation of matched scatterers is employed to estimate image registration parameters. To mitigate the glint effects of scatterers, the random sample consensus (RANSAC) algorithm is applied for parameter estimation. By extracting global spatial information, the constructed feature curves exhibit greater stability and reliability. Additionally, using multiple dominant scatterers ensures adaptability to low signal-to-noise (SNR) ratio conditions. The effectiveness of the method is validated through both simulated and natural ISAR image sequences. Comparative performance results with traditional image registration methods, such as the SIFT, SURF and SIFT+SURF algorithms, are also included.

1. Introduction

ISAR has been used for space target monitoring for many years, which can provide useful information about space targets such as shape structure, body attitude, and working mode, etc. [1,2,3,4]. To improve ISAR image quality and exploit more complete information, image fusion [5] and the three-dimensional (3D) target reconstruction technique [6,7,8] have been extensively investigated, which usually involve jointly processing radar data of the same target taken at different times from different viewpoints and/or by different sensors. Image registration is a crucial step in integrating or comparing different measurements in these applications. For example, different images need to be rotated to a common coordinate system prior to high resolution image formation in the image fusion; the same scatterers from different images need to be correctly registered to compare position differences or phase differences in 3D reconstruction. According to the imaging mechanism of radar systems, an ideal ISAR image can be represented by a two-dimensional (2D) point cloud, which is significantly different from an optical image [9]. Low SNR, sparse scatterer distribution, scatterer defocusing, scatterer glinting, and secondary reflections of targets in ISAR images impact the performance of image registration algorithms. Therefore, image registration methods should be reevaluated and modified according to the characteristics of ISAR images, even if they work well in matching synthetic aperture radar (SAR) images [10], optical images [11], and medical images [12].

ISAR image registration algorithms based on joint translational compensation have been proposed to resolve the translation mismatch problem in the signal domain [13,14,15,16,17,18]. These algorithms utilize the relative rotation of the line of sight (LOS) to compensate for image shifts. In [13], the relationship between time-varying rotation parameters and the compensation phase is derived from the observation geometry in a multiple antenna configuration, such as an L-shaped structure. Several studies have focused on image registration algorithms for non-cooperative space targets, where the relative rotation of target is unknown. In [13,14,15], the rotation trajectory is estimated by tracking dominant scatterers in one-dimensional range profiles along the cross direction. To solve the problem of cross-term phase interference, the authors in [16] proposed measuring the phases of joint-channel phase difference profiles generated by raw echo data from different receiving antennas. Reference [17] combined joint range alignment and joint Doppler centroid tracking to register multichannel ISAR images. In [18], a joint wave path difference compensation algorithm was developed to enable precise image registration for InISAR imaging with SFB-SA signals.

Registration techniques processed in the image domain have also been developed, which can avoid the impact of fluctuating quality in range profiles, especially in low SNR scenarios or when there are few dominant scatterers in the range profiles. Image alignment was achieved using cross-correlation techniques by iteratively searching for the translation and rotation parameters [19]. However, the drawback of this method is its high computational complexity. Since cross-correlation can be calculated in the frequency domain based on the Fourier theorem, computational efficiency can be improved by using the fast Fourier transform (FFT). The authors of [20] utilized phase correlation in the frequency domain to estimate the relative translation offset between two ISAR images. For image rotation parameter estimation, [21] decomposed the image rotation using a three-step FFT-based convolution interpolation process to improve computational efficiency, and [20] applied the polar mapping approach, which turns the rotation of an image into a translation. Since this type of method uses the intensity information of the image, the content of the image has a significant influence on the registration result. In [22], a registration method was proposed by tracking and matching scatterers based on the Kanade–Lucas–Tomasi (KLT) algorithm. The KLT algorithm makes use of spatial intensity information to directly search for the position that yields the best match. The rapid advancement of deep learning has led to significant progress in SAR image registration [23,24]. Deep learning’s feature extraction and matching capabilities have greatly improved image registration accuracy and efficiency. However, there has been limited research on the application of deep learning methods for ISAR image registration, which requires further investigation.

Feature-based methods are powerful tools in remote sensing, medical imaging, and computer vision applications. By focusing on distinctive features within images, these methods can effectively align images despite changes in viewpoint or partial overlaps. Researchers have explored feature-based registration techniques in the field of ISAR image registration [25,26,27,28,29,30,31]. For instance, line features [25] and rectangular features [26] are applied for the image registration of continuously changing ISAR image sequences. The extraction of these advanced features primarily relies on the target’s contours. For a discrete ISAR image, image preprocessing techniques, such as morphological operations, are necessary to produce continuous contours. Reference [27] proposed an image registration method by integrating SIFT and a local correlation matching method. In [28], a combination of scale-invariant feature transform (SIFT) and speeded up robust features (SURF) was employed to align ISAR image sequences and achieve large-angle imaging of ISAR targets. While effective, SIFT and SURF are optimized for images with distinct, well-defined features such as corners, edges, and blobs, but they are sensitive to image noise [32,33]. A well-focused ISAR image can be taken as a composition of many 2Ds sinc functions, with slight blur, located at the projected positions of the 3D scatterers in the ISAR imaging plane [1]. The lack of rich textures in many ISAR images poses challenges, necessitating alternative feature detection methods to achieve optimal results.

Sparsely distributed scatterers resemble point clouds more than color patches. Thus, using scatterers rather than other features to describe an ISAR image is more appropriate. Additionally, using structural information instead of boundary or gradient information in ISAR image registration is more reasonable.

In this work, we propose utilizing the center points of the sinc envelope as feature points (FPs) and constructing a spatial-based feature descriptor by employing the characteristics of distance, angle, and SNR of scatterers. This descriptor exploits essential information from ISAR images and is expected to provide a more stable performance in parameter estimation for image registration.

The innovations of this article are listed as follows:

(1)

A spatial feature of dominant scatterers is proposed to describe the spatial information of ISAR images, offering a more fundamental approach compared to boundary-based or gradient-based features in ISAR image registration.

(2)

The proposed method projects the two-dimensional spatial structure of the target into a one-dimensional feature vector, translating image rotation into a shift of the feature vector, which can be easily matched through circular cross-correlation.

(3)

The traditional random sample consensus (RANSAC) method was enhanced with a consistency check for the rough rotation angle, allowing closer integration with the proposed features and thereby improving the accuracy and stability of mapping parameter estimation, even when the percentage of outliers approaches 50%.

The remainder of this paper is organized as follows: Section 2 analyzes the transform model of ISAR images from different viewpoints. Section 3 introduces the proposed ISAR image registration method, including spatial feature model construction, feature matching, and transform model estimation. Section 4 presents simulated and experimental results to validate the efficacy of the proposed method. Finally, Section 5 provides the discussion, and Section 6 presents the conclusions.

2. Transform Model of ISAR Images

The conventional ISAR imaging methodology generates two-dimensional representations of targets by projecting the reflectivity function onto an imaging projection plane (IPP) [34]. By analyzing the relationship between imaging coordinate systems from two distinct viewpoints, an ISAR image registration model can be derived. This section will analyze the transformation model for registering two ISAR images obtained from different arcs by a ground-based radar imaging system.

Figure 1 illustrates the imaging geometry of a spatial target. An ISAR image sequence is generated for the entire visible orbit. Without a loss of generality, we assume that $t_n$ and $t_{n+1}$ represent the imaging moments for the n-th and n+1-th images in the sequence, with $\mathrm{\Delta }t$ being the time interval between these moments. Let $T_n\left(i_n,j_n,k_n\right)$ and $T_{n+1}\left(i_{n+1},j_{n+1},k_{n+1}\right)$ denoting the imaging coordinate systems corresponding to these moments. $j_n$ and $j_{n+1}$ are the range axes aligned with the imaging LOS. $i_n$ and $i_{n+1}$ are the Doppler axes perpendicular to the range axis and the effective rotation vector. $k_n$ and $k_{n+1}$ are the normal vectors of the IPP, consistent with the direction of the effective rotation vector.

During target observation, LOS rotates with the target’s movement, while the target adjusts its attitude to maintain three-axis stability or accomplish specific missions. This results in relative rotation between the radar and the target. In the imaging model described in [9], a relative LOS vector is defined within the target’s body coordinate system, assuming that the target’s body coordinate system aligns with the instantaneous imaging system $T_n$. The imaging LOS vector $j_{n+1}$ can be expressed as

$$j_{n+1}\left(\mathrm{\Delta }t\right)={\left[\cos \beta \left(\mathrm{\Delta }t\right)\sin \alpha \left(\mathrm{\Delta }t\right),\cos \beta \left(\mathrm{\Delta }t\right)\cos \alpha \left(\mathrm{\Delta }t\right),\sin \beta \left(\mathrm{\Delta }t\right)\right]}^{\mathrm{T}}$$

(1)

where $\alpha$ is the azimuth angle of the LOS, $\beta$ is the elevation angle of the LOS, and the superscript T represents the transpose of the matrix.

By differentiating and normalizing LOS vector, the Doppler vector is derived, as shown in Equation (2). This allows for the determination of the instantaneous IPP at time $t_{n+1}$. The normal vector of the IPP is obtained by taking the cross-product of the range axis and the Doppler axis, as shown in Equation (3).

$$i_{n+1}=\frac{d{j}_{n+1}}{d\mathrm{\Delta }t}/\left|\frac{d{j}_{n+1}}{d\mathrm{\Delta }t}\right|$$

(2)

$$k_{n+1}=i_{n+1}\times j_{n+1}$$

(3)

The transformation matrix between the imaging coordinate systems at these two different time instances is

$$\begin{array}{lll}{R}_{t_{n+1}}^{t_n}& =& \left[i_{n+1},j_{n+1},k_{n+1}\right]\\ & =& \left[\begin{array}{ccc}\frac{-\dot{\beta }\sin \beta \sin \alpha +\dot{\alpha }\cos \beta \cos \alpha }{\sqrt{{\dot{\beta }}^2+{\left(\dot{\alpha }\cos \beta \right)}^2}}& \cos \beta \sin \alpha & \frac{-\dot{\beta }\cos \alpha -\dot{\alpha }\cos \beta \sin \beta \sin \alpha }{\sqrt{{\dot{\beta }}^2+{\left(\dot{\alpha }\cos \beta \right)}^2}}\\ \frac{-\dot{\beta }\sin \beta \sin \alpha -\dot{\alpha }\cos \beta \cos \alpha }{\sqrt{{\dot{\beta }}^2+{\left(\dot{\alpha }\cos \beta \right)}^2}}& \cos \beta \cos \alpha & \frac{\dot{\beta }\sin \alpha -\dot{\alpha }\cos \beta \sin \beta \cos \alpha }{\sqrt{{\dot{\beta }}^2+{\left(\dot{\alpha }\cos \beta \right)}^2}}\\ \frac{\dot{\beta }\cos \beta }{\sqrt{{\dot{\beta }}^2+{\left(\dot{\alpha }\cos \beta \right)}^2}}& \sin \beta & \frac{\dot{\alpha }\cos \beta }{\sqrt{{\dot{\beta }}^2+{\left(\dot{\alpha }\cos \beta \right)}^2}}\end{array}\right]\end{array}$$

(4)

where $\dot{\alpha }$ and $\dot{\beta }$ are the derivatives of $\alpha$ and $\beta$, respectively.

Suppose P is an arbitrary scatterer on the target. Its coordinates in the imaging coordinate system $T_n$ are $p_n={[x_n,y_n,z_n]}^{\mathrm{T}}$ and $p_{n+1}={[x_{n+1},y_{n+1},z_{n+1}]}^{\mathrm{T}}$ in the imaging coordinate system $T_{n+1}$. Based on the analysis above, the coordinate transformation relationship for the scatterer P in these two imaging systems is $p_n=R_{t_{n+1}}^{t_n}{p}_{n+1}$.

In ISAR imaging, target motion compensation techniques, such as envelope alignment and phase focusing, are employed to compensate for target motion and enhance image quality. However, the self-focusing procedure in translational motion compensation causes misalignment of ISAR images. We denote the induced offsets in ISAR processing as $q={\left[\mathrm{\Delta }x,\mathrm{\Delta }y,0\right]}^{\mathrm{T}}$, where $\mathrm{\Delta }x$ and $\mathrm{\Delta }y$ represent the shifts in the cross-range and range dimensions of the image, respectively. The registration model can be represented as follows:

$$\left[\begin{array}{l}{x}_n\\ y_n\end{array}\right]=\left[\begin{array}{ccc}\frac{-\dot{\beta }\sin \beta \sin \alpha +\dot{\alpha }\cos \beta \cos \alpha }{\sqrt{{\dot{\beta }}^2+{\left(\dot{\alpha }\cos \beta \right)}^2}}& \cos \beta \sin \alpha & \frac{-\dot{\beta }\cos \alpha -\dot{\alpha }\cos \beta \sin \beta \sin \alpha }{\sqrt{{\dot{\beta }}^2+{\left(\dot{\alpha }\cos \beta \right)}^2}}\\ \frac{-\dot{\beta }\sin \beta \cos \alpha -\dot{\alpha }\cos \beta \sin \alpha }{\sqrt{{\dot{\beta }}^2+{\left(\dot{\alpha }\cos \beta \right)}^2}}& \cos \beta \cos \alpha & \frac{\dot{\beta }\sin \alpha -\dot{\alpha }\cos \beta \sin \beta \cos \alpha }{\sqrt{{\dot{\beta }}^2+{\left(\dot{\alpha }\cos \beta \right)}^2}}\end{array}\right]\left[\begin{array}{l}{x}_{n+1}\\ y_{n+1}\\ z_{n+1}\end{array}\right]+\left[\begin{array}{l}\mathrm{\Delta }x\\ \mathrm{\Delta }y\end{array}\right]$$

(5)

From Equation (5), it can be noted that the registration model is related with the $\alpha$ and $\beta$. $\alpha$ is always the primary rotational component of LOS, providing the foundation for achieving high resolution in the azimuth direction. In contrast, $\beta$ represents the spatial variation of the IPP, resulting in a height-dependent displacement.

Let ${\left|z_{n+1}\right|}_{\max }$ denote the target’s maximum height. The cross-range displacement $\mathrm{\Delta }{x}_z$ satisfies $\mathrm{\Delta }{x}_z\approx \left|z_{n+1}\right|\tan \beta \sin \alpha \le \beta {\left|z_{n+1}\right|}_{\max }$, and the range displacement $\mathrm{\Delta }{y}_z$ satisfies $\mathrm{\Delta }{y}_z\approx \left|z_{n+1}\right|\sin \beta \cos \alpha \le \beta {\left|z_{n+1}\right|}_{\max }$. According to [35], when the coordinate error is less than one-eighth of a resolution unit, its impact on image correlation and the interferometric phase can be neglected. Thus, the constraint can be derived as $\beta \le 0.125 \min \left({\rho }_r,{\rho }_a\right)/{\left|z_{n+1}\right|}_{\max }$, where ${\rho }_r$ and ${\rho }_a$ are the range and cross-range resolution, respectively. Assuming ${\left|z_{n+1}\right|}_{\max }$ is 10 m, with a bandwidth of 1 GHz and a resolution of 0.15 m, the permissible angular deviation is $\beta \le 0.1074°$.

For slow-rotating space targets, changes in IPP occur more gradually and do not shift rapidly over short periods. By limiting the interval between arc segments, the variation in IPP can be ignored. Section 4.1 provides a detailed analysis of $\alpha$ and $\beta$. Both $\alpha$ and $\beta$ are influenced by various factors, such as the satellite’s orbital parameters, the radar station’s position, and the satellite’s attitude. Notably, $\alpha$ is much larger than $\beta$, particularly when the satellite is passing through the zenith. When the interval between two arcs is less than 10 s, the influence of $\beta$ and $\dot{\beta }$ on image registration can be ignored.

Based on the previous analysis, it can be assumed that $\cos \beta \approx 1$, $\sin \beta \approx \beta$, and $\dot{\beta }\approx 0$. Therefore, the mathematical registration model between two ISAR images from adjacent imaging arcs can be simplified as a rigid transformation that includes both translation and rotation. This can be expressed as

$$\left[\begin{array}{l}{x}_n\\ y_n\\ 1\end{array}\right]=\left[\begin{array}{ccc}\cos \alpha & \sin \alpha & \mathrm{\Delta }x\\ -\sin \alpha & \cos \alpha & \mathrm{\Delta }y\\ 0& 0& 1\end{array}\right]\left[\begin{array}{l}{x}_{n+1}\\ y_{n+1}\\ 1\end{array}\right]$$

(6)

Through the registration of consecutive images, sequential image registration is ultimately employed to achieve long interval image registration.

3. Method

The image registration of ISAR images is the process of overlaying different ISAR images obtained from different viewpoints [11]. The principle of image registration involves extracting and matching representative features of ISAR images to establish mapping parameters between these images. In this section, a feature-based image registration method is proposed by exploiting the concordance of spatial structures of multiple scatterers in various ISAR images.

3.1. Feature Detection

Unlike optical images, ISAR images are not composed of color blocks. Instead, they are a set of slightly blurred 2D sinc functions due to system errors and imperfections in ISAR imaging processing. Assuming that the radar transmits a broadband signal with bandwidth $B$ and center frequency $f_0$, $T_{obs}$ is the coherent processing time under the n-th viewing. After imaging processing, the obtained ISAR image can be expressed as:

$$I_n\left(\widehat{t},f_d\right)\approx {\displaystyle \sum _{i=1}^K{\rho }_{i}B{T}_{obs}\mathrm{sinc}\left[B\left(\widehat{t}-\frac{2r_n^i}{c}\right)\right]\mathrm{sinc}\left[T_{obs}\left(f_d-f_n^i\right)\right]\exp \left(-j\frac{4\pi f_0}{c}{j}_n^{\mathrm{T}}{p}_i\right)}$$

(7)

where $\widehat{t}$ is the fast time, $f_d$ is the Doppler frequency, $K$ represents the number of scatterers on the target, ${\rho }_i$ is the scattering coefficient, $c$ is the velocity of light, and $r_n^i$ and $f_n^i$ denote the 2D location of scatterer P_i in the ISAR image, while $p_i$ is its 3D coordinate in the body coordinate system. For in-orbit space targets, the cumulative rotation angle can be estimated using high-precision orbital measurements and attitude information. This ensures that two adjacent ISAR images have similar resolution. After cross-range scaling, the images are further adjusted to the same size through scaling and cropping.

Traditional salient structures and features, such as prominent regions, lines, region corners, line intersections, and points on curves with high curvature, used for optical image registration, are based on variations in color regions of images. In ISAR imaging, many edges form between dominant scatterers and the background. However, these edges are unstable due to noise, defocusing, and sidelobe interference from adjacent scatterers. Consequently, the performance of methods that rely on edge inflection points or corners as FPs degrades for ISAR image registration. This instability worsens when scatterers are sparsely distributed within well-focused ISAR images. Therefore, it is more appropriate to use the peak information of scatterers rather than their boundary or gradient information for ISAR image registration. In addition to amplitude and phase, the spatial structure of dominant scatterers is a key identifying feature in ISAR images. This paper extracts scatterers with high SNR as FPs and analyzes their relative distribution to construct a feature descriptor.

3.2. Spatial Feature Descriptor

A feature descriptor is a vector representation of an image feature, encapsulating essential information about the feature’s characteristics in a format useful for further processing or analysis. When a reference scatterer is selected for each set of FPs, each FP relative to this reference FP can be modeled as a distance vector. Considering counterclockwise as the positive direction, the angle-distance distribution of these distance vectors can describe the spatial relationships of these FPs.

Figure 2 illustrates a schematic diagram of the reference and sensed images used in image registration. The sets of FPs in the reference and sensed images are denoted as $A=\left\{A_1,A_2,A_3,....,A_{N-1},A_N\right\}$ and $B=\left\{B_1,B_2,B_3,....,B_{M-1},B_M\right\}$, where $N$ and $M$ correspond to the number of dominant scatterers extracted from each image.

Let FP $A_1$ be the reference. We can denote the vector $H^{ori}\left(A_1\right)$ as the distance vectors of FPs relative to the reference $A_1$ in the set.

$$H^{ori}\left(A_1\right)=\left[\overrightarrow{A_1{A}_1},\overrightarrow{A_1{A}_2},\dots ,\overrightarrow{A_1{A}_i},...,\overrightarrow{A_1{A}_N}\right]=\left[0,r_{A_1}^2{e}^{j{\zeta }_{A_1}^2},\dots r_{A_1}^i{e}^{j{\zeta }_{A_1}^i},...,r_{A_1}^N{e}^{j{\zeta }_{A_1}^N}\right]$$

(8)

where $r_{A_1}^i$ denotes the modulus of the vector $\overrightarrow{A_1{A}_i}$ and ${\zeta }_{A_1}^i$ denotes the angle.

FP $A_1$ and FP $B_8$ are the projections of the same scatterer on the IPPs of the reference and sensed images, respectively. Similarly, the distance vectors of the FPs set relative to the reference $B_8$ in the sensed image are

$$H^{ori}\left(B_8\right)=\left[\overrightarrow{B_8{B}_1},\overrightarrow{B_8{B}_2},\dots ,\overrightarrow{B_8{B}_i},...,\overrightarrow{B_8{B}_M}\right]=\left[r_{B_8}^1{e}^{j{\zeta }_{B_8}^1},r_{B_8}^2{e}^{j{\zeta }_{B_8}^2},\dots r_{B_8}^i{e}^{j{\zeta }_{B_8}^i},...,r_{B_8}^M{e}^{j{\zeta }_{B_8}^M}\right]$$

(9)

where $r_{B_8}^i$ denotes the modulus of the vector $\overrightarrow{B_8{B}_i}$, and ${\zeta }_{B_8}^i$ denotes the angle.

The indices of the same scatterer in the two sets can differ significantly, and even the total number of FPs may not be the same. Elements with the same index in $H^{ori}\left(A_1\right)$ and $H^{ori}\left(B_8\right)$ may be different, namely $r_{A_1}^i{e}^{j{\zeta }_{A_1}^i}\ne r_{B_8}^i{e}^{j{\zeta }_{B_8}^i}$. Therefore, computing the conjugate cross-correlation of $H^{ori}\left(A_1\right)$ and $H^{ori}\left(B_8\right)$ may lead to failure in feature matching. To resolve this problem, the distance vectors of the FPs should be sorted and padded before feature matching, and then we have

$${\mathrm{H}}_{A_1}={\mathrm{P}}_{\mathrm{sz}}\left\{\underset{{\zeta }_{A_1}^i\uparrow }{\mathrm{sort}}\left(H^{ori}(A_1)\right)\right\}$$

(10)

$${\mathrm{H}}_{B_8}={\mathrm{P}}_{\mathrm{sz}}\left\{\underset{{\zeta }_{B_8}^i\uparrow }{\mathrm{sort}}\left(H^{ori}(B_8)\right)\right\}$$

(11)

where $\mathrm{sort}(\cdot )$ is the sorting function that rearranges the distance vectors of FPs in increasing order according to the vector angle, ${\zeta }_{A_1}^i$ or ${\zeta }_{B_8}^i$. ${\mathrm{P}}_{\mathrm{sz}}(\cdot )$ is the process of scatterer superposition and zero-padding in the generation of $H_{A_1}$ and $H_{B_8}$. After this process, scatterers with the same angle will have the same index in $H_{A_1}$ and $H_{B_8}$, enabling the use of cross-correlations of $H_{A_1}$ and $H_{B_8}$ to describe their similarity. In practical applications, the correlation of such discrete features is highly sensitive to angular errors. Therefore, an improved continuous feature is considered.

Noise affects the estimation accuracy of scatterer coordinates, which subsequently impacts the angular precision of the corresponding distance vectors of FPs. According to [36], the coordinate errors in range and cross-range are independent Gaussian random variables. Under the small error approximation, the angular error is a linear combination of the coordinate errors and can be approximated as Gaussian noise. Therefore, the angle value in the feature descriptor can be replaced with a parameterized Gaussian function to achieve an angular-error-tolerant feature. Finally, the continuous feature descriptor of the scatterer $A_1$ is expressed as the accumulation of N-1 Gaussian functions, which is

$${\mathrm{H}}_{A_1}^{\mathrm{G}}\left(\zeta \right)={\displaystyle \sum _{i=1}^N\frac{1}{{\sigma }_{A_1}^i\sqrt{2\pi }}{e}^{-\frac{{\left(\zeta -{\zeta }_{A_1}^i\right)}^2}{2{\left({\sigma }_{A_1}^i\right)}^2}}}$$

(12)

where the range of variable $\zeta$ is from −$\pi$ to $\pi$, ${\zeta }_{A_1}^i$ is the mean of the Gaussian function, and ${\sigma }_{A_1}^i$ is the standard deviation, which is calculated based on the SNR and the modulus of the distance vectors of FPs. According to [36], the 2D coordinate errors between two scatterers in the ISAR image are expressed as

$$\mathrm{\Delta }{x}_{A_1}^i=\left(\frac{1}{\sqrt{SN{R}_{A_1}}}+\frac{1}{\sqrt{SN{R}_{A_i}}}\right)\frac{\sqrt{3}}{\pi }{\rho }_a\hspace{1em}\hspace{1em}\mathrm{\Delta }{y}_{A_1}^i=\left(\frac{1}{\sqrt{SN{R}_{A_1}}}+\frac{1}{\sqrt{SN{R}_{A_i}}}\right)\frac{\sqrt{3}}{\pi }{\rho }_r$$

(13)

where $\left(\mathrm{\Delta }{x}_{A_1}^i,\mathrm{\Delta }{y}_{A_1}^i\right)$ denote the cross-range error and range error, respectively, and $SN{R}_{A_1}$ and $SN{R}_{A_i}$ represent the peak SNR of scatterer $A_1$ and $A_i$.

Figure 3 shows the schematic diagram of angle error $\mathrm{\Delta }{\zeta }_{A_1}^i$. The standard deviation of ${\zeta }_{A_1}^i$ can be derived as

$${\sigma }_{A_1}^i=f\left(SN{R}_{A_i},SN{R}_{A_1},r_{A_1}^i\right)=\arcsin \left(\sqrt{\mathrm{\Delta }{x}_{A_{1,i}}^2+\mathrm{\Delta }{y}_{A_{1,i}}^2}/r_{A_1}^i\right)$$

(14)

We present an example that utilizes this new feature descriptor. We assume that each line segment in Figure 2 is 5 m. The lengths and angles of each vector from $A_1$ to other points are listed in

Table 1. The length and angle of $\overrightarrow{A_1{A}_i}$.

Scatterer Number	1	2	3	4	5	6	7	8	9	10
$r_{A_1}^i$ (m)	0	5.00	7.07	11.18	14.14	18.02	15.81	15.00	7.50	11.28
${\zeta }_{A_1}^i$ (degree)	0	90.00	45.00	63.43	45.00	33.69	18.43	0	0	12.80

. Given that the SNR of each scatterer is 10 dB and the resolutions, ${\rho }_r$ and ${\rho }_a$, are 1 m, we can calculate the variance ${\sigma }_{A_1}^i$ of the angle ${\zeta }_{A_1}^i$ using Equations (13) and (14). The feature of scatterer $A_1$ is then obtained using Equation (12). Figure 4 illustrates the feature curves of scatterer $A_1$ along with its 10 components. It is worth noting that the modulus of the distance vector and the amplitude of the FPs are incorporated into the feature descriptor through ${\sigma }_{A_1}^i$. As shown, a longer modulus corresponds to a taller and narrower Gaussian waveform, while a shorter modulus results in a shorter and wider Gaussian waveform.

In Figure 5, the feature curve of $A_1$ is compared with those of $B_4$ and $B_8$. It can be observed that the $A_1$ curve and the curve $B_4$ are quite different, whereas the $A_1$ curve and the $B_8$ curve show a stronger correlation after an angle shift, despite some outlier pairings. The proposed method projects the two-dimensional structure into a one-dimensional vector and maps the image rotation into a translation of the one-dimensional feature vector.

3.3. Feature Matching

After obtaining the feature descriptors of the FPs, the next step is feature matching. Image rotation causes a shift in the feature curves. Therefore, this paper employs the circular cross-correlation method for feature matching [27,37]. The circular cross-correlation of scatterer $A_i$ and scatterer $B_j$ is

$${\chi }_{A_i{B}_j}\left(\zeta \right)={\displaystyle \int {\mathrm{H}}_{A_i}^{\mathrm{G}}\left(\phi \right){\mathrm{H}}_{B_j}^{\mathrm{G}}\left(\phi -\zeta \right)d\phi }$$

(15)

Performing circular cross-correlation on all the FPs yields a cross-correlation matrix ${\gamma }_{N\times M}$.

$${\gamma }_{N\times M}=\left[\begin{array}{cccc}{\gamma }_{11}& {\gamma }_{12}& \cdots & {\gamma }_{1M}\\ {\gamma }_{21}& {\gamma }_{22}& \cdots & {\gamma }_{2M}\\ \cdots & \cdots & \cdots & \cdots \\ {\gamma }_{N1}& {\gamma }_{N2}& \cdots & {\gamma }_{NM}\end{array}\right]$$

(16)

where ${\gamma }_{ij}=\max ({\chi }_{A_i{B}_j})$. In this circular cross-correlation matrix, the row number represents the index of the scatterer in the reference image, and the column number represents the index of the scatterer in the sensed image. Furthermore, the circular cross-correlation can be obtained via the FFT, which greatly improves computational efficiency.

For scatterer index $i$ in the reference image, its matching scatterer index $\widehat{j}$ in the sensed image can be obtained by searching for the maximum value in the $i$-th row. Let $\mathrm{argmax}(\cdot )$ stand for the argument of the maxima, it can be expressed as

$$\widehat{j}=\underset{j}{\mathrm{argmax}}\left({\gamma }_{ij}\right),\left(j=1,2,\dots ,M\right)$$

(17)

Feature matching is typically achieved by searching for the maximum value in each row. However, false matches are inevitable in image feature matching due to factors such as symmetrical structures, repetitive patterns, and noise. Additionally, viewpoint changes and occlusion contribute to these errors. Estimating registration parameters accurately is challenging in the presence of outliers because they introduce noise and bias, leading to incorrect estimates. To ensure precise registration, it is crucial to identify and remove outliers before estimating the registration parameters.

3.4. Two-Step Transform Model Estimation Method

In this subsection, we provide a two-step transform model estimation process, which includes outlier detection, removal, and parameter estimation.

The first step uses a statistical method for outlier filtering. Since the angular shift that corresponds to the maximum circular cross-correlation is close to the image’s rotation angle, a histogram analysis is applied to the angle shifts of all matching pairs. This histogram-based method automatically identifies incorrect matches, such as mirror mismatches caused by symmetrical structures. The angle intervals can be set based on a rough estimate of the image rotation. By adjusting the histogram’s central angle, most correct matches can be grouped within one interval, while matches outside this range are considered outliers. However, this method is less effective at excluding redundant adjacent scatterer pairs with similar rotation angles.

In the second step, a model-based approach identifies adjacent scatterer pairs using the RANSAC algorithm. The RANSAC algorithm fits a model to the data and identifies redundant outliers based on their deviation from the model. The parameters of the transformation model are estimated from the residual inliers.

RANSAC is an iterative method used to estimate parameters of a mathematical model from a set of observed data containing outliers [38]. According to the RANSAC algorithm, a subset of matched scatterer pairs is randomly selected, and the related parameters are estimated by solving the transform model equation. This allows redundant matches and mismatches to be excluded through consistency checks. It iterates through two steps: generating a hypothesis from random samples and verifying it against the data [39]. The most probable hypothesis, supported by the most inlier candidates, is chosen during the iterations. The RANSAC algorithm can be formulated by

$$\left[\mathrm{\Delta }\widehat{x},\mathrm{\Delta }\widehat{y},\widehat{\alpha }\right]=\underset{\mathrm{\Delta }x,\mathrm{\Delta }y,\alpha }{\mathrm{argmin}}\left\{{\displaystyle \sum _{i=1}^{N_e}\mathrm{Loss}\left(\mathrm{Err}\left(A_i,B_{\widehat{j}};\mathrm{\Delta }x,\mathrm{\Delta }y,\alpha \right)\right)}\right\}$$

(18)

where $\mathrm{\Delta }x$, $\mathrm{\Delta }y$, $\alpha$ are the registration parameters as defined in Section 2, $\mathrm{\Delta }\widehat{x}$, $\mathrm{\Delta }\widehat{y}$, $\widehat{\alpha }$ are their estimated values, $A_i$ and $B_{\widehat{j}}$ are scatterer pairs for matching the reference image and sensed image, respectively, $N_e$ represents the number of matched scatterer pairs after the initial outlier exclusion, $\underset{\mathrm{\Delta }x,\mathrm{\Delta }y,\alpha }{\mathrm{argmin}}(\cdot )$ is the arguments of the minima, $\mathrm{Err}(\cdot )$ is the error function, and $\mathrm{Loss}(\cdot )$ is the loss function. The explicit definitions of these functions and parameter estimation procedure are provided below.

Assuming L scatterers are randomly selected, L must be at least three due to the need to estimate three parameters. The choice of L should balance model accuracy and algorithm efficiency. By substituting the scatterers’ coordinates into the transformation model in Equation (6), it can be derived

$$A_{2L\times 7}{h}_{7\times 1}=0$$

(19)

where $h={\left[\begin{array}{ccccccc}\cos \alpha & \sin \alpha & \mathrm{\Delta }x& -\sin \alpha & \cos \alpha & \mathrm{\Delta }y& 1\end{array}\right]}^T$ is the column vector related to parameters to be solved, $A$ is a Matrix related to scatterer’s coordinates, and the 2l-1th row and 2lth row are

$$A_{\left(\begin{array}{l}2l-1,:\\ \hspace{1em}2l,:\end{array}\right)}=\left[\begin{array}{ccccccc}{x}_{n+1}^l& y_{n+1}^l& 1& 0& 0& 0& -x_n^l\\ 0& 0& 0& x_{n+1}^l& y_{n+1}^l& 1& -y_n^l\end{array}\right]$$

(20)

where $x_n^l$ and $y_n^l$ denote the scatterer’s coordinates in the n-th ISAR image, while $x_{n+1}^l$ and $y_{n+1}^l$ denote the coordinates of the corresponding scatterer in the n+1-th ISAR image.

Through SVD decomposition of matrix $A$, we obtain $A_{2L\times 7}=U_{2L\times 2L}{\Sigma }_{2L\times 7}{V}_{7\times 7}^T$. The singular vector $v_{\min }$ corresponding to the smallest singular value is the non-zero solution that minimizes $Ah$ [40]. Thus, $h\approx v_{\min }$. The mapping function parameters between different ISAR images can be obtained as follows:

$$\mathrm{\Delta }\widehat{x}=v_{\min }\left(3\right)\hspace{1em}\mathrm{\Delta }\widehat{y}=v_{\min }\left(6\right)\hspace{1em}\widehat{\alpha }=\arctan \left[v_{\min }\left(2\right)/v_{\min }\left(1\right)\right]$$

(21)

The coordinates of scatterers in the sensed image can be recalculated based on the estimation results of this randomly selected subset of matched scatterer pairs. Let $H$ represent the estimated transformation matrix, which is shown in Equation (6), then the new coordinates of scatterer ${\tilde{B}}_{\widehat{j}}$ can be obtained as ${\tilde{B}}_{\widehat{j}}=H{B}_{\widehat{j}}$. The Euclidean distance between scatterers is used as the error function $\mathrm{Err}(\cdot )$, and the distance error $e$ is expressed as

$$e=\mathrm{Err}\left(A_i,{\tilde{B}}_j\right)=\sqrt{{\left(x_n^{A_i}-x_{n+1}^{{\tilde{B}}_{\widehat{j}}}\right)}^2+{\left(y_n^{A_i}-y_{n+1}^{{\tilde{B}}_{\widehat{j}}}\right)}^2}$$

(22)

where $x_n^{A_i}$ and $y_n^{A_i}$ denote the coordinates of scatterer $A_i$, and $x_{n+1}^{{\tilde{B}}_{\widehat{j}}}$ and $y_{n+1}^{{\tilde{B}}_{\widehat{j}}}$ denote the coordinates of scatterer ${\tilde{B}}_{\widehat{j}}$.

The subset of matched scatterer pairs is considered an inlier candidate, when its error $e$ is within a predefined threshold, formulated as

$$\mathrm{Loss}\left(e\right)=\left\{\begin{array}{l}0,e<TH\\ 1,\mathrm{otherwise}\end{array}\right.$$

(23)

This threshold, $TH$, depends on the image resolution and SNR and can be set to $3{\sigma }_{\mathrm{dis}}$ or $5{\sigma }_{\mathrm{dis}}$. ${\sigma }_{\mathrm{dis}}$ is the standard deviation of the distance between the same scatterers in the two images when it is ideally registered, and it can be calculated by Equation (13).

Once the outliers are identified, they can be eliminated from the dataset before estimating the registration parameters. The remaining scatterers can be taken as inliers and used to estimate the registration parameters using least squares estimation (LSE).

3.5. Flowchart of the ProposedMethod

The entire procedure of the proposed ISAR image registration method can be summarized as feature detection, feature matching, transform model estimation, and image transformation. The flowchart is shown in Figure 6.

The specific processing steps are as follows.

Step 1: Analyze the noise power in the noise region to set the threshold for scatterer extraction and extract dominant scatterers as FPs in each image using the CLEAN technology.

Step 2: Calculate the distance vectors between any two scatterers within the image and construct the feature descriptor for each FP based on shape structures.

Step 3: Calculate the circular cross-correlation matrix and the angle shift of all matching scatterer pairs via FFT and then obtain the feature matching results.

Step 4: Discard false matched pairs by a consistency check of the estimated rotation angle.

Step 5: Discard redundant adjacently outliers and estimate the parameters of transform model by RANSAC.

Step 6: Image registration based on the estimated parameters.

4. Experiments and Results

In this section, several experiments based on both simulated and measured ISAR images are conducted to validate the effectiveness of the proposed method. Additionally, the well-known feature registration algorithms, SIFT, SURF and SIFT+SURF, are compared with the proposed method. The real ISAR images used in this section are sourced from a report published by the German FGAN Lab in 2012 (@Fraunhofer FHR) [21]. This report provides an ISAR image sequence of the US space shuttle, obtained with the TIRA system after long arc imaging.

4.1. Azimuth and Elevation Angles of LOS

In this subsection, the azimuth and elevation angles for different visible passes, imaging times, and time intervals are simulated. In this simulation, the radar station is set in Beijing (39.9° N, 116.4° E, 88 m), and the satellite orbit is the TIANHE space station. The two-line orbital elements (TLE) of the TIANHE space station are listed in

Table 2. Two-line orbital elements (TLE) of the TIANHE space station.

1 48274U 21035A 23059.87944444 00018297 00000-0 21063-3 0 9991

2 48274 41.4746 96.4827 0003993 340.8847 14.0436 15.61544460104882

. The observation time for the first visible pass is from 2:22:07 to 2:32:31 on 1 March 2023, and for the second visible pass, from 3:59:00 to 4:09:26 on 1 March 2023.

Figure 7a,b illustrate the variations in azimuth and elevation angles, while Figure 7c,d depict the derivation of these angles. According to the definition in [9], the center moment of coherent processing interval is defined as the imaging moment. In Figure 7, the time values on the horizontal axis represent the imaging moment of the reference image, while the imaging moment of the sensed image is offset by a time interval $\mathrm{\Delta }t$ with respect to the reference moment. During the observation, the LOS form a curved surface in the target’s body coordinate system, which expands more in the horizontal plane than the vertical change. As a result, the azimuth angle variation is significantly larger than the elevation angle variation, particularly near the zenith-passing time. When the arc segment interval is less than 10 s, with alpha less than 10° and beta less than 0.05°, and a resolution of 0.15 m, the displacement for a target at a height of 10 m is 0.0086 m. This displacement, being less than one-tenth of a resolution unit, is negligible.

4.2. Experiment with Simulated Images

This section presents the experimental results based on simulated ISAR images. The scatterer model contains 80 dominant scatterers, extracted from a real ISAR image of the US space shuttle given in Figure 8a. It is assumed that there is a 20% outlier rate between the reference image and the sensed image to simulate phenomena such as angle glint and occlusion effects in ISAR imaging. So, scatterers are randomly divided into three groups according to the number of outliers: set1 (54 scatterers), set2 (13 scatterers), and set3 (13 scatterers), as shown in Figure 8a. The scatterers in set1 and set2 are used to generate the reference image shown in Figure 8b, and the scatterers in set1 and set3 are used to generate the sensed image shown in Figure 8c. The sensed image has been affine-transformed according to the parameters listed in

Table 3. Estimation results of the mapping parameters.

Method	$\mathbf{\Delta }\mathit{x}$ (Pixel)	$\mathbf{\Delta }\mathit{y}$ (Pixel)	$\mathit{\alpha }$ (Degree)	$\mathbf{Error} \mathbf{\Delta }\mathit{x}$	$\mathbf{Error} \mathbf{\Delta }\mathit{y}$	$\mathbf{Error} \mathit{\alpha }$
Ground Truth	−10	20	−15	-	-	-
SIFT	−9.5278	20.1025	−14.8710	0.4721	0.1025	0.1284
SURF	−9.6035	19.6803	−14.4706	0.3965	−0.3197	0.5294
SIFT + SURF	−10.2273	20.0843	−14.9511	−0.2273	0.0843	0.0489
Proposed method	−10.0911	20.0119	−14.9668	−0.0911	0.0119	0.0332

.

The proposed image registration algorithm is validated using the sensed and reference images from Figure 8. Complex Gaussian noise is added to each image to simulate a noise environment with a 24 dB SNR. After extracting FPs, the feature descriptor of each point is constructed according to Equation (12), and the circular cross-correlation of the feature descriptors between any two reference FPs is calculated.

Figure 9a shows the maximum circular cross-correlation matrix, where the brightness of each cell indicates the cross-correlation of two feature curves, implying the reliability of the matching. It is observed that a few bright points are sparsely distributed in the maximum circular cross-correlation matrix, which can be taken as candidates for matched scatterer pairs. Finally, scatterer pairs are obtained by searching for the maximum value row by row, and these are taken as the coarse matches before the consistency check. The histogram of shift feature curves for matched pairs is shown in Figure 9b, with an interval of 20° between adjacent bins. It is observed that a tall column in the histogram is contributed by the correct matches and adjacent matches, while the scattered small values correspond to the false matches.

Figure 10 shows the feature matching results: red lines mark the mismatched pairs identified from the histogram, green lines mark the redundant adjacent matches identified during the RANSAC process, and yellow lines mark the correctly matched pairs used to estimate the mapping parameters. As per outlier detection, there are 8 false matches, 5 redundant adjacent matches and 54 correct matches.

Figure 11a depicts the superposition of the reference image and the sensed image before registration, which shows an obvious displacement. Figure 11b depicts the superposition after registration, where the two images are overlapped perfectly. Furthermore, the estimation results listed in Table 3 are compared with those of classical image registration methods, such as the SIFT algorithm, the SURF algorithm and the SIFT+SURF algorithm. To ensure the objectivity and fairness of the various experiments, all four methods used the RANSAC algorithm during the parameter estimation phase, randomly selecting four scatterer pairs per iteration with a total of 10,000 iterations. It can be seen that the estimation error of the proposed method is much smaller.

4.3. Experiment with Real Images

In this subsection, two different real ISAR images of a US Space Shuttle are used to test the effectiveness of the proposed method. The reference and sensed images are shown in Figure 12a,b, respectively. We sequentially extracted 200 dominate scatterers according to their amplitude from each ISAR image.

The matched pairs are shown in Figure 13. There are 65 pairs of false matches marked with red lines, 112 pairs of redundant adjacently matches marked with green lines, and 23 pairs of best matches marked with yellow lines.

The sensed image is transformed based on the estimated parameters, and the ISAR images before and after registration are compared in Figure 14. Figure 14a shows the superposition of the two ISAR images before registration, and Figure 14b shows the registration results. It can be observed that the main body is perfectly overlapped, although the tails caused by secondary reflection are displaced. Finally, the proposed method is compared with the SIFT, SURF, and SIFT + SURF methods. The estimation errors and the normalized image correlation of the two images after registration are listed in

Table 4. Estimation results of the mapping parameters and normalized image correlation.

Method	Estimated Value			Normalized Image Correlation
	$\mathbf{\Delta }\mathit{x}$ (Pixel)	$\mathbf{\Delta }\mathit{y}$ (Pixel)	$\mathit{\alpha }$ (Degree)	Whole Image	Filtered Image
SIFT	3.0767	4.0650	−12.4428	0.6911	0.7136
SURF	1.2862	5.2708	−16.7880	0.6902	0.7181
SIFT+SURF	5.2214	6.0782	−16.4055	0.6907	0.7196
Proposed method	3.6828	4.8772	−14.3980	0.7360	0.7658

. The filtered image, with its tail removed, ensures that the results of image correlation are not affected by the secondary reflections in ISAR images.

4.4. Robustness Analysis

To further evaluate the robustness of the proposed image registration algorithm, we conducted 1000 Monte Carlo simulations under varying SNR levels and different outlier ratios. In the robustness versus SNR experiments, noise was added to the original ISAR image, and the SNR was adjusted from 10 dB to 30 dB in 2 dB increments using the scatterer model and mapping parameters described in Section 4.2. The mapping parameter estimation results and the image registration success rates for the three methods were compared.

Successful image registration is evaluated based on a preset estimation error threshold. The threshold for range and cross-range shifts is set to half a resolution unit, while the estimation error for the rotation angle is set at 1.6°, ensuring that the shift of the farthest scatterer remains within half a resolution unit. Figure 15 shows the success probabilities of the four image registration methods at different SNR levels. The results indicate that the proposed method achieves the highest success rate under low SNR conditions, whereas the SIFT method has the lowest success rate. The mean error (ME) and root mean squared error (RMSE) for range shift, cross-range shift, and rotation are depicted in Figure 16. Although the success rates of the SURF and SIFT + SURF methods are comparable to the proposed method when the SNR exceeds 18 dB, Figure 16 shows that the proposed method provides higher accuracy in parameter estimation, as indicated by the lower mean and variance of the registration parameters. The presence of noise reduces the number of detectable scatterers, introduces more outliers, and increases the localization error, consequently lowering both registration accuracy and method stability.

Similarly, the impact of outliers is analyzed using 1000 Monte Carlo simulations at 24 dB SNR. The ratio of outliers is varied from 0% to 50% in steps of 5%. Figure 17 illustrates the success probabilities of the four registration methods at different outlier ratios. It can be noted that the success rate of the SIFT+SURF method is lower than that of the SURF method when the proportion of outliers is high. Although the combination of SIFT and SURF increases the number of FPs, the mismatch pairs introduced by SIFT make it more difficult for RANSAC to find inliers when the same number of iterations is used, resulting in reduced registration accuracy. Figure 18 depicts the estimation ME and RMSE of the transformation parameters. It can be observed that the proposed method maintains its performance even as the outlier ratio increases to 50%.

4.5. Computational Complexity Analysis

In this subsection, the computational complexity of the proposed method is analyzed. It is evaluated for four core stages: feature point extraction, feature descriptor construction, feature matching, and model parameter estimation. Assuming the image size is $N_x\times N_y$, containing $K_s$ scatterers, the computational complexity of each stage is denoted as $C_d$, $C_f$, $C_m$ and $C_e$ respectively.

(1)

At the feature point extraction stage, the sinc-based CLEAN method is used to extract scatterers. Assuming $Z_s$ is the number of pixels in the segmented area, the computational complexity is approximately $C_d\approx O\left(K_s{N}_x{N}_y+10K_s{Z}_s\right)$ [41].

(2)

At the feature descriptor construction stage, the computation is primarily focused on calculating the vector angles and magnitudes as well as the summation of Gaussian functions. Assuming the length of the Gaussian function is $G_s$, the computational complexity is approximately $C_f\approx O\left(K_s^2{G}_s\right)$.

(3)

Feature matching is implemented using the circular cross-correlation method. Assuming the feature vector length is $F_s$, the computational complexity using the conventional time-domain correlation is approximately $C_m\approx O\left(0.5K_s^2{F}_s^2\right)$. By applying FFT for frequency-domain processing, the complexity can be reduced to approximately $C_m\approx O\left(0.5K_s^2{F}_s+0.5K_s{F}_s\log F_s+0.25K_s^2{F}_s\log F_s\right)$.

(4)

Model parameter estimation is performed iteratively. Each iteration involves singular value decomposition of the observation matrix, transformation of feature point coordinates, and distance loss calculation. Assuming the observation matrix size is $A_m\times A_n$ and the number of scatterers after consistency checking is $K_{s2}$, the computational complexity for $I$ iteration is approximately $C_e\approx O\left(A_m{A}_n^{2}I+A_m^2{A}_{n}I+11K_{s2}I\right)$.

From the above analysis, it can be concluded that the total computational complexity is approximately $C_T=2C_d+2C_f+C_m+C_e$.

The computational complexities of different methods are compared through simulations.

Table 5. Related parameters.

Parameter Description	Symbol	Value
Image Size	$N_x\times N_y$	$487\times 487$
Number of FPs	$K_i;K_u;K_s;K_{s2}$	168; 198; 200; 135
Feature Length	$F_i;F_u;F_s$	128; 128; 1024
Scales	$S_i$$;S_u$	5; 5
Octaves	$V_i$	3
Neighboring Region Size	$Z_i$$;Z_u$$;Z_s$	16; 20; 16
Gaussian Kernel Size	$G_i$$;G_s$	16; 512
Observation matrix Size	$A_m\times A_n$	$8\times 7$
Iteration Number	$I$	10,000

Notes: Subscript i, u and s represents the SIFT, SURF, and proposed method, respectively.

lists the relevant parameters and the specific values used in Section 4.3.

Table 6. The computational complexity formulas in the four stages.

Method	SIFT	SURF	Proposed Method
$C_d$	$O\left(N_x{N}_y{G}_i^2{S}_i{V}_i\right)$	$O\left(3N_x{N}_y{S}_u\right)$	$O\left(K_s{N}_x{N}_y+10K_s{Z}_s\right)$
$C_f$	$O\left[K_i{Z}_i^2\left(2+F_i\right)\right]$	$O\left[K_u{Z}_u^2\left(1+F_u\right)\right]$	$O\left(K_s^2{G}_s\right)$
$C_m$	$O\left(K_i^2{F}_i\right)$	$O\left(K_u^2{F}_u\right)$	$\begin{array}{l}O\left(0.5K_s^2{F}_s+0.5K_s{F}_s\log F_s\right.\\ \hspace{1em}\left.+0.25K_s^2{F}_s\log F_s\right)\end{array}$
$C_e$	$O\left(A_m{A}_n^{2}I+A_m^2{A}_{n}I+11K_{i}I\right)$	$O\left(A_m{A}_n^{2}I+A_m^2{A}_{n}I+11K_{u}I\right)$	$O\left(A_m{A}_n^{2}I+A_m^2{A}_{n}I+11K_{s2}I\right)$

presents the computational complexity formulas for each method at various stages. Figure 19 shows the simulation results of the computational complexity using different methods.

The SIFT method exhibits relatively high computational complexity, primarily because it requires generating a multi-scale Gaussian pyramid during the keypoint extraction stage. The SURF method effectively reduces computational complexity by employing integral images, Hessian matrix approximation, simplified feature descriptors, and accelerated orientation computation, making it the least complex of the four methods. The combination of the SIFT and SURF methods retains relatively high computational complexity. The computational complexity of the proposed method lies between that of SIFT, SIFT+SURF and SURF. In the proposed method, due to the use of longer feature descriptors and the circular cross-correlation method for matching, the feature matching stage contributes the second most to the overall computational complexity. However, with the introduction of preliminary outlier filtering, the complexity of the proposed method in the parameter estimation stage is relatively lower compared to the other three methods.

5. Discussion

In the above research, it was demonstrated that using the spatial distribution of dominant scatterers for registration is a more effective approach in ISAR image registration. For sparse and discrete ISAR images, the effectiveness of the SIFT and SURF methods decreases. This is because SIFT and SURF identify inflection points, corners, and textured areas as FPs, but in ISAR images, the edges formed between scatterers and the background are often unstable and unreliable, leading these methods to detect numerous invalid FPs. Additionally, SIFT and SURF rely on local gradient information around FPs for feature description, which restricts the amount of information that can be extracted from ISAR images. The proposed method uses the distance and angular information between all scatterers, adopting a global approach. Therefore, it is more effective for ISAR registration.

The proposed method encounters limitations in scenarios with significant IPP variations, such as target maneuvers, large viewing angle differences between distant stations, or images captured from widely separated arc segments at the same location. Variations in IPP result in projection discrepancies and image decorrelation. Improving radar image registration under such conditions remains a challenging but essential task and will be addressed in future research.

6. Conclusions

This paper introduces a novel ISAR image registration algorithm based on the spatial distribution of dominant scatterers. The method projects two-dimensional images into one-dimensional feature vectors by utilizing the distance vectors between scatterers. Image rotation is mapped to a shift in the feature vector, facilitating feature matching through circular cross-correlation. To mitigate error impact, a continuous Gaussian-based feature descriptor is employed. Finally, consistency checks and the RANSAC algorithm are combined to estimate registration parameters, enhancing efficiency by filtering out scatterer pairs with significant errors. Simulation results show that the proposed algorithm delivers more accurate parameter estimates than SIFT, SURF, and SIFT+SURF across various SNR levels. Additionally, it demonstrates greater stability in estimating transformation parameters, even when outliers constitute up to 50%. Real ISAR image registration results further confirm the accuracy and effectiveness of the approach.