BACKGROUND OF THE INVENTION

[ 0001 ] The present invention relates generally to a method and apparatus for processing space-borne radar backscatter data. In particular, the method and apparatus process filtered and terrain corrected backscatter generated by satellite-borne synthetic aperture radar (SAR) apparatus to detect the deforestation or degradation of forest vegetation.

[ 0002 ] It is desirable to measure deforestation and forest degradation in order to monitor: i) individual forest carbon projects; ii) individual protected forest areas; iii) forest estate belonging to a sovereign state, including both protected and unprotected areas; and/or iv) forest reserves belonging to a private enterprise. Reduction in forest coverage has a significant impact on the pace of global warming and reduces bio-diversity.

[ 0003 ] Satellite-based optical remote sensing data is commonly used for monitoring changes in forest area: examples include optical date captured by MODIS and Landsat satellites. However, optical data suffers from: a) cloud (and smoke) obscuring images of forests at a given point in time; and b) not being able to readily detect forest degradation.

[ 0004 ] In the case of cloud obscuration, it is feasible that circumstances can result in series of images of certain areas of forest (in particular, forests in the tropics) that are never cloud-free in optical imagery. When areas are obscured, even temporarily, it is not possible to determine whether there has been deforestation at a site. If a patch of forest is obscured for two years, say, that may be sufficient time for forest to be degraded, and then for the canopy gap to close back, and therefore for that change not to be detected at all.

[0005] Forests are complex three dimensional environments with vegetation at different heights relative to the underlying terrain and the leaf canopy presented by the tallest trees. In addition to the possibility of secondary regrowth (after disturbance removing some trees), it is also possible for lower tiers of vegetation to be significantly degraded while the leaf canopy remains unchanged. Degradation such as this is not easily detected as the spectral signature in optical data from degraded forest is not easily distinguishable from the intact forest.

[0006] It is possible to estimate aboveground biomass (AGB) from LiDAR data, see for instance EP2163846. While the data from LiDAR is laser data which penetrates the forest canopy and can give information on mid-tier and lower tier vegetation; there are significant obstacles to the provision of effective AGB estimation (not least of these obstacles being that there is no current LiDAR satellite operational) . It is furthermore considered impractical to attempt to provide global continuous coverage using LiDAR data from satellite platforms, and thus the use of LiDAR data for estimating rates of deforestation and degradation is limited to users who can afford to task LiDAR data collection from aircraft.

[0007] MODIS data is freely available and has been used to determine the extent of fire damage, for example, in geographical areas otherwise difficult to survey (see CN103150569) . Landsat data has also widely been used to map deforestation, and is now also free. Recent publications have released free-to-use maps of deforestation from Landsat, though with limitations remaining due to cloud cover and likely frequent errors of omission of small-scale deforestation, and no attempt to map forest degradation.

[0008] Radar data can overcome some of these problems. Due to the wavelengths of electromagnetic energy emitted by satellite radar units, the detecting radiation does not interact with cloud and smoke particles, hence the sensor can 'see through' that cloud and smoke.

[0009] Long wavelength radar bands, like IEEE L-band (with a wavelength of between 15cm and 30cm) , IEEE P-band or Very High Frequency (VHF) , are most clearly suitable for measuring forest structure and forest degradation. Where radar units use longer wavelengths (such as the L-band radar at ~23cm wavelength as used on the Japanese ALOS-PALSAR satellite) , the energy interacts with the limbs and trunks of trees. The energy reflected to the sensor (called "backscatter" ) can then be related to the amount of biomass in the forest. It can also then detect forest degradation, as the quantity of limbs and trunks in the forest are reduced through logging and fires. The use of L-band radar data to monitor forests is discussed in a series of published articles by the authors of the present document. However, long time series of L-band radar data do not exist, and there is no current (or near/mid-term) prospect of free or low-cost L-band radar data.

[0010] Currently, the longest wavelength SAR sensor ever used in a satellite operates in the L-band: operation at P- band is however contemplated, a P-band mission being planned for 2020 (BIOMASS) .

[0011] It has been known to use even longer wavelengths, for instance VHF, to measure biomass (discussed in patent documents US2010225531 & US5886662, for instance) . VHF radar is known to be good at mapping aboveground biomass and changes in biomass. These deployments use aircraft as a platform: the extremely long wavelength of VHF means no satellite mission could be designed using conventional technology.

[0012] Existing techniques are thus not fit for many purposes: in the case of optical approaches there are limitations due to cloud cover and the ability to detect forest degradation (as opposed to deforestation) ; with L-band radar approaches there is a lack of free or low-cost data, and a lack of data availability, with no operational satellite series planned.

[0013] It would therefore be desirable to provide an improved and more nuanced monitoring of deforestation and/or forest degradation.

SUMMARY OF THE DISCLOSURE

[0014] In one aspect of the present disclosure, there is provided a method for processing radar backscatter data from forested terrain, the method comprising: obtaining a plurality of time-separated rasters of radar backscatter data from a remote sensing platform; masking said rasters to exclude confounding variables; filtering said rasters to reduce the effect of radar speckle; comparing the pixel values of the data in respective rasters, where any given pixel corresponds to substantially the same geographic location; for each geographic location, determining whether the difference in pixel values between rasters exceeds a predetermined threshold value; and generating a representation of changes in forest coverage for each geographic location in accordance with the outcome of said determination. [0015] In certain embodiments, the remote sensing platform may be a satellite.

[0016] In certain embodiments, the remote sensing platform may include synthetic aperture radar, SAR, equipment.

[0017] In certain embodiments, the radar backscatter data is generated using the SAR equipment and the SAR equipment is sensitive to radar in the C-band.

[0018] In certain embodiments, the method may further comprise: estimating the intensity of disturbance using the magnitude of the change between radar scenes or other ancillary data. Here the intensity of disturbance may be categorized into at least two classes.

[0019] Examples of such classes of disturbance may include ( 'degraded' and 'deforested' ) , further categorical classes may be used (e.g. 'high', Medium' and 'low' intensity change) .

Alternatively a continuous variable may be derived estimating the percentage of canopy cover lost or total carbon stock changes due to the disturbance.

[0020] In certain embodiments, the method may further comprise: estimating a level of confidence associated with the detected change predicted to have occurred at a particular date .

[0021] Here the estimated level of confidence may be categorized into at least two classes. Alternatively the estimated level of confidence may be evaluated on a continuous scale. In either case, the level of confidence may be determined using the radar backscatter and/or other ancillary data . [0022] The estimate of confidence may be expressed either categorically (e.g. ^> high' , Medium' and lovi' confidence), or continuously, for example a percentage confidence.

[0023] In a further aspect of the present disclosure, there is provided an apparatus for processing radar backscatter data, the apparatus comprising a processor which is operable to execute a method in accordance with any one of the above- described aspects.

[0024] Another aspect of the present disclosure provides a computer program comprising instructions arranged, when executed, to implement a method in accordance with any one of the above-described aspects. A further aspect provides machine-readable storage storing such a program.

[0025] Various further aspects and embodiments of the present disclosure are provided in the accompanying independent and dependent claims.

[0026] It is an aim of certain embodiments of the present disclosure to solve, mitigate or obviate, at least partly, at least one of the problems and/or disadvantages associated with the prior art. Certain embodiments aim to provide at least one of the advantages described below.

BRIEF DESCRIPTION OF THE DRAWINGS

[0027] The invention, together with objects and advantages thereof, may best be understood by reference to the following description of preferred embodiments together with the accompanying drawings in which: [0028] Figure 1 illustrates the detection of backscatter from a typical multi-tier forest using SAR equipment on board a satellite platform;

[0029] Figure 2 illustrates certain components of one type of SAR equipment;

[0030] Figure 3 illustrates functional components of a satellite platform including SAR equipment such as that illustrated in Figure 2;

[0031] Figures 4A and 4B illustrates the basic principle of the SAR technique and shows the change in range to a given geographical feature over time from a moving platform;

[0032] Figures 5A, 5B and 5C chart respective phases in the processing of radar backscatter data in accordance with certain embodiments of the present disclosure;

[0033] Figure 6 shows three "raw" rasters of backscatter data before and after masking;

[0034] Figure 7 shows three masked rasters of backscatter data before and after filtering;

[0035] Figure 8 shows the identification of change pixels at successive time frames, resulting in a cumulative output plot of change pixels in accordance with certain embodiments of the present disclosure; and

[0036] Figure 9 charts an alternative to the first phase in the processing of radar backscatter data in accordance with certain alternative embodiments of the present disclosure.

DETAILED DESCRIPTION OF PREFERRED EMBODIMENTS

[0037] The detailed description set forth below in connection with the appended drawings is intended as a description of presently preferred embodiments of the invention, and is not intended to represent the only forms in which the present invention may be practised. It is to be understood that the same or equivalent functions may be accomplished by different embodiments that are intended to be encompassed within the scope of the invention. In the drawings, like numerals are used to indicate like elements throughout. Furthermore, terms "comprises," "comprising," or any other variation thereof, are intended to cover a nonexclusive inclusion, such that module, circuit, device components, structures and method steps that comprises a list of elements or steps does not include only those elements but may include other elements or steps not expressly listed or inherent to such module, circuit, device components or steps. An element or step proceeded by "comprises ..." does not, without more constraints, preclude the existence of additional identical elements or steps that comprises the element or step .

[0038] Throughout the description and claims of this specification, the singular encompasses the plural unless the context otherwise requires. In particular, where the indefinite article is used, the specification is to be understood as contemplating plurality as well as singularity, unless the context requires otherwise.

[0039] It will be appreciated that features and aspects of the present disclosure described above in relation to the first and other aspects of the invention are equally applicable to, and may be combined with, embodiments of the invention according to the different aspects of the invention as appropriate, and not just in the specific combinations described above. Furthermore features of the dependent claims may be combined with features of the independent claims in combinations other than those explicitly set out in the claims . [ 0040 ] To measure changes in forest area over large areas of land, it has become accepted to use space-borne radar backscatter data. Radar has the capability to penetrate through tree canopy and cloud cover.

[ 0041 ] A common definition of forest is an area greater than 1 ha in size consisting of trees with a potential to grow at least 5 meters in height, with a canopy cover of at least 30 %.

[ 0042 ] A forest is said to have been "degraded" if human action reduces its canopy cover but it remains above 30 %, and therefore is still forest. A forest is "deforested" if it is modified by humans so it no longer meets the definition of forest. The change in canopy cover as a forest is degraded or deforested is thought to be particularly evident in radar data, since theoretical models of radar backscatter predict substantial decreases in backscatter as size and density of trees decrease.

[ 0043 ] When compared to the optical methods of detecting changes in state of forest areas (using spectral decomposition to identify leaf colour and density), radar has advantages. Radar backscatter responds to density, size, orientation and water content of scattering elements on the surface of the target region; which means that it can be more sensitive to changes in the quantity or characteristics of the woody material in the target region.

[ 0044 ] The shorter wavelength radiation employed in IEEE C- band (~6cm wavelength, frequencies ~ 6GHz) radar, has been discounted for use in monitoring forest change over time because has been presumed to interact too weakly with features of forest biomass to provide an effective metric of the vegetation . [ 0045 ] However, the present authors have discovered that it is possible to use time series (also referred to as "stacks") of data from shorter wavelength radar to detect deforestation and forest degradation.

[ 0046 ] These stacks of shorter wavelength radar (e.g. C- band radar) data are processed into estimates of deforestation and forest degradation at a given time period.

[ 0047 ] The approach works by:

• determining intact forest areas

• excluding confounding variables in the data set (for example correcting for the presence of slopes; etc.)

• calculating statistics on the properties of the C band radar data over time

• calculating statistics on the properties of the radar data over the areas of forest over time; and

• integrating these statistics dynamically to determine the change in the radar signal over time and the change in radar data forest over time in order to estimate where and when deforestation and degradation have occurred.

[ 0048 ] Crucially, these statistics account for a) natural variation in forest, and b) changes in radar backscatter readings which occur due to other seasonal factors (e.g. the presence of drier vegetation and soil in dry season) .

[ 0049 ] As a consequence of the above, there is provided:

• a method for processing shorter wavelength radar data to detect deforestation and forest degradation more accurately and for larger sets of data than has previously been possible: the unique filtering and data processing scheme developed here allows genuine changes to be distinguished more readily from random radar noise (and "speckle" arising from coherency) .

• the method also deals with change dynamically and automatically, in that it accounts for changes in the forest and radar signal which are independent of actual changes in the degree of forest degradation/deforestation itself .

• Overall, the algorithm can provide automated estimates of deforestation at multiple instances through the year, unobscured by cloud or seasonal variation in forest,

[0050] Certain embodiments use a type of radar data which has hitherto not been exploited for forest observation as it has been assumed to be inappropriate, and has not been widely available. Shorter wavelength radar data is nevertheless available: for instance, the Sentinel 1 satellite is currently beginning to produce large volumes of C-band data spanning geographical regions across the world, including the tropical forest regions where there is the particular demand for monitoring. Sentinel-1 is planned to image areas of interest to the detection of deforestation, such as the tropics, around four times each year.

[0051] Figure 1 illustrates a satellite 120 with SAR equipment on board passing a typical multi-tier forest 110.

[0052] One example of such equipment is the phased array- type L-band synthetic aperture radar (PALSAR) sensor of the Advanced Land Observing Satellite (ALOS) . This sensor has the capability of collecting cross-polarised (i.e. HV, horizontal- send vertical-receive) data in addition to like-polarised (or "co-polarised"), horizontal-send horizontal-receive (HH) data. The cross-polarised data can only be from scattering elements that change the polarisation of incoming electromagnetic radiation: a complex three dimensional object such as a tree produces a strong response in the cross-polarised data whereas moisture in soil will not. Such polarization shifts in the scattered radar signal may thus be used feature analysis, discriminating between features that appear similar in data of a single polarization. Prior art methods for quantifying changes in forest state (e.g. forest degradation) have used this cross-polarised data to distinguish vegetation from underlying terrain.

[0053] It is noted that for the purposes of the present discussion horizontal polarization indicates that the electromagnetic radiation is polarized with an electric field component varying in a direction parallel with the earth' s surface, whereas vertical polarization indicates that the electromagnetic radiation is polarized with an electric field component varying in a direction normal to the earth' s surface. These orientations are illustrated as "V" and "H" in Figure 1.

[0054] Figure 1 also shows schematically the transmission of radar pulses by the SAR equipment, the backscatter of certain pulses from a typical multi-tier forest 110 and the detection of those backscatter pulses using SAR equipment on board the satellite platform 120.

[0055] Figure 2 illustrates certain components of one type of SAR equipment. Here the SAR equipment 200 includes a pulse transmitter unit 230 and an antenna 210. The SAR equipment 200 in Figure 2 further includes a circulator 250, which allows pulse signals generated by the transmitter unit 230 to reach the antenna, the transmitter unit 230 being electrically connected to the antenna 210. The generated pulse signals supply energy to the antenna 210 and the antenna 210 sends pulse trains of corresponding electromagnetic radiation at a known frequency in a direction and polarization defined by the antenna type.

[0056] A small proportion of the radiation pulse echoes from a target in the direction of the antenna 210. This "backscattered" radiation (often simply termed "backscatter" ) is the returning portion of a previously transmitted radar pulse .

[0057] An antenna (often the same antenna that transmitted the radiation pulse, as is the case illustrated in Figure 2) may then be configured to receive the backscattered radiation and to transfer the power in the returned radiation, via the circulator 250, into a received electrical signal. This received signal is then transmitted via an amplifier 240 and received at a receiver unit 220, where it is sampled, stored and/or displayed. The receiver unit 220 and the transmitter unit 230 may be collocated, especially in cases where the antenna 210 is shared between transmit and receive paths.

[0058] Certain radar arrangements allow the polarization of the antenna, when configured to create an electromagnetic radiation pulse, to be set differently from the antenna, when configured to receive the return radiation pulse. In the preceding discussion, the ALOS PALSAR sensor transmits horizontal polarized electromagnetic radiation while receiving the returned radiation as both horizontal and vertical polarized electromagnetic radiation. In the case of the standard mode of the Sentinel 1 satellite, the radiation is transmitted in the vertical polarization and received in both vertical and horizontal (W and VH) .

[0059] In Figure 2, a further optional antenna 212 is illustrated. Antenna 212 may be of different antenna type to antenna 210 and may thus transmit electromagnetic radiation in a polarization different from antenna 210 (horizontal, H say rather than vertical, V, polarization) . An optional switching component 260 may ensure that pulse signals generated by the transmitter unit 230 may be directed to reach the further antenna 212. The generated pulse signals supply energy to the antenna 212 and the antenna 212 sends pulse trains of corresponding electromagnetic radiation at a known frequency in a direction and polarization defined by the further antenna type .

[0060] The same further antenna 212 may also be configured to receive backscattered radiation and to transfer the power in the returned radiation, via a further circulator 252, into a further received electrical signal. This further received signal is then transmitted via a further amplifier 242 and received at a further receiver unit 222, where it is sampled, stored and/or displayed.

[0061] The resulting stored data therefore corresponds to at least one polarization band (i.e. W, VH, HV and/or HH) depending upon the various characteristics of the transmit and receive antennas. Where data from more than one band is available, a further data set corresponding to total power

(TP) may also be stored.

[0062] Figure 3 illustrates additional functional components of a satellite platform 300. The SAR equipment 200 of Figure 2 is incorporated in a radar module 340. The satellite platform 300 also includes: a processor unit 330, which may include one or more processors and/or digital signal processors; a memory unit 310 configured to store at least temporarily data output by the one or more receiver units 220, 222; and a communication module 320 configured to communicate data from the satellite platform to a ground-based control centre . [0063] Backscatter data is gathered in ^x slant range' , as SAR is a sidewise looking sensor type. Antennas for SAR (such as antenna 210 and further antenna 212 in Figure 2) are directed at a predetermined angle from the nadir (i.e. the point vertically below the satellite platform) . This is illustrated in Figure 4A: the nadir at time To is denoted 404 and the antenna is directed at an off-nadir angle to give a footprint 410. Processing, including a Digital Elevation Model (DEM) is therefore required to convert this signal into a conventional-looking raster image, mapped onto a standard map projection. The same target is shown being measured at other time points, Ti, T ₂ and T ₃ .

[0064] Figure 4B illustrates the basic principle of SAR. As the satellite platform 120 proceeds along its path, the Doppler effect means that pulse trains emitted in the direction of a target object 420 at a time i-At before closest approach are blue-shifted relative to a pulse train that arrives when the satellite approaches the target most closely

(at time Ti) and red-shifted relative to the Ti pulse train after that time (i.e. at time Ti+At) . As a result, the same target object will appear in backscatter data sets at different frequency offsets from the closest approach data set. The backscatter data can be processed to compensate for this effect resulting in a data set that would otherwise require a larger antenna "aperture" - thus the name "synthetic aperture" .

[0065] As Figure 4B shows, the footprint 410 of the SAR beam moves with the platform giving a data set corresponding to a strip 450 (or "swath") of the Earth's surface. The arrangement of the antenna at an off-nadir angle means that the nadir path 440 is offset from this strip 450, the offset being termed the "range".

[ 0066 ] Data from the strip 450 is stored as an array of timestamped samples, often termed "rasters".

[ 0067 ] Figure 5A to 5C show certain processing operations applied to C-band SAR backscatter data in accordance with certain embodiments of the present disclosure.

[ 0068 ] Raw data from C-Band backscatter measurements is processed, in a number of phases, to generate a data set that accurately models changes in forest coverage. In a first phase (illustrated in Figure 5A) , the radar backscatter data rasters are stacked and masked (to excluding variables associated with effects other than backscatter from vegetation; such as non-forest areas, anthropogenic features or steep terrain) . The "cleaned up" data stacks are filtered to correct for statistical abnormalities in the time domain (illustrated in Figure 5B) . Many of the statistical abnormalities are due to intrinsic properties of the radar technique used - random scatter in radar data may be due to the synthetic aperture speckle arising from coherency and other noise. By observing the same geographical area on more than one occasion it is possible to correct for such random scatter, and distinguish between random changes and sudden, actual changes in backscatter caused by degradation or deforestation. The data for pixels corresponding to individual geographic locations in the corrected, clean data stacks are then categorized to determine whether vegetation at that location has suffered degradation during the time period spanned by the data (Figure 5C) .

[ 0069 ] The first phase (illustrated in Figure 5A) comprises inputting a stack 510 (i.e. an ordered array of rasters) of co-registered and terrain corrected C-band radar backscatter data (in this instance, both W and VH radar backscatter data) and digital elevation model (DEM) data (in metres) .

[0070] In respect of the DEM data, a DEM layer is extracted from the input data (i.e. "parsed") - operation S512. The DEM layer is then used to create a model of the slope of the underlying terrain - operation S514. A mask is prepared to mask out pixels of the data stack that correspond to gradients that exceed a threshold angle (e.g. gradient > 5 degrees) operation S516.

[0071] In respect of the radar backscatter data, metadata is retrieved - operation S522 - so that the backscatter data can be associated with a given location.

[0072] An elevation masked stack is created by applying the DEM layer mask to the radar backscatter data (the W data, say) - operation S524.

[0073] Where values in the stack fall below a threshold value (e.g. <1000 units) or represent an error code (e.g. "- inf") , that value is replaced by a null value (i.e. "removed") - operation S526.

[0074] The resulting stack data 520 (i.e. VV stack data) is thus cleaned of data from regions having steep inclines (since the SAR response to vegetation cover on such terrain is likely to be atypical; and steep slopes will generate no SAR data) and other anomalous data (such as might arise from areas of human habitation or open surface water) .

[0075] Further masking may be applied at this stage to limit the area of coverage to an area of interest of a party requesting the deforestation survey data. In certain embodiments, the entire data set is processed now that slope and erroneous data has been masked out (this previous stage being somewhat more computationally expensive) . Once processed the processed data stacks are conveniently archived to facilitate future access for data over this geographic region in the time frame covered by the data.

[0076] A forest cover mask is then created - operation S532 by comparing representative portions of the data (for example, successive 3 pixel by 3 pixel blocks of data in each raster) to a threshold value (which may be derived empirically and may also be confirmed experimentally - operation S540) . In the discussion of the following embodiment, it is assumed that the stack data being processed is W stack data and the rasters are W rasters, where the first W raster is denoted WO, the second Wl and so on. The reader will readily appreciate that data from an HH stack and/or a VH stack (i.e. cross-polarised backscatter data) could be used in place of the W stack and/or to augment the information in the W stack. One instance of such an alternative embodiment is described below in respect of Figure 9.

[0077] In certain alternative embodiments, HH or VH stack data is used in place of W stack data and references to the stack data being W stack data (and rasters being W rasters) in the discussion of Figures 5A to 5C may be substituted by references to HH or VH stack data and HH or VH rasters respectively .

[0078] In certain embodiments, a majority value window is used to allow a degree of approximation - operation S534. For the exemplary case of a three pixel by three pixel window, the presence of above-threshold values in seven or more pixels may be deemed as indicating that all nine pixels should be marked as "forest" in the forest cover mask. [0079] The resulting stack of W data 530 is thus corrected (by masking) for anomalous data, regions of steep gradients and for non-forest regions (such as rivers, exposed rock, savannah, or open grassland) . See Figure 6.

[0080] In the next phase, illustrated in Figure 5B, the "cleaned up" data stacks 530 are then temporally filtered. A separate stack of data is created in the same data structure as the W data stack using a nine by nine pixel window (for example) to determine the mean backscatter data - operation S536. Clearly, the number and indeed the arrangement of pixels in a window is not limited to a nine by nine pixel square, as illustrated in Figure 7, the reader will readily appreciate that any predetermined arrangement of pixels may be used .

[0081] The data for each of the original W rasters 510 (which have been preserved through the masking phase) are then divided by the respective mean W raster from the mean backscatter data stack - operation S538: the result being a mean ratio (i.e. original W raster / mean W raster) S542.

[0082] The original W raster 510 is then corrected (by multiplication by a weighting proportional to the mean ratio) - operation S544.

[0083] The standard deviations 550 of the corrected VV rasters 530 are calculated (in sequential groups of three rasters i.e. SD(W0, Wl, W2 ) ; then SD(W1,W2,W3) etc.) - operation S546.

[0084] Separately the mean difference 560 between neighbouring W rasters is calculated sequentially (i.e. mean (Wl) -mean (W2) , mean (W2 ) -mean (W3 ) etc.) . This provides a normalisation between rasters - operation S548. [0085] In a third phase (illustrated in Figure 5C) , data in the corrected rasters is categorized to determine whether vegetation at that location has suffered degradation during the time period spanned by the data.

[0086] The pixel difference between original neighbouring W rasters is calculated sequentially (i.e. W2-W3, W3-W4, etc.) operation S552. See Figure 8.

[0087] The pixel differences correspond to changes of backscatter data for respective individual geographic locations between the times of the respective W rasters.

[0088] A pixel is then classified as a " forest_change" pixel - operation S558 - when the pixel difference both: a) exceeds a threshold value multiplied by the standard deviation 550 for a sequential group of three rasters (e.g. for the two preceding rasters WO and Wl and the present raster W2 ) - tested at operation S554; and b) exceeds the mean difference 560 between corrected sequential W rasters - tested at operation S556. Otherwise the pixel is marked as "no change" operation S562.

[0089] The threshold value in operation S554 may be determined empirically - operation S580.

[0090] In certain embodiments, a further "changes" raster is generated to represent whether any given the pixel has changed over the time period of the stack. If more than one change is detected at a pixel over the period, that pixel is nevertheless represented as a changed pixel (values > 1 are represented as =1) - operation S564.

[0091] An output raster 570 is then generated. This output raster has a plurality of layers (i.e. "multiband": in the field of remote sensing, each layer may be referred to as a "band") . The multiband output raster provides one change period per band; the final band showing all changes detected over the time period spanned by the original data.

[0092] In addition, a metadata file is generated. The metadata includes a respective datestamp for each band number of the output raster.

[0093] In other words, certain embodiments provide snapshots of the cumulative changes in deforestation or degradation for a predefined geographic region over time.

[0094] An alternative embodiment would output a single band, with each pixel coded to the date at which a change was detected, with a separate metadata file stating the date to which each pixel code corresponds.

[0095] In certain embodiments, a positive classification as a " forest_change" pixel as a result of operations S554 and S556 may be used to provide a more detailed indication of the degree of change in forest coverage. To this end the single threshold in operation S554 may be replaced by more than one threshold corresponding to boundaries between categories of forest reduction - operation S566.

[0096] In a further aspect of the invention, additional operations are applied in the provision of a cleaned input stack of backscatter data. These may include operations that specifically exclude pixels whose backscatter characteristics imply the presence of bodies of water and/or buildings. Such operations facilitate the exclusion of rivers, lakes and other water bodies as well as forest areas that are within or near cities or towns, such as in urban parks etc.

[0097] Figure 9 illustrates a variation on the operations of Figure 5A. As in Figure 5A, rasters of raw data from C- Band backscatter measurements are stacked and masked (to excluding variables associated with effects other than backscatter from vegetation; such as non-forest areas, anthropogenic features or steep terrain) . The later phases of the algorithm in this variant, in which the "cleaned up" data stacks are filtered to correct for statistical abnormalities in the time domain and the data are categorized to determine whether vegetation at that location has suffered degradation during the time period spanned by the data, correspond to the operations illustrated in Figures 5B and 5C and are therefore omitted for brevity.

[0098] As in Figure 5A, a stack 910 of co-registered and terrain corrected C-band radar backscatter data (in this instance, both W and VH radar backscatter data) and DEM data is input.

[0099] The treatment of the DEM data in Figure 9 is identical to that in the algorithm of Figure 5A. A DEM layer is extracted from the input data ("parsed") - operation S912. The DEM layer is then used to create a model of the slope of the underlying terrain - operation S914. A mask is prepared to mask out pixels of the data stack that correspond to gradients that exceed a threshold angle (e.g. gradient > 5 degrees) - operation S916.

[00100] In respect of the radar backscatter data, metadata is retrieved - operation S922 - so that the backscatter data can be associated with a given location.

[00101] An elevation masked stack is created by applying the DEM layer mask to the radar backscatter data (the W data, as illustrated here, for instance) - operation S924.

[00102] Where values in the stack fall below a threshold value (e.g. <1000 units) or represent an error code (e.g. "- inf") , that value is replaced by a null value (i.e. "removed") - operation S926.

[00103] Where cross-polarised (say VH) backscatter data is obtained from the input stack S925, and where pixels in this cross-polarised data are determined to have backscatter characteristics consistent with the presence of a body of water those pixels are masked from the like-polarised (here, W) stack, operation S928. Water bodies present a relatively smooth surface and as a result incident radar radiation experiences low depolarisation and the overall power of the backscatter radiation is also low because a greater proportion of incident radiation is reflected in a direction away from the receive antenna.

[00104] The resulting stack data 920 (i.e. W stack data) is thus cleaned of data from regions having steep inclines (since vegetation cover on such terrain is likely to be atypical) and other anomalous data (such as might arise from areas of human habitation or open surface water) .

[00105] It will be appreciated that embodiments can be realized in the form of hardware, software or a combination of hardware and software. Any such software may be stored in the form of volatile or non-volatile storage, for example a storage device like a ROM, whether erasable or rewritable or not, or in the form of memory, for example RAM, memory chips, device or integrated circuits or on an optically or magnetically readable medium, for example a CD, DVD, magnetic disk or magnetic tape or the like. It will be appreciated that the storage devices and storage media are embodiments of machine-readable storage that are suitable for storing a program or programs comprising instructions that, when executed, implement embodiments of the present invention. [00106] Accordingly, embodiments provide a program comprising code for implementing apparatus or a method as claimed in any one of the claims of this specification and a machine-readable storage storing such a program. Still further, such programs may be conveyed electronically via any medium, for example a communication signal carried over a wired or wireless connection and embodiments suitably encompass the same.

[00107] Features, integers or characteristics described in conjunction with a particular aspect, embodiment or example of the invention are to be understood to be applicable to any other aspect, embodiment or example described herein unless incompatible therewith.

[00108] It will be also be appreciated that, throughout the description and claims of this specification, language in the general form of "X for Y" (where Y is some action, activity or step and X is some means for carrying out that action, activity or step) encompasses means X adapted or arranged specifically, but not exclusively, to do Y.

[00109] The description of the preferred embodiments of the present invention has been presented for purposes of illustration and description, but is not intended to be exhaustive or to limit the invention to the forms disclosed. It will be appreciated by those skilled in the art that changes could be made to the embodiments described above without departing from the broad inventive concept thereof. It is understood, therefore, that this invention is not limited to the particular embodiment disclosed, but covers modifications within the scope of the present invention as defined by the appended claims.