Remote Sensing for Non‐Technical Survey

Abstract

This chapter describes the research activities of the Royal Military Academy on remote sensing applied to mine action. Remote sensing can be used to detect specific features that could lead to the suspicion of the presence, or absence, of mines. Work on the automatic detection of trenches and craters is presented here. Land cover can be extracted and is quite useful to help mine action. We present here a classification method based on Gabor filters. The relief of a region helps analysts to understand where mines could have been laid. Methods to be a digital terrain model from a digital surface model are explained. The special case of multi‐spectral classification is also addressed in this chapter. Discussion about data fusion is also given. Hyper‐spectral data are also addressed with a change detection method. Synthetic aperture radar data and its fusion with optical data have been studied. Radar interferometry and polarimetry are also addressed.

1. Problem statement

Any sufficiently advanced technology is indistinguishable from magic.

—Arthur C. Clarke

Experience shows that a lot of effort is sometimes spent to clear areas that are not actually mines [1]. It is therefore paramount to have a good assessment of the areas that are contaminated so that clearance is done on actually mined areas. Remote sensing can be used for that. It provides information on large areas that surveyors cannot enter because of the suspicion of mine contamination.

In the early days of research, the main objective was to find ways to detect mines by airborne survey by developing the right sensors. Nowadays, the goal is to use remote sensing to collect the right data about an area, to merge it with information from other sources, for instance expert knowledge and then improve awareness on the suspicious areas to update our assessment of their contamination.

This chapter describes some techniques developed to use remote sensing to help this activity.

A key objective is the detection of elements in the remote sensing data that may indicate the possible presence of mines. Trenches are such indicators. Section 2 describes the detection of long trenches on visible images.

Another example of the detection of indicators is the detection of craters as in Section 3.

The texture of an area as seen in airborne data may be an indication of the land use and therefore may help understand the contamination. A good example is the use of Gabor filters as in Section 4.

Remote sensing can also be used to provide a three‐dimensional (3D) representation of an area. This information proves to be invaluable to photo‐interpreters when they want to understand a certain area. A method to provide this surface, or terrain, model is in Section 5.

Classification of the land use may benefit from the multi‐spectral sensors as described in Section 6.

Section 7 explains how to fuse the data from several sensors to obtain a common and improved understanding of an area.

Hyper‐spectral sensors can bring a lot of added information for instance for change detection as explained in Section 8.

Radars are less often used in mine action but prove to be quite useful. A good example can be found in Section 9.

In order to help the mine action community to use remote‐sensing data, a geoinformation platform has been created and is available on the TIRAMISU website: http://www.fp7‐tiramisu.eu/.

2. Long trench detection (Vinciane Lacroix)

Using aerial photographs for individual mine recognition has been proposed in Humanitarian Demining since 1998 (e.g., see Ref. [2]) although with little success. Maathuis then introduced the concept of indirect minefield indicators in Ref. [3]. This concept has been further developed in the SMART project (e.g., see Ref. [4]) and used in AIDSS, an operational system designed to find indicators of mine presence (IMP) and indicators of mine absence (IMA) (e.g., see Ref. [5]). In the scope of the EU FP7 TIRAMISU project, new methods to detect and map IMAs and IMPs on aerial and satellite images were developed and tested [6–9]. The choice of indicators results from an analysis involving photo‐interpreters and possibly former members of the military who took part in the conflict that caused the contamination.

In the case of military conflicts, trenches are examples of IMPs. The first step for detecting an IMP/IMA consists in translating the indicator in terms of basic image features. As trenches can be detected thanks to their long straight linear shadows, an automatic dark line detection tool is proposed for their extraction.

The low occurrence of trenches in the huge amount of aerial photos collected during a flight over a suspected hazardous area (SHA) makes the detection task overwhelming for a photo‐interpreter. Figure 1 gives an idea of the amount of data collected during a typical campaign. The line in Figure 1 shows the path of a plane, which performed an aerial campaign over the SHAs of Bihać, Bosnia and Herzegovina, in 2010.

Figure 1.

Example of a partial aerial campaign over Bihać region in Bosnia and Herzegovina. Left: the partial fight path. Right: a stripe made of 17 colour photographs taken during the flight.

Five stripes of 17 colour photographs of size 4288 × 2848 were taken using Nikon D90 camera. On the right of the figure, the corresponding photographs of one of the airplane tracks are overlaid.

Various scenarios using the dark line detector were proposed to extract some suspicious photographs from the whole set. The result of such a scenario is shown in Figure 2. Note that the area is covered by forest producing tree shadows generating many other dark lines.

Figure 2.

Extraction of the most suspicious image thanks to the long trench detection tool: Image #9 has been identified as having the longest dark linear structures.

The trench detector based on dark line detection is described in Ref. [8] and summarized hereafter.

The success of line detection relies on (i) a line filter producing the contrast of the line and its direction at each pixel, (ii) an efficient non‐maximum suppression to extract the line axis and (iii) an appropriate linking of line elements. In order to use the line detection process for trench detection, a set of relevant local properties enabling to discriminate a trench from other elements should then be computed.

The line filter is based on the gradient line detector, which exploits the change in gradient direction at each side of a line [10]. Each pixel p is considered as a potential line element: in an eight‐neighbourhood, dp, the dot product of the gradient of the intensity at pixels arranged symmetrically around p is computed (see Scheme 1); in the presence of a line, dp is negative for some of the four pairs, and |dp|, its absolute value, is the highest for the pair lying in the direction perpendicular to the line (quantized direction). A more precise direction is provided by the difference of the gradient vectors G located at the pixels belonging to this pair. The direction of each of these gradient vectors enables to discriminate between dark and bright lines (see Scheme 1 (centre) and (right)), while the square root of |dp| provides quantitative information related to the contrast in the neighbourhood of p.

Scheme 1.

Local neighbourhood at p (left); dark line at p (centre); bright line at p (right).

A line filter providing a vector field is thus available over the image, based on the above computation. Non‐maxima suppression is then used to obtain the line median axis. Local line width is then computed by associating the line axis with its borders; the gradient norm at these borders enables to compute the local line contrast, while linking axis elements enables to obtain various average attributes over line segments: contrast, direction and width.

Red, green, blue and possibly panchromatic channels may be available for trench detection. As linear shadows should be detected, the channel offering the highest dark line contrast is kept for line extraction. This will most probably be the red channel. The IMP (i.e., trenches) should then be translated into image features according to image resolution and scene characteristics. In this case, the minimum length (in pixels) of the line segments and the minimum average width (in pixels) should be provided. Several scenarios were considered, all of them analysing the series of images taken during the aerial campaign displayed in Figure 1. In the first one, information about potential trench orientation is provided thanks to old annotated scanned maps, focusing the detection on a specific orientation. In another one, the detection in a specific orientation will be avoided; this could be an option in areas where tree shadows may generate a lot of dark lines. In the most general scenario, no constraint on the orientation is provided. In all cases, line features satisfying the constraints are provided as a list vectors ranked based on the average line contrast. The vectors may be superimposed on the image or imported in a geographic information system. The most suspicious image is the one providing the larger total length of valid segments.

When looking at the results, the photo‐interpreter may consider the image as suspicious and note it for further visual analysis, object‐based image analysis [9] or further processing such as ortho‐photo production.

The result of launching the dark line detector with a minimum length of 180 pixels and a minimum width of five pixels (the image resolution is about 10 cm) with no orientation constraint provides only three suspected photographs, all of which contain trenches. Part of the most suspected one is displayed together with the superimposed detected linear objects.

3. Crater detection (Vinciane Lacroix)

3.1. Introduction

In the scope of the TIRAMISU project, image‐processing tools to help non‐technical survey were developed and tested (see Section 2). In this framework, the Cambodian Mine Action Centre (CMAC) expressed the need for having the mapping of craters in the eastern part of Cambodia as they might provide indication on the presence of unexploded ordnance (UXO). The presence of UXOs resulting from the US bombing during the late 1960s and 1970s is still preventing the use of the land in Cambodia. When dropped, the bombs produced craters that may still exist today. Many of them are filled with water so that they appear as circular objects on satellite images.

Similar work was made by Hatfield‐Consultants [11] for Laos. The authors used historic Corona satellite images; they computed differences between the original image and its smoothed version and used these differences in an unsupervised K‐means fuzzy classifier.

We rather used data acquired by the WorldView‐2 (WV2) instrument. We used circularity of the craters with the assumption of water or bare soil inside.

This initial circle detection method applied to the panchromatic image is described in Ref. [12] and has been validated as a circle detection technique on various sets of images used for circle detection [13] and on artificial images made of controlled shapes (circles, ellipse, squares and triangles). The method called CGC (constrained gradient for circle) has been enhanced and further applied to the pan‐sharpened Band 8 (covering near‐infrared (NIR) range from 860 to 900 nm) and to various normalized difference indices. In order to validate the method as a crater detection tool, since ground truth was not available, a visual crater detection has been performed by an independent photo‐interpreter at IGEAT (Institut de Gestion de l’Environnement et d’Aménagement du Territoire). A summary of the whole process is provided in this section.

3.2. Circle detection overview

The state of the art in circle detection is provided in Ref. [12]. The most existing approaches use the fact that some basic elements (pixels, edge elements or ‘edgels’, or connected segments) are part of a circumference so that they can be combined to generate a centre‐radius pair hypothesis. To our knowledge, none of these methods directly uses a centre‐radius hypothesis. We designed a ‘Roundness Index’ filter providing at each pixel a radius and a contrast corresponding to the most contrasted circle having the current pixel as centre. The current design of the filter is such that a unique circle is assigned to a pixel (i.e., if the pixel is the centre of several circles, only the most contrasted one will be kept). The filtering involves a computation on pre‐defined digital circles of radius spanning the radius range to detect (see Figure 3, left). The Roundness Index computation described in the next section is based on the average flux and circulation of the gradient of the intensity along the tested circle and a gradient angle compatibility test. Indeed, in the case of a perfect circle, all gradient directions of pixels located on (and nearby) the circumference should point along the normal―that is, the line joining the current point to the centre of the circle; in the case of an imperfect circle, a small difference α should be tolerated. If the test is positive for a minimum proportion of pixels on the considered circumference, the pixel is considered as a potential circle centre candidate, and the radius providing the highest contrast is kept. The local maxima of the contrast are then selected. A final test on the variation of the gradient on the identified digital circle is performed in order to keep circular shapes only. A reasonable range for the tolerated gradient difference α is [0.06, 0.14] radian. And a reasonable range for the fraction of pixels satisfying the gradient angle compatibility is [0.7, 1]. The method should not be used for too small radius (i.e., above or equal to 6 pixels) because the digitalization generates errors on the gradient direction.

Figure 3.

Left: Various digital circles of odd radius (from 7 to 11) generated by Bresenham's circle algorithm. Right: local computation at $p$ located on the digital circle of radius 7, involving projection of $G$, the gradient of the intensity, along two perpendicular directions (normal and tangent to the circle).

3.3. Round index computation

The computation of the Gaussian gradient $G$ and angle $A$ is a prerequisite to the filtering process. Let $r_{\text{0}}$ and $r_n$ be, respectively, the minimum and maximum radius of the circles to detect. Let $C={\displaystyle \left\{C_{\text{0}},\dots ,C_k,\dots ,C_n\right\}}$ be the list of digital circles located at the origin, corresponding to the increasing radii $r_{\text{0}},\dots ,r_k,\dots ,r_n$. In order to initialize the process, the following information about each $C_k$ is stored: $Z$ the number of pixels in the circle and, for each of these pixels, the integer Cartesian coordinates $i$ and $j$, and the angle in polar coordinates $\theta$. Bresenham's circle algorithm may be used to compute $i$ and $j$. Finally, let $\alpha$ be the tolerance angle, that is, the maximum angle difference between the gradient and the normal to the assumed circle, and let $b$ represent the minimum portion of pixels satisfying the gradient angle compatibility test.

Two rasters, $R$, the local radius, and $F$, the filtered image providing the ‘Roundness Index’, are initialized. The image is scanned; each pixel $c$ is considered as a potential circle centre candidate; each digitized circle in $C$ is translated at $c$; the average flux $\overline{\left|G_n\right|}$ (absolute value of the projection of $G$ along the normal) and average circulation $\overline{\left|G_t\right|}$ (absolute value of the projection of G along the tangent) over all pixels of the translated circle is computed (see an example in Figure 3, right). If $\overline{\left|G_n\right|}>\overline{\left|G_t\right|}$ and if the gradient angle compatibility test are satisfied for at least $bz$ pixels, $c$ is accepted as a circle centre candidate. The compatibility test at pixel $P$ located on circle $C_j$ is $\left|\theta -A_p\right|<\alpha$, where $A_p$ represents the gradient angle at $p$ translated by $c$ and $\theta$ is the polar angle of $P$ in $C_j$. The maximum average flux $\overline{\left|G_n\right|}$ computed over all digital circles and the radius $r$ of the circle providing that maximum values are stored in $R$ and $F$, respectively. Bright circles are discriminated from dark circles based on the sign of the gradient projection (bright circles are such that the gradient is pointing towards the centre), so that the average flux is computed for dark and bright circles separately. For practical purpose, dark circles are stored as negative values. The result of the filtering process holds in these two rasters $R$ and $F$.

3.4. Crater detection

The main assumption underlying the crater detection is that craters are circular. Moreover, we assume that they are either filled with water or made of bare soil. This is probably too restrictive but limits the false‐positive rate especially in the presence of trees whose crowns are circular. Therefore, we have identified Band 8 covering near‐infrared (NIR) range from 860 to 900 nm as ideal band to detect dark circles. Water has indeed a low reflectivity in this part of the spectrum. A threshold at 0.6 on the normalized difference vegetation index (NDVI) computed from WV2 data has been added: only circles having their average NDVI equal to or lower than threshold are considered as potential craters. NDVI is computed on WV2 as the normalized difference between Band 6 $B_{\text{6}}$ and Band 5 $$, that is, $\frac{B_{\text{6}}-B_{\text{5}}}{B_{\text{6}}+B_5}$. All potential normalized differences have also been considered as a basis for dark circle detection.

We purchased WV2 data taken on 26 November 2011, covering 100 km² of the eastern part of Cambodia near the border with Vietnam, in a rural zone (Choam Kravien), with a spatial resolution of about 0.5 and 2 m for panchromatic band and multi‐spectral band, respectively, and a pixel depth at acquisition of 11 bits. The dark circle detection has been performed on the data after ortho‐rectification and atmospheric correction.

3.5. Validation

In order to validate the crater detection algorithm, as ground truth was not available, a photo‐interpreter analysed five zones of 4000 × 4000 pixels, each one covering 4 km², based on a coloured composition (red = Band 7, green = Band 5, blue = Band 3). The photointerpreter provided for each detected crater the centre, the radius, the land cover, the circularity expressed in percentage (100% representing perfect circles) and a confidence measurement also expressed in percentage. She found 334 craters over the five zones (20 km²). In the following, ‘visual crater’ and ‘CGC circle’ denote a crater identified by the photo‐interpreter and a circle detected using the automatic method, respectively. Table 1 summarizes some statistics derived from the visual crater detection over the five zones.

Confidence	#	Circularity	#	Radius	#	Land cover	#
[75, 100]	102	[80, 100]	32	Larger than 12	12	water	187
[50, 75]	109	[60, 80]	122	[9, 12]	80	Bare soil	143
[25, 50]	86	[40, 60]	109	[6, 9]	171	Vegetation	2
[0, 25]	37	[20, 40]	58	[3, 6]	69	Unknown	2
		[0, 20]	13	Lower than 3	2
Total	334		334		334		334

Table 1.

Statistics on visual crater detection: confidence, circularity, radius and land cover.

The following remarks can be drawn from Table 1:

• As far as the land cover is concerned, the assumption of bare soil or water inside crater makes sense as only four visual craters over 334 (less than 1%) have another land cover according to the photo‐interpreter.

• The circularity assumption (i.e., assuming that the crater is almost a circle) is far from true: 152 visual craters only (45%) have a circularity coefficient larger than 60.

• The radius of the visual crater can be very small: 71 (21%) have a radius lower than 6 pixels, which are the minimum radius recommended for our circle detection technique.

• Only 211 visual craters (63%) are detected with a confidence larger than 50.

Given these remarks, it seems a priori difficult to detect the craters based on the CGC method: less than 1/3 (33%) satisfies the required assumption of a radius larger or equal to 6 and circularity larger or equal to 70%. The percentage drops to 28% if one considers only the craters with a confidence larger or equal to 75%. Nevertheless, the comparison between the visual craters and the CGC circles has been performed. In order to obtain quantitative results, precision $p$ and recall $r$ have been computed over the five zones; $p$ and $r$ are defined as $p=\frac{\text{TP}}{\text{TP}+\text{FP}}$ and $r=\frac{\text{TP}}{\text{TP}+\text{FN}}$ where $\text{TP}$, $\text{FP}$ and $\text{FN}$ are the number of true positive, false positive and false negative, respectively. A false positive is a CGC circle, which is not identified as a visual crater. A true positive is a CGC circle, which corresponds to a visual crater. A false negative is a visual crater that has been missed by the automatic method. Precision is also referred as the positive‐predictive value while recall is also called true‐positive rate, sensitivity or hit rate. These definitions are better understood in view of Table 2, where it is assumed that visual detection corresponds to ‘Ground Truth’.

		Ground Truth
		Crater	No Crater	Sum	Ratio
Detection	Circle	TP	FP	Detected	P = TP/(TP + FP)
	No Circle	FN	TN
	Sum	Existing
	Ratio	r = TP/(TP + FN)

Table 2.

Precision p and recall r definitions.

The circle detection has been performed on Band 8 (with NDVI of >0.6), and on various normalized difference indices: ${\text{ND}}_{\text{83}}$, ${\text{ND}}_{\text{84}}$, ${\text{ND}}_{\text{85}}$, ${\text{ND}}_{\text{73}}$, ${\text{ND}}_{\text{74}}$, ${\text{ND}}_{\text{75}}$ and ${\text{ND}}_{\text{65}}$, where ${\text{ND}}_{ij}$ refers to the normalized difference between $B_i$ and $B_j$, that is,${\text{ND}}_{ij}=\frac{B_i-B_j}{B_i+B_j}$. A true positive is counted if the centre of the CGC circle is falling inside a visual crater. Fusion using logical operator ‘OR’ has been considered. Several precisions and recalls (summed up on Figure 4) have been computed in different ways:

1. All visual craters have been considered with an equal weight.

2. Visual craters (TP and FN) are counted according to the confidence assigned by the photo‐interpreter: the weight assigned is 25th of the confidence so that a crater of confidence 100 has a weight of 4 compared to a crater with the lowest confidence having a weight of 1. False positives have an equal weight of 1.

3. Only visual craters having a radius larger than or equal to 6 and circularity $c\ge \text{70}$ are equally considered. Nevertheless, CGC circles corresponding to visual craters having radius or circularity outside these bounds are not considered as FP.

4. Only visual craters having a radius larger than or equal to 6 and circularity $c\ge \text{70}$ are considered and weighted according to their confidence. Similarly, CGC circles corresponding to visual craters having a radius or circularity outside these bounds are not considered as FP.

Figure 4.

Precision and recall of the CGC automatic crater detection method.

Selected results among these four sets of experiments are presented in Figure 5. The selection is based on the ‘best precision’, the ‘best recall’ and the ‘best average of precision and recall’, taking each considered image independently and fusion (OR) of all of possible results.

Figure 5.

Examples of FP superimposed on part of Band 8 (left) and part of ND₇₄ (right). Bright and dark circles represent CGC circles detected on band 8 and ND₇₄ respectively.

As expected, results are less good if all the visual craters are considered in the computation, independently of their circularity and radius. Indeed, the CGC method assumes a minimum circularity of the shape and should not be used for detecting too small circles. As long as the craters satisfy the working hypothesis, the detection of craters showing much evidence is very good (more than 80% for the best fusion result).

Weighting the craters according to the confidence given by the photo‐interpreter always improves the results, which means that the CGC circles mainly correspond to visual craters with high confidence.

In general, recall is better than precision, which means that they are more FP than FN, which could be interpreted as CGC method generating too many FP. This statement should be qualified; indeed, there are 116 FP for the experiment providing the best precision, which corresponds to a density of 5.8 FP per km² or 1.45 FP for 10⁸ pixels, which is quite reasonable compared to the 334 visual craters corresponding to a density of 16.7 craters per km²or 4.17 visual craters for 10⁸ pixels. For the experiment providing the best recall, they were 738 FP (36.9 per km² or 9.22 for 10⁸ pixels) which is still acceptable.

The distribution of FP is not uniform: in particular, one zone presents a lot of false positives, which are generated either by the shadow of small‐squared houses (on Band 8) or by the house themselves (on ${\text{ND}}_{\text{74}}$) as can be seen in Figure 5.

The best precision is always obtained for the ND of Band 7 and Band 4―which means it is the best image to consider if low false‐positive rate is searched for. CGC method applied to that image provided 161 FP and 116 TP.

The best recall is a fusion of results obtained on different images. If all visual craters are considered, the fusion involves Band 8, ${\text{ND}}_{\text{83}}$, ${\text{ND}}_{\text{84}}$, ${\text{ND}}_{\text{75}}$ and ${\text{ND}}_{\text{65}}$. On the other hand, if craters outside the domain of CGC are filtered out, Band 8 fused with ${\text{ND}}_{\text{84}}$ and ${\text{ND}}_{\text{75}}$ provides the best results. That experiment provided 738 FP and 223 TP.

The CGC circles have been presented to the photo‐interpreter. While some circles (FP) have been easily identified as house shadow or house themselves (see Figure 5), others were recognized by the photo‐interpreter as potential craters that could have been missed. Indeed, this visual inspection is quite tedious and is prone to errors. Other examples of TP, FP and FN are shown in Figure 6.

Figure 6.

Examples of TP, FP and FN over part of Band 8 (left) and ND₇₄ (right).

3.6. Conclusions

This study showed that, according to a photo‐interpreter, the naive idea of a crater being seen as perfect circle on satellite image is not true; indeed, only 45% of the craters visually detected can be considered as such. Nevertheless, the CGC method used over Band 8 and fused with various normalized difference indices enables to detect at least half of them, mainly the ones having a radius larger than 6 pixels, and the ones for which the photo‐interpreter gives a larger confidence, all this with a reasonable amount of false positives. The fusion method providing the best recall should probably be used in the current context as it should be more important not missing a crater than having more false positives.

The performances would be enhanced if the method could be adapted to small circle detection, for example, by making the detection parameter dependant on the radius.

4. Segmentation with Gabor filters (Mahamadou Idrissa)

4.1. Introduction

In satellite image interpretation, classification is the operation by which an operator would like to detect different kinds of region like forest, urban zone, waterways, and so on. As the scene is a set of points (pixels) with intensity in grey‐scale values in several bands, most methods use these grey‐scale values to determine the kind of terrain. However, a single ground cover usually occupies a number of pixels with some variability in their grey‐scale values. A more satisfactory interpretation of the scene should thus include textural aspects of regions.

In general, texture is characterized by invariance of certain local attributes that are periodically or quasi‐periodically distributed over a region. There are many approaches for analysing image texture. Haralick [14] proposed a set of features (energy, entropy, maximum probability, correlation, etc.) based on grey level co‐occurrence matrices (GLCMs). Some statistical techniques use Markov Random Field models to characterize textures [15]. The spectral approach [16, 17] to texture analysis is referred to as the multi‐channel filtering approach. Textures are characterized by their responses to filters, each one being selective in spatial frequency and orientation.

We present a spectral approach to extract texture features. The textured input image is decomposed into feature images using a bank of Gabor filters. These feature images are used to form feature vectors and each of them corresponds to one dimension of the feature space. Then, we present the Fuzzy c‐means‐clustering algorithm used for unsupervised classification of the input pixels based on their associated feature vector. This method considers clusters as fuzzy sets, while membership function measures the possibility that each feature vector belongs to a cluster. At last, we present methods for evaluating how well different textures are separated in feature space, as well as measuring classification performance. In most applications, the number of classes is unknown. Here, we propose a method for choosing the best number of classes and we apply it to a real texture representation problem.

4.2. Gabor filters

Gabor filters perform a local Fourier analysis thanks to sine and cosine functions modulated by a Gaussian window. In the complex space, these filters are defined as follows:

$G{\displaystyle \left(x,y|X,Y,k_x,k_y\right)}=e^{\frac{-({(x-X)}^2+{(y-Y)}^2)}{2{\sigma }^2}}\text{.}{e}^{j(k_{x}x+k_{y}y)}$E1

where x,y represent the spatial coordinates while represent the frequency coordinates. X and Y are the spatial localization of the Gaussian window.

As signal is discrete, two simplifications are proposed [18]. The first one makes use of short‐time Fourier transform (STFT) enabling a fixed window size independently of the filter frequency. The second introduces binomial window as approximation of the Gaussian. The basis functions of this decomposition are as follows:

$\begin{array}{l}{C}_{k,l}(n,m)=W(n,m)\times \text{cos2}\pi (\frac{kn}{N+\text{1}}+\frac{lm}{M+\text{1}})\\ S_{k,l}(n,m)=W(n,m)\times \text{sin2}\pi (\frac{kn}{N+\text{1}}+\frac{lm}{M+\text{1}})\end{array}$E2

where $W^{\text{2}}(n,m)=\frac{\text{1}}{{\text{2}}^{(N+M)}}{C}_N^{(\frac{N}{\text{2}}+n)}{C}_M^{(\frac{M}{\text{2}}+m)}$ is the (N + 1) × (M + 1) binomial window. The coefficients are given by

$C_j^i={\displaystyle \left\{\begin{array}{l}\frac{j!}{i!(j-i)!}\text{ 0}\le i\le j\hfill \\ \text{0 elsewhere}\hfill \end{array}\right\}}$E3

The indexes $n=-\frac{N}{\text{2}},\mathrm{\dots },\frac{N}{\text{2}}$ and $m=-\frac{M}{\text{2}},\mathrm{\dots },\frac{M}{\text{2}}$ are the window spatial coordinates with N and M even integers. The filters selectivity in frequency (expressed in number of cycles per pixel) and orientation (in radian) is derived from the following equation:

$\begin{array}{l}f=\text{2π}\sqrt{{\left(\frac{k}{N+\text{1}}\right)}^{\text{2}}+{\left(\frac{l}{M+\text{1}}\right)}^{\text{2}}}\\ \theta =\text{arctan}\frac{k}{N+\text{1}}\frac{M+\text{1}}{l}\end{array}$E4

with $k=-\frac{K}{\text{2}},\mathrm{\dots },\frac{K}{\text{2}}$, $l=-\frac{L}{\text{2}},\mathrm{\dots },\frac{L}{\text{2}}$ and K ≤ N, L ≤ M. Figure 7 shows a subset of filter functions (their corresponding sine and cosine terms) for a 31 × 31 binomial window with K = 4 and L = 4, while Table 3 shows their corresponding sine and cosine terms.

Figure 7.

Filter bank for N = M = 30 (grey = 0, white = positive values, black = negative values).

C0,0	C0,1	S0,1	C0,2	S0,2
C1,0	C1,1	S1,1	C1,2	S1,2
S1,0	C−1,1	S−1,1	C−1,2	S−1,2
C2,0	C2,1	S2,1	C2,2	S2,2
S2,0	C−2,1	S−2,1	C−2,2	S−2,2

Table 3.

Filters legend: respectively, for the sine and cosine terms.

A set of feature images is obtained by convolving the image I with each filter providing a local ‘Energy’ given by

$E_{k,l}(n,m)={{\displaystyle \left[(C_{k,l}\hspace{1em}I)(n,m)\right]}}^{\text{2}}+{{\displaystyle \left[(S_{k,l}\hspace{1em}I)(n,m)\right]}}^{\text{2}}$E5

Here, ⊗ denotes the convolution and k,l the frequency‐orientation coordinates.

4.3. Classification

Many algorithms have been developed for supervised and unsupervised classification. In supervised classification, training sets are needed whereas unsupervised classification classifies images automatically by finding clusters in the feature space. One of the unsupervised data‐clustering methods is the hard k‐means‐clustering algorithm. It assigns each sample (feature vector) to one and only one of the clusters. This method assumes that boundaries between clusters are well defined. The model does not reflect real data where samples are in general mixed. In fuzzy c‐means algorithm [19], unlike the k‐means algorithm, samples belong to all clusters with various membership degrees. Membership degree is defined as a function of the distance between sample and the ith‐cluster centre. The fuzzy c‐means may be seen as a generalization of the k‐means method with membership values in the interval [0,1] rather than in the set {0,1}.

4.4. Classification validity

Classification of data should be of high quality, that is, all samples should have a large membership degree for at least one cluster. This problem is related to how many classes there are in the data. In fuzzy c‐means algorithm, the number of clusters is required though in many applications this information is unknown. A method for measuring performance is needed to compare the goodness of different classification results. Gath has defined an ‘optimal partition’ of the data into subgroups as follows [20]:

1. Clear separation between the resulting clusters.

2. Minimum volume of the clusters.

3. Maximum number of data points concentrated around the cluster centre.

In literature, this is known as cluster validity problem; there exists a wide variety of cluster validity parameters. Gath has proposed the partition density PD more related to the geometry of the dataset. The method accounts for variability in cluster shapes and the number of data points in each of the subsets:

$PD=\frac{S}{FHV}$E6

where $s={\displaystyle \sum {}_{i=1}^C{\displaystyle \sum {}_{j=1}^n}{u}_{ij}}$ is the fuzzy cardinality of dataset and $F_{HV}={\displaystyle \sum {}_{i=1}^C[\det {\left({\displaystyle \sum {}_i}\right)}^{1/2}}$ the fuzzy hyper‐volume of dataset.

This parameter is only sensitive to well‐compact substructures in the dataset and does not take the overlapping between them into account.

For this work, we have defined a cluster validity parameter, which combines Gath partition density PD (Eq. (6)) and a resemblance function r defined by Bloch in Ref. [21]. This function r is used as separation measurement between clusters

$r=\frac{f(A\cap B)}{f(A\cup B)}$E7

where A and B are two fuzzy sets; f(X) is the fuzzy cardinality of a set X. The intersection A ∩ B is the largest fuzzy set, which is contained in both A and B; the union is the smallest fuzzy set containing both A and B. The new definition of Compactness and Separation coefficient is

$nCS=\frac{PD}{\text{max}(r)}$E8

The factor max(r) corresponds to the maximum resemblance value measured between two clusters. This normalized factor always takes a value close to zero as soon as the clusters are disjoint; the higher resemblance value is close to one. It can also be interpreted as a punishing factor applied to the partition density PD to avoid the over‐segmentation problem.

With this definition, the procedure to select the best number of classes is to cluster the dataset for many different numbers of classes; the best one should give the maximum nCS index.

4.5. Experiments

We now apply this texture classification method to a grey‐scale image (Figure 8) taken in November 1998 during an airborne minefields survey in Mozambique [22]. In the surveyed area, minefields are homogeneous regions with no human activity inside and surrounded by agricultural fields. Since the vegetation in these minefields is different from the vegetation outside, texture‐based classification is a good approach for the interpretation of this image. Unfortunately, exact class maps or texture models are not available. Thus, the use of the proposed cluster validity parameter will help us to find the optimal class number, which maximizes the final classification accuracy.

Figure 8.

Aerial image of a minefield area in Mozambique.

Gabor filters are computed with an 11 × 11 binomial window and five discrete values of frequency coordinates k, l {−2,−1,0,1,2} (empirically determined). We have used only 9 features selected from the 13 feature images. The feature selection scheme is based on a simple method, which consists in sorting the feature images based on their amount of energy and pick up as many features as needed to achieve 95% of the total energy.

The optimal number of clusters that we find for this image is 5 (Figure 9), and the classification result is presented in Figure 10. As we can see in Figure 11, the dark area corresponding to the suspected minefield is determined and the other types of vegetation as well.

Figure 9.

Compactness‐separation function (nCS) for the aerial image.

Figure 10.

Classification in five classes.

Figure 11.

Minefield region classification.

4.6. Conclusions

A fuzzy‐clustering approach to textured image classification has been presented. The texture features are extracted using a set of Gabor filters with different frequencies and orientations.

The fuzzy c‐means algorithm has been successfully used for discriminating different types of textured image, but the drawback is that one has to specify the number of clusters. We thus discussed the use of cluster validity parameters. A modification of Compactness and Separation validity function is proposed as index to estimate the optimal number of texture categories, which have proven adequate for texture discrimination.

5. From DSM to DTM (Charles Beumier)

The bare ground is of interest for applications such as flood‐risk evaluation, environment planning or terrestrial navigation. In the specific case of minefield delineation, a terrain model is a precious piece of information to spot where camps could have been settled and hence where mines could have possibly been laid or buried.

In a mine context, the earth surface model is ideally captured by a remote‐sensing method, resulting in a fast and cost‐effective solution if the area of interest may be flown over. LIDAR and photogrammetry are the two classic options to deliver a dense digital surface model (DSM) with high resolution. Since it is taken from the air, this model contains elevated objects such as trees and buildings. These usually complicate the automatic extraction of features of the bare soil useful for detection.

By contrast, a digital terrain model (DTM) is a representation of the bare ground surface. It is advantageously derived from the DSM by automatic or semi‐automatic procedures. Most of them filter elevated points and fill filtered areas by interpolation between remaining points.

In the specific case of the TRAMISU project, gigantic DSMs (one or several giga‐pixels) were produced with a spatial resolution of 0.3 m. Fortunately, a rougher level of detail for the desired DTM is sufficient and we opted for 3‐m resolution, leading to tractable images in terms of speed and memory.

We implemented a DTM from DSM approach based on the sate‐of‐the‐art solution ‘lasground’ distributed by rapidlasso Gmbh in the LASTOOLS suite, and inspired from the work of Ref. [23] about LIDAR point filtering. The method filters out recursively locally elevated points from a point cloud and returns the reduced cloud normally made of ground points.

We first converted the DSM raster to a point cloud after low‐pass filtering and regularly sub‐sampling the raster to get a reasonable file size. The lasground programme was then run to deliver a reduced point cloud, normally made of ground points. The main parameter (‘step’) was set to 10 to accommodate for rather mountainous/forest scenes. Finally, a raster was created from the resulting point cloud and filled in by interpolation thanks to the las2dem of LASTOOLS.

The project TIRAMISU successfully employed the produced DTM. By subtracting it to the DSM, obtaining the so‐called normalized DSM (nDSM), the vegetation and buildings were easily localized. Figure 12 shows perspective views for one example of extracted DTM and the corresponding nDSM.

Figure 12.

Digital surface model and digital terrain model.

6. Multi‐spectral classification (Nada Milisavljevic)

Daedalus, the multi‐spectral scanner used in the project SMART [24], records the electromagnetic spectrum in 12 channels ranging from visible blue to thermal infrared, with the spatial resolution of 1 m. The appearance of objects in the different Daedalus channels may vary widely, depending on the molecular reflection and absorption characteristics of the matter of the objects.

Generally speaking, multi‐spectral remote sensing takes benefit from the fact that the way that different types of material or of land cover absorb and reflect solar radiation depends on the wavelength. This means that any type of land cover or material has its own multi‐spectral signature. This signature can be represented as a plot of the fraction of reflected radiation in function of the incident wavelength. For instance, the reflectance of bare soil increases with the increasing wavelength. On the contrary, clear water generally has a low reflectance, but it reaches the maximum at the blue region of the spectrum and decreases with the increasing wavelength. Finally, the reflectance of vegetation is minimal in blue and red parts of the spectrum; the reflectance increases as the wavelength increases towards the near‐infrared part of the spectrum, and becomes significantly higher than in the visible part. Apart from these general tendencies of soil, water and vegetation, differences exist in the spectral signature (i.e., shape of the reflectance spectrum) between different types of the same class. In other words, the exact shape of the signature of bare soil is a function of the type of soil, its moisture, composition, and so on. The main idea behind our work with the Deadalus data is based on the fact illustrated by these examples. Namely, in order to perform land‐cover classification, we analyse the shape of the multi‐spectral signatures and then assign pixels to classes based on their similar signatures.

There are two main types of multi‐spectral classification. In supervised classification (applied when the knowledge about the data and types of land‐cover classes is good), classes are related to known features on the ground. In unsupervised classification, the pixels are analysed in an automatic way and divided in spectrally distinct classes (i.e., accumulations of points in a multidimensional feature space). The unsupervised classification is used when we have a low knowledge about the data before classification, and the classification result is given to the image analyst, who should then give meaning to the obtained classes.

In our case, there is no reliable information on the minefield indicators (as a matter of fact, one goal of our work is to determine useful minefield indicators). In addition, we have to be careful, taking into account the context, and avoid any misclassification that could have dangerous consequences. Consequently, our plan is to proceed as follows:

• as a first step, to perform an unsupervised classification which should provide us with classes of different land cover;

• then, based on the available ground‐truth information, we attach meaning to typical and wide‐spread classes of land cover (such as forest, water, grass and cornfield);

• finally, the remaining regions that are not labelled could be related to land‐cover anomalies, which means that they are possibly dangerous and should be analysed further.

Note that in order not to discard, misclassify or miss any mined area, the first two steps have to be performed slowly and carefully. Another reason for being cautious lies in the fact that the mine density is often low in mine‐affected regions of Croatia. Furthermore, combining these results with data coming from synthetic aperture radar (SAR) (Section 7), we can improve the reliability of the obtained results.

During the application of this unsupervised classification, we take advantage of another piece of information related to the multi‐spectral signature analysis—the information provided by comparing the reflectance in the red part (r) and in the near‐infrared of the spectrum (n), the so‐called vegetation index. NDVI, that is, the normalized difference vegetation index is calculated as $\frac{n-r}{n+r}$ and

• is below 0 for snow and water;

• is close to 0 for bare soil and rocks;

• is between 0 and 1/3 in case of lower vegetation;

• increases above 1/3 with the vitality and density of the plant canopy (e.g., for fresh vegetation and forests).

Water and forests are not potential minefield indicators. Thus, before we apply a land‐cover classification, we can mask areas with NDVI below 0 and above 1/3, taking into account that there is no fresh biomass in several year‐old pioneer vegetations of the mine‐suspected fields.

Thus, the way our unsupervised classification proceeds can be summarized in the following way. As a first step, the 12 Daedalus channels are analysed pixel per pixel. A first analysed pixel is assigned to a first class, and its multi‐spectral signature (i.e., its values p_1,1, p_1,2,…, p_1,12in the 12 channels) is stored as the representative of that class, noted as c_1,1, c_1,2,…, c_1,12. The following pixel, which has values p_2,1, p_2,2,…, p_2,12in the 12 Daedalus channels, is compared with the representative of the first class based on their multi‐spectral distance t, as $t=\sqrt{{\displaystyle {\sum }_{i=\text{1}}^{\text{12}}{{\displaystyle \left(p_{\text{2,}i}-c_{\text{2,}i}\right)}}^{\text{2}}}}$.

The analysed pixel is grouped in the first class if t is lower than or equal to a tolerance set in advance (k), and the representative multi‐spectral signature of that class is then updated (replaced by averaging the signature of the prior existing pixel in the class and the newly found pixel). If t is greater than k, the analysed pixel is the founder of a new class, which means that its values initialize the multi‐spectral representative of that class. For all pixels having NDVI > 1 or NDVI < 0, a separate rejection class is formed. This procedure is repeated so that each pixel signature is compared with all existing classes. Each pixel is always grouped to the class having the minimum distance t or it creates a new class if t is greater than k. Once all pixels are analysed, we obtain a first output of our unsupervised classification, as well as the final multi‐spectral signature representative for each of the classes. At this step, there are two possibilities: our method can either end or it can go on for a new iteration, where the results of the last iteration are used as a priori classes. As a matter of fact, in the first iteration, the process of finding the multi‐spectral signatures of the classes consists in updating the class representatives each time a new class member is found, which means that misclassifications for pixels belonging to similar vegetation types are possible. After the first iteration, the class representatives are fixed, so they can be used to reclassify the pixels and correct most misclassifications (if not all). Furthermore, in such a way, a strong influence of the choice of the tolerance, k, that influences the number of classes is significantly diminished.

The described method is illustrated for the mine‐suspected area shown in Figure 13. Figure 14 contains results of the method for k = 19 (20 classes) and k = 22 (10 resulting classes), while results obtained for a larger area are given in Figure 15. If the number of classes increases, both the level of details and noise (i.e., isolated pixels that are grouped with other classes than their surroundings) increase as well, where the former could result in identifying possible anomalies in land cover. Although the resulting classified data are relatively noisy, we do not smooth the image in order to preserve possible minefield indicators.

Figure 13.

Daedalus geocoded images of a mine-suspected area in Croatia. Left: visible, right: infrared.

Figure 14.

Results of an unsupervised classification of Daedalus data for the test area in Croatia: 10 classes (left), 20 classes (right).

Figure 15.

Left: RGB composite of geocoded Daedalus channels of a larger area in Croatia; r,g,b = ch5, ch3, ch2. Right: result of an unsupervised classification of Daedalus data, 10 classes.

7. Multisensor data fusion (Nada Milisavljevic)

As shown in detail in Ref. [25], in order to find indicators of mine absence and mine presence as well as to provide image analysts with adequate tools to analyse and interpret potentially mined zones during the process of area reduction, various pieces of information obtained using Daedalus and airborne full‐polarimetric SAR, together with the existing context information (all integrated in a geographical information system), are classified and then combined.

As far as the data fusion module is concerned, we have worked on several numerical fusion methods. The proposed methods are to a large extent original thus representing a result of the project SMART. Results are obtained on three test sites using the three most promising methods: fusion using belief functions [26] and having a global discounting factor, fusion using belief functions and using confusion matrices for a more specific discounting, and fusion using fuzzy logic. Then, in the next step, we have shown how the results obtained on this basic level can be improved by introducing additional knowledge in the fusion process (e.g., knowledge on the width of the roads, on the existence of rivers, etc.). Finally, as the third step, a spatial regularization at a regional level further improves the results, based on the idea that isolated pixels of one class quite unlikely appear in another class. Note that the results are at least as good as those obtained for each class using the best classifier for that particular class. In other words, they are globally better than any separate input detector or classifier, which proves the improvement that is brought by fusion. The user can influence the choice of the classifiers as well as the choice of some parameters (some supervision is still required in the choice of parameters for the fuzzy fusion approach in particular).

As an illustration, Figure 16 contains results with fuzzy fusion. For more detailed information, please refer to Ref. [25].

Figure 16.

Results with fuzzy logic fusion: basic version (left), after knowledge inclusion and spatial regularization (right).

8. Hyper‐spectral change detection (Michal Shimoni)

8.1. Introduction

Hyper‐spectral imaging collects information from the visible to the microwave regions of the electromagnetic spectrum by taking a series of images over a range of narrow and continuous wavelengths. This results in the creation of a three‐dimensional hyper‐spectral imaging cube, where the x‐ and y‐axes present the spatial dimension of the image and the z‐axis denotes the spectrum as dependent on the number of bands. For the visible to short‐wave infrared (SWIR) wavelengths, a continuous reflectance spectrum for each pixel is formed, whereas for the mid‐wave infrared (MWIR) to long‐wave infrared (LWIR) wavelengths, an emittance spectrum is created. These wave ranges are, respectively, named in the literature as the reflective and the emittance domains.

Hyper‐spectral imagery has been applied with success to many security and defence applications including the detection of buried mines and occluded improvised explosive devices (IEDs) [27]. In the reflective domain, spectral‐based methods were applied for the detection of surface landmines using anomaly detection [28], matched filtering [29], the constrained energy minimization (CEM) [30] and the discrete wavelet transform (DWT) [31]. The presence of buried landmines was investigated using the changes in the spectrum of the covered vegetation or soil upon burying the devices [32–37].

In the radiative domain, the detection of disturbed soil exploits the ‘reststrahlen effect’, in which the emissivity absorption of a particular soil is changed within the 8–10‐μm spectral range [38]. However, the reststrahlen effect has been shown to be quite variable depending on the geographic location and the soil composition, and it requires the use of higher‐order decision fusion (i.e., fusion of thermal, textural and spectral information) to achieve target/clutter discrimination from a single thermal hyper‐spectral image [34]. Our study tries to overcome this limitation by applying spectral‐based change detection of buried IED using temporal thermal hyper‐spectral scenes. Our study assesses the detection of small buried objects using multi‐temporal thermal hyper‐spectral images and by applying two change detection methods based on multivariate statistical algorithms: the Chronochrome (CC) and the class‐conditional change detector (QCC).

8.2. Dataset

A measurement campaign took place on the 21 October 2009 in the location of a quarry in the periphery of Quebec City (Canada).

The TELOPS imager (Hyper‐Cam) [39] was located at a distance of 2.5 m away from the area of interest. The sensor scanned the area site looking downwards at an angle of 26° and with an instantaneous field of view of 1.4 μrad (Figure 17). The spectral images were collected using 660 bands in the wavelengths 800–1350 cm⁻¹. A magnifying telescope was installed for improving the spatial resolution. The detector array was set to 220 × 200 pixels, resulting in an FOV of 17.6° by 16°.

Figure 17.

View of the experimental layout showing the Hyper‐Cam field of view, buried targets, and undisturbed area.

Four aluminium plates consisting of 6.5 × 8 × 0.25 inches (Figures 17 and 18) were buried in a sandy soil and were framed with several rocks (diameter of 10–25 cm). The plates were buried at different depths: 2, 5, 10 and 15 cm (Figure 18) the day prior to the experiment on 20 October 2009, between 12h00 and 13h00. Part of the site was undisturbed as indicated at the upper part of Figure 18.

Figure 18.

Target site showing the location of the buried targets.

Images were collected from 6h40 am on the 21 October 2009 in nine sequences consisted with the following:

• 50 acquisitions of a blackbody at a given temperature;

• 100 acquisitions of the target site;

• 50 acquisitions of a blackbody at a second given temperature.

Table 4 summarizes the environmental conditions in the nine measured sequences.

Sequence	Start time	Ambient temp °C	Humidity %	Ground temp (E = 0.98)
IED#1	06:40	1.3	80	−7.3
IED#2	07:55	1.8	79	−2.5
IED#3	09:10	3.1	78	7
IED#4	10:50	6.5	57	8.5
IED#5	00:00	9.9	45	16
IED#6	13:16	9.1	45	12–14
IED#7	14:24	9.3	44	9.5–11
IED#8	15:34	7.9	44	3.5–4.5
IED#9	17:24	4.4	54	−1.3

Table 4.

Experiment log.

8.3. Change detection methods

For automatic target detection of the buried objects, two change detection methods were applied: the Chronochrome (CC) and the class‐conditional CC (QCC).

8.3.1. Chronochrome (CC)

For the two hyper‐spectral matrices $x$ and $y$, the diagonal matrices $T_x$ and $T_y$ are written as follows [40]:

$\begin{array}{l}{x}_i=T_x{\rho }_i+d_x\hfill \\ y_i=T_y{\rho }_i+d_y\hfill \end{array}$E9

where ${\rho }_i$ is the spectral reflectance and $d_x$ and $d_y$ are the changes between the observations. The spatial position index on the vector is dropped for notation convenience:

$\begin{array}{l}x=T_{xy}y+d_{xy}\hfill \\ T_{xy}=T_x{T}_{y^{-\text{1}}}\hfill \\ d_{xy}=d_x-T_{xy}{d}_y\hfill \\ \widehat{x}={\widehat{T}}_{xy}y+{\widehat{d}}_{xy}\hfill \end{array}$E10

And the change residual image ($\delta$) is defined as follows:

$\delta =x-\widehat{x}-{\displaystyle \left({\widehat{T}}_{xy}y+{\widehat{d}}_{xy}\right)}$E11

Using second‐order statistics, the transformation parameters ${\widehat{T}}_{xy}$ and ${\widehat{d}}_{xy}$ can be estimated using the mean vectors m_x and m_y and the covariance matrices C_x and C_y. Following, the covariance matrices are diagonalized:

$\begin{array}{l}{C}_x=V_x{D}_x{V}_{x^T}\hfill \\ C_y=V_y{D}_y{V}_{y^T}\hfill \end{array}$E12

where $V_x$ and $V_y$ are the eigenvector matrices, and $D_x$ and $D_y$ are the diagonalized covariance matrices. The Chronochrome (CC) change detection method is written as

$\begin{array}{l}{\widehat{T}}_{xy}^{(CC)}=C_{xy}{C}_{yy}^{-\text{1}}\hfill \\ {\widehat{d}}_{xy}^{(CC)}=m_x-{\widehat{T}}_{xy}^{(CC)}{m}_y\hfill \end{array}$E13

8.3.2. Class‐conditional CC (QCC)

The QCC [41, 42] is an image with a normal mixture model that allows the transformation parameters ${\widehat{T}}_{xy}$ and ${\widehat{d}}_{xy}$ to vary between spectral classes. In QCC, each spectrum (x) is defined by the class index q (where q = 1, 2, …, Q) with a prior probability P(q) to belong to each respective class. In this method, a class‐conditional probability function p(x|q) is assigned to the transformation parameters in Eq. (2):

$\widehat{x}|q={\widehat{T}}_{xy}|q\text{.}y+{\widehat{d}}_{xy}|q$E14

And the change residual image (δ₁) is defined as

${\delta }_{\text{1}}=x-\widehat{x}|q=x-{\displaystyle \left({\widehat{T}}_{xy}\right|q\text{.}y+{\widehat{d}}_{xy}\left|q\right)}$E15

As a classification method, we applied the stochastic expectation maximization (SEM) [43] after principal component analysis on the reference image (x). The SEM is a quadratic‐clustering method that addresses the bias against overlapping class‐conditional probability density functions and uses a Monte‐Carlo class assignment. Because in thermal imagery, the temperature variations will cause a large non‐linear radiance offset [43], we selected a stochastic model rather than a linear‐mixing model. We assume that local thermal variation occurs due to changes in the upwelling radiance and the down‐welling illumination from nearby object emissions.

8.4. Results

The nine sequenced acquisitions were resulted in 72 pairs of images to which the CC was applied. From Figure 19, one can learn that the obtained changes are spectral‐based and not thermal‐based. Apart from the illumination and the shadow effects, the obtained changes did not vary along the day. Nevertheless, because the thermal spectra are influenced by the down‐welling illumination, clear changes were recorded only from sequence #3 (i.e., 09h10 am). When the up‐spectra increase in the afternoon, the spectral changes between the disturbed soil and its surrounding are reduced.

Figure 19.

Several CC change detection results.

In Figure 19, the disturbed soil is clearly detected using the sequenced images #3, #4 and #5 (i.e., image x = #3). It is also clearly detected using the short temporal sequences (i.e., x = 3; y = 4) as the long temporal sequences (i.e., x = 1; y = 5).

Further, the detection of the buried objects was applied using the QCC method. Following the results acquired using the CC method, the transformation parameter‐based spectral classes were obtained only using the hyper‐spectral images x = #3, #4, #5. We implemented the QCC with Q = QSEM (i.e., the maximum number of classes Q may be obtained using the expectation maximization). In our case study, the Q = 5 classes: three different types of soils (included the disturbed soil), stones and shadows (Figure 20). The spectral doublet features of quartz that are located at 1176 and 1123 cm⁻¹ (8.5 and 8.9 μm) are clearly identified in the spectra of the undisturbed soils shown in Figure 21. On the other hand, the spectral features of the quartz were suppressed in the recently disturbed soil.

Figure 20.

SEM results using Q = 5.

Figure 21.

Spectra of the three detected soils using SEM.

Figure 22 presents the results of the QCC. In comparison to the CC method, the disturbed soil is clearly identified in the different sequence images when using the QCC. Using the two applied methods, the changes were clearly detected using the short temporal sequences as the long temporal sequences.

Figure 22.

The QCC results.

By merging the QCC results from different sequences (Figure 22) and after applying a threshold filter, the detection of the disturbed soil became significant.

8.5. Conclusions

This research found that spectral‐based change detection is a valuable method for the detection of buried IEDs under disturbed soil condition. The change detection methods (Chronochrome (CC)) and the class‐conditional (QCC) were able to detect changes even using short temporal sequences. The efficiency of these methods increases using pair of images from large daily temporal differences (e.g., morning and evening). A thermal hyper‐spectral sensor with very high spectral resolution (as the TELOPS Hyper‐Cam) was found to be a sensitive tool for the detection of small changes in the upper layer of the soil that were caused by the buried small objects.

9. Synthetic aperture radar (Dirk Borghys and Michal Shimoni)

9.1. Introduction

Remote sensing by synthetic aperture radar (SAR) systems provides information that is complementary to information available from optical images. Moreover, in areas under risk where rapid land‐cover mapping and suspected mined area reduction are required, the advantages of SAR offering cloud penetration and day/night acquisition are evident in comparison to optical data.

The current chapter summarizes results obtained by using different SAR modalities for helping in mined‐area reduction.

All the results presented in this chapter were obtained by automatic image processing. They could be further improved by introducing human intervention in a semi‐automatic approach.

9.2. Properties of SAR images

The aim of this paragraph is to provide a very brief overview of the principal properties of SAR images. For a more detailed and rigorous explanation, the interested reader is referred to the plethora of books, websites and online tutorials on SAR, for example, [44–48].

In SAR imaging, an electromagnetic wave is emitted towards the earth surface where it is scattered by the elements on that surface. The energy scattered into the direction of the sensor is then measured and converted into an image. Usually, the emitter and the receiver are located on the same platform (airplane, satellite or UAV). This is called mono‐static SAR, and in this case, the radar backscattering is measured.

The backscattered energy depends on properties of the scattering surface (geometrical properties, dielectric constant, etc.) as well as on the characteristics of the SAR system (wavelength, incidence angle, used polarization, imaging mode, etc.). SAR is operated in the microwave at different frequencies (Table 5). In general, each frequency has its advantage and the longer the wavelength, the higher the penetration depth of the SAR energy.

Band	Wavelength (cm)	Frequency (GHz)	Main applications
X	2.4–3.8	12.5–8	Urban
C	3.8–7.5	8–4	Land cover, soil moisture and ice
S	7.5–15	4–2	Security, manmade targets
L	15–30	2–1	Forest and soil moisture
P	30–100	1–0.3	Tomography

Table 5.

Wavelength/frequency bands used in SAR remote sensing.

The geometrical properties determine the type of scattering that occurs on the surface. For a given SAR system, the type of scattering and the dielectric constant of the scattering surface (determined by the material type and humidity) determine the amount of backscattered energy and thus the intensity of the image. For a given type of land cover, the type of backscattering depends on the frequency of the SAR system.

In a SAR system, the emitted and received radar waves are polarized. Usually, the SAR system emits either horizontally (H) or vertically polarized waves (V), and the receiver is able to capture also either H or V waves, or both. This leads to four possible polarization configurations denoted HH, HV, VH or VV. A fully polarimetric system (PolSAR system) can emit/receive these four modes simultaneously. In PolSAR, the tight relation between the physical properties of natural media and their polarimetric behaviour leads to highly descriptive results that can be interpreted by analysing underlying scattering mechanisms [49, 50].

Different SAR images acquired from a slightly different location can, under certain conditions, be combined into a set of interferometric couples. SAR interferometry provides information about the 3D structure of the terrain and can be used for generating digital elevation models (DEMs) [51]. Interferometry can also be used for detecting subtle changes by a method called coherent change detection (CCD). A complete chapter of this book is devoted to the subject.

The combination of interferometric and polarimetric information (PolInSAR) provides information on scattering type and 3D structure of the image objects [52, 53]. PolInSAR data provide information about the structure and the complexity of the observed objects. When utilized concurrently, these different capabilities allow substantial improvements in land‐cover determination [54, 55].

Table 6 gives an overview of past, current and planned civilian SAR remote‐sensing satellites used for applications over land. The table shows that most satellites operate in C‐ or X‐band, that since 2007 metric resolution satellite SAR images are available and that several SAR systems offer polarimetric capabilities. The column labelled ‘Min. SR’ gives the best spatial resolution available from that satellite. Usually, this best spatial resolution is obtained only for single‐band spotlight mode images. The best obtainable spatial resolution in full‐polarimetric mode is usually three times less good. The images from some of the satellites mentioned in the table are freely available (PALSAR, Sentinel‐1). Furthermore, in an attempt to have a larger community benefit from satellite SAR images, space agencies freely offer SAR toolboxes that include most of the relevant processing tools and are able to process most of the available SAR satellite images [56–59].

Satellite	Lifetime	Agency	Band	Polarization	Min. SR
ERS‐1	1991–2000	ESA (EU)	C	VV	25 m
JERS	1992–1998	NASDA (JP)	L	HH	25 m
ERS‐2	1995–2012	ESA (EU)	C	VV	25 m
RADARSAT‐1	1995–	CSA (CA)	C	HH	10 m
ENVISAT‐ASAR	2002–2012	ESA (EU)	C	HH/HV/VV	15 m
ALOS‐PALSAR	2006–	JAXA (JP)	L	PolSAR	10 m
TerraSAR‐X	2007–	DLR (GE)	X	PolSAR	1 m
Tandem‐X	2008–	DLR (GE)	C	PolSAR	1 m
Cosmo‐Skymed constellation	2007–	eGEOS (IT)	X	PolSAR	1 m
RADARSAT‐2	2008–	MDA (CA)	C	PolSAR	1 m
Sentinel‐1	2014–	ESA (EU)	C	Dual‐Pol	5 m
NovaSAR‐S	2015–2022	SSTL‐EADS (UK)	S	PolSAR	6 m
Radarsat constellation	2018 (planned)	CSA (CA)	C	PolSAR	3 m

Table 6.

Overview of civilian SAR remote‐sensing satellites.

9.3. Problem statement, dataset and region of interest

The SAR dataset used in this chapter was acquired in the frame of the SMART project [4]. SMART was dedicated to help mined area reduction by providing methods and integrated tools to help the use of airborne and space data.

In the project, a large dataset of airborne images was acquired over four test areas in Croatia. The flight campaign took place in August 2001 and was performed by the German Aerospace Agency (DLR). Metric resolution multi‐spectral images (Daedalus line scanner with 12 bands in VIS‐LWIR), very high resolution (5‐cm) panchromatic visual images and metric resolution multi‐frequency SAR data from the E‐SAR system of the DLR were acquired. The SAR data consist of PolSAR data in L‐ and P‐band as well as HH‐ and VV‐polarized data in X‐ and C‐band (Table 7).

E‐SAR parameter	X‐band	L‐band	P‐band
Central frequency	9.6 GHz	1.3 GHz	450 MHz
Polarization	VV	HH‐VV‐HV‐VH	HH‐VV‐HV‐VH
Nr of looks	1	1	1
Spatial resolution	1.5 m	2 m	4 m
Altitude above sea level	3160 m	3160 m	3160 m
Radiometric resolution	<2 dB	<2 dB	<2 dB
Incidence angle	50°	50°	55°
Azimuth beam width	17°	18°	30°
Elevation beam width	30°	35°	60°

Table 7.

Characteristics of the E‐SAR data (©DLR) acquired in the SMART project.

The land‐cover mapping serves for ‘suspected area’ reduction and not for mine detection. The end users defined nine land‐cover classes of interest: residential areas, roads, forests, wheat and corn fields, pastures, abandoned areas, bare soil, and rivers. These remote‐sensing data were mainly used for identifying and detecting indicators of mine presence or mine absence. A list of such indicators was compiled in the project and the different available data were investigated for their adequacy for detecting the defined indicators. Table 8 shows the indicators for which SAR data were found to improve the detection capability. The table shows that a combination of image classification and detection of linear features contributes to the detection of most indicators. For detecting power supply poles, a dedicated detector was developed by DLR, which exploits the specific effect of dihedral scattering on PolSAR images. Note that this is an example where objects much smaller than the resolution cell can be detected using SAR image processing. The relative 3D structure of the terrain can be easily extracted from interferometric image couples [51].

Class	Indicators of mine presence	Methods/tools
C1	Cultivated land (A)	Classification
C2	Asphalted roads (A)	Classification and detection of linear features
C3	Agricultural areas that are no longer in use (P)	Classification and change detection
C4	Edges of forests (P)	Classification
C5	Trenches and man‐made embankments (P)	Detection of linear features
		Interactive enhancement and extraction
C6	Soft edges of hardtop roads (P)	Classification and detection of linear features
C7	Crossroads, especially crossings of main roads with tracks no longer in use (P)	Detection of linear features, classification
C8	River banks (P)	Detection of linear features, classification
C9	Power supply poles (P)	Dedicated detector based on polarimetric features
C10	Hilltops and elevated plateaus (P)	DEM generation SAR interferometry

Table 8.

Indicators of mine presence (P) or absence (A) for which SAR data can contribute to the detection.

The remainder of the current chapter will discuss how the different SAR modalities contribute to increased classification accuracy for the other indicators. The research area is located in Glinska Poljana (Croatia), which is a post‐war land mine‐affected zone (Figure 23). The Croatian Mine Action Centre found the SAR to be the most valuable remote‐sensing source due to the high frequency of its acquisition and its independency on cloud covers. Training and validation datasets have been constructed using the collected ground‐truth data. Due to the risk incurred in mine‐suspected areas, the ground‐truth survey could be applied in only a small part of the research zone and its peripheries. Figure 23 presents the locations of the training and the validation sets superimposed on the L‐band SAR image.

Figure 23.

Training and validation dataset locations superimposed on the L‐band E‐SAR.

Based on the confusion matrix determined using the validation set, four statistics are calculated for the validation: global Kappa coefficient (κ) [60], overall accuracy (OA), producer's accuracy (PA) and user's accuracy (UA) [61].

9.4. Fusion of single‐channel SAR and optical data

For investigating the potential of ‘future’ satellite remote‐sensing data for mine‐area reduction‐related land‐cover classification, the available airborne data were used for simulating these satellite data. The Daedalus images were sub‐sampled to the spectral bands of a typical optical satellite (Pleiades). A four‐channel optical image was thus simulated (B1: 0.45–0.52 µm, B2: 0.52–0.60 µm, B3: 0.63–0.69 µm, B4: 0.76–0.90 µm). For the SAR data, VV‐polarized X‐band data were used.

For the land‐cover classification from the single‐channel SAR data, a feature‐based supervised classification method was developed and the information available from the single‐channel X‐band SAR image was augmented by combining intensity information with spatial information. Figure 24 presents a schematic overview of the developed approach.

Figure 24.

Processing scheme for land‐cover classification from single‐channel SAR data.

The introduced spatial information consisted of a combination of textural features derived from grey level co‐occurrence matrices (GLCMs) [62] and the results from a detector of bright and dark lines [63]. The line detector was applied after speckle reduction [64]. Feature‐level fusion based on stepwise logistic regression (LR) [65] was applied to obtain the classification map. The stepwise LR implicitly performs a feature selection for each class to be detected [66].

Figure 25 shows the speckle‐reduced X‐band VV‐polarized SAR amplitude image of part of the region of interest. The results of the two line detectors and of the extraction of textural parameters are shown in Figures 26 and 27, respectively. Figure 28 shows the classification results obtained from the SAR image. From the classification results, it appears that the method based on only a single channel allows detecting the location of forests and edges as well as large parts of the roads and some of the buildings. It also provides an idea of the location of abandoned land and fields in use with vegetation; however, the distinction between them is difficult because the two classes are mixed. Nevertheless, using only the single‐band SAR image the distinction between water and radar shadows is very difficult, because both surface types result in a very low backscattering. In Table 9, the global classification accuracy as well as the user and producer accuracy per class is presented.

Figure 25.

Speckle‐reduced VV‐polarized, X‐band E‐SAR image (©DLR).

Figure 26.

Results of the dark and bright line detectors.

Figure 27.

Images of the various textural features.

Figure 28.

Classification map obtained from the single‐channel X‐band SAR image.

Class name	SAR		Optical		Fusion
	PA	UA	PA	UA	PA	UA
Roads	0.47	0.41	0.86	0.28	0.86	0.30
Water	0.28	0.05	0.93	0.90	0.94	0.91
Residential areas	0.08	0.98	0.32	0.17	0.41	0.30
Abandoned land	0.66	0.78	0.74	0.82	0.81	0.89
Fields with vegetation	0.60	0.51	0.74	0.81	0.81	0.80
Trees and shrubs	0.58	0.55	0.88	0.83	0.88	0.86
Total accuracy	0.59		0.76		0.82

Table 9.

Classification results obtained from SAR processing, optical image processing and fusion.

For building a classification map from the optical data, a minimum distance classifier was applied. Decision‐level fusion was used for combining the classification from SAR and optical data. Classification results obtained from the optical data and by fusion of optical and SAR are also presented in Table 9. The table shows that while the optical image leads to better classification results than the single‐channel SAR image, the fusion of both results improves further the classification accuracy. This is particularly the case for the class ‘abandoned land’, which is a very important indicator for possible mine presence.

9.5. PolSAR and PolInSAR

In the past recent years, several polarimetric decomposition methods have been developed for various applications. Each application combines the PolSAR or PolInSAR information in order to characterize a certain type of scattering process which was determined through complementary information as experience, physical grounds [67] or by systematic selection and reduction process [68]. In very complex scenes, it is useful to exploit the discriminative power offered by the combination of a great number of these features.

At RMA, we investigated the complementarity of different frequencies and the fusion of PolSAR and PolInSAR data for land‐cover classification in mine‐covered areas. Several PolSAR and PolInSAR features are extracted in the imaging process, each combining amplitude, phase and correlation information in order to extract specific characteristics of the scene. For land‐cover classification, two levels of fusion are applied: feature and decision levels. In the feature‐level fusion, a logistic regression (LR) is used for feature selection and for fusing the selected features by optimizing the log‐likelihood function. The obtained probability images are then combined using a soft‐decision fusion method, the neural network (NN), in order to obtain the final classification results.

9.5.1. Feature extraction

A polarimetric SAR system measures the electro‐magnetic signal and its polarization state that are backscattered from the scene. The interaction of the transmitted wave with a scattering object transforms its polarization. Discrimination of different types of mechanisms is possible using the polarimetric signatures that are strongly depended on the scattering process. It has been demonstrated that the inclusion of SAR polarimetry can lead to a significant improvement in the quality of the data analysis in comparison to conventional single‐polarization SAR data [69]. In our study, for describing different properties of the land‐cover objects and for the fusion processes, a total of 72 different PolSAR and PolInSAR features were extracted. Nine PolSAR features and nine PolInSAR features are presented in Figures 29 and 30, respectively. The figures demonstrate the variation and the inherent complementarity between the different features. For a detailed description of the used feature set the authors refer to and the references therein, see Ref. [70].

Figure 29.

PolSAR features dataset extracted using different decomposition methods and by calculating the PolSAR coherences.

Figure 30.

PolInSAR features dataset extracted using different decomposition methods and by calculating the PolInSAR optimal coherences.

9.5.2. Fusion methodology

In this study, four independent sources of information: L‐band PolSAR, L‐band PolInSAR, P‐band PolSAR and P‐band PolInSAR, are considered. The image processing resulted in the extraction of a large features set: 25 L‐band PolSAR features, 25 P‐band PolSAR features, 13 L‐band PolInSAR features and 13 P‐band PolInSAR features [70].

The end users identified nine different land‐cover classes using a ground‐truth survey: roads, residential areas, wheat crop fields, corn crop fields, forests, rivers, bare soils, and abandoned areas (no trees or shrubs). The processing of the classification found two fundamental technical problems. The first consists in identifying the parameters that capture the information from the radar signals and that allow to distinguish the different classes of land cover. The second consists in selecting a technique that fuses these features in an efficient manner that will allow distinguishing the land‐cover classes.

Figure 31 presents the selected methodological scheme. The neural network, (NN) which is a non‐linear method, is selected as a classifier because linear classifiers do not provide satisfying results in the fusion of data from different sources [71, 72]. However, the NN cannot be applied directly on the original 76 features because the training set should be larger than available. To overcome this shortcoming, we separated the fusion process for two levels. In the feature‐level fusion, a logistic regression (LR) [65] is used for combining the different features extracted from the PolSAR and the PolInSAR data. The LR is used as intelligent feature reduction method. The LR is also a linear classifier that is applied separately on each information source and on each of the nine classes. The results are nine probability images for each information source. The obtained probability images are then fused using a multi‐layer perceptron neural network with one hidden layer [73, 74] to yield the final classification.

Figure 31.

The fusion process scheme.

9.5.3. Results

As mentioned before, only the probability images derived by the LR using different combinations of SAR features are used as input to the NN method.

The NN results per class using only L‐band or only the P‐band E‐SAR dataset are, respectively, presented in Tables 10 and 11. Table 10 shows that high UA and PA results using only the L‐PolSAR or only the L‐PolInSAR dataset were obtained for the classes ‘forests’, ‘wheat fields’ and ‘abandoned areas’. For the other classes, due to the low PA and UA results and due to the lack of pattern in the results, it is difficult to identify which feature set is more significant, the PolSAR or the PolInSAR. The highest PA (98.87 and 91.10%) and UA (98.71 and 88.38%) results are obtained for the abandoned area class using the L‐band PolSAR and PolInSAR probability datasets and the NN classifier.

L SAR/Land cover class	L‐PolSAR		L‐PolInSAR
	PA (%)	UA (%)	PA (%)	UA (%)
Residence	31.88	33.55	1.25	4.44
Road	6.28	7.96	24.37	24.25
River	23.92	31.82	13.49	32.89
Forest	88.00	72.33	94.58	71.67
Bare soil	41.00	50.00	67.55	85.60
Abandoned area	98.87	98.71	91.10	88.38
Wheat	87.45	55.18	94.40	77.62
Corn	30.32	16.05	43.15	36.01
Pasture	0.00	0.00	19.43	17.28

Table 10.

The NN results per class using E‐SAR L‐Band datasets.

P SAR/Land cover class	P‐PolSAR		P‐PolInSAR
	PA (%)	UA (%)	PA (%)	UA (%)
Residence	95.62	52.85	15.94	17.11
Road	36.93	37.69	14.32	12.98
River	21.94	40.00	7.01	16.18
Forest	83.13	62.66	38.00	65.73
Bare soil	82.76	49.57	39.15	17.97
Abandoned area	99.84	82.27	89.16	62.05
Wheat	51.35	76.88	91.51	40.0
Corn	0.00	0.00	0.0	0.0
Pasture	0.26	5.88	0.0	0.0

Table 11.

The NN results per class using E‐SAR P‐band.

High UA and PA results are obtained for the ‘abandoned area’ class by using only the P‐PolSAR or only the P‐PolInSAR dataset (Table 11). For the other classes, lower UA and PA results are obtained using P‐PolInSAR probability sets than those obtained using the P‐band PolSAR probability sets. Using the NN classifier, similar trend was obtained for the ‘abandoned area’ class with the highest PA (99.84 and 89.16%) and UA (82.87 and 62.05%) results and using P‐band PolSAR and PolInSAR probability datasets, respectively.

The ‘pasture’ class obtained very low PA and UA results in both L‐band and P‐band datasets. Analysis of the confusion matrices shows that the ‘pasture’ class was mainly confused with bare soil. As expected, the signal of the L‐ and the P‐frequencies penetrates the pasture areas that are covered with very short scattered grass, as in bare soil.

A comparison between the results was obtained using single dataset (Tables 10 and 11) and the results after fusion (Table 12) show that process using the NN significantly improves the results of the land‐cover classification.

Land cover class	NN classification
	The best fusion set	PA (%)	UA (%)
Residences	All‐SAR	90.31	68.81
Roads	All‐SAR	41.84	62.67
River	All‐SAR	61.64	81.73
Forest	All‐SAR	91.60	88.08
Bare soil	All‐SAR	89.32	79.92
Abandoned areas	All‐SAR	99.51	99.35
Wheat	All‐SAR	99.03	97.53
Corn	L‐PolSAR + L‐PolInSAR	65.31	66.27
Pastures	L‐PolInSAR + P‐PolSAR	31.35	29.58

Table 12.

The highest PA and UA results per class obtained using the NN classifier.

The best results are derived using the ‘All‐SAR’ probability datasets for seven of the nine land‐cover classes (Table 12). The highest PA and UA results were obtained for all the classes using the fused PolSAR and PolInSAR datasets. These results emphasize again the complementarities of the different SAR frequency datasets and of the PolSAR and PolInSAR information. The results of the ‘abandoned area’ class are the highest with PA (99.51%) and UA (99.35%). The results for this class have a high importance for the Croatian mine action (i.e., end users) because it considered as suspected area. Moreover, it was demonstrated that using a pixel‐wise classification based on optical data or a single SAR dataset, it is difficult to distinguish this class. Applying the LR and multinomial logistic regression (MNLR) methods only on the PolSAR data resulted with lower user accuracy (UA) of 80% and producer accuracy (PA) of 66% for the ‘abandoned areas’ in ‘Glinska Poljana’ [66], than the results were presented in this study. Using the PolSAR and multi‐spectral data and a Dempster‐Shafer (DS) method, the ‘abandoned areas’ are classified with PA of 84% and using a fuzzy method with PA of 89% [25].

After the fusion process, the road and river classes are still confused. However, those classes can be easily distinguished by a photo‐interpreter. Table 13 presents the overall accuracies and the Kappa values obtained using the NN for the land‐cover classifications. The results show that the classification accuracy improves for the complete land‐cover classification, using fused datasets from different SAR frequencies. A kappa of 0.809 and the highest accuracy (84.00%) are obtained using the All‐SAR dataset. The results emphasize again the complementarities of the PolSAR and PolInSAR information and of data acquired in different frequencies.

LR probability dataset	NN
	Overall accuracy (%)	Kappa
L‐PolSAR	52.13	0.437
P‐PolSAR	59.72	0.522
L+P‐PolSAR	67.48	0.618
L‐PolInSAR	61.04	0.539
P‐PolInSAR	37.69	0.522
L+P‐PolInSAR	66.16	0.601
L‐PolSAR + L‐PolInSAR	63.87	0.576
P‐PolSAR + P‐PolInSAR	65.55	0.596
L‐PolSAR + P‐PolInSAR	63.22	0.569
L‐PolInSAR + P‐PolSAR	69.96	0.646
All‐SAR	84.00	0.809

Table 13.

Overall accuracy for land‐cover classification obtained using the NN classifier.

Figure 32 presents the land‐cover classification obtained by implementing the All‐SAR dataset and using the NN classifier. Many false alarms are present in the class ‘road’. Please note that we did not select ‘shadow’ as class because we did not visualize this phenomenon in the original image scene. However, common artefacts as shadow and layover are also part of the SAR images and future improvement should include this class.

Figure 32.

NN land‐cover classification of Glinska Poljana obtained using All‐SAR dataset.

9.6. Conclusion

The results presented here show that a feature‐based fusion method for classifying SAR image data is a valuable approach. When only a single‐channel SAR image is available, a rough land‐cover map can be extracted by combining the available intensity information with spatial information such as line detection results and textural features. Fusion with a classification obtained from optical images does provide a significant improvement in classification accuracy with respect to using optical data alone.

Features from different SAR frequencies were found to be complementary and adequate for land‐cover classification. It was also found that PolInSAR features are complementary to the PolSAR features and essential for producing an accurate classification of complex scenes composed of very diverse land‐cover classes. Fully polarized and interferometric SAR data proved to produce valuable remote‐sensing information and can be used to obtain accurate information for areas under danger of mine presence.

The current availability of high‐resolution SAR image data from various satellites on one hand and the SAR‐processing toolboxes freely available on the other hand makes it easier to apply an approach as described in this chapter. Nevertheless, an efficient exploitation of the full information richness offered by SAR images and its complementarity to information from optical data still requires acquiring a good understanding of the SAR image phenomenology.

10. Conclusion: the future of remote sensing in mine action

The information collected by remote sensing can help mine action in many other ways that survey a suspicious area to decide where to send the clearance teams. In some mined countries, the availability of up‐to‐date maps can be poor. Remote sensing can be used to create recent maps to help the planning of mine action operations. An analysis of the relief and the land cover can help understand the road network and decide the best route to access a given location. Deciding which areas to clear first is often a difficult decision. Remote sensing can provide input regarding the risk of the presence of mines or the vulnerability of local population. In case of major disasters such as flooding, remote sensing can be used to assess the change in the mine contamination.

For decades, airplanes and satellites were the only vectors to obtain remote‐sensing data. Drones offer now new opportunities. With the reduction of size and costs of cameras, they change drastically the way remote sensing can be used. If the autonomy of the smaller drones is still limited, they still allow for fast deployment, which can help reactivity of mine action organizations.

With drones offering a good solution for small‐range survey and satellite for large areas, the question now is about the relevance of conventional aircraft. Have they still a place between drones and satellites? U.S. Air Force is currently training more pilots of drones than pilots of fighter aircraft. Will aircraft slowly phase out for remote sensing?

Acknowledgments

The authors of this chapter would like to thank the Belgian Ministry of Defence and the European Union.

Section 2: The aerial campaign was supported by financial aid of the US State Department, through the International Trust Fund for Demining and Mine Victims Assistance (ITF), and kindly shared by HCR Centre for testing, development and training Ltd.

Section 3: The research leading to these information and results has received funding from the European Community's Seventh Framework Programme (FP7/2007–2013) under grant agreement n° 284747 (TIRAMISU). Noveltis and IGEAT, two partners of TIRAMISU, contributed to this crater detection study. Noveltis pre‐processed the data (georeferencing and atmospheric correction). IGEAT made the visual interpretation and provided valuable information for the data preparation and results discussion.

Sections 6, 7 and 9: The authors thank all the researchers of the project SMART, co‐funded by the European Commission (IST‐2000‐25044) and involving the following partners: CROMAC, DLR, ENST, IXL, Renaissance ASBL (RMA), RST, TRASYS, ULB‐IGEAT and Zeppelin.