SAR Image Change Detection Algorithm Based on Different Empirical Mode Decomposition

Show more

1. Introduction

Synthetic aperture radar (SAR) imaging has a typical characteristic of obtaining data in all-weather and all-time, so it has been widely used. But SAR is a microwave imaging remote sensing, so the received electromagnetic echo signal is a typical nonlinear and non-stationary signal, and there are some difficulties with traditional methods to process SAR images. Although the wavelet transform and multi-scale geometric analysis theory are very good modern signal processing methods, they are related to the choice of parameters and basic functions in signal processing process. Different parameters and basic functions will bring some influences on the final results. Empirical mode decomposition (EMD) is a kind of completely data-driven and adaptively multi-scale transform theory. It does not need to choose the basic function, but it can adaptively adjust according to the characteristics of the signal itself. Therefore, the EMD theory is better than the wavelet transform and the multi-scale geometric analysis theory; in other words, it can make up for the shortcomings of these typical multi-scale transform theories. EMD can not only deal with linear and stationary signals, but also can deal with nonlinear and non-stationary signals. Especially in non-stationary and nonlinear signal processing, it has great advantages and has been widely used in many fields [1] - [6] .

The definition of empirical mode decomposition and the related theory were proposed by Huang in 1998 [7] . EMD theory can decompose the signal into a set of linear intrinsic mode function (IMF) according to the characteristics of the signal itself, and IMF has its physically meaning. Since the EMD theory has been proposed, although it is not very mature, the theory has been improved and applied, and many fruitful results has been achieved. For example, aiming at the aliasing phenomenon in EMD decomposition, Wu et al. proposed the theory of ensemble empirical mode decomposition (EEMD) in 2009 [8] . It can not only solve the phenomenon of mode aliasing, but also can restrain the phenomenon of boundary point flying caused by too few points. The principle of EEMD is to use the characteristic that EMD has the approximate two enter filter banks to Gauss white noise. Due to the addition of Gauss noise in signal processing process in EEMD, some residual noises are introduced. In order to overcome this phenomenon, Yeh et al. proposed the complementary ensemble empirical mode decomposition (CEEMD) theory [9] , which can effectively eliminate the noise, and can make the further improvement to the one-dimensional EMD decomposition theory. When the one-dimensional EMD processes the two-di- mensional image signal, it not only does not consider the pixel spatial relationship, but also the calculation speed is very slow. So Billing et al. proposed the theory of bidimensional empirical mode decomposition (BEMD) [10] . Because BEMD directly decomposes the two-dimensional image signal, considering the spatial relationship of pixel, it has been widely applied in remote sensing image processing [11] [12] . Based on the BEMD development, Rehman and Mandic proposed the trivariate empirical mode decomposition (TEMD) [13] and the multivariate empirical mode decomposition (MEMD) [14] in 2010.

The SAR image is the mapping from the feature space to the image space. SAR receives the electromagnetic echo is scattered from the ground objects, and then after a series of complex processing, then the SAR image is obtained. The echo signal received by SAR is a typical non-stationary and nonlinear signal; meanwhile, the SAR imaging mechanism is the coherent imaging, so the SAR image contains a lot of speckle noise. However, SAR is a microwave imaging which can obtain target area data under all-weather and all-time, and is very suitable for dynamic monitoring and change information acquisition. According to the characteristics of SAR images, after studying the principle of EEMD and BEMD, this paper proposed a new fusion change detection algorithm for SAR images based on different empirical mode decompositions, and it is referred to FCD-EMD algorithm. The feasibility of FCD-EMD algorithm is tested with SAR image data, and good experimental results are obtained.

2. Empirical Mode Decomposition Theory

2.1. EEMD Theory

Using EMD theory to decompose an input signal, we can obtain a finite number of intrinsic mode function (IMF) and a residual quantity. Whether the input signal is linear, nonlinear or stationary, non-stationary, the decomposed IMF features are linear and stationary.

For an input signal $x\left(n\right)$ , after processing of EMD decomposition, it can be expressed as the sum of the modal characteristic component ${C}_{i}\left(n\right)$ and the residua $R\left(n\right)$ , and definition is given in Equation (1).

$x\left(n\right)={\displaystyle {\sum}_{i=1}^{N}{C}_{i}\left(n\right)}+R\left(n\right)$ (1)

Here
${C}_{i}\left(n\right)$ is the i^{th} IMF component, N is the total number of IMF components.

Although EMD is a fully adaptive data-driven and multi-scale decomposition theory, it is easy to produce the aliasing when there is interference signal in the signal to be decomposed. EEMD adds white noise into the original signal, uses the characteristic that the mean value of white noise is zero, and makes the ensemble average to the IMF feature components obtained by EMD decomposition, which can well eliminate the aliasing. The realization process is as follows.

1) Generate N stripe white noise which is as same as the equal length with the input original signal $x\left(n\right)$ , the white noise is added into $x\left(n\right)$ , and the signal which is added by the noise is ${x}_{i}\left(n\right)$ .

${x}_{i}\left(n\right)=x\left(n\right)+{h}_{i}\left(n\right),\text{\hspace{0.17em}}i=1,2,\cdots ,N$ (2)

Here
${x}_{i}\left(n\right)$ the signal which is added by white noise after the i^{th} time,
${h}_{i}\left(n\right)$ is the white noise that is added by the i^{th} time.

2) All
${x}_{i}\left(n\right)$ is decomposed by EMD, respectively. Finally, M number of
${C}_{ij}\left(n\right)$ and one Residual
${R}_{i}\left(n\right)$ is got. Here
${C}_{ij}\left(n\right)$ denotes the j^{th} IMF component obtain by adding white noise after the i^{th} time,
$i=1,2,\cdots ,N$ ,
$j=1,2,\cdots ,M$ .

3) ${C}_{ij}\left(n\right)$ and ${R}_{i}\left(n\right)$ perform the ensemble average processing operation, respectively. When N is large enough, the sum of the IMF components which has been added by white noise will tend to zero. The results of the ensemble average are shown as following.

${C}_{j}\left(n\right)=\frac{1}{N}{\displaystyle {\sum}_{i=1}^{N}{C}_{ij}\left(n\right)}$ (3)

$R\left(n\right)=\frac{1}{N}{\displaystyle {\sum}_{i=1}^{N}{R}_{i}\left(n\right)}$ (4)

2.2. BEMD Theory

The decomposition process and idea of BEMD is the same as that of EMD decomposition method. A complex two-dimensional signal is decomposed into a number of bidimensional intrinsic mode function (BIMF) components ${C}_{i}^{\text{B}}\left(m,n\right)$ and a residual ${R}^{\text{B}}\left(m,n\right)$ by the screening method. Every BIMF component is independent of each other, including the local feature signals of the original signal of different frequencies and different spatial.

Next the two-dimensional image signal is as an example to introduce the screening process of BEMD method. $I\left(m,n\right)$ denotes the two-dimensional image signal that will be processed, ${I}_{\mathrm{max}}\left(m,n\right)$ denotes the maximum envelope surface, ${I}_{\mathrm{min}}\left(m,n\right)$ denotes the minimum envelope surface, and the screening steps are as follows.

1) Finding the local maxima and local minima of the curved surface in an image.

2) Using the interpolation method perform the interpolation operation on the extreme points, and the maximum and minimum envelope surfaces ${I}_{\mathrm{max}}\left(m,n\right)$ and ${I}_{\mathrm{min}}\left(m,n\right)$ are obtained.

Computing the mean of the envelope surface ${I}_{\text{mean}}\left(m,n\right)$ , and ${I}_{\text{mean}}\left(m,n\right)=\left[{I}_{\mathrm{max}}\left(m,n\right)+{I}_{\mathrm{min}}\left(m,n\right)\right]/2$ ; extracting high frequency detail information, i.e. $d\left(m,n\right)=I\left(m,n\right)-{I}_{\text{mean}}\left(m,n\right)$ .

4) Determine whether $d\left(m,n\right)$ satisfies the two conditions of BIMF. If it does not satisfy, sign $I\left(m,n\right)=d\left(m,n\right)$ , repeat steps form step (1) to step (4), until the conditions is satisfied. At this time, $d\left(m,n\right)$ is a BIMF component, sign ${C}_{1}^{\text{B}}\left(m,n\right)$ , ${C}_{1}^{\text{B}}\left(m,n\right)=d\left(m,n\right)$ .

5) Sign ${R}_{1}^{\text{B}}\left(m,n\right)=I\left(m,n\right)-{C}_{1}^{\text{B}}\left(m,n\right)$ , and it is as a new signal to be decomposed, repeat the operations from step (1) to step (4) until the second BIMF component ${C}_{2}^{\text{B}}\left(m,n\right)$ is obtained. Here the residual is ${R}_{2}^{\text{B}}\left(m,n\right)={R}_{1}^{\text{B}}\left(m,n\right)-{C}_{2}^{\text{B}}\left(m,n\right)$ . Repeat the above steps, until the remainder ${R}_{i}^{\text{B}}\left(m,n\right)$ has only one extreme point surface, and no longer can any BIMF components be extracted from it.

The image $I\left(m,n\right)$ can be expressed as a sum of several BIMF components and the residual, i.e. Equation (5).

$I\left(m,n\right)={\displaystyle {\sum}_{i=1}^{M}{C}_{i}^{\text{B}}\left(m,n\right)}+{R}^{\text{B}}\left(m,n\right)$ (5)

Here $M$ is the number of BIMF feature components.

In the step (4) to determine whether the screening process is over, which see if it satisfies the two conditions in the definition of BIMF. The standard deviation (SD) of the two screening results is usually used to determine it. And the definition is given as following.

$\text{SD}={\displaystyle {\sum}_{i=0}^{M}\left[\frac{{\left|{d}_{i-1}\left(m,n\right)-{d}_{i}\left(m,n\right)\right|}^{2}}{{d}_{i-1}^{2}\left(m,n\right)}\right]}$ (6)

Here the value of SD is generally between 0.2 and 0.3, and
${d}_{i}\left(m,n\right)$ denotes the i^{th} residual.

3. FCD-EMD Algorithm Description

When the EMD method is used to process an image, if one wants to extract the much IMF feature information, the original image should have the better direction information. For the SAR images with artificial targets, the direction information is usually obvious. BEMD method is suitable for the processing of SAR images, not only can get abundant BIMF feature information, but also keep the good spatial information of image pixels. Therefore, if the IMF feature and BIMF feature perform the fusion processing, it not only can get more abundant spatial information and direction information, but also is very helpful for target detection and change information obtaining. According to that the SAR imaging mechanism is sensitive to the azimuth; this paper proposes a new SAR image target change detection algorithm based on two different types of EMD decomposition theory, i.e. the FCD-EMD algorithm and its flow chart frame is shown in Figure 1. It can be seen from Figure 1; the proposed change detection algorithm consists of two parts. One is the change detection based on EEMD, which uses EEMD to obtain the IMF feature images and to get the change information; another part is to use BEMD obtaining the BIMF feature images and to get the change information. Finally, the two parts of the change detection results are fused and the final change information of the target area is got. The change detection process based on EEMD and BEMD is similar, therefore, the following detailed description of the change detection based on EEMD is introduced. The change detection based on BEMD can refer to the implementation of the EEMD method.

1) Input different temporal SAR images, i.e. input SAR images ${I}_{1}$ and ${I}_{2}$ which were obtained from the same target area but at different time ${t}_{1}$ and ${t}_{2}$ .

2) SAR images ${I}_{1}$ and ${I}_{2}$ performs the EEMD decomposition operations in horizontal and vertical direction, respectively.

3) The IMF feature images of horizontal and vertical directions on the same scale are fused to obtain the total IMF feature image for every scale, which are the parts including the frequency information and the detailed information.

4) Select the appropriate IMF feature components. In FCD-EMD algorithm, the first to third scale feature images are selected to perform the next fusion processing operation.

5) The selected IMF feature images are fused to obtain the fused feature images FIMF1 and FIMF2, respectively.

6) The corresponding pixels of feature images of FIMF1 and FIMF2 perform the subtraction operation, and the feature difference image FIMF-D1 will be got.

Figure 1. The block diagram of FCD-EMD algorithm.

7) The expectation maximization (EM) algorithm is used to generate an unsupervised detection threshold ${T}_{1}$ in the feature difference image FIMF-D1. EM algorithm is a maximum likelihood estimation method for incomplete data problems, it does not need any external data or prior knowledge, and the detailed implementation process can refer to the document [15] .

8) Obtain the result image ${D}_{1}$ of change detection with EEMD. Using the obtained detection threshold ${T}_{1}$ performs the detection and determination operations in the feature difference image. If ${D}_{1}\left(m,n\right)\ge {T}_{1}$ , this shows that the change happens; If ${D}_{1}\left(m,n\right)<{T}_{1}$ , this shows that there is no changes. Here ${D}_{1}\left(m,n\right)$ denotes the gray value of the feature difference image at $\left(m,n\right)$ pixel.

9) Obtain the result image ${D}_{2}$ of change detection with BEMD. The change detection process based on BEMD is similar to EEMD change detection, and it also includes the above several steps. The difference is that BEMD directly decomposes the SAR image, and the two-dimensional feature images BIMF1 and BIMF2 are obtained from the SAR image ${I}_{1}$ and ${I}_{2}$ , respectively. Then the feature selection and fusion processing is performed to generate the feature difference image and the detection threshold; finally, the change detection result ${D}_{2}$ is obtained. After obtaining ${D}_{1}$ and ${D}_{2}$ , they perform the fusion operation, and the change information between SAR images at ${t}_{1}$ and ${t}_{2}$ can be obtained.

10) Fuse the results obtained by the two different methods, and the fusion rule is given by

$D\left(m,n\right)=\alpha \times {D}_{1}\left(m,n\right)+\beta \times {D}_{2}\left(m,n\right)$ (7)

Here $D\left(m,n\right)$ denotes the final result, $D\left(m,n\right)$ is the result with EEMD and ${D}_{2}\left(m,n\right)$ is the result with BEMD. The parameters $\alpha $ and $\beta $ are their weight value, respectively, and the value ranges of them is $\left[0,1\right]$ , and here $\alpha =\beta =0.5$ .

4. Experimental Results and Analysis

In order to verify the feasibility of the proposed algorithm, the simulation experiments and the verification test of measured data are carried out. The experimental SAR image data came from the Sandia National Laboratories [16] and the Canadian RADARSAT satellite. For SAR image data from the Sandia, part of them performs the simulation experiments, which are shown in Figure 2. The targets are tanks and the background is bushes in the SAR images. Here Figure 2(A) and Figure 2(D) are the original SAR images. Assume that Figure 2(A) is the SAR image before change, i.e. the SAR image ${I}_{1}$ at time ${t}_{1}$ ; Figure 2(D) is the image after change, that is, the SAR image ${I}_{2}$ at time ${t}_{2}$ . They are IMF features and BIMF features of Figure 2(A) in Figure 2(B) and Figure 2(C), respectively. Similarly, Figure 2(E) and Figure 2(F) are IMF features and BIMF features of Figure 2(D), respectively. What Figure 3 shows is the change detection results obtained by different empirical mode decomposition methods. But there are some changes, that is, the images which are used to extract the change information are no longer the original SAR images, but they are the feature images obtained by empirical mode decomposition. The change detection results are shown in Figure 3(a), which is based on one-dimensional EEMD. At this time, the feature images at ${t}_{1}$ and ${t}_{2}$ time are shown in Figure 2(B) and Figure 2(E), respectively. In Figure 3, the change shown in Figure 3(A) is that the SAR image at time ${t}_{1}$ is relative to the SAR image at time ${t}_{2}$ , i.e. the area of information weakening; the change of ${t}_{2}$ relative to ${t}_{1}$ is shown in Figure 3(B), i.e. the information enhancing area. Figure 3(C) are all change results, including both enhancement and weakening, namely the fusion of Figure 3(A) and Figure 3(B). The change detection results obtained by BEMD are shown in Figure 3(b). The image at ${t}_{1}$ is the feature image shown by Figure 2(C), and the image at ${t}_{2}$ is the feature image in Figure 2(F). The result of the FCD-EMD algorithm is shown in Figure 3(c). In fact, it is the fusion result of Figure 3(a) and Figure 3(b). It can be seen in Figure 3 that for those SAR images which is more obvious or their direction information is abundant, these detection methods all can get the target feature information and also can detect the change information, and the focus of each method is different. IMF features in EEMD

Figure 2. Simulation of SAR images and their IMF features. ((A), (D) is original SAR images; (B), (E) is IMF feature images; (C), (F) is BIMF feature images).

(a) Change detection results of EEMD

(b) Change detection results of BEMD

(c) Change detection results of FCD-EMD

Figure 3. Change detection results of different methods (t_{2} relative to t_{1}, (A) weakened area, (B) enhanced area, (C) fusion of (A) and (B)).

focus on the acquisition of direction information and geometric detail information; BIMF features in BEMD focus on the acquisition of the spatial information of the target (the correlation of neighboring regions). Through the fusion, the formation of complementary advantages, one can get a better change detection effect, which is shown in Figure 3. At the same time, it also shows that the fusion algorithm of two different EMD methods can be used to detect the change information for different temporal SAR images, and it is feasible and has wide application prospect and popularization value.

Figure 4 shows the change detection results which are obtained after the fusion of different IMF and BIMF features. The change information shown in Figure 4(A) is the SAR image at time ${t}_{1}$ relative to ${t}_{2}$ . Contrary, what Figure 4(B) shows is ${t}_{2}$ relative to ${t}_{1}$ . Figure 4(C) is the sum of Figure 4(A) and Figure 4(B). Figure 3 and Figure 4 show that the sequence of fusion of different methods will affect the final change detection results. The change detection results in Figure 3 are better than that in Figure 4, so the FCD-EMD algorithm proposed in this paper is the first to carry out the change detection and then to perform fusion operations.

The SAR images in Figure 5 are from the Canadian RADARSART satellite, the imaging area is Bengbu area, Anhui province, China. Here, the imaging time of the SAR image in Figure 5(A) was in July 2001, assuming that it is the SAR image at time ${t}_{1}$ , that is, it is the SAR image before change; the SAR image in Figure 5(D) obtained in July 2005, and assume that it is the SAR image at time ${t}_{2}$ , i.e. the SAR image after change. Figure 5(B) and Figure 5(C) are the IMF and BIMF feature information image of Figure 5(A), respectively; similarly, Figure 5(E) and Figure 5(F) are IMF and BIMF feature information image of Figure 5(D).

Figure 4. Results of at first carry out fusion and then detect change information (t_{2} relative to t_{1}, (A) weakened area, (B) enhanced area, (C) fusion of (A) and (B)).

Figure 5. Original SAR images and their IMF and BIMF feature images. ((A), (D) is original images; (B), (E) is IMF feature images; (C), (F) is BIMF feature images).

Figure 6. Change detection results of SAR images with different methods.

The results of change detection in Figure 5 are shown in Figure 6. What Figure 6(a) shows is obtained by Figure 5(B) and Figure 5(E), the results of Figure 5(C) and Figure 5(F) are shown in Figure 6(b), and the fused results of Figure 6(a) and Figure 6(b) are shown in Figure 6(c). Figure 6 further shows that the fusion FCD-EMD algorithm by different EMD theories can be used for real SAR image change detection, and it is feasible and effective. Because the IMF feature information is extracted by row or column, therefore, as long as there is azimuth information, it can acquire them; in other words, EEMD is more sensitive to the directional information, as can be seen in Figure 6(a). In this experiment, the change information based on the IMF feature is more than that of the BIMF feature, as shown in Figure 6(a) and Figure 6(b). However, there is more false alarm information in Figure 6(a). The reason is that IMF affected by speckle noise is greater than BIMF. In essence, another obvious feature of the FCD-EMD algorithm is that it can greatly reduce the effect of speckle noise. Because speckle noise mainly concentrates in the high-frequency part, when one selects the intrinsic mode function feature, the IMF feature at first scale may be removed without affecting the detection results, and the speckle noise is reduced. This merit is also reflected in Figure 5(C) and Figure 5(F).

5. Conclusion

Empirical mode decomposition is a fully adaptive data-driven signal processing theory. Although there are many kinds of decomposition methods, each method has its advantages and disadvantages. EEMD can obtain very good directional information, and BEMD can get the spatial information of the pixel. The proposed FCD-EMD algorithm fuses their advantages, holding good directional and spatial characteristics; at the same time, the new method can reduce the influence of speckle noise on obtaining change information with multi-temporal SAR images. So more accurate information will be obtained by the FCD-EMD algorithm that that of the single EMD method. The experimental results of simulation experiments and the measured SAR images show that the FCD-EMD algorithm can obtain good effects. This shows that the new proposed algorithm is feasible. In the future work, we will make a deep study of the quantitative analysis of the method.

Acknowledgements

The work was supported by Natural Science Foundation of China (No. 61379031, 41574008), Natural Science Basic Research Plan in Shaanxi Province of China (No. 2016JM6052) and Special Foundation for Special Talents of Xijing University (No. XJ17T04).

References

[1] Hao, H., Wang, H.L. and Rehman, N.U. (2017) A Joint Framework for Multivariate Signal Denoising Using Multivariate Empirical Mode Decomposition. Signal Processing, 135, 263-273.

https://doi.org/10.1016/j.sigpro.2017.01.022

[2] Camarena-Martinez, D., Valtierra-Rodriguez, M., Perez-Ramirez, C.A., et al. (2016) Novel Down Sampling Empirical Mode Decomposition Approach for Power Quality Analysis. IEEE Transactions on Industrial Electronics, 63, 2369-2378.

https://doi.org/10.1109/TIE.2015.2506619

[3] Mishra, V.K., Bajaj, V., Kumar, A., et al. (2016) Analysis of ALS and Normal EMG Signals Based on Empirical Mode Decomposition. IET Science, Measurement & Technology, 10, 963-971.

https://doi.org/10.1049/iet-smt.2016.0208

[4] Ren, Y., Suganthan, P.N. and Srikanth, N. (2016) A Novel Empirical Mode Decomposition with Support Vector Regression for Wind Speed Forecasting. IEEE Transactions on Neural Networks and Learning Systems, 27, 1793-1798.

https://doi.org/10.1109/TNNLS.2014.2351391

[5] Liu, G., Li, L., Gong, H., et al. (2016) Multisource Remote Sensing Imagery Fusion Scheme Based on Bidimensional Empirical Mode Decomposition and Its Application to the Extraction of Bamboo Forest. Remote Sensing, 9, 1-19.

https://doi.org/10.1080/01431161.2016.1253899

[6] Yang, M.D., Huang, K.S., Yang, Y.F., et al. (2016) Hyperspectral Image Classification Using Fast and Adaptive Bidimensional Empirical Mode Decomposition with Minimum Noise Fraction. IEEE Transactions on Geoscience and Remote Sensing Letters, 13, 1950-1954.

https://doi.org/10.1109/LGRS.2016.2618930

[7] Huang, N.E., Shen, Z., Long, S.R., et al. (1998) The Empirical Mode Decomposition and the Hilbert Spectrum for Nonlinear and Non-Stationary Time Series Analysis. Proceedings of the Royacl Society of London. Series A: Mathematical, Physical and Engineering Sciences, 454, 903-995.

https://doi.org/10.1098/rspa.1998.0193

[8] Wu, Z. and Huang, N.E. (2009) Ensemble Empirical Mode Decomposition: A Noise-Assisted Data Analysis Method. Advances in Adaptive Data Analysis, 1, 1-41.

https://doi.org/10.1142/S1793536909000047

[9] Yeh, J.R., Shieh, J.S. and Huang, N.E. (2010) Complementary Ensemble Empirical Mode Decomposition: A Novel Noise Enhanced Data Analysis Method. Advances in Adaptive Data Analysis, 2, 135-156.

https://doi.org/10.1142/S1793536910000422

[10] Billing, G., Flandrin, P., Goncalves, P., et al. (2007) Bivariate Empirical Bode Decomposition. IEEE Transactions on Signal Processing Letters, 14, 936-939.

https://doi.org/10.1109/LSP.2007.904710

[11] He, Z., Wang, Q., Shen, Y., et al. (2013) Multivariate Gray Model-Based BEMD for Hyperspectral Image Classification. IEEE Transactions on Instrumentation and Measurement, 62, 889-904.

https://doi.org/10.1109/TIM.2013.2246917

[12] Chen, Z., Luo, S., Xie, T., et al. (2014) A Navel Infrared Small Target Detection Method Based on BEMD and Local Inverse Entropy. Infrared Physics & Technology, 66, 114-124.

https://doi.org/10.1016/j.infrared.2014.05.013

[13] Rehman, N. and Mandic, D.P. (2010) Empirical Mode Decomposition for Trivariate Signals. IEEE Transactions on Signal Processing, 58, 1059-1068.

https://doi.org/10.1109/TSP.2009.2033730

[14] Rehman, N. and Mandic, D.P. (2010) Multivariate Empirical Mode Decomposition. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Science, 466, 1291-1302.

https://doi.org/10.1098/rspa.2009.0502

[15] Moon, T.K. (1996) The Expectation-Maximization Algorithm. Signal Processing, 6, 47-60.

https://doi.org/10.1109/79.543975

[16] Chavez, J. (1997) Pathfinder Airborne ISR and Synthetic Aperture Radar (SAR) Systems.

http://www.sandia.gov/radar/index.html