Trace Coherence: A New Operator for Polarimetric and Interferometric SAR images

Quadratic forms play an important role in the development of several polarimetric and interferometric synthetic aperture radar (Pol-InSAR) methodologies, which are very powerful tools for earth observation. This paper investigates integrals of Pol-InSAR operators based on quadratic forms, with special interest on the Pol-InSAR coherence. A new operator, namely Trace Coherence, is introduced, which provides an approximation for the center of mass of the coherence region (CoRe). The latter is the locus of points on the polar plot containing all the possible coherence values. Such center of mass can be calculated as the integral of Pol-InSAR coherences over the scattering mechanisms (SMs). The trace coherence provides synthetic information regarding the partial target as one single entity. Therefore, it provides a representation, which is not dependent on the selection of one specific polarization channel. It may find application in change detection (e.g., coherent change detection and differential DEM), classification (e.g., building structure parameters), and modeling (e.g., for the retrieval of forest height). In calculating the integral of the Pol-InSAR coherences, an approximate trace coherence expression is derived and shown to improve the calculation speed by several orders of magnitude. The trace coherence approximation is investigated using Monte Carlo simulations and validated ESA (DLR) L-band quad-polarimetric data acquired during the AGRISAR 2006 campaign. The result of the analysis using simulated and real data is that the average error in approximating the integral of the coherence region is 0.025 in magnitude and 3° in phase (in scenarios with sufficiently high coherence).

In calculating the integral of the Pol-InSAR coherences, an approximate Trace Coherence expression is derived and shown to improve the calculation speed by several orders of magnitude.
The Trace Coherence approximation is investigated using Monte Carlo simulations and validated ESA (DLR) L-band quad-polarimetric data acquired during the AGRISAR 2006 campaign. The result of the analysis using simulated and real data is that the average error in approximating the integral of the Coherence Region is 0.025 in magnitude and 3 degree in phase (in scenarios with sufficiently high coherence).

I. INTRODUCTION
1 Synthetic Aperture Radars (SAR) are powerful sensors able to acquire high resolution im-2 ages of scene reflectivity at microwave frequencies [1]. Such products are complementary 3 to optical images and have the advantage of measuring with almost any weather conditions 4 and at night time. Also, microwaves can penetrate some class of targets providing informa-5 tion on the internal target structure [2]. For those reasons, SAR has been largely used for 6 stationary target detection and recognition, multi-pass target change detection and retrieval 7 of biophysical parameters (specially related to vegetation). Finally, the use of multiple po-8 larimetric channels (i.e. 2 or 3 channels) or multiple flight passes (i.e. baselines) increases 9 substantially the amount of observables allowing the development of more powerful method-10 ologies [3], [4]. The combination of polarization and interferometry is often referred to as 11 Pol-InSAR [5], [6]. In the last decades a large variety of PolSAR and Pol-InSAR method-12 ologies were proposed. A very short list of applications are retrieval of parameters [7], [8], 13 [9], detection/classification [10], [11], [12], [13] and change detection [14], [15], [16], [17]. 14 In the following, a very brief introduction to Pol-InSAR is provided with the purpose of 15 presenting the mathematical formalism exploited in the rest of the paper. eficial. In this paper, a parametrization proposed by Cloude  . The latter constrains the covariance matrix to 60 be 3 × 3. Given a partial target, the power backscattered by a specific SM can be calculated 61 considering the quadratic form of the covariance matrix T with an appropriate projection 62 vector ω that represents the SM: where Θ represents the support of the projection vector which is a unitary complex sphere.

70
S is equal to the surface of such sphere and the integral is divided by S because we are not 71 interested in the size of the support. In other words, we want that a unitary function (i.e.

72
T equals to the identity matrix) provides a unitary integral. S can assume different values 73 depending on the dimension of the space in which ω lives (i.e. dual-or quad-pol data).

74
Interestingly, it is not necessary to know the exact value of S for the following derivation.

75
The final solution of the integral is:  [23].

87
The center of mass of the CoRe (i.e. the average of all the coherence points) depends on 88 the density of points inside the loci. This is defined as the integral: In this paper, a new operator named Trace Coherence is introduced: Motivated by the previous result, we hypothesize that γ tr can approximate the γ integral: It is clear, that the result obtained with a single quadratic form cannot be extended straight-92 forwardly, since the coherence operator is nonlinear. On the other hand, it is possible to 93 prove that the previous equality holds when the matrices have some specific structures.
To investigate this, the eigenvector basis of the first covariance matrix T 11 can be used to 95 represent the space. Following the previous nomenclature, the integral can be rewritten as: The U ij matrices multiplying T 11 are not written because the product is traceless.

97
Looking at Eq. 8, it is possible to tell that the equality holds in the following situations:

98
Proof 1: If the three matrices are equal but differ only by a scaling factor, the eigenvectors of 100 T 11 are able to vanish the off diagonal terms of T 22 and T 12 as well. Additionally, decorrelation.

116
This condition forces T 12 to be diagonal (i.e. the polarimetric channels are not cor-117 related independently on the interferometric information). In this situation, the off 118 diagonal terms of all the matrices T 11 , T 22 and T 12 vanish independently on the basis 119 used. Moreover, A, B and C will be always equal to 1 3 , since each element of any 120 orthonormal set will contain one third of the total matrix energy. Therefore, A, B, 121 and C simplify and they leave the integral equal to γ tr .

122
In all the other situations, it is not possible to prove mathematically that the integral is 123 equal to γ tr and therefore it has to be considered an approximation. Interestingly, Proof 124 2 and Proof 3 coincide with the boundary conditions for the polarimetric behavior of par-125 tial targets (completely polarized and de-polarized, respectively). Therefore, we may hope 126 that intermediate situations will have similar behavior. In order to test the approximation, 127 simulated and real data are exploited in the following.

129
Before proceeding with tests, it is valuable to spend few words explaining some advan-130 tages of using γ tr . From a general point of view, the main advantage of using γ T r compared 131 to a single scattering mechanism solution (i.e. single channel, optimum polarization, ex-132 tremes of CoRe) is that γ T r represents a synthetic information about the CoRe intrinsically 133 based on the idea that the observed target is partial and therefore composed by several SM.

134
With synthetic information it is meant a quantity that is able to combine, compact and syn-135 thesize a larger amount of information, which would otherwise need many more numbers.  In this section, a sensitivity analysis is performed by means of a Monte Carlo method.

196
Some preliminary results on such analysis can also be found in [31].

198
The simulations performed in this work assume the scattering vectors to be Gaussian. The

199
impact of different types of texture could be investigated in the future.

200
The simulations were performed as follows:

201
(1) A Monte Carlo method is used to generate N realizations of scattering vectors drawn by a 202 3D-Complex Gaussian distribution. This is performed twice (one for each acquisition).

203
Therefore, we generate two sets of "white random vectors": k w 1 (i) andk phase difference is selected. This allows to model the shape of the coherence loci on the polar plot. The way this is done is by generating k w where a is a 3 dimensional real vector with each element included in the interval [0, 1], and • is the Hadamard (or element-wise) product. a and b contain information regarding 211 respectively interferometric correlation and phase differences of the three simulated SM.

212
(3) The white random vectors are colored using two asymptotic covariance matricesṪ 11 and 213Ṫ 22 (representing the partial targets observed in the two acquisitions):Ṫ The simulated Pol-InSAR matrices are calculated by averaging N realizations of the 216 outer product of the simulated vectors: where . N is the finite average of N realizations. Since the simulation pro-218 vides random variables that are close to be independent, the value of ENL can be ap-219 proximated by N. In the real scenario, the realizations would be neighbor pixels and the 220 average would be done by a spatial filter (e.g. a boxcar).

221
(5) Points 1 to 4 are repeated K times to evaluate statistics. They represent K experiments.

222
In other words, a set of K covariance matrices is produced: T 11 (k), T 22 (k) and T 12 (k), 223 with k = 1, ..., K. Each of these realizations is slightly different due to speckle. iment the CoRe is evaluate using L points. A block diagram of the simulation procedure is 231 presented in Figure 1.

232
For the sake of simplicity, we concentrate on simulations with the assumption of ESM: should clarify that the MCI is a very good approximation of the integrals (we believe better 260 than γ tr as shown in next section), but it is not immune from errors. Therefore, part of the 261 errors that we estimate in next sections could be related to the MCI. It was shown that the equality between γ tr and the integral of γ holds when the entropy of 264 both covariance matrices is either 0 or 1. In this section we would like to put this under test.

265
The simulator was used to generate CoRe with triangular shapes. This is for the sake of In this section, we concentrate the tests on H = 0.5, since this seems to be the worst 296 scenario in terms of entropy. 297 We want to analyze the effect of using different SM to simulate the partial target. The appears to be rather similar, with a triangular shape and points concentrated in the upper 308 right corner. The variation between the different CoRe is only due to speckle, because they 309 use the same underlying covariance matrix to generate the Monte Carlo realizations. Figure   310 6.a presents the error ∆ as defined previously. Interestingly, this is always around 0.03.

311
The fluctuation can be explained because each experiment has slightly different covariance 312 matrices (due to speckle) that therefore generate slightly different CoRe. It can be inferred 313 that the error is independent of the specific SM. Abstracting this result, we could say that the can also be easily proofed using the definition of γ tr and noticing that the Trace is basis The magnitude of the error is presented in Figure 8.

334
It is possible to observe that the error depends on R. In particular, it is minimum around 335 R = 0.5, that is, when the moving SM is more aligned with the other two. The error increases 336 when the loci are more stretched and the point density is less uniform. Fortunately, even in the worst conditions, the error seems to be in mean around 0.04.

338
A.5 Worst case scenario 339 We want to devise an experiment that is the most challenging for the approximation.

340
Please note, the simulations performed in this section are likely to be unrealistic and they 341 only purpose is to gain understanding of the approximation. Therefore, we are not suggest-342 ing that the peculiar shapes presented in this section could be observed in real data. 343 We want to create a point density that is largely unbalanced inside the CoRe. A way to between γ tr (k) and the peak (i.e. or mode) of the CoRe P eak(k). Please note, the word 367 "mode" may be an abuse of notation because it requires interpreting the CoRe as a random 368 process. This is true when the loci are obtained by a Monte Carlo method, but in general it is more proper to talk about point density as deterministic 3-D surfaces on the polar plot.

370
In order to test this idea, the point density is estimated with a histogram and the peak 371 location is determined. Care was taken to have the histogram bin size small enough to 372 accurately capture the peak location, but not too small to produce a jagged histogram. Once 373 the peak location is determined, the differences between P eak(k) and γ tr (k) or M CI(k)   Previous sections adopted the ESM assumption. In this section, we want to gain some 413 understanding regarding the approximation when the ESM hypothesis is not fulfilled. Figure   414 13 displays the results of two tests.

415
(1) The SM of the partial targets stay the same. The entropy of the first target is H I = 0.5 416 and the entropy of the second target is varied between 0 and 1 (H II ∈ [0, 1]).

417
(2) The entropy of T 11 and T 22 is 0.5. The dominant α of the first target is α I 1 = 0 and the 418 dominant α of the second target is varied between 0 and 90 degrees (α II 1 ∈ [0, 90]).

419
As expected, the approximation is dependent on the specific selection of the two partial speckle) is larger. This is because the structure of the T 12 matrix becomes more variable.

425
As a final remark, it is important to keep in mind that the way non-ESM targets are simu-426 lated can impact strongly the analysis. In this work, we set that the projection of the second 427 partial target over the first maintains the same correlation. This is to say that the second par-428 tial target is obtained by the first one plus an additive component (that is clearly uncorrelated 429 with the first target). This therefore does not cover the case when the first target is substituted 430 by a completely new second target. The latter scenario will show a much larger decorrelation 431 and we may expect that the CoRe will cluster around the zero in a uniform way. This should 432 improve the approximation.

433
A.9 Summary of simulations 434 It was observed that the approximation depends on the CoRe shape and point density.

435
Further experiments showed that γ tr is significantly biased toward the peak of the density 436 at a level that it could be possible that γ tr represents the peak. Unfortunately, proving this 437 property is not trivial, unless the peak and mean have the same location. This happens in the 438 following situation: (1) ESM hypothesis and single decorrelation mechanism: the CoRe is a circle with a sym-440 metric density

441
(2) The entropy is unitary: even if the distribution is rather smooth, the peak and the mean 442 would both be the middle point between the three scattering mechanisms.

443
(3) The entropy is null: in this case, the CoRe collapses to a single point 444 The previous conditions are the same in which it was possible to proof that γ tr is equal to 445 the integral.

446
As a final remark, in some applications, knowing the location of the peak may be even

455
Some preliminary test of γ tr on the AGRISAR dataset can be found in [32].

456
The main parameters of the acquisitions exploited in this work are shown in Table I.
457 Figure 14 shows the RGB Pauli images of a portion of the entire scene that will be used 458 as initial test area. The two acquisitions considered here were carried out the 13 th of June 459 and the 5 th of July and they have a nominal baseline equal to zero. However, in the exploited 460 data, there is still a residual baseline. In the future some test will also be focused on baselines 461 largely different from zero. The scene presents several agricultural fields with some buildings  (farms). The color coding of the Pauli RGB is Red: 1 √ 2 |HH − V V | 2 , Green: It can be observed that some of the fields have experienced a change between the two 465 acquisitions, while others appear to be rather stationary. Also, the image contains bright 466 point targets that allow to test a variety of entropy values. while the even-bounce (or horizontal dihedral) seems to be the weakest for most of the fields.

474
Also, the cross-polarization channel (which is often associated with volume scattering) is 475 stronger than the dihedral scattering (since there is a volume component), but it is lower than 476 the surface scattering, since it suffers more from volume decorrelation. Additionally, all the 477 SM which present a low backscattering (e.g. HV channel on bare ground) will suffer from 478 noise decorrelation.

479
The images for the wrapped interferometric phases are reported in figure 16. This final section is dedicated to estimate the approximation error. Figure 18 shows the 509 CoRe for three generic points in the image covering winter wheat and field grass. The latter 510 are just a few representatives of the shapes that we can encounter in this dataset. Due to 511 the fact that we are often in the condition of polarimetric non-stationarity (i.e. the ESM 512 hypothesis is often not fulfilled), the loci can assume shapes that differ from triangles or 513 ellipses (even presenting regions that have a non-convex shape, but still connected). In order to have a more quantitative comparison for the quality of the approximation the 530 scene is sampled with a grid of 50 pixels width in range and azimuth and the resulting 1200 531 pixels are used to extract statistics. Figure 19.a presents the magnitude of the difference 532 between γ tr and MCI, while 19.b presents the same difference as complex numbers on a 533 polar plot. Interestingly, the difference seems to be rather contained with values on average 534 around 0.02. Moreover, the distribution of the errors on the polar plot is quite homogeneous 535 which suggests that there are no biases.

536
To investigate these last points some histograms are shown in Figure 20. Additionally, 537  The final test concerns the estimation of the difference between the peak and γ tr or MCI.

557
Unfortunately, the analysis could not be accomplished successfully because we did not man-558 age to produce reliable estimates of the peak location. The plots for this analysis are not 559 presented here because unreliable, however it is still possible to observe that γ tr appears to 560 be generally closer to the peak. As discussed previously the problem in estimating the peak 561 is that the density can have rather small derivatives for values of entropy higher than 0.5 and 562 a Monte Carlo search with histograms fails.

564
The main goal of this paper was to propose a new operator for Polarimetric SAR Inter-565 ferometry (Pol-InSAR), namely the Trace Coherence γ tr . This operator is an approximation 566 of the center of mass of the coherence region (CoRe), that can be formally evaluated as the 567 integral of the Pol-InSAR coherence γ over all the projection vectors.

568
The following mathematical proofs were given:

569
(1) The integral of a quadratic form is equal to the third part of its trace.

570
(2) The integral of γ is equal to the new operator Trace Coherence γ tr in the following .

575
The approximation was tested using Monte Carlo simulations and real data. γ tr represents 576 a good approximation in several situations. In particular, the error depends on the density of 577 points in the CoRe. The accuracy of the approximation degrades as the CoRe size increases 578 and the skewness of the CoRe increases. This is because γ tr appears to be closer to the 579 location of the density peak that can differ from the mean. As a consequence, γ tr can also 580 be used as an approximation of the peak location when this is very different from the mean.

581
The latter is an interesting feature since it is particularly hard to retrieve the peak location 582 when the polarimetric entropy is higher than 0.5.

583
The tests on real data (AGRISAR 2006, DLR) showed an average error of approximately 584 0.025 in magnitude and less than 3 degrees in phase; however, the average phase error can 585 increase to as much as 7 degrees for low values of coherence (around 0.2).

586
As a future work, a larger dataset with available ground measurements will be used to 587 validate specific applications in the context of coherent change detection and classification.

588
Specifically, a set of indexes will be designed and tested against different forest types to 589 understand if different forest structures can be discriminated by this synthetic information.

591
Here, the proof of the integral of a single quadratic form is equal to the third part of the 592 matrix Trace is provided.

593
The integral can be rewritten as: whereT = T T race [T ] and Ω = ω ω * T . The latter passage was obtained calling the property of cyclic permutation of the Trace (i.e. the first ω is moved after the second ω).

595
In our case, T is Hermitian and therefore it can be diagonalized. Without loss of generality, we can use the eigenvector basis u i with i = 1, 2, 3 to represent any vector in the space. The integral variable ω can be written as a linear combination of the eigenvector basis. Therefore, |a| 2 + |b| 2 + |c| 2 = 1, |a| 2 = A, |b| 2 = B, |c| 2 = C, It can be easily proven that the matrix Ω can be decomposed in the sum of three U i matrices plus the sum of matrices with zero trace (i.e. in the eigenvector basis they only have offdiagonal elements). Ω = AU 1 + BU 2 + CU 3 +  The three integrals are computed on the projections of the vector ω over the eigenvector basis that varies when ω is changed. Moreover, the three integrals have the same value, since each of the components will cover the same volume of space while ω is swept over the entire unitary complex sphere. Since they are squared values they vary inside a (non-negative) cube of unitary side. To conclude, the three integrals have to sum to one (i.e. the volume is unitary) and they have to be equal, therefore each of the integrals has to be equal to one third. The sum of the three eigenvector matrices in the eigenvector basis is clearly the identity 601 matrix I (i.e. they are the standard basis), therefore the solution of the integral is: ACKNOWLEDGEMENT 603 The AGRISAR2006 data were acquired by the E-SAR airborne system of DLR and they 604 were provide by ESA.

605
The author would like to thank Simon Zwieback (ETH Zurich) for the engaging conversa- for the support of the work that was carried out at ETH Zurich.