Issue 
Natl Sci Open
Volume 3, Number 5, 2024
Special Topic: Microwave Vision and SAR 3D Imaging



Article Number  20230082  
Number of page(s)  10  
Section  Information Sciences  
DOI  https://doi.org/10.1360/nso/20230082  
Published online  29 April 2024 
RESEARCH ARTICLE
Segmentationaided phase unwrapping for 3D reconstruction with airborne arrayInSAR images
The Key Laboratory for Information Science of Electromagnetic Waves (MoE), Fudan University, Shanghai 200433, China
^{*} Corresponding authors (emails: fm_hu@fudan.edu.cn (Fengming Hu); fengxu@fudan.edu.cn (Feng Xu))
Received:
11
December
2023
Revised:
25
March
2024
Accepted:
26
April
2024
Unmanned aerial vehicle (UAV) array InSAR is a new type of singleflight 3D SAR imaging system with the advantages of high coherence and resolution. However, due to the low altitude of the platform, the elevation ambiguity of the system is smaller than the maximal terrain elevation. Since the spatial phase unwrapping is conducted based on the assumption of phase continuity, the inappropriate ambiguity will cause significant unwrapping errors. In this work, a 3D phase unwrapping algorithm assisted by image segmentation is proposed to address this issue. The Markov random field (MRF) is utilized for image segmentation. The optimal spatial phase unwrapping network is achieved based on the segmentation results. The asymptotic 3D phase unwrapping is adopted to get the refined 3D reconstruction. Results based on the real airborne arrayInSAR data show that the proposed method effectively improves the elevation ambiguity.
Key words: 3D reconstruction / image segmentation / phase unwrapping / arrayInSAR
© The Author(s) 2024. Published by Science Press and EDP Sciences.
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
INTRODUCTION
Synthetic aperture radar (SAR) is an imaging radar, working on spaceborne or airborne platforms, has the capabilities of highresolution, dayandnight, and weatherindependent imaging, which is widely used in various applications[1]. The targets in 2D radar imaging are mixed along the crossrange direction. Using two SAR images with a slight different viewing angle, the interferograms caused by phase differences can be obtained[2]. The interferometric phase can be used for terrain elevation mapping, denoted as SAR interferometry (InSAR)[3]. The spatial phase unwrapping is the main step of InSAR with the assumption of phase continuity[4]. However, this assumption often fails if the elevation difference is larger than the height ambiguity. The multibaseline InSAR is an efficient way to solve this issue by the baseline diversity[5] and achieving precise height measurements of the terrain[6]. However, conventional spaceborne SAR system acquires repeatpass images using a single antenna, which requires a long period to accumulate multibaseline InSAR (MBInSAR) observations. To solve this problem, the array InSAR system with a linear antenna array uses the multiinput multioutput (MIMO) technique to acquire a stack of multibaseline images in a single flight, which significantly improves the practical capability of 3D reconstruction[7].
In the practical application, the linear antenna array with a large number of channels will incur high costs. Moreover, the array InSAR system with more channels will increase the system and processing complexity because of the crosschannel calibration and decoupling. To reduce the number of channels, an asymptotic 3D phase unwrapping (PU) is proposed to achieve a reliable 3D reconstruction using only threechannel array InSAR images[7]. This type of 3D PU assumes that the acquired MB acquisitions include both short and long baseline interferograms, following the successful unwrapping criteria. If the systematical height ambiguity is less than the actual height difference, the unwrapping errors in the initialization by the short baseline interferogram will bias the final results. The singlepass arrayInSAR interferograms are with high coherence for both artificial and natural targets, which increases the probability of phase jumping.
A potential improvement for reliable phase unwrapping is separate the artificial targets to get a refined unwrapping network, denoted as segmentationaided phase unwrapping. However, different from the optical images, it is more difficult extract the information from the SAR images. For target objects in SAR images, the information may overlap, which can lead to the inability to accurately extract the phase information of the target in the image, and cause difficulties for PU.
Image segmentation algorithms can be theoretically divided into two categories. The first category includes traditional segmentation algorithms such as threshold segmentation[8], edge detection[9], and clustering algorithms[10]. The second category includes segmentation algorithms based on specific theories, such as graph theorybased segmentation algorithms[11] and the Markov random fieldbased segmentation algorithm[12], etc.
In this work, we propose a segmentationaided phase unwrapping (PU) to avoid the unwrapping errors, leading to a better 3D reconstruction. First, we iteratively segment the SAR image using a Markov random field model[13]. Then, the unwrapping network can be refined using the segmented image. Finally, the experiments are conducted using real airborne arrayInSAR data.
The remainder of this paper is organized as follows. We introduce the image segmentation and the asymptotic 3D PU in Section "Methodology". Real data description from Ref. [14], results of Image segmentation and on 3D reconstruction are shown in Section "Results and analysis". The conclusions are presented in Section "Conclusions".
METHODOLOGY
The flowchart of the proposed segmentationaided 3D reconstruction is shown in Figure 1. First, the SAR amplitude image is used for the segmentation. Setting the number of categories for classification and the number of iterations, Kmeans clustering is used for the initialization [15]. Then, the posterior probability of pixels is obtained by combining Markov random field[12] and Bayesian formulation [16]. The labels of all pixels are updated iteratively. The segmented image is used to refine the unwrapping network using the classified labels. Finally, 3D reconstruction of airborne array InSAR images is conducted using the refined network.
Figure 1 Flowchart of the proposed segmentationaided 3D reconstruction. 
Image segmentation
Let X_{s} represent the amplitude value of pixels in the SAR image, and Y_{s} represent the classification label of each pixel in the SAR image. For a binary classification question, Y_{s} belongs to {1, 2}. For preprocessing of the image, according to the amplitude of each pixel X_{s}, using Kmeans clustering algorithm[15] to obtain the classification label of each pixel X_{s}.
The image segmentation problem can be considered as the problem of estimating the optimal label ${\widehat{Y}}_{s}$ with the given amplitude X_{s} by the ExpectationMaximization algorithm (EM)[17], which is defined as follows:${\widehat{Y}}_{s}=\underset{{Y}_{s}}{\mathrm{arg}\mathrm{max}}P\left({Y}_{s}{X}_{s}\right)\mathrm{.}$(1)According to Bayesian formula P(Y_{s}X_{s}) as follows[16]:$P\left({Y}_{s}{X}_{s}\right)=\frac{P\left({X}_{s}{Y}_{S}\right)\cdot P\left({Y}_{s}\right)}{P\left({X}_{s}\right)}\mathrm{.}$(2)The P(X_{s}) can be obtained through observation. However, the P(Y_{s}) and P(XsY_{S}) need the following method to solve it.
The P(Y_{s}) can be solved by analyzing the relationship between pixels of SAR images. The Markov random field[12] is an efficient way to model the spatial relationships between pixels. Regard the SAR image as a collection of discrete pixels, for each pixel in SAR images, the relationship between pixels can be described by a neighborhood system[18], which is shown in Figure 2. Conventional neighborhood system includes two types. One is the firstorder neighborhood system, where pixels are surrounded by four surrounding neighborhoods, as shown in Figure 2(a). The other one is a secondorder neighborhood system, where pixels are surrounded by eight surrounding neighborhoods, as shown in Figure 2(b).
Figure 2 Neighborhood system. (a) Firstorder neighborhood system; (b) secondorder neighborhood system. 
For preprocessed images by Kmeans clustering algorithm, using a secondorder neighborhood system to analyze the relationship between pixels. In a secondorder neighborhood system, each pixel has the label Y_{s}(Y_{s} belongs to {1, 2}). By counting the number of pixels with label values Y_{s} equal to 1 and Y_{s} equal to 2 in the neighborhood system, the prior probability of the central pixel in the neighborhood system can be obtained, as follows:$P\left({Y}_{s}=1\right)=\frac{{N}_{1}}{N}$(3)$P\left({Y}_{s}=2\right)=\frac{{N}_{2}}{N}$(4)where N_{1} and N_{2} represent the number of pixels with Y_{s} equals to 1 and Y_{s} equals to 2 in the neighborhood, respectively. N represents the number of pixels in the neighborhood, excluding the center pixel.
Assuming that each pixel in the SAR image is statistically independent and satisfies a Gaussian distribution. Divide all pixels with Y_{s} equal to 1 and Y_{s} equal to 2 in the SAR image, and extract information on all pixels X_{s} for each class to construct a function, the P(X_{s}Y_{S}) can be solved as follows[19]:$P\left({X}_{S}{Y}_{S}\right)=\frac{1}{\sqrt{2\pi \sigma}}\cdot {e}^{\frac{{\left(x\mu \right)}^{2}}{2{\sigma}^{2}}}\mathrm{.}$(5)
Through this approach, the P(X_{s}Y_{S}=1) and P(X_{s}Y_{S}=2) are obtained, and then the posterior probability P(Y_{s}X_{s}) is calculated for each pixel in the image. Using the maximum P(Y_{s}X_{s}) of each pixel for image iteration until get the optimal label ${\widehat{Y}}_{s}$, and the image segmentation part is completed.
3D PU
Following the conventional InSAR process, the flat earth effect can be well estimated using the SAR system’s geometric parameters. After removing this phase, 3D PU can be applied to the flattened phase. Note that the performance of 3D PU depends on the probability of phase discontinuity. In the multibaseline SAR interferograms, such an assumption is evaluated using the relationship between the interferograms with varying baseline lengths. Inappropriate baseline diversity will increase the probability of phase discontinuity, leading to unwrapping errors. Here, the segmented image is used to refine the processing by considering the labels during the phase unwrapping. Only pixels with the same labels will be unwrapped in the same network.
The spatial PU is the integration of the phase gradient in the spatial domain. Its reliability depends on the selection of the phase difference. However, since the phase difference is always smaller than π, it is impossible to check whether the arc satisfies the PC assumption using a single baseline interferogram. Alternatively, using the smoothing criteria with a pair of interferogram [7], the reliability of the spatial PU can be improved.
During the process of 3D PU, the spatial PU is conducted on the interferogram with the shortest baseline, and the initial height estimation can be obtained as follows [20]:${\widehat{h}}_{\text{S}}=\frac{\lambda}{4\text{\pi}}\frac{R\mathrm{sin}{\theta}_{\text{inc}}}{{B}_{\perp \mathrm{,}\text{S}}}{\varphi}_{\text{S}}\mathrm{,}$(6)where φ_{S} denotes the unwrapped phase of the interferogram with the shortest baseline. R is the slant range. B_{perp} denotes the equivalent perpendicular baseline. θ_{inc} is the local incident angle and λ is the radar wavelength.
Using the initial height, the spatial PU of the long baseline interferogram follows the interferogram pair PU approach [7]. Supposing that the unwrapped phase difference between two scatterers p and q is $\Delta {\varphi}_{p\mathrm{,}q}$, denoted as an arc, the pseudo unwrapping phase and the corresponding ambiguity of this arc can be estimated as follows:$\Delta {\varphi}_{p\mathrm{,}q\mathrm{,}\text{pseudo}\mathrm{,}\text{L}}=\frac{4\text{\pi}}{\lambda}\frac{{B}_{\perp \mathrm{,}\text{L}}}{R\mathrm{sin}{\theta}_{\text{inc}}}\Delta {\widehat{h}}_{p\mathrm{,}q\mathrm{,}\text{S}}=\frac{{B}_{\perp \mathrm{,}\text{L}}}{{B}_{\perp \mathrm{,}\text{S}}}\Delta {\varphi}_{p\mathrm{,}q\mathrm{,}\text{S}}\mathrm{,}$(7)${\widehat{n}}_{p\mathrm{,}q\mathrm{,}\text{L}\mathrm{,0}}=\frac{\Delta {\varphi}_{p\mathrm{,}q\mathrm{,}\text{pseudo,L}}\Delta {\phi}_{p\mathrm{,}q\mathrm{,}\text{L}}}{2\text{\pi}}\mathrm{.}$(8)${\widehat{n}}_{p\mathrm{,}q\mathrm{,}\text{L}}=\text{round}\left\{{\widehat{n}}_{p\mathrm{,}q\mathrm{,}\text{L}\mathrm{,0}}\right\}\mathrm{,}$(9)where $\Delta {\phi}_{p\mathrm{,}q\mathrm{,}\text{L}}$ denotes the unwrapped phase difference between two scatterers p and q in a long baseline interferogram. round{·} denotes the rounding operation. Whether a long baseline interferogram can be successfully unwrapped depends on the uncertainty of the estimated phase ambiguity, denoted as SU criteria, which is defined as follows:${\widehat{n}}_{p\mathrm{,}q\mathrm{,}\text{L}\mathrm{,0}}{\widehat{n}}_{p\mathrm{,}q\mathrm{,}\text{L,true}}<\frac{1}{2}\mathrm{.}$(10)
In practical scenarios, the aforementioned equation can rewritten as follows:$\frac{{B}_{\perp \mathrm{,}\text{L}}^{2}}{{B}_{\perp \mathrm{,}\text{S}}^{2}}{\sigma}_{{\varphi}_{\text{S}}}^{2}+{\sigma}_{{\varphi}_{\text{L}}}^{2}<{\left(\frac{\pi}{{u}_{a}}\right)}^{2}\mathrm{,}$(11)where ${\sigma}_{{\varphi}_{\text{S}}}^{2}$ and ${\sigma}_{{\varphi}_{\text{L}}}^{2}$ denote the phase variances of the short baseline interferogram and the long baseline interferogram, respectively. Note that the multibaseline interferograms satisfied this SU criteria and can be regarded as the optimal baseline design.
If ${\widehat{n}}_{p\mathrm{,}q\mathrm{,}\text{L}\mathrm{,0}}$ satisfies the SU criteria, the unwrapped phase of the arc can be calculated as$\Delta {\varphi}_{p\mathrm{,}q\mathrm{,}\text{L}}=\Delta {\phi}_{p\mathrm{,}q\mathrm{,}\text{L}}+2\text{\pi}\cdot {\widehat{n}}_{p\mathrm{,}q\mathrm{,}\text{L}}\mathrm{.}$(12)
Using the same reference point during the spatial PU of the short baseline interferogram, the unwrapped phase of the long baseline interferogram is obtained by the spatial PU denoted as φ_{L}. Note that only arcs satisfying the SU criteria are used in the spatial PU.
With the unwrapped phase φ_{L}, the final height can be estimated using the follows [20]:${\widehat{h}}_{\text{L}}=\frac{\lambda}{4}\frac{R\text{}\mathrm{sin}\text{}{\theta}_{\text{inc}}}{{B}_{\perp \mathrm{,}\text{L}}}{\varphi}_{L}\mathrm{.}$(13)
RESULTS AND ANALYSIS
In this section, we show the results obtained by the proposed methodology using the real airborne arrayInSAR images[14] and the main parameters of the SAR system are list in Table 1.
Main parameters of the system
The following is a brief description of the experimental region, along with the results of image segmentation. Finally, the segmentation results are used for phase unwrapping and achieving 3D reconstruction.
Experimental region
The dataset covers the main region around the Lingang Business Building in Tianjin, China. The optical images are shown in Figure 3. The building consists of two 15story towers and a twostory podium, and covers an area of 22905.6 m^{2}, with a height of 69.2 m and a spacing of 57.6 m. The podium is 12.4 m high, 105.6 m long, and 72.4 m wide. The structure of the building is representative, and the surroundings are relatively open[14].
Figure 3 Optical images and photos of the investigated area (Source: Google Earth, Date: December 8, 2023). 
Results of image segmentation
Using the proposed MRF segmentation approach, Figure 4 shows the original image and the corresponding segmentation results.
Figure 4 Original airborne arrayInSAR image and segmented images. (a) Original airborne arrayInSAR image; (b) image of segmentation based on threshold; (c) image of segmentation based on MRF. 
The Figure 4(b) shows a simple segmentation by the threshold, which can roughly separate buildings and background areas. However, the threshold segmentation does not fully utilize the relationships between pixels in the image and thus shows many natural targets in the final segmentation.
Figure 4(c) shows the segmentation by the proposed algorithm, which can better utilize the relationships between pixels in the image to get a more accurate result of segmentation. Finally, the image is divided into two categories: the white mask is buildings, while the dark one denotes the low reflection areas such as trees and grasslands.
In the conventional unwrapping network, the arcs are generated using the neighboring principle. The rapid phase change between the building and the ground may lead to unwrapping errors. Using the segmented image, the generation of the unwrapping network can be optimized by removing arcs formed by different types of scatterers, as shown in Figure 5. Figure 5(a) shows the conventional unwrapping network, the connection between building and road will induce the probability of rapid phase change, i.e., unwrapping errors. In the proposed optimal network, the individual network within each class will significantly reduce the unwrapping errors.
Figure 5 Delaunay triangulation of a local image. (a) Connection of buildings and roads; (b) separation of buildings and roads. 
Results of 3D reconstruction
Since the baseline design of this dataset satisfied the SU criteria, it is possible to achieve a good success rate of phase unwrapping using the asymptotic 3D phase unwrapping method. Figure 6 shows the 3D reconstruction results by conventional 3D PU and segmented imageassisted 3D PU, respectively. In Figure 6(a), we present the results of a 3D reconstruction process conducted without the use of image segmentation techniques. It is evident that the basic reconstruction of the experimental area has been completed, with rough depictions of the heights and geographical positions of the buildings. The experimental area is divided into two parts during the reconstruction process: the highrise building area and the surrounding lowrise buildings. By analyzing the highrise section in the experimental area, we compared the color of the reconstructed top part of the highrise section with the color bar on the right side of the image and eventually concluded that the height of the highrise building is approximately between 60 to 70 m, which is consistent with the actual height of the highrise buildings, demonstrating a certain level of accuracy in the reconstruction results. However, it should be noted that in the actual experimental area, there are two highrise buildings. However, in Figure 6(a), only one of them is reconstructed, while the other one is not. This situation may be the result of interference between ground and building information. To address this limitation, image segmentation was employed as an aid in Figure 6(b) for 3D reconstruction. This improvement not only preserved the precise reconstruction of building height as shown in Figure 6(a), but also overcame the limitation of incomplete reconstruction of the experimental area in Figure 6(a), significantly enhancing the accuracy of the 3D reconstruction results.
Figure 6 Results of 3D reconstruction. (a) Without the image segmentation; (b) with the image segmentation. 
CONCLUSIONS
This study proposes an image segmentationaided 3D reconstruction for 3D imaging of array InSAR data. By using MRF and Bayes, the segmentation of SAR image is conducted with both good accuracy and computing efficiency. The segmentation algorithm fully considers the relationship between pixels in the image and their surrounding pixels, and can better utilize spatial correlation and statistical features in the image, resulting in quite good segmentation results. Under the constraint of segmentation results, 3D phase unwrapping can avoid phase jumps and effectively solve the problem of elevation ambiguity. The method of segmentationaided phase unwrapping for 3D reconstruction fully utilizes image characteristics to achieve reliable phase unwrapping, which effectively improves the practical application of unmanned aerial vehicle (UAV) array InSAR systems. This work only considers isolated buildings. Future work will focus on more precise segmentation for more complex 3D reconstruction with height ambiguity.
Data availability
The original data are available from the corresponding authors upon reasonable request.
Author contributions
X.Z. conducted the research, analyzed the data, wrote the manuscript, and revised it. F.H. conceived and designed the research, collected and interpreted data, wrote the manuscript, and revised it. F.X. conceived the idea, provided expertise in data analysis, and critically revised the manuscript for important intellectual content.
Conflict of interest
The authors declare no conflict of interest.
References
 Moreira A, PratsIraola P, Younis M, et al. A tutorial on synthetic aperture radar. IEEE Geosci Remote Sens Mag 2013; 1: 643. [Article] [Google Scholar]
 Massonnet D, Feigl KL. Radar interferometry and its application to changes in the Earth’s surface. Rev Geophys 1998; 36: 441500. [Article] [Google Scholar]
 Xing M, Pascazio V, Yu H. A Special Issue on Synthetic Aperture Radar Interferometry [From the Guest Editors]. IEEE Geosci Remote Sens Mag 2020; 8: 67. [Article] [Google Scholar]
 Rosen PA, Hensley S, Joughin IR, et al. Synthetic aperture radar interferometry. Proc IEEE 2000; 88: 333382. [Article] [Google Scholar]
 Yu H, Xing M, Yuan Z. Baseline design for multibaseline InSAR system: A review. IEEE J Miniat Air Space Syst 2021; 2: 1724. [Article] [Google Scholar]
 Zebker HA, Goldstein RM. Topographic mapping from interferometric synthetic aperture radar observations. J Geophys Res 1986; 91: 49934999. [Article] [Google Scholar]
 Hu FM, Wang F, Yu HW, et al. Asymptotic 3D phase unwrapping for very sparse airborne array insar images. IEEE Trans Geosci Remote Sens 2022; 60: 115. [Google Scholar]
 Cai H, Yang Z, Cao X, et al. A new iterative triclass thresholding technique in image segmentation. IEEE Trans Image Process 2014; 23: 10381046. [Article] [Google Scholar]
 Ziou D, Tabbone S. Edge detection techniques—An overview. Распознавание образов и анализ изображен/Pattern Recognition and Image Analysis: Advances in Mathematical Theory and Applications 1998; 8: 537559. [Google Scholar]
 Pappas TN, Jayant NS. An adaptive clustering algorithm for image segmentation. In: Proceedings of the International Conference on Acoustics, Speech, and Signal Processing. IEEE, 1989. 16671670. [Google Scholar]
 Abinash S, Pattnaik S. An experimental study of graph theorybased image segmentation techniques to improve graph cut algorithm. AIP Conf Proc 2022; 2481: 020001. [NASA ADS] [CrossRef] [Google Scholar]
 Omati M, Sahebi MR. Markov random fieldbased segmentation algorithm for detection of land cover changes using uninhabited aerial vehicle synthetic aperture radar polarimetric images. Int J Comput Inf Eng 2017; 11: 873877. [Google Scholar]
 Chen X, Zheng C, Yao H, et al. Image segmentation using a unified Markov random field model. IET Image Process 2017; 11: 860869. [Article] [Google Scholar]
 Qiu XL, Jiao ZK, Yang ZL. Key technology and preliminary progress of microwave vision 3D SAR experimental system. J Radars 2022; 11: 119. [Google Scholar]
 Likas A, Vlassis N, Verbeek JJ. The global kmeans clustering algorithm. Pattern Recognit 2003; 36: 451461. [Article] [Google Scholar]
 Zhang H. The optimality of naive Bayes. In: Proceedings of the Seventeenth International Florida Artificial Intelligence Research Society Conference (FLAIRS 2004). 2004. [Google Scholar]
 Moon TK. The expectationmaximization algorithm. IEEE Signal Process Mag 1996; 13: 4760. [Article] [Google Scholar]
 Zhong P. Image segmentation based on Markov random fields with adaptive neighborhood systems. Opt Eng 2006; 45: 097202. [Article] [Google Scholar]
 Ribeiro MI. Gaussian probability density functions: Properties and error characterization. Institute for Systems and Robotics, Lisboa, 2004. [Google Scholar]
 Hanssen RF. Radar Interferometry: Data Interpretation and Error Analysis. Dissertation for Doctoral Degree. Delft: Delft University of Technology, 2001. [Google Scholar]
All Tables
All Figures
Figure 1 Flowchart of the proposed segmentationaided 3D reconstruction. 

In the text 
Figure 2 Neighborhood system. (a) Firstorder neighborhood system; (b) secondorder neighborhood system. 

In the text 
Figure 3 Optical images and photos of the investigated area (Source: Google Earth, Date: December 8, 2023). 

In the text 
Figure 4 Original airborne arrayInSAR image and segmented images. (a) Original airborne arrayInSAR image; (b) image of segmentation based on threshold; (c) image of segmentation based on MRF. 

In the text 
Figure 5 Delaunay triangulation of a local image. (a) Connection of buildings and roads; (b) separation of buildings and roads. 

In the text 
Figure 6 Results of 3D reconstruction. (a) Without the image segmentation; (b) with the image segmentation. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.