{'fraction_non_alphanumeric': 0.04762918796138558, 'fraction_numerical': 0.013250047321597578, 'mean_word_length': 4.084817131857555, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Graph Signal Processing for Heterogeneous Change Detection Part II: Spectral Domain Analysis

This is the second part of the paper that provides a new strategy for the heterogeneous change detection (HCD) problem, that is, solving HCD from the perspective of graph signal processing (GSP). We construct a graph to represent the structure of each image, and treat each image as a graph signal defined on the graph. In this way, we can convert the HCD problem into a comparison of responses of signals on systems defined on the graphs. In the part I, the changes are measured by comparing the structure difference between the graphs from the vertex domain. In this part II, we analyze the GSP for HCD from the spectral domain. We first analyze the spectral properties of the different images on the same graph, and show that their spectra exhibit commonalities and dissimilarities. Specially, it is the change that leads to the dissimilarities of their spectra. Then, we propose a regression model for the HCD, which decomposes the source signal into the regressed signal and changed signal, and requires the regressed signal have the same spectral property as the target signal on the same graph. With the help of graph spectral analysis, the proposed regression model is flexible and scalable. Experiments conducted on seven real data sets show the effectiveness of the proposed method.


I. INTRODUCTION
In the heterogeneous change detection (HCD), because the multi-temporal images are acquired by different sensors and show quite different appearances and characteristics, the core of HCD is to find the connection between heterogeneous images to make them comparable. As introduced in the part I [1], there are two main challenges in HCD under the unsupervised condition: first, how to find a robust and universal connection (transformation) between heterogeneous images, which is stable even if the HCD scene is very complex; second, how to alleviate the negative influence of the unknown changed samples in the transformation.
In the part I of this paper [1], we have introduced the idea of solving the HCD from the perspective of graph signal processing (GSP). We construct a graph for each image to capture the structure information, and then treat each image as the signal on the graph. The changes between images caused by the event produce two outcomes: first, changes in the graph; second, changes in the graph signal. In this way, we can analyze the responses of different signals (i.e., X and Y) on different systems (i.e., H (S t1 ) and H (S t2 )) defined on the graphs to detect the changes: 1) we compute the changes by measuring the structure difference between the graphs from the vertex domain (in part I), which is calculated by comparing the output signals of the same signal on different graph filters, i.e., the differences between H (S t1 ) Y and H (S t2 ) Y, or H (S t1 ) X and H (S t2 ) X; 2) we can compute the changes in graph signals between X and Y by comparing the spectral properties of different signals on the same graph filter, i.e., the difference between H (S t1 ) X and H (S t1 ) Y, or H (S t2 ) X and H (S t2 ) Y in the spectral domain (in this part II), as shown in Fig. 1.
To avoid the leakage of heterogeneous data, in the measuring of structure difference of part I [1], we weaken the effect of difference in the original pixel values on image structure, and focus more on the changes in the connectivity between vertices. Therefore, we choose the graph shift operator S as the normalized average weighting matrix W avg , or normalized random walk matrix P and Laplacian matrix L rw to measure the change level, for example, d x i = (W avg t2 X) i − (W avg t1 X) i . Similarly, directly comparing H (S t1 ) X and H (S t1 ) Y will also cause the leakage of heterogeneous data. However, Setting S = W avg at this time does not solve the problem. Alternatively, we analyze the graph signal on the spectral domain, which can eliminate the influence of heterogeneous data from different domains.
Nevertheless, spectral domain analysis also presents another challenge: how to get the change of signal from the difference of spectral domain. Since the goal of HCD is to find the region where the change occurred during the event, it is directly corresponding to the region represented by the vertex in the graph. Therefore, the change measurement of d x i = (W avg t2 X) i − (W avg t1 X) i in part I that finds the changes in the vertex domain can directly output the HCD results: the changed region. As we known, in the classical signal processing, the Fourier coeffcients are the integration or summation form of the signal in the time domain. Similarly, the graph frequency coeffcients are also the summation form of the graph signal in the vertex domain. Therefore, directly comparing the spectra of two graph signals in the spectral domain is not possible to find the changed vertex (i.e., changed region).
To address the challenge introduced by the spectral analysis, we use a signal decomposition method: we decompose the source signal Y into the regressed signal Z and the changed signal ∆ as Y = Z + ∆, and require the regressed signal Fig. 1. Two strategies for the HCD from GSP perspective: calculating the structure difference by letting the same graph signal pass through different graph filters (part I); and calculating the signal difference by letting different graph signals pass through the same graph filter (part II).
Z have the same spectral property as the target signal X on the graph G t1 . With this decomposition model, we can obtain the changed signal ∆, which can be directly used to detect the areas of change. On the other hand, if we treat Z as the regression image, because the regressed signal Z is separated from Y, thenẐ andŶ belong to the same domain. And because Z and X have the same spectral property, their structure is consistent. Therefore, this signal decomposition model can also be regarded as an image regression method.

A. Related image regression works
It is a natural and intuitive idea to use an image regression approach to solve HCD problems, which transforms one image to the domain of the other image, and thereby converts the HCD into homogeneous change detection (CD) by comparing the regressed image and the target domain image. The homogeneous pixel transformation (HPT) method builds up the mappings between pixel values of heterogeneous images with a K-nearest neighbors (KNN) based kernel regression [2], which is supervised by the unchanged samples. To avoid the reliance on labeled data, an affinity matrix difference (AMD) based image regression method (AMD-IR) is proposed [3], which first employs the AMD to pick probably unchanged pixels as the pseudo-training set, and then uses different regression methods to complete the image translation. Some deep learning based regression methods have also been proposed for HCD, such as the coupling translation network (CPTN) [4], deep translation based CD network (DTCDN) [5], conditional generative adversarial network (cGAN) [6], AMD based X-Net and adversarial cyclic encoder network (ACE-Net) [7], cycle-consistent adversarial networks (CycleGANs) based method [8], and the image style transfer-based method (IST) [9].
The proposed signal decomposition based regression method has three attractive features when comparing with the above regression-based HCD methods: first, different from the common regression-based methods that aim to learn a luminance transformation to map one image to the domain of the other image, i.e., the pixel-to-pixel transformation, the proposed method is a structure regression based on the spectral analysis, which is robust to the noise and complex HCD conditions. Second, the unknown changed samples will affect the unsupervised regression process, whether manually constructed or learned regression functions. Since it is impossible to know in advance which pixels are changed, these methods can only train the regression function by selecting a pseudo-training set (e.g., AMD-IR [3], X-Net and ACE-Net [7]) or by treating all samples (pixels) as unchanged and then correcting the regression function step by step with an iterative process to refine the training samples (e.g., cGAN [6] and IST [9]), which is a coarse-to-fine process. However, the signal decomposition model, Y = Z + ∆, directly outputs both the regression signal Z and changed signal ∆, which means that the changed samples are removed from the original Y, so the negative influence of the changed samples on the regression process is reduced. Third, the proposed signal decomposition based regression method is very intuitive and interpretable by using the graph spectral analysis, and it is also very flexible by designing graph filters with different spectral characteristics.

B. Contributions
The main contributions of this part II are as follows.
• We analyze the spectral properties of different images on the defined graphs, and illustrate the connection between changes in heterogeneous images and differences in their spectral properties. • We propose a signal decomposition based regression method for HCD based on the spectral analysis, which decomposes the source signal into the regressed signal and changed signal, and constrains the spectral properties of the regressed signal and target signal to be consistent. • We analyze the vertex domain filtering (VDF) based HCD method proposed in part I in the spectral domain. We also give some discussions about the proposed regression method, and show that the proposed method can also be extended to other applications. • Experimental results on seven data sets demonstrate the effectiveness of the proposed method by comparing with some state-of-the-art (SOTA) methods (source code will be made available at https://github.com/yulisun/HCD-GSPpartII).

C. Outline
The remainder of this part is structured as follows: Section II describes the related basics of GSP. Section III analyzes an N × 1 column vector of ones the spectral property of the heterogeneous images on different graphs. Section IV proposes the signal decomposition based regression method based on the spectral domain analysis. Section V gives some discussions of the proposed method. Section VI presents the experimental results. Finally, Section VII concludes this paper and mentions the future work. For convenience, Table I lists some important notations used in the rest of this part. One notation should be noted: X i represents the i-th row of matrix X in this paper, not the i-th column as commonly used.

A. spectral filtering on graphs
Let G = {V, E, W} be a graph, and f = [f 1 , · · · , f N ] T be the signal on the graph G. Different from the Fourier domain analysis for classic signal processing, the spectral representation of graph signals employs the eigenspectra (or simply "spectra" hereafter) of the graph shift operator S given by [10] where U is an orthonormal matrix of the eigenvectors u k in its columns and Λ is a diagonal matrix of the corresponding eigenvalues λ k . The graph shift operator can be chosen as the adjacent matrix A, the random walk matrix P, or the Laplacian matrix L, based on different applications. For the undirected graph, we have S = S T , and U −1 = U T . The graph Fourier transform (GFT) e f of a graph signal f is then defined as The inverse graph Fourier transform (IGFT) is defined as Consider a shift-invariant system defined as the polynomials in the graph shift operator S of the form [10]- [12] H where S 0 = I, and h 0 , h 1 , · · · , h M are system coefficients.
The output f out of the system (4) with the input signal f in is is the transfer function of the system. Based on (5), we have e f out = U −1 f out = H (Λ) e f in , which is the spectral domain filtering of the graph signal.
Definition 2. (Spectral ordering of the weighting matrix) [12], [15]: if we sort the spectra of weighting matrix W of the graph as γ 1 ≥ γ 2 ≥ · · · ≥ γ N , then γ 1 represents the lowest frequency and γ N represents the highest frequency.
The definitions of frequency are induced by the energy of signal change (i.e., total variation), that is, we call the frequency components with smaller variations as low frequencies and call the frequency components with higher variations as high frequencies [15]. Specially, the Definition 1 is based on the 2-Dirichlet form, i.e., a quadratic function as Definition 2 is based on the 1-Dirichlet form total variation (TV) on a graph with signal f as where the normalized weighting matrix W norm = W/γ max with γ max = max |γ k | , k = 1, · · · , N .

IMAGES ON DIFFERENT GRAPHS
Given two heterogeneous images ofX andŶ acquired by different sensors over the same region at different times (before and after an event), we first construct a KNN graph to represent the structure of each image as in part I [1]: G t1 = {V t1 , E t1 , W t1 } and G t2 = {V t2 , E t2 , W t2 } for pre-and postevent images, respectively. We set X = {X 1 , X 2 , · · · , X N } and Y = {Y 1 , Y 2 , · · · , Y N } to be the graph signals with X i and Y i being the feature vectors on the i-th vertex.
In this part, we mainly consider the case of forward regression that transforms the pre-event image to the domain of the post-event image with Y = Z + ∆, and do not consider the backward regression of transforming the post-event image to the domain of pre-event image, which is a similar process to the former. Therefore, for the sake of simplicity in defining notation, we use the eigenvalue decompositions of L t1 and W t1 in the latter part of this paper as L t1 = UΛU −1 and W t1 = VΓV −1 in descending order of eigenvalues,

A. Low-pass property
As the frequency can be indicated by the rate of change between the vertices on the edges as illustrated by the TV L and TV W , we can intuitively conclude that X is the lowpass signal on G t1 , that is, the total variations of TV Lt1 (X) and TV Wt1 (X) are very small. This can be observed in the construction process of the KNN graph, where the i-th vertex and j-th vertex are connected if and only if i ∈ N x j (i.e., X i belongs to the KNN of X j or X j belongs to the KNN the X i ), which means that X i and X j are very similar.
Remark 1. The graph signal X is an approximate low-pass signal on the KNN graph G t1 , and the graph signal Y is an approximate low-pass signal on the KNN graph G t2 .
Specifically, Remark 1 can also be demonstrated by using the Definitions 1 and 2, that is, the high frequency component of e X is almost zero. Substitute L t1 and X into (6), we have where corresponding to the large λ k should also be very small. According to the Definition 1 that large λ k represents the high frequency, we have that the high frequency component in e X is very small. We define the normalized symmetrized weighting matrix as where the diagonal matrix Φ is defined such that W1 = 1 and W T 1 = 1. W sym is a symmetric doubly stochastic matrix, and it can be obtained by applying the Sinkhorn-Knopp balancing algorithm [16]- [18] to iteratively normalizes the rows and columns of W until convergence. For the weighting matrix W sym , we have 1 = γ 1 ≥ γ 2 ≥ · · · γ N ≥ 0, which follows from the property of symmetric doubly stochastic matrix. We have whereX = V −1 X. Due to the conservation of energy, we also have ∥X∥ Spectral properties of graph signals X and Y on the KNN graph G t1 , where X, Y and G t1 are constructed from the pre-and post-event images in Fig. 3.
corresponding to the small γ k is also very small. According to the Definitions 1 and 2 that small γ k represents the high frequency, we have that the high frequency component is almost zero inX. Figure 2(a) and 2(c) show the illustration of this Remark 1.

B. High-pass property
Next, we consider the spectral property of Y on the graph G t1 . As shown in Fig. 2, we find that different from the e X (oȓ X) on G t1 , the high frequency component of e Y (orY) on G t1 is not equal to zero, which seems like the noisy component in the classic signal processing.
We decompose the Y into the regressed signal Z and changed signal ∆ as Y = Z + ∆, where Z is the assumed unchanged signal that represents the translated X in the domain Y.
Remark 2. The regressed signal Z is the low-pass signal on the KNN graph G t1 . The high frequency component of e Y (oȓ Y) on G t1 is introduced by the changes caused by the event.
First, we have that the structure of Z is as same as X, that is, if X i and X j represent the same kind (or different kinds) of object, then Z i and Z j also represent the same kind (or different kinds) of object. Then, the similarity relationships of Z i and Z j is the same as X i and X j . Therefore, based on the definition of frequency or by using the TV Lt1 (Z) or TV Wt1 (Z), we have that the regressed Z have the same spectral property as the X on the graph G t1 : both Z and X are the approximate low-pass signals on G t1 .
Second, Taking the GFT of the Y on G t1 with Laplacian matrix L t1 , we have Because the high frequency component of e Z is approximately equal to zero, we have e Y k ≈ e ∆ k for high frequency λ k , that is, the high frequency component of e Y in graph G t2 is introduced by the changes ∆ caused by the event, as illustrated by the Fig. 2(b) and 2(d). Meanwhile, it should be noted that Fig. 3. Regression and changed images generated by the low-pass and high-pass filters defined on the KNN graph G t1 , respectively. By comparing (a2), (a3) with (a1), we can find that the originalX is very similar to the reconstructedX ′ andX ′′ , with the low frequency components of e X L andX L , respectively. For the regressedŶ ′ (b2) andŶ ′′ (b3) with the e Y L andY L , it can be found that their structures are consistent with the pre-event image, however the appearances are similar to that of the post-event image. For the generated difference images of (c2) and (c3), we can see that the high frequency components of e Y H andY H can be used to detect the changes.
the low frequency component of e ∆ k is not always equal to zero. For example, there will also be cases where adjacent vertices change at the same time, as shown in Fig. 3(c).

C. Regression with the low frequency components
If we divide the e Y into the low frequency and high frequency parts as e Y = e Y L + e Y H by using an ideal lowpass filter with cut-off eigenvalue λ kc , i.e., With this low-pass filter, we can also divide e Z and e ∆ as e If we ignore the low frequency component of e ∆ L , we can obtain an approximate Z ′ and ∆ ′ as: Figure 3 shows the approximate regression image and the changed image.

IV. REGRESSION BASED MODEL
From the decomposition model Y = Z+∆, we have that it requires the regressed Z have the same spectral property (i.e., low-pass) as X in the KNN graph G t1 , which can be used as a constraint for the Z. On the other hand, for the changed signal ∆, we have a change prior based sparse constraint for change detection problem, which is based on the fact that only a small part of the area changed and most areas remain unchanged during the event in practice. Therefore, we can obtain a decomposition based regression model for HCD problem where g (Z) represents the spectral constraint for Z, and f (∆) represents the prior regularization for ∆, and α > 0 is a balancing parameter. First, for the spectral constraint g (Z), we have different choices. For example, we can constrain e Z H ≈ 0 by using the penalties: represents the high-frequency transformation matrix for L t1 ; 2) are the ideal high-pass filters with the cut-off eigenvalues λ kc and γ kc , respectively.
Although the penalty with ideal high-pass filter is intuitive, it has two drawbacks: first, it requires the eigenvalue decomposition for the graph shift operator (e.g., L t1 or W t1 ), which requires a high computational complexity for large scale graph; second, it requires the selection of cut-off eigenvalue λ kc or γ kc , i.e., the cut-off frequency of the high-pass filter.
Based on the fact that the high frequency component of e Z almost equal to zero, we have that the value of is very small: first, for the small k, although λ k ̸ = 0, however, e Z k → 0; second, for the large k, although e Z k ̸ = 0, however, λ k → 0 with λ N = 0. Therefore, we can set g (z) as Second, for the sparsity regularization of f (∆), it means that ∆ only exits on a small part of vertices and remains zero on other vertices. Therefore, it requires that the number of non-zero rows in ∆, i.e., ∥∆∥ 2,0 , is very small. To meet this requirement, f (∆) can be chosen in different forms, such as the ℓ 2,0 -norm [19], the ℓ 2,1 -norm [20], the ℓ 2,p -norm with p ∈ (0, 1) [21], and the ℓ ∞,1 -norm [22], [23].
Combining the spectral constraint g (Z), and the sparsity regularization f (∆), we have the decomposition based regression model as follows

A. Optimization
By using the alternating direction method of multipliers (ADMM), the augmented Lagrangian function of (16) can be written as where R ∈ R N ×My is a Lagrange multiplier, and µ > 0 is a penalty parameter. The minimization problem of (17) can be solved by the alternating direction method, which iteratively updates one variable at a time and fixes the others.
Z-subproblem. Given the current points (Z t , ∆ t , R t ) at the t-th iteration, the minimization of (17) with respect to Z can be formulated as It can be solved by taking the first-order derivative of objective function to zero, then Z can be updated by where I N ∈ R N ×N represents an identity matrix. ∆-subproblem. Give the fixed points � Z t+1 , ∆ t , R t , the minimization of (17) with respect to ∆ can be formulated as which can be solved by the proximal operator as Input: Signal Y, graph G t1 , parameters of α, µ, ξ 0 . Initialize: Set ∆ 0 , R 0 = 0, and calculate H (L t1 ). Repeat: 1: Update Z according to (19). 2: Update ∆ according to (23), (24), or (25). 3: Update R according to (26). Until stopping criterion is met. Output: The regressed signal Z and changed signal ∆.
with Q t+1 = Y − Z t+1 + R t µ , and the proximal operator defined as Depending on different regularization forms of f (∆), we have different closed-form solutions for updating ∆ t+1 .
If we have known in prior the size of changed regions, i.e., the row sparsity level of ∆, we can construct a forced constraint as f (∆) = 0, if ∥∆∥ 2,0 ≤ τ ∞, otherwise , and then ∆ t+1 can be updated by the hard thresholding operator as where p τ is the top τ values' indices vector of Q t+1 i 2 |i = 1, · · · , N S with descending order. Multiplier updating. Finally, with the fixed points � Z t+1 , ∆ t+1 , R t , the Lagrangian multiplier can be updated as The procedure of solving the problem (16) is summarized in Algorithm 1 of Table II. The algorithm terminates when the maximal number of iterations is reached or the relative difference between two iteration results ξ t+1 < ξ 0 , where

B. DI and CM calculation
Once the regressed signal is computed from Algorithm 1, the regression imageẐ can be obtained by extracting the pixel value in Z when G t2 is a patch-wise graph or extracting the mean features in Z when G t2 is a superpixel-wise graph. With the changed signal ∆ output by Algorithm 1, we can obtain the DI as DI m,n = P Cy c=1 (ŷ (m, n, c) −ẑ (m, n, c)) 2 1/2 for patch-wise graph and DI m,n = ∥∆ i ∥ 2 , (m, n) ∈ Ω i for superpixel-wise graph. Once the DI is obtained, the binary CM solution can be regarded as an image segmentation problem by using thresholding method or clustering method, such as the Otsu threshold method [25], K-means clustering [26], fuzzy cmeans (FCM) clustering [27], or the Markov random field (MRF) based segmentation method [28].

A. Why the VDF-HCD in part I works
In the part I of this paper, we measure the change level by using the vertex domain filtering as As Y is a low-pass signal on G t2 , we have L t2 Y ≈ 0. Therefore, d y = L t1 Y − L t2 Y mainly contains information about ∆, which can be used to measure the changes.
Similarly, if we choose H (S) = W avg , we have H (Γ) = Γ, which gives a large value for low frequency (i.e., large eigenvalue as Definition 2), then it acts as a low-pass filter for the graph signal. As Y is the low-pass signal on G t2 , we have W avg t2 Y ≈ Y. With the GFT (W avg t1 Y) = Γ Z +∆ , the high-frequency components ofY (i.e.,∆) are reduced by H (Λ). Then W avg t1 Y can be regarded as an approximate Z. Therefore, d y = W avg t1 Y − W avg t2 Y ≈ Z − Y can be used to measure the changes. In fact, when we choose L = I N − W avg , we have H (Γ) = I N − H (Λ).

B. How to choose the filter H (S)
For the choices of the H (S) for vertex and spectral domain methods, we need to construct a low-pass filter for H (A) in part I, and a high-pass filter for H (L) in part I and part II based on the above analysis, Suppose Ψ (Λ) is the desired graph transfer function of a filter, we need to use the H (L) = P M m=1 h m L m to approximate the filter. Two methods can be used for the design of this spectral domain filter [10], [14]: first, using the least-squares approximation of h T being the vector of system coefficients to be estimated, Θ λ being the Vandermonde matrix form of the eigenvalues λ k , and ψ = [Ψ (λ 1 ), · · · , Ψ (λ N )] T being the diagonal vector of Ψ (Λ). Second, using the polynomial approximation of M degree, e.g., the Chebyshev polynomial series, such as the filters used in Fig. 4(c) of part I [1].

C. KNN graph construction
From the Fig. 3, we can find that the requirements for the KNN graph is that: first, G t1 can represent the structure of the image; second, X is the low-pass signal on the graph G t1 . These correspond to the two challenges of the KNN graph: the choice of K and the weighting metric.
First, a very small K is not appropriate, e.g., K = 1, each vertex is connected with is nearest neighbor. In this case, although the X is definitely the low-pass signal on the graph G t1 , but the G t1 is not robust. In this case, each vertex can only get information from its nearest neighbor, which will cause many unnecessary disconnected subgraphs in G t1 . Besides, this type of graph can not adequately characterize the structure ofX, which also means that perhaps other images that differ fromX will also be low-pass signal in G t1 . On the other hand, a very large K is also not appropriate, e.g., K = N , the complete graph, i.e., each vertex accepts information from all other vertices whether they are similar or not. In this case, we can easily find that X is not a low-pass signal on the graph, or it has a large pass band at least.
Here, we recommend the adaptive probability graph with the model where dist x i,j = ∥X i − X j ∥ 2 2 , and the regularization of R (W t1 , β) aims to make the graph smooth and avoids the trivial solution. For example, if we ignore this R (W t1 , β), then each vertex only connects with its nearest vertex with probability 1.
We present two different regularizations for W: 1) the weighted ℓ 2 -norm as used in [28], which uses the parameters β i to determine the number of neighbors k i of each vertex based on kselection strategy. Therefore, we can find that G t1 is a datadependent KNN graph with adaptive neighbors selection and adaptive weight calculation.
2) The entropy regularizer of i,j similar as [29], which attempts to maximize the information entropy of W t1 and uses the parameters β i to adjust the distribution of weights. In particular, we can find that G t1 constructed by (28) is very suitable for the proposed regression model (16): the regularization term � X T L t1 X in the objective function of graph construction model (28) is consistent with the penalty term Tr � Z T H (L t1 ) Z in the regression model.

D. Extended graphs
In this paper, we construct the KNN graph for each image to capture the structure information. An important reason why the KNN graph was chosen is that it can distinguish between the changed and unchanged signals in the spectral domain: the original X and regressed Z are low-pass signals in the KNN graph G t1 , while the high-frequency components of e Y are introduced by the changed signal ∆. Similarly, if we can Fig. 4. Regression image (b) and difference image (c) generated by the bandstop and band-pass filters defined on the KFN graph G ′ t1 respectively, where X, Y and G ′ t1 are constructed from the pre-and post-event images in Fig. 3. The high-frequency component of e X in the graph G ′ t1 , i.e., e X 4901:5000 is large. The approximate regressed image can be obtained with partial low and high frequency components, i.e., e Y Bs , and the difference image can be calculated with partial mid-frequency components, i.e., e Y Bp .
find such a graph G ′ that the changed and unchanged signals can be discriminated in the spectral domain with G ′ , then more regularization terms g (Z) can be added to the regression model (14).

If we construct a K-farthest neighbors (KFN) graph
if and only if dist x i,j is among the K-largest elements of the vector dist x i or the vector dist x j . In this way, the change rate of signal X between the vertices on the edges of G ′ t1 is very rapid, that is, the total variation of TV L ′ t1 (X) is very large. Therefore, the high-frequency component of e X in the graph G ′ t1 is also large, so as to the regression signal Z, as shown in Fig. 4. To exploit this property, we can add a KFN graph induced repulsive regularization term g (14), which requires neighboring nodes of Z connected by the KFN graph G ′ t1 to be further apart (because they do not represent the same type of object). Figure 4 shows an example of the usage of KFN in HCD. Similarly, other types of graph properties can be extended in the proposed GSP based HCD framework.

E. Framework of GSP-SDA based CD
Here, we give a general framework for CD problem (with both homogeneous and heterogeneous multi-temporal remote sensing images) based on the GSP with spectral domain analysis (SDA).
Step 1. Construct the graph and graph signals.
Step 2. Chose the penalties of g (Z) and f (∆) to obtain the regression model of (14).
Step 3. Solve the minimization problem to obtain the changed signal.
Step 4. Segment the DI to obtain the final CM. In the step 1, it requires that the graph signals are distinguishable in the spectral domain of the graph. In the step 2, the g (Z) is corresponding to the spectral property of the signal, such as the low-pass property of the KNN graph; the f (∆) can also be other prior knowledge based penalties in addition to the sparsity penalty, such as the low-rank, smoothness, and some known labels of the semi-supervised (or supervised) CD problem. In the step 4, there are many existing methods that can be used for this binary classification problem.

VI. EXPERIMENTS
In this section, we analyze the performance of the proposed SDA-HCD. We first provide a a brief description of the experimental setting, then present the experimental results on different HCD data sets, and finally some detailed discussions are made.

A. Experimental setting
We use the same heterogeneous data sets (Datasets #1-#7 as listed in Table III of part I) and evaluation metrics as in part I [1], e.g., the receiver operating characteristic (ROC) and precision-recall (PR) curves, the areas under the ROC curve (AUR) and PR curve (AUP), overall accuracy (OA), Kappa coefficient (Kc), and F1-measure (Fm). For the SDA-HCD, we construct the superpixel-wise graph with N = 10 4 , use the mean, median, and variance values of superpixel to construct the graph signals of X and Y, and choose the adaptive model (28) with weighted ℓ 2 -norm to construct the graph G t1 . We choose ℓ 2,1 -norm of f (∆) = ∥∆∥ 2,1 and set the balancing parameter α = 0.05 for the regression model (14). We choose the H (L t1 ) = L t1 + h 2 L 2 t1 + h 3 L 3 t1 with h 2 = h 3 = 1 for the results in section VI-B, and leave the discussion of H (L t1 ) in section VI-C.

B. Experimental results
We compare the proposed method with following regression based HCD methods on all the heterogeneous data sets.
HPT [2]. A supervised kernel regression method that estimates the mapping pixels using known unchanged pixels estimated mapping pixels using known unchanged pixels. In the experiments, we use 40% of the unchanged pixels as training samples. AMD-IR [3]. An unsupervised regression method uses the AMD to pick probably unchanged pixels as the pseudotraining set.   [28] 0 PSGM [30]. An unsupervised method that learns a selfexpression based patch similarity graph matrix (PSGM) for one image, then computes the regression image by multiplying the other image with this learned matrix.
SCASC [28]. An unsupervised method that uses the regression model with sparse constrained adaptive structure consistency (SCASC), which is similar to the proposed SDA-HCD but without using the high-order information and spectral analysis.
1) Image regression: In the first experiment, we verify the effectiveness of the SDA based image regression model. Figure 5 shows the regression images generated by different methods on all the heterogeneous data sets, which transforms the pre-event image to the domain of post-event image. Intuitively, these methods basically complete the image translation. That is, in the unchanged area, the regression image and the post-event image are similar, while they appear very different in the changed area, which can be used to detect the changes. Visually, we can see that SDA-HCD performs better on Datasets #1, #3 and #7; PSGM performs better on Dataset #5 and worse on Dataset #7; HPT performs better on Datasets #4 and #6. Also, we can find that the superpixel-based regression method have a block smoothing effect, such as the regression images of SCASC and SDA-HCD on Datasets #4 and #6.
2) Difference image: In order to evaluate the ability of the algorithm to measure change, we show the DIs generated by different methods in Fig. 6, and plot the corresponding ROC and PR curves in Fig. 7. The AUR of ROC curves and the AUP of PR curves are listed in Table III. In the DIs of Datasets #1 and #5, PSGM performs best, SDA-HCD comes second, as shown in Fig. 6(b)-(f) and Fig. 7. Meanwhile, we can also find that the DIs generated by PSGM, SCASC and SDA-HCD in Fig. 6 are sparse (e.g., especially in the DIs of Datasets #1, #2 and #7), which is due to the prior sparsity based regularization used in these methods. In the DIs of Datasets #4 and #6, HPT performs better than other methods. However, these results are attributed to the fact that HPT employs a large amount of labeled unchanged samples (40%). SDA-HCD can obtain better DIs on the Datasets #3 and #7 than other comparison methods, as demonstrated by the ROC and PR curves in Fig. 6. Comparing SCASC and SDA-HCD, we can find that they produce similar DIs, because they both belong to the structured graph based regression methods; however, there are differences in the performance of their DIs in Fig. 7 and Table III, i.e., SDA-HCD outperforms SCASC, which is due to the use of high-order information and spectral analysis in SDA-HCD. In general, as can be seen from the results of Figs. 6 and 7 and Table III, the DIs of SDA-HCD can well distinguish between changed and unchanged areas, and SDA-HCD gains the highest average AUP of 0.588, which is 4.4% higher than the second ranked SCASC.
3) Change maps: Fig. 6(h)-(k) show the CMs of different methods on all the evaluated data sets. Intuitively, the results generated by these comparison methods can generally reflect main information of changes. To be specific, the CMs provided by HPT, AMD-IR are affected by the salt-and-pepper noise as show in Fig. 6(h) and 6(i), especially on the Dataset #3, which is because their pixel-to-pixel regression process is not robust enough to cope with complex and variable HCD conditions. As a result, a large number of unchanged pixels are misclassified into changed ones, e.g., the CMs of HPT on Datasets #3 and #5, the CMs of AMD-IR on Datasets #1, #3 and #7. The PSGM achieves relatively good performances on some data sets, such as Datasets #1 and #5, but also has higher FP on some data sets, such as Dataset #7. SCASC can achieve satisfactory CMs with less FP on Datasets #2, #4 and #7, thanks to the exploitation of graph based robust structure consistency between heterogeneous images. Nevertheless, due to the fact that SCASC ignores the higher-order neighborhood information hidden in the graphs, making its performance not as good as SDA-HCD. For example, interpreted in detail, the results of SDA-HCD are more consistent with the ground truth, especially in some regions marked by the yellow boxes on the Datasets #1 and #7 in Fig. 6(j) and 6(k). Table IV reports the quantitative evaluation results of comparison methods on the data sets. Clearly, the quantitative results are consistent with the visual analysis of Fig. 6. As can be seen from Table IV, the SDA-HCD outperforms others on Datasets #4, #5 and #7, and obtain very competitive performance on other data sets. The average OA, Kc and Fm obtained by SDA-HCD on all the evaluated data sets are about 0.955, 0.660, and 0.683, respectively. These scores are higher than other methods, e.g., the average Kc and Fm are improved by 1.7% and 1.6% respectively compared to the second ranked SCASC. This demonstrates SDA based method can efficiently   improve the HCD performance. Finally, in order to further compare the performance of the proposed graph signal processing based HCD methods, i.e., the VDF-HCD in part I [1] and SDA-HCD in part II, the results obtained by some representative and SOTA methods [31]- [47] are summarized in Table V, except for M3CD [48], FPMS [49], CICM [50], NPSG [51] and IRG-McS [52], HPT [2], AMD-IR [3], PSGM [30], and SCASC [28], which have been compared in detail in this paper. Among these comparison approaches, DFR-MT [31], CACFL [32], AFL-DSR, DCCAE [37], DCCA [38], DPFL [43], LT-FT [44], SSL [45], ACE-Net [7], X-Net [7], SCCN [46], and ASDNN [47] are deep learning based methods. For the sake of fairness, we directly quote the results of the corresponding data sets in their original published papers in Table V. Because the datasets used in each paper are not identical, Table V is not aligned. As can be seen in Table V, the proposed graph signal processing based methods (VDF-HCD and SDA-HCD) consistently yields better or very competitive accuracy across different data sets by comparing with these SOTA approaches, which again demonstrates the effectiveness of the proposed GSP perspective for HCD problem.

C. Discussions
1) The choices of H (L t1 ): As analyzed in Section IV, the proposed decomposition based regression model (14) requires the spectral constraint g (Z) to penalize high frequency components of e Z H , such that the regressed Z have the same spectral property (low-pass) as X in the KNN graph G t1 , i.e., e Z H ≈ 0.
If we decompose H (L t1 ) = H λ. This also intuitively explains why SDA-HCD, which uses higher-order information and spectral analysis, performs better than the SCASC, as illustrated by the results of Figs. 5 to 7 and Tables III and IV.
In Table VI, we list the quantitative results (average scores on all the evaluated data sets) of SDA-HCD with different H (L t1 ) , including L t1 , L t1 +L 2 t1 , L t1 +L 2 t1 +L 3 t1 , and L t1 + L 2 t1 +L 3 t1 +L 4 t1 . With respect to the first-order H (L t1 ) = L t1 , high-order H (L t1 ) improves by 4.9% and 2.6% on average AUP and Fm in Table VI, respectively. By contrast, the fourthorder H (L t1 ) deteriorates the performance a little bit on the average scores, which is perhaps caused by the excessively narrow bandwidth of the transfer function.
2) Computational analysis: The main computational complexity of the proposed method is concentrating on the graph construction and image regression. For the former, calculating the distance matrix between all the patches or superpixels requires O � M x N 2 /2 , and sorting the distance matrix by column to construct the adaptive graph (28) with weighted ℓ 2 -norm R (W t1 , β) requires O � N 2 log N [28]. For the image regression model of Algorithm 1, updating Z with (19) requires O � N 3 , updating ∆ requires O (M y N ) by using the closed-form proximal operator (21), and updating R with (26) requires O (M y N ). Although the complexity of Algorithm 1 is very high, which requires O � N 3 for each iteration, two acceleration strategies are available to improve the efficiency of the Algorithm 1 as introduced in [28]: computing the matrix inversion of (2H (L t1 ) + µI N ) −1 off-line in advance or solving the Z-subproblem (18) with preconditioned conjugate gradient (PCG) method.
Table VII reports the computational time of each process of SDA-HCD with different N on Datasets #1 and #3, which is performed in MATLAB 2016a running on a Windows desktop with Intel Core i7-8700K CPU. In Table VII, t pre , t gc , and t ir represent the computational times spent in the pre-processing, graph construction, and image regression respectively, t total represents the total running time. As can be seen in Table VII, the running time of SDA-HCD is mainly determined by the graph scale N , and image regression is the most timeconsuming process in SDA-HCD.

VII. CONCLUSION
In this part, we analyze the GSP for HCD from the spectral domain. We first show the spectral properties of the heterogeneous images on the same graph, and illustrate that it is the changes between images that cause the differences in their spectral properties. Based on our finds, we propose a signal decomposition based regression model for HCD, which decomposes the source signal into the regressed signal and changed signal, and constrains the spectral properties of the regressed signal and target signal to be consistent. Experimental results on seven heterogeneous data sets demonstrate the effectiveness of the proposed SDA-HCD.
In this paper, we provide a new perspective for HCD problem and propose two methods with GSP, i.e., VDF-HCD from the vertex domain and SDA-HCD from the spectral domain. However, we only consider the simple graph, linear filtering and graph Fourier transform in this paper. Our future work includes studying hypergraph, graph wavelet analysis and graph neural networks on HCD problem.