2-D Rayleigh Autoregressive Moving Average Model for SAR Image Modeling

Two-dimensional (2-D) autoregressive moving average (ARMA) models are commonly applied to describe real-world image data, usually assuming Gaussian or symmetric noise. However, real-world data often present non-Gaussian signals, with asymmetrical distributions and strictly positive values. In particular, SAR images are known to be well characterized by the Rayleigh distribution. In this context, the ARMA model tailored for 2-D Rayleigh-distributed data is introduced -- the 2-D RARMA model. The 2-D RARMA model is derived and conditional likelihood inferences are discussed. The proposed model was submitted to extensive Monte Carlo simulations to evaluate the performance of the conditional maximum likelihood estimators. Moreover, in the context of SAR image processing, two comprehensive numerical experiments were performed comparing anomaly detection and image modeling results of the proposed model with traditional 2-D ARMA models and competing methods in the literature.

An alternative approach for modeling such type of data is the use of the Rayleigh distribution [Bayer et al., 2020a].
This model is considered in signal and image processing [Zanetti et al., 2015, Sumaiya andKumari, 2018], being relevant in the context of synthetic aperture radar (SAR) image modeling, due to its good characterization of image pixel amplitude values [Yue et al., 2019, Kuruoglu and Zerubia, 2004, Jackson and Moses, 2009, Kuttikkad and Chellappa, 2000. In particular, this distribution is known to well fit SAR data homogeneous regions Quegan, 2004, Yue et al., 2021].
The main justification for supposing the Rayleigh distribution as a model for the amplitude of SAR data is that the assumption of a large number of reflectors in an observed image allows one to invoke the central limit theorem, according to which the distribution of the real and complex parts of the received signal are independent and approaches the Gaussian distribution with zero mean and constant variance [Kuruoglu andZerubia, 2004, Yue et al., 2021]. Thus, under these assumptions, the amplitude values of complex SAR data are exactly Rayleigh distributed.
Frequently, the SAR image modeling is performed assuming constant parameters [Kuruoglu andZerubia, 2004, Jackson andMoses, 2009]. In cases where this assumption is not suitable, an alternative is to use a regression model, where each observation has one specific estimated mean [Wang and Ouchi, 2008, McCullagh and Nelder, 1989, Palm et al., 2019. Motivated by these SAR image characteristics, Palm et al. [2019] proposed a regression model based on the Rayleigh distribution for SAR image modeling. However, image pixels usually present a resolution spatial dependence [Jackson and Moses, 2009, Yan et al., 2018, and consequently, a 2-D model can be used as a venue for addressing such a problem.
The use of one-or two-dimensional dynamical models are often employed in image applications. For example, in Almeida-Junior and Nascimento [2021], the authors proposed a new 1-D ARMA model considering the G 0 distribution to estimate the intensity values of SAR image pixels. Bayer et al. [2020a] derived an 1-D Rayleigh-based dynamical model useful to land-use classification in SAR images. In Bustos et al. [2009b], a 2-D Gaussian ARMA model was applied in image filtering schemes. Additionally, the literature presents several studies for 1-D dynamical models based on different distributions, useful in a multitude of scientific applications; see., e.g., Benjamin et al. [2003], Bayer et al. [2017], Rocha and Cribari-Neto [2009], Möller and Weiß [2020], Palm et al. [2021].
However, to the best of our knowledge, a two-dimensional ARMA model assuming the Rayleigh distribution is not present in the literature and this paper aims at proposing a first treatment on the topic. Our goal is two-fold. First, we derive a two-dimensional ARMA model for non-Gaussian spatial autocorrelated images, where the observed signal is asymmetric and strictly positive. For the proposed model, we introduce parameter estimation, large data record inference, and the quantile residuals. Second, we propose an image modeling tool based on the derived spatial model estimated parameters and an anomaly detector for non-Gaussian SAR images. Considering control charts of the proposed model residuals, the introduced detection scheme measures the deviations of an observed pixel value from its estimated mean value. The residual-based control charts have already been employed in remote sensing data change and anomaly detection [Bayer et al., 2020b, Brooks et al., 2013.
Anomaly detection is a popular field in signal processing, machine learning, and statistics [Talagala et al., 2020, Kadri et al., 2016, Quatrini et al., 2020, Kwon and Nasrabadi, 2005. In particular, anomaly detection in noisy image is explored for quality control purposes in different manufacturing applications, such as composites, steel, and textile production [Yan et al., 2018]. The anomaly detection problem can be addressed considering different techniques, depending on the way that anomalies are defined. For example, different types of expected outputs or input data have particular problem formulation and need to be addressed through specific techniques [Talagala et al., 2020]. On the other hand, our proposal is based on a simple residual analyze considering control chart with a fixed theoretical threshold, which is defined based on the residual distribution. This paper is organized as follows. In Section 2, we describe the proposed spatial model, provide conditional maximum likelihood estimators, and present a hypothesis testing methodology. Section 3 details the introduced image modeling tool and proposed an anomaly detection algorithm. Section 4 presents Monte Carlo simulations and two empirical analyses of the derived detector applied to SAR images. Section 5 brings final remarks and concludes the paper.

Mathematical Setup
Recently, a regression model [Palm et al., 2019] and an 1-D ARMA model [Bayer et al., 2020a] based on the Rayleigh distribution have been proposed. The Rayleigh ARMA (RARMA) model introduced in Bayer et al. [2020a] relates the mean of an one-dimensional discrete-time signal to a linear predictor through a strictly monotonic, twice differentiable link function g(·), where g : R + → R. The goal of this section is to extend the 1-D model presented in Bayer et al. [2020a] and introduce to the 2-D case.

The Model
The proposed 2-D Rayleigh autoregressive and moving average model, hereafter referred to as 2-D RARMA, is defined according to where the two-dimensional autoregressive operator Φ(z 1 ,z 2 ) and moving average operator Θ(z 1 ,z 2 ) are furnished, respectively, by and β ∈ R is an intercept; the quantities z i 1 g(y [n,m] As suggested in Basu and Reinsel [1993] for an unilateral spatial ARMA, we assume φ (0,0) = θ (0,0) = 0. Replacing the . . .
The Broyden-Fletcher-Goldfarb-Shanno (BFGS) method [Press et al., 1992] with analytic first derivatives was adopted as the nonlinear optimization algorithm [Nocedal and Wright, 1999] to solve (3). The BFGS method was selected due to its superior performance for non-linear optimization [Mittelhammer et al., 2000]. The initial values for the constant (β) and the autoregressive (φ) parameters were derived from the ordinary least squares estimate associated to the linear regression model, which has a closed matrix form. The response vector is and the covariate matrix is given by As suggested in Bayer et al. [2020a], we set θ = 0 as initial values.

Large Data Record Inference
Under some usual regularity conditions, the CMLE is consistent and asymptotically normally distributed [Andersen, 1970].
Thus, in large data records, we have that where d −→ represents convergence in distribution and N(0,I −1 (γ)) denotes the multivariate Gaussian distribution with null mean and covariance matrix I −1 (γ). The conditional Fisher information matrix, I(γ), is discussed in detail in the Appendix.
To derive a hypothesis testing methodology tailored for the 2-D RARMA model parameters, the likelihood-based detection theory [Pawitan, 2001, Kay, 1998] can be considered. Let γ be partitioned in a parameter vector of interest γ I of dimension ν, and a vector of nuisance parameters γ J of dimension [(p + 1) 2 + (q + 1) 2 − 1] − ν [Kay, 1998]. In addi-tion, H 0 : γ I = γ I 0 is the hypothesis of interest and H 1 : γ I = γ I 0 the alternative hypothesis, where γ I 0 is a fixed column vector of dimension ν. The Wald statistic is given by [Kay, 1998] where ) ⊤ is the CMLE under H 1 and I −1 ( γ 1 ) γ I γ I is a partition limited to the estimates of interest of the estimated conditional Fisher information matrix. Under H 0 , the test statistic, T W , asymptotically follows the chi-square distribution with ν degrees of freedom, χ 2 ν [Kay, 1998]. The hypothesis test consists of comparing the computed value of T W with a threshold value, ǫ, which is obtained based on the χ 2 ν distribution and the desired probability of false alarm [Kay, 1998].
To test the overall significance of a fitted model, we considered the following hypotheses where γ ⋆ = (φ ⊤ ,θ ⊤ ) ⊤ . Using the Wald test described above, we reject H 0 when T W > ǫ. In this situation, γ ⋆ = 0, indicating that at least some of the autoregressive and moving average parameters are nonzero and the spatial correlation among the pixels is significant.

Image Modeling and Anomaly Detection
In this section, we propose an image modeling and an anomaly detection tool based on the proposed 2-D RARMA model. We also discuss model selection strategies. For such, we introduce the estimated values of µ [n,m] and present the residuals of the 2-D RARMA model.

Image Modeling
The modeled image is obtained by applying the estimated values of µ [n,m], µ[n,m], in the 2-D RARMA(p,q) model structure, given by (1), and evaluating it at γ. Therefore, the fitted signal is given by where n = w + 1, w + 2,... , N and m = w + 1, w + 2,... ,M. Hence, similar to the 1-D model, the image border is not included in the modeling process, since the resulting fitted image has (N − w) × (M − w) pixels.

Anomaly Detector
Residuals can be useful for performing a diagnostic analysis of the fitted model and can be defined as a function of the observed and predicted values [Kedem and Fokianos, 2005]. We employ the quantile residuals [Dunn and Smyth, 1996], defined as where Φ −1 denotes the standard normal quantile function. When the model is correctly fitted, for construction, the quantile residuals are approximately Gaussian distributed with zero mean and unit variance, i.e., r[n,m] ∼ N (0,1) [Dunn and Smyth, 1996]. The residual analysis for other classes of 1-D non-Gaussian dynamical models is discussed in Kedem and Fokianos [2005] and Dunn and Smyth [1996].
Large values of r [n,m] can be interpreted as anomaly changes in the image behavior. To capture such variations of the residual values, we adopted the use of control charts. Since the residuals have approximately unitary variance, the control chart detects an image change if the residual value is outside the control limit ±L. We adopted L = 3, since it is expected that residuals are randomly distributed around zero and inside the interval [−3,3], about 99.7% of the observations, since Brooks et al., 2013. If the residual value is outside this range, the analyzed pixel is understood to differ from the expected behavior according to the 2-D RARMA model fitted in the region of interest and, consequently, some anomaly change might have occurred.
Notice that the proposed model relies on neighboring pixels from the northwest direction, as shown in Figure 1. Thus, to take into account the other directions in an omnidirectional manner, thus ensuring that all surrounding pixels are considered, the 2-D RARMA fitting is also applied to the 90 • , 180 • , and 270 • rotated region of interest to capture information from the versions of the southwest, southeast, and northeast directions. Results are combined according to the morphological union of the resulting binary images. Search for the best model for each rotation is computationally expensive.
Consequently, considering a trade-off between simplicity and efficacy of the anomaly detection method, we suggest using the same model order for the four directions.
To further increase the performance of the proposed detector, a post-processing step using mathematical morphological operations, such as erosion, dilation, opening, and closing operations, can be considered [Edmond et al., 2000, Gonzalez et al., 2009. Such operations aim at (i) removing small spurious pixel groups which are regarded as noise and (ii) preventing the splitting of the interest targets into multiple substructures [Gonzalez et al., 2009]. The resulting data is the detected image. The proposed ground type change detection method is summarized in Algorithm 1.
Algorithm 1 Anomaly detection method based on the 2-D RARMA(p,q) model Require: Interest image X input Ensure: Anomaly detection image X detected 1) Select region of interest X selected ⊂ X input which anomaly detection is to be tested against.
3) For each resulting fitted image, compute residuals r k [n,m] relative to X input .

4) Obtain four binary images as follows
5) Compute binary image from the following pixel-wise Boolean union:X ← 3 k=0X k .
6) Apply morphological operators as a final post-processing step: X detected ← post-processing(X).

Model Selection
To perform the model selection, we considered the three-stage iterative Box-Jenkins methodology [Box et al., 2008], which is based on the following steps: identification, estimation, and statistical model checking. For the first step, we suggest to use the Akaike's (AIC) [Akaike, 1974] and Schwartz's (SIC) [Schwarz et al., 1978] information criteria to define some rough boundaries on the choice of p and q orders. The AIC and SIC are given by where κ = (p + 1) 2 + (q + 1) 2 − 1 is the number of estimated parameters in the fitted model. Lower AIC and SIC values are related to more suitable models. The estimation step is done by the conditional maximum likelihood method, as explored in Section 2.3.
Regarding the model checking step, we suggest to consider the following approach: (i) test the overall significance of

Numerical Results
In this section, we aim at evaluating the CMLE of the 2-D RARMA model parameters and assessing the performance of the proposed image modeling and anomaly detector. For such, the proposed analyses were performed in the context of SAR image processing. We performed three numerical experiments: (i) a simulated data analysis to assess the proposed estimators and (ii) two computations aiming at anomaly detection based on actual SAR images.

Analysis with Simulated Data
Rayleigh distributed data y [n,m] were generated by the inversion method with mean given by (1)

Analysis with Real Data
In this section, we perform two anomaly detection experiments considering two different amplitude SAR images. For that, we use the methodology proposed in Section 3.2. We do not have access to the complex data of the tested SAR images, and consequently, it is not possible to check if the real and complex parts approaches the Gaussian distribution with zero mean and constant variance, as discussed in [Kuruoglu andZerubia, 2004, Yue et al., 2021]. However, evidence of the usefulness of the Rayleigh distribution for describing image pixel amplitude values is found in Kuruoglu and Zerubia [2004], Yue et al. [2021]. Moreover, the results obtained in this paper suggest that the 2-D RARMA model provides good anomaly detection and modeling for SAR data.

CARABAS II
The SAR image considered in this experiment was collected by the CARABAS II system described in the previous subsection. As reported in Ulander et al. [2005], the spatial resolution of CARABAS II is about 3 m in azimuth and range. The ground scene of the selected image is dominated by pine forest, a lake, and 25 military vehicles [Lundberg et al., 2006].
The forest and lake regions characterize most of the image area and they follow a homogeneous pattern. The military vehicles deployed in the SAR scene [Lundberg et al., 2006] introduce more representative behavior changing when compared to the forest and lake regions. For instance, the sample mean value of forest, lake and military vehicles areas are about 0.1267, 0.1148, and 0.2863, respectively, i.e., the sample mean value of military vehicles region is roughly three times the forest and lake regions. Pixels related to power line areas show similar amplitude values with the targets, and consequently, are strongly associated to the false alarm detection, as discussed in Lundberg et al. [2006]. The considered SAR image in this experiment is shown in Figure 3.
The model selection was based on the steps described in Section 3.3. We employed the search space restricted to models with (p,q) ∈ {0,1,2} and the size of the considered region of interest was N = M ∈ {10,20,40,80}. As a result, we obtained the following optimal parameters: (i) p = q = 1 and (ii) N = M = 80. The selected region of interest in this section was forest. For the post-processing step, we followed the methodology defined in Lundberg et al. [2006], where we considered two morphology operations, namely, an opening operation followed by a dilation. The opening uses a 3 × 3 pixels square structuring element, whose size is determined by the system resolution; the dilation considers a 7 × 7 pixels structuring element, which is linked to the approximate size of the military vehicles. Table 3 shows the estimated parameters and standard error (SE) for each rotated image, as described in the second step of the Algorithm 1. The overall significance Wald test p-value can be found in Table 3, showing that the spatial autocorrelation is significant for a probability of false alarm equal to 0.05. Figure 4 shows the residual versus index charts of the 2-D RARMA(1,1) models. As expected, the residuals are randomly distributed and present values close to zero.
For comparison purposes, we also obtained the detection results based on the 2-D ARMA(1,1) model. Detection results for both 2-D RARMA(1,1) and 2-D ARMA(1,1) models can be found in Figures 5(a) and 5(b), respectively. For a better visualization, Figure 5 shows the zoomed images in the region where ground type changes were detected. The proposed method detected 24 military vehicles and five false alarms. In contrast, the 2-D ARMA(1,1) model can only detect 16 military vehicles and two false alarms.
We also compared the proposed methodology with three different competing approaches: (i) constant false alarm rate    To further compare the image modeling performance of the evaluated models, we computed the MSE and MAPE of all pixels. Table 4

San Francisco Bay
The considered SAR image in this section is the San Francisco Bay, obtained by the AIRSAR sensor at band L with four looks . As reported in Safaee and Sahebi [2019], the San Francisco Bay image is multilooked by a factor of approximately 10 m in azimuth and range. Figure 6 shows To perform the ground type change detection in the San Francisco SAR image, we set the following parameters in Algorithm 1 adopting the same methodology described in the previous subsection. As discussed in Yue et al. [2019], each statistical model is suitable for one type of specific scenario of terrain surface. In particular, the Rayleigh distribution is adequate to represent homogeneous areas [Oliver and Quegan, 2004], such as pasture, crops, and sea ground type ; consequently, the region of interest in this section is linked to the ocean area. Hence, non-ocean regions (forest and urban areas) are expected to trigger a detection, suggesting an anomaly change. Considering the methodology described in Section 3.3, we have that (i) p = 1 and q = 0; and (ii) N = M = 40. In the post-processing step, we considered two mathematical morphology steps, namely, closing and opening operations. The dilation considered in both steps used a 10 × 10 pixels square structuring element, whose size is determined by the system resolution and the size of the evaluated areas.
Because q = 0, we have that θ = 0, and therefore, γ ⋆ = (φ ⊤ ,0 ⊤ ) ⊤ in (5). The estimated parameters and SE for each rotated image, as described in the second step of the Algorithm 1, are shown in Table 5. The p-values of the Wald test are also reported in Table 5, indicating that the spatial autocorrelation is significant for a probability of false alarm equal to 0.05. Regarding the diagnostic analysis of the fitted model, as displayed in Figure 7, the model residuals are randomly distributed and present values close to zero.
The anomaly detection results can be found in Figure 8. Detection results were compared to the ones based on the 2-D ARMA(1,0). The results for the detectors considering the 2-D RARMA and 2-D ARMA can be found in Figures 8(a) and 8(b), respectively. Both evaluated detectors identified the difference among the ocean ground type from the urban and forest areas in the tested SAR image. We also compared the proposed methodology with the ones presented in Gomez et al.     MSE and MAPE measures, as can be verified in Table 6. The elapsed time execution for the 2-D RARMA model was equal to 3.29 seconds, while for the Gaussian-based model, we had an elapsed time of 2.76 seconds, i.e., difference between both models is less than one second in this experiment.

Conclusions
In this paper, we proposed the 2-D RARMA model for SAR image modeling. We introduced an inference approach for the model parameters, the conditional Fisher information matrix, and the asymptotic properties of the CMLE. Monte Carlo simulations were used to evaluate the performance of the CMLE. The proposed model was applied for modeling and anomaly detection in SAR images, showing competitive results when compared to 2-D ARMA models. Moreover, although the proposed approach requires much less information when compared to Lundberg et Gomez et al. [2017], for the San Francisco image. The proposed model is presented as a suitable tool for image modeling and anomaly detection in the context of Rayleigh distributed data, in general, and SAR image processing, in particular.