UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Application of the Getis statistic to the monitoring of riparian zones using mult-temporal radarsat images Eichel, Frank Herbert 2004

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-ubc_2005-0035.pdf [ 32.34MB ]
Metadata
JSON: 831-1.0075005.json
JSON-LD: 831-1.0075005-ld.json
RDF/XML (Pretty): 831-1.0075005-rdf.xml
RDF/JSON: 831-1.0075005-rdf.json
Turtle: 831-1.0075005-turtle.txt
N-Triples: 831-1.0075005-rdf-ntriples.txt
Original Record: 831-1.0075005-source.json
Full Text
831-1.0075005-fulltext.txt
Citation
831-1.0075005.ris

Full Text

APPLICATION OF THE GETIS STATISTIC TO THE MONITORING OF RIPARIAN ZONES USING MULT-TEMPORAL RADARSAT IMAGES by FRANK HERBERT EICHEL B.S.F., The University of British Columbia, 1979 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE in THE FACULTY OF GRADUATE STUDIES (Forestry) THE UNIVERSITY OF BRITISH COLUMBIA December 2004 © Frank Herbert Eichel, 2004 A B S T R A C T Monitoring the condition of a large number of dispersed riparian management areas on northern Vancouver Island is hampered by poor ground access, the high cost of flying, and a lack of resources. Furthermore, monitoring them with high spatial resolution satellites, such as IKONOS, QuickBird, and Terra, is frustrated by persistent cloud cover. A possible alternative is to use radar imagery to facilitate monitoring. Satellite-based radar sensors, including Canada's RADARSAT, generate their own illumination and consequently they are able to acquire imagery at any time of the day or night. The imagery can also be acquired during stormy weather since the radar satellite's illumination is not obstructed by clouds. Radar backscattering from forested areas typically exhibits low spatial autocorrelation because these areas are structurally quite heterogeneous. On the other hand, freshly logged cut-blocks appear to be comparatively homogeneous on radar images. Consequently backscatter from these areas should exhibit higher spatial autocorrelation. Due to similar characteristics exhibited by windthrown areas within riparian management areas and along cut-block boundaries, it was postulated that the Getis statistic could be applied to multi-temporal R A D A R S A T images to detect both the increased backscatter and higher spatial autocorrelation associated with windthrow damage in these areas. The areas that were examined in this thesis are located on northern Vancouver Island in a triangular area between Port Hardy, Port McNeill, and Port Alice. The area is entirely forested and is under active management for timber production. Depending on location, existing stands of mature timber are comprised of one or more species of western hemlock, western red cedar, and Sitka spruce. Regenerated areas are dominated by western hemlock and western red cedar with wetter areas being reforested with Sitka spruce. ii Four high spatial resolution R A D A R S A T images provided the subscenes used in this study. These images, with a pixel spacing of 3.125 metres and a resolution of approximately 8 metres, were acquired in Fine 2 Mode in December 1996, August 1997, November 1998, and March 1999 on identical ascending orbits under similar rainy conditions. To better understand the behaviour of the Getis statistic with radar data, it was first used to detect increased levels of backscatter resulting from clearcutting. It was then applied to riparian zones to determine if it could detect increased backscatter resulting from windthrow damage. Each subscene was converted to a series of Getis value images by passing five kernels ranging from 3x3 to 11x11 in size over the image. A single Getis value representing the highest local spatial autocorrelation was selected for each pixel from the five Getis value images and then written to a maximum Getis (MaxGetis) value image. The associated MaxGetis Distance images indicated that, for the data used in this research, clearcut areas were just as heterogeneous as forested areas. Therefore, without an improvement in data homogeneity, it was not possible to determine the extent of disturbances such as clearcutting and windthrow on the premise that they exhibited high spatial autocorrelation over more extensive areas than adjacent forest. However, the characteristics of the MaxGetis image were such that the difference between areas of high and low backscatter was significantly enhanced and consequently multi-temporal composites of MaxGetis images were found to be especially useful for visualizing areas exhibiting high backscatter. Although thresholding MaxGetis difference images provided a means for determining where significant increases in backscatter had occurred, its reliability for detecting fresh windthrow damage or any other structural change in the landscape for that matter was thrown into doubt by the high number of apparent false alarms appearing in areas where no disturbance was known to have occurred. These false alarms arose as a result of a 95% thresholding level being applied to the difference image to isolate significant increases in MaxGetis values. Raising the threshold to 99% substantially reduced the number of false alarms although a substantial number remained. This led to the conclusion that the nature of the data used in this research does not allow iii thresholded MaxGetis difference images to be used as a reliable means of detecting significant increases in backscatter due to some form of disturbance. This research also determined that backscatter emanating from windthrow areas is considerably weaker than it is from freshly logged areas, thus making it difficult to distinguish windthrow damage from forest. Consequently, it is suggested that further research is required to determine if other types of radar data, such as multi-look and fully polarimetric data, are more suitable for detecting windthrow damage. iv T A B L E OF CONTENTS A B S T R A C T :. ii TABLE OF C O N T E N T S v LIST OF T A B L E S vii LIST OF T A B L E S - APPENDIX A vii LIST OF F IGURES viii A C K N O W L E D G E M E N T S . xi Chapter 1 Introduction 1 1.1 Background... 1 1.2 Monitoring Windthrow 3 Chapter 2 Literature Review 10 2.1 Change Detection Techniques 10 2.1.1 Multi-Temporal Composites 13 2.1.2 Differencing 14 2.1.3 Ratio Images 15 2.1.4 Backscatter Change Detection 15 2.1.5 Classification 16 2.1.6 Coherence 17 2.1.7 Summary 18^ 2.2 Spatial Autocorrelation 19 2.3 Developing a Change Detection Technique using the Getis Statistic 22 Chapter 3 Methods 23 3.1 Location and Climate 23 3.2 Data Preparation 24 3.3 Registering the Datasets 28 Chapter 4 Evaluating Change Detection Approaches 29 4.1 Comparing Change Detection Approaches 29 4.2 Change Detection Using Established Techniques 32 4.2.1 Multi-temporal Composites 32 4.2.2 Image Differencing 34 4.2.3 Image Ratioing 36 v 4.2.4 Summary '. 41 4.3 Change Detection using the Getis Statistic 43 4.3.1 Getis Value Images 44 4.3.2 Maximum Getis Distance Images 68 4.3.3 Speckle Filtered Data 72 4.3.4 Change in the Next Period 75 Chapter 5 Testing the MaxGetis Differencing Technique 80 5.1 Strip 31 a 80 5.2 Strip 16 83 5.3 Edge Blowdown 87 5.4 Edge Blowdown at M5307 90 5.5 Other Areas 94 Chapter 6 Discussion 97 Chapter 7 Conclusion 104 7.1 Suggestions for Further Research 104 7.2 Summary of Results 106 References 109 Appendix A 115 vi LIST OF T A B L E S Table 1: Date and weather at time of acquisition of R A D A R S A T data used in this research 27 Table 2: Example showing which Getis Value is chosen as the MaxGetis Value 59 LIST OF T A B L E S - APPENDIX A Table A-1 . Feature-based threshold results for datasets processed by the C H D E T algorithm 116 Table A-2. Feature-based threshold results and summary statistics for Getis format data 117 Table A-3. Proportion of pixels in riparian zones and blowdown areas indicating significant change on MaxGetis difference images 118 LIST OF FIGURES Figure 1: Riparian zones are required alongside streams to protect aquatic and streamside ecosystems (Source: Riparian Management Area Guidebook, 1995, B.C. Ministry of Forests) 2 Figure 2: Portion of R A D A R S A T F2 image acquired on March 24, 1999 7 Figure 3: The riparian zone outlined in red on this image acquired on December 4, 1996 cuts through a cut-block whose boundary is indicated in yellow. Recent windthrow in the middle and lower end of the zone is revealed by the brighter pixels 8 Figure 4: Progressive windthrow damage incurred by a riparian zone in cut-block M5307 near Rupert Inlet. These R A D A R S A T F2 images were acquired on December 4, 1996 (top left), August 25, 1997 (top right), November 24, 1998 (bottom left) and March 24, 1999 (bottom right) 11 Figure 5: Map of study area showing footprint of the 4 R A D A R S A T images 24 Figure 6: 1995 aerial photograph showing the approximate boundary of the cut-block as taken from a logging completion map 30 Figure 7: 1998 (left) and 1997 images of the area in and around cut-block S501 31 Figure 8: Histograms of the 1997 and 1998 (3° images 32 Figure 9: Multi-temporal composite image showing areas of high backscatter in 1998 (red), 1997 (green), and 1996 (blue) 33 Figure 10: Difference image created from the 1997 and 1998 J3° data in intensity format and the corresponding histogram on the right. A normal curve has been overlaid on the histogram to show that the data in the image on the left is not normally distributed 35 Figure 11: Logarithmically scaled ratio image of the 1997 and 1998 amplitude data 36 Figure 12: Logarithmically scaled ratio image filtered with a 7 x 7 Kuan M M S E filter 37 Figure 13: Bitmaps corresponding to +1 standard deviation overlaid on the logarithmically scaled ratio image 39 Figure 14: Frequency distribution of the logarithmically scaled ratio image with the normal curve superimposed on it 40 Figure 15: Positive threshold bitmaps corresponding to a P F A of 5% overlaid on the logarithmically scaled ratio image. The bitmaps on the right are unthinned while those on the left have been thinned by setting the minimum number of neighbours to 3 in the C H D E T algorithm 41 Figure 16: Getis value images generated by a 3x3 kernel using the 1998 (left) and 1997 (right) radar brightness data 46 Figure 17: Getis value images generated by a 5x5 kernel using the 1998 (left) and 1997 (right) radar brightness data 46 Figure 18: Getis value images generated by a 7x7 kernel using the 1998 (left) and 1997 (right) radar brightness data 47 Figure 19: Getis value images generated by a 9x9 kernel using the 1998 (left) and 1997 (right) radar brightness data 47 Figure 20: Getis value images generated by an 11x11 kernel using'the 1998 (left) and 1997 (right) radar brightness data 48 Figure 21: Histograms of the 1997 and 1998 Getis value images generated by a 3x3 kernel 49 Figure 22: Profile of pixel values on the 1997 radar brightness image (black), 3x3 Getis value image (red), and 11x11 Getis value image (blue) 50 viii Figure 23: Profile of grey (pixel) values on 1998 radar brightness image (black), 3x3 Getis value image (red), and 11x11 Getis value image (blue) 51 Figure 24: Getis value difference images generated from the 3x3 data (left) and the 11x11 data (right).. 53 Figure 25: Frequency distributions of the 3x3 (left) and 11x11 Getis difference images with a normal curve overlaid on them in black 54 Figure 26: 3x3 (top) and 11x11 (bottom) Getis value difference images with 95% threshold bitmaps overlaid on them (left). Images on the right show the bitmap coverages after 5x5 mode filtering.... 55 Figure 27: 7x7 Getis value difference image with 95% threshold bitmaps on the left. On the right are the same bitmaps mode filtered by a 5x5 kernel 57 Figure 28: MaxGetis images generated from the 1998 (left) and 1997 source data 60 Figure 29. Histograms of the 1997 and 1998 MaxGetis images 61 Figure 30: 1997 (3° image (left) and the MaxGetis image of the same, area 62 Figure 31: 1998 (3° image (left) and the MaxGetis image of the same area 62 Figure 32: 1998-1997 MaxGetis difference image and corresponding histogram 63 Figure 33: 95% threshold applied to the 1998-1997 MaxGetis difference image (left). On the right are the same bitmaps mode filtered with a 5x5 kernel 64 Figure 34: Multi-temporal composites of MaxGetis and radar brightness images 65 Figure 35. 1998 (left) and 1997 Maximum Getis Distance images 69 Figure 36. 1998 (left) and 1997 coefficient of variation images derived from the M G D images 70 Figure 37: M G D values over a portion of the 1998 (left) and 1997 M G D images 71 Figure 38: Net changes in maximum Getis distance from 1997 to 1998 72 Figure 39: MaxGetis difference image derived from speckle filtered data (left) with 95% threshold bitmaps overlaid on the right 74 Figure 40: MaxGetis threshold bitmaps from Figure 39 mode filtered with a 5x5 kernel 74 Figure 41: Difference image (left) generated from 1997 and 1998 11x11 Getis value images and overlaid with the 95% threshold bitmaps on the right. The 11x11 Getis value images were generated from speckle filtered data '. 75 Figure 42: 1999 (left) and 1998 (3° image of the cut block area 77 Figure 43: 1999 (left) and 1998 MaxGetis images of the cut-block area 77 Figure 44: Multi-temporal composite of the 1999 and 1998 MaxGetis images of the cut-block area (left) arid the 1999-98 difference image on the right 78 Figure 45: Multi-temporal composite of the 1999, 1997, and 1996 (3° images (left) and the equivalent as a MaxGetis multi-temporal composite 81 Figure 46: 1999-96 MaxGetis difference image (left) overlaid with 95% threshold bitmaps. On the right, the same difference image has been overlaid with the mode filtered 1999-96 bitmaps in red and the mode filtered'bitmaps from the 1997-96, 1998-97, and 1999-98 difference images as cyan, magenta, and yellow polygons 82 Figure 47: Aerial photo showing location and condition of Strip 16 in 1995 83 Figure 48: August 1997 LANDSAT image (left) and November 2000 A S T E R image confirming the decimation of Strip 16 84 Figure 49: Multi-temporal composite of the 1999, 1997, and 1996 p° images 85 ix Figure 50: 1999-96 MaxGetis difference image with mode filtered 95% threshold bitmaps overlaid in red. The bitmaps have been overlaid on the aerial photo on the right 86 Figure 51: 1997 LANDSAT image showing location of new cut-block and road right-of-way leading to it. 87 Figure 52: Multi-temporal composite of the R A D A R S A T B° images showing the same blowdown damage in the northeast corner of the cut-block 88 Figure 53. 1997-96 MaxGetis difference image of the blowdown area overlaid with mode filtered 95% threshold bitmaps 89 Figure 54: 1999-96 MaxGetis difference image (left) and 1998-97 MaxGetis difference image with mode filtered 95% threshold bitmaps 90 Figure 55: Multi-temporal composites of the 1999, 1997, and 1996 MaxGetis (left) and p° (right) images of the northeast corner of cut-block M5307 91 Figure 56: The 1999-1997-1996 MaxGetis multi-temporal composite overlaid with threshold bitmap vectors derived from various MaxGetis difference images to highlight areas of significant change. On the right are the same bitmap vectors after 5x5 mode filtering 92 Figure 57: The 1999-1998-1997 MaxGetis multi-temporal composite overlaid with the 1999-97 (white), 1999-98 (cyan), and 1998-97 (black) threshold bitmap polygons to highlight areas of significant change on the 1999 image. On the right are the same bitmap polygons after 5x5 mode filtering.... 93 x A C K N O W L E D G E M E N T S I am indebted to my principle advisor, Peter Murtha, for all of the data and field observations used in my research, as well as for his patience and encouragement as I explored various avenues and reworked the data. My gratitude also goes to Mike Wulder for advice regarding the application of the Getis statistic in this research, to Edmond Nezry of PR IVATEERS N.V. who speckle filtered the data and to Peter Marshall and Don Leckie, my academic advisor and research advisor respectively, for their advice, support and encouragement. Mike Wulder and Don Leckie are both with the Pacific Forestry Centre in Victoria. The R A D A R S A T data, provided to Professor Murtha under A D R O Project 384 and used in this research, are copyrighted by the Canadian Space Agency. xi Chapter 1 Introduction 1.1 Background Under the British Columbia Forest Practices Code, a strip of vegetation must remain alongside any stream that passes through or beside an area where timber has been harvested. This so-called riparian management area consists of a riparian management zone that is up to 50 metres wide. It may also include an additional 20-30 metre wide riparian reserve zone if the stream is in a community watershed or if it provides habitat for fish. Figure 1 illustrates the position of the two zones in a situation where the stream forms the boundary for a cut-block on the right. Partial cutting is permitted in the riparian management zone; however, no harvesting of any kind is permitted in the riparian reserve zone (BCMOF and BC Environment, 1995). 1 Riparian Management Area Ripa.'ion Riparian Reserve Manac-Kmenl Riparian Maragemfint Area Figure 1: Riparian zones are required alongside streams to protect aquatic and streamside ecosystems (Source: Riparian Management Area Guidebook, 1995, B.C. Ministry of Forests). The function of these zones is to protect and sustain streamside and aquatic ecosystems. Undisturbed vegetation alongside streams is essential for maintaining water quality and preventing downstream sedimentation of spawning areas. Both zones provide food to a diverse range of resident animals and birds, while also serving as wildlife corridors. The riparian reserve zone in particular provides essential food and nutrients to fish and other aquatic inhabitants of the stream. On the coast of B.C., high winds and heavy rains are typical of winter storms. During such storms, individual trees in fully stocked stands rely on others around them to protect them from windthrow. However, once this protection is removed by timber harvesting, trees in riparian zones and along cut-block boundaries that become 2 exposed to high winds are at risk of windthrow, especially if the soil is unstable or saturated. The resulting damage can either be catastrophic or progressive. In either case, the protective function of the riparian zone may be seriously impaired or lost altogether. 1.2 Monitoring Windthrow The problem of windthrow in riparian zones and along cut-block boundaries is known to be quite significant on northern Vancouver Island (e.g., Mitchell, 1995, Mitchell et al., 2001). Although little can be done to forestall windthrow in existing riparian zones, there are two reasons for continuing to monitor their condition: • An understanding of soil types, seasonal weather patterns, prevailing winds and the influence that topography has on the latter can be integrated with data on the distribution and severity of windthrow damage in riparian zones. This information can be used to design cut-blocks and associated riparian zones so that the zones are more resistant to windthrow. • The data can be used by licensees and the Ministry of Forests to identify windthrown timber for possible salvage that could otherwise be lost to decay or fire. This reduces the potential fire hazard and prevents insect infestations from spreading from downed timber to nearby stands. Continuing research on the problem and the attendant development of practical guidelines to minimize windthrow damage in riparian zones established in the future in this area would be facilitated if the condition of the existing zones could be assessed at frequent intervals. However, frequent bad weather and the requirement to deactivate roads once they are no longer needed to access timber hampers monitoring the condition of riparian zones from the air and on the ground. Persistent cloud cover also makes it difficult to acquire aerial photography and imagery from satellites equipped with sensors operating in the visible portion of the electromagnetic spectrum. Despite 3 the ability to acquire 170 km by 183 km images of the area every 16 days at a very low cost (USD 600 per scene), only 14% of the LANDSAT images acquired over a period of 6 years have been cloud-free and usable for monitoring (Peter Murtha, personal communication)1. Similar constraints would prevent the acquisition of imagery by higher spatial resolution satellites such as IKONOS, QuickBird, and T E R R A that also rely on incident sunlight. A possible solution to the monitoring dilemma is offered by RADARSAT, a radar satellite that is capable of imaging the same area from the same orbit every 24 days (RADARSAT International, 2000) 2. Unlike LANDSAT 7 which relies on incident light, R A D A R S A T provides its own illumination in the ' C microwave band (at 5.3 Ghz) and is able to acquire an image day or night under any weather conditions. It also has a higher resolution (~8 metres) than LANDSAT 7 (15 metres panchromatic, 30 metres multi-spectral) and thus is potentially capable of detecting change at a smaller scale than LANDSAT. However, of greater significance is that radar data provide the ability to readily detect structural (physical) changes in the scene due to the nature of the electro-magnetic energy transmitted by the sensor and its interaction with scene elements. Radar is actually an acronym for RAdio Detection And Ranging. Imaging radars are a specific type of radar that transmits electromagnetic waves as pulses with specific characteristics at an oblique angle towards a target area, such as the earth's surface, and receives the reflected energy or backscatter returning from that same area. Digital signal processing of the backscatter data uses the round-trip time and incidence angle of the transmitted electromagnetic pulses and their characteristics to construct an image (Henderson and Lewis, 1998; R A D A R S A T International, 1995). Variations in the strength or magnitude of backscatter are dependent on surface roughness, topography, inherent reflectivity of a target, surface moisture and the 1 Professor Emeritus, Faculty of Forestry, University of British Columbia. 2 R A D A R S A T has since been renamed RADARSAT-1 once development of RADARSAT -2 was announced, 4 moisture content of the target. Smooth water surfaces act as specular reflectors in that the electromagnetic wave is reflected in the direction opposite to that from whence it came. As the water surface roughens, more of the energy is reflected back towards the sensor. Interaction with a vegetation canopy is considerably more complex as a pulse may reflect off the surface of leaves and branches many times before exiting the canopy altogether. Some of the reflected energy returns to the sensor thus producing variations in tone in the resulting image. The degree of penetration into the canopy depends on the wavelength of the electromagnetic energy; the longer the wavelength (i.e. the lower the frequency of the electromagnetic wave), the deeper the penetration. C-band radars, of which R A D A R S A T is one, have limited canopy penetration. Surface moisture and vegetation moisture content also determine the strength of backscatter with wet vegetation returning more energy that dry vegetation. It is also important to note that man-made objects such as buildings are very strong reflectors and consequently they typically appear as bright objects (Henderson and Lewis, 1998; R A D A R S A T International, 1995). Topography also has a significant influence on the proportion of the incident energy returning to the sensor. Slopes facing the sensor will reflect more energy than slopes facing away from the sensor. Depending on the incidence angle, so-called radar shadows can be induced in slopes facing away from the sensor if they are sufficiently steep. Shadows can also by cast by trees, such as along the west boundary of some the cut-blocks in the images used in this research (Henderson and Lewis, 1998; R A D A R S A T International, 1995). RADARSAT. launched in November 1995, is Canada's first earth observation satellite. It has an electrically steerable antenna which allows it to collect data at incidence angles ranging from 10° to 60°. Depending on the level of detail and area of coverage required, backscatter data can be collected in one of seven beam modes: Fine, Standard, Wide, ScanSAR Narrow, ScanSAR Wide, Extended High, and Extended Low. Each beam mode has one or more beam positions defined by a near and far range incidence angle which in turn determines the swath width. The resolution of the data 5 collected by each beam mode ranges from less than 10 metres for Fine mode to 100 metres for ScanSAR Wide mode. The base resolution is determined by the incidence angle and pulse duration. Signal processing techniques make use of the changing characteristics of the pulse over its duration to obtain a higher spatial resolution image than what would otherwise be obtainable from the pulse duration alone (RADARSAT International, 1995; 2000). The single-look R A D A R S A T Fine 2 mode data used in this research was acquired for A D R O (Application Development and Research Opportunity) Project 384 at an incidence angle range of 39-42°. The resolution of the data is 7.9 metres in the range direction (perpendicular to the satellite's orbital track) and 8.4 metres in the azimuth (along-track) direction. Four-look Standard 4 mode data were also obtained much less frequently for the same project at an incidence angle range of 36-42°. The resolution of these data are much lower than the Fine mode data at 23 metres in the range direction and 27 metres in the azimuth direction (RADARSAT International, 2000; Murtha, 2000c). The Fine 2 mode data were chosen for this research for their higher spatial resolution and therefore greater opportunity to capture elevated levels of backscatter from small scale windthrow events. Figure 2 shows part of a R A D A R S A T Fine 2 mode (F2) image acquired on March 24, 1999. The image was acquired on an ascending orbit with the illumination coming from the left. Forested areas are highly textured and dark grey in colour. Recently harvested cut-blocks, soggy from frequent and lengthy rainstorms occurring before and during the acquisition of this image, appear somewhat brighter. It is speculated that the slash and stumps, in combination with the exposed micro-topography, are responsible for a greater proportion of the incident energy being reflected back to the sensor thereby making the cut-blocks appear brighter than the surrounding forest. Murtha (2000b) has also reported a strong correlation between scene reflectivity and the type of landform upon which cut-blocks are situated. Cut-blocks situated on till-derived landforms are typically darker than those located on lacustrine sediments, even when very wet. In 6 other environments with less slash remaining after logging, cut-blocks are darker than adjacent forest areas, even under wet conditions (Leckie, 1997; Leckie etal., 1998). Figure 2: Portion of RADARSAT F2 image acquired on March 24,1999. Riparian zones located inside very wet cut-blocks typically appear on radar imagery as narrow, dark ribbons. The trees within the zone appear dark as a result of canopy volume scattering and radar shadows cast by the trees. Damage in them is revealed by brighter pixels in the areas where windthrow has occurred (Murtha, 1997) (Figure 3). 7 Figure 3: The riparian zone outlined in red on this image acquired on December 4,1996 cuts through a cut-block whose boundary is indicated in yellow. Recent windthrow in the middle and lower end of the zone is revealed by the brighter pixels. The objective of this thesis is to examine the effectiveness of a spatial autocorrelation statistic, called the Getis statistic, for detecting and defining the extent of windthrow damage in riparian management areas and along cut-block boundaries. Its performance will first be compared to that of existing methods for detecting structural changes in the landscape associated with the appearance of a freshly logged cut-block on a radar image acquired in November 1998. It is postulated that a patch of windthrown timber exhibits a signature similar to that of a freshly logged cut-block and consequently the Getis statistic should provide similar results. The remainder of this thesis is laid out as follows: • Chapter 2 reviews existing methods used for change detection and briefly describes the application of the Getis statistic to multispectral and passive microwave data. • Chapter 3 describes the study area, and the data and their preparation for analysis. 8 • Chapter 4 demonstrates the application of existing change detection methods and the Getis statistic for detecting change associated with the appearance of a new cut-block on a 1998 radar image. • Chapter 5 illustrates the application of the MaxGetis differencing technique to the detection of windthrow in riparian zones and along cut-block boundaries. • Chapter 6 discusses the techniques used and the issues that influence the success of the Getis statistic for change detection with the data used in this research. • Chapter 7 summarizes the results and makes suggestions for further research. 9 Chapter 2 Literature Review 2.1 Change Detection Techniques Monitoring the condition of riparian zones with the aid of multi-temporal radar data is essentially a change detection problem. Windthrow damage in these zones is revealed by the appearance of bright pixels on a radar image acquired after the event. Figure 4 illustrates progressive windthrow damage in the riparian zone outlined in Figure 3. The bright pixels within the riparian zone on the December 1996 image suggest windthrow damage has already occurred at the very top, middle, and bottom of the zone. Subsequent images show a progressive increase in the number of bright pixels in the riparian zone. That this is due to progressive windthrow damage has been independently confirmed in the field (Peter Murtha, personal communication)3. 3 Professor Emeritus, Faculty of Forestry, University of British Columbia 10 Figure 4: Progressive windthrow damage incurred by a riparian zone in cut-block M5307 near Rupert Inlet. These RADARSAT F2 images were acquired on December 4, 1996 (top left), August 25, 1997 (top right), November 24, 1998 (bottom left) and March 24, 1999 (bottom right). Since windthrow is usually progressive, assessment of the severity of the problem in a particular riparian zone is best achieved using a series of images. A wide range of techniques have been developed for detecting change between two or more imaging dates. These include (Khorram, 1999; Jensen, 1996): 1. Change Detection Using Write Function Memory Insertion. This method utilizes individual bands of data acquired on 2 or 3 acquisition dates to create a two- or 11 three-colour composite. Colour additive theory is then used to discern the location and extent of changes occurring between one date and the next. 2. Image Algebra Change Detection. Change is identified between two dates of imagery by using simple image algebra, either by band ratioing or subtraction, followed by thresholding, to identify areas of change. 3. Multi-date Composite Image Change Detection. Identical bands from two dates of imagery are combined into a single composite dataset. Principal components analysis or classification (either supervised or unsupervised) is then employed to identify clusters of pixels that can be correlated with change occurring between the two dates. 4. Post-Classification Comparison Change Detection. Images acquired on two dates are classified according to the same rules and the two classified images are then compared using a GIS matrix analysis technique. 5. Change Detection Using a Binary Change Mask Applied to Date 2. This method involves classifying the "before change" image (Date 1). A band from each of the Date 1 and Date 2 images are brought together in a new dataset and image algebra is used to locate areas of change. The area of change is then used as a mask on the Date 2 image and the area under the mask is then classified. This method has the advantage that it identifies the location of the change and identifies the cover type for each area of change. 6. Manual On-Screen Digitization of Change. Change can also be identified by "heads-up" digitizing with both the "before change" and "after change" images displayed side by side on a computer monitor. All of these techniques were initially developed for detecting change on multispectral images. Despite the high level of inherent noise in radar data, variants of these techniques have been successfully adapted for detecting change with synthetic aperture radar data. These techniques are described in more detail in the following sections. 12 2.1.1 Multi-Temporal Composites The simplest means of detecting change between imaging dates is to construct a three colour red-green-blue (RGB) or multi-temporal composite. Three images acquired on three different dates are assigned in descending order to the red, green, and blue colour guns of a video monitor. Colour additive theory is then used to identify the interval during which the increase in backscatter associated with change has occurred. The duration and timing of the increase in backscatter determines the predominant colour of each individual pixel in the composite. If the pixel is bright on only one image, its colour will be the colour assigned to that particular image in the composite. Areas that are predominantly red indicate that the increase in backscatter occurred between the two most recent imaging dates. Areas that are predominantly yellow indicate that such an increase was sustained over the last two imaging dates. Areas that are bright on the images assigned to the red and blue colour guns will appear magenta while those bright on the images assigned to the blue and green colour guns will appear cyan. Areas bright on all three imaging dates will appear white. Decreases in backscatter can occur as a result of vegetation blocking the reflective surface. This technique has been used extensively to detect windthrow damage in riparian zones on northern Vancouver Island (Murtha, 1998a, 1998b, 1998c, 2000a, 2000b; Murtha and Mitchell, 1998). Both Radarsat Fine 2 mode and Standard 4 mode have been used, with the former being more effective for detecting damage. Data acquired during both the summer when the associated cut-blocks are relatively dry and the winter when the cut-blocks are relatively wet provided the best combination for detecting damage (Murtha, 1998b). The relative brightness of cut-blocks imaged on the same date has since been found to depend a great deal on the landform upon which the cut-block sits (Murtha, 1998b; Murtha, 2000c). A plausible link between landform and susceptibility to windthrow damage has also been established (Murtha and Mitchell, 1998). 13 Multi-temporal composites have also been used to detect landscape level change associated with widespread fires on Borneo. Pre- and post-fire multi-temporal composites were assembled for mapping fire damage using two speckle filtered images (pre-fire image assigned to the green gun, post fire image assigned to the blue gun) and channel 2 of a principal components analysis assigned to the red gun. Speckle is the 'noiselike' characteristic of images produced by coherent imaging systems such as the synthetic aperture radar system used to collect the R A D A R S A T data discussed in this paper. It is a natural outcome of the incident electromagnetic wave interacting with point and distributed targets in the scene. Distributed targets, such as forested areas, are comprised of a great many discrete point scatterers. As electromagnetic waves interact with a distributed target, each scatterer generates a backscattering wave that exhibits a phase and amplitude that differs from that of the incident wave. The total returned modulation of the incident wave is the summation of these individual backscattered waves less any attenuation caused by propagation and scattering in the troposphere as the wave travels back to the receiving antenna. The constructive and destructive interference inherent in the summation of backscattered waves of different amplitude and phase within a resolution cell results in the variable brightness between adjacent pixels that is characteristic of speckle (Oliver and Quegan, 1998). The extent of the damage could be mapped after the pre- and post-fire composites had been brought into a GIS environment. Although the images were filtered to remove the effects of speckle, this method preserved the geomorphology and texture information at a relatively high spatial resolution (Ruecker and Siegert, 2000; Siegert et al., 1999; Siegert and Nakayama, 2000; Siegert and Rucker, 1999). Multi-temporal composites have also been used to locate deforested areas in Brazil (Kux et al., 1998). 2.1.2 Differencing A very common practice used with multispectral imagery is to simply subtract an older image from a more recent image after the two images are registered. However, in a study involving single-look complex data in intensity format, Rignot and van Zyl (1993) found that changes will not be detected to the same degree in high intensity regions as in low intensity regions. This is due to the nature of the probability distribution function 14 of the difference image which is sensitive to the relative size of the difference between the mean intensity values of the two images as well as to the magnitude of the mean intensity level in the reference image. They also found that, because the variance of the distribution of the difference increases with intensity level, the probability of error in detecting change is higher in high intensity regions than in low intensity regions. In their view, both of these problems were sufficient to rule out the application of the difference method to radar data. 2.1.3 Ratio Images An alternative to the difference method is the ratio method. The ratio method has a number of advantages over the difference method, the most significant one being that it is unaffected by calibration errors which are typically the same in repeat-pass imagery (Rignot and van Zyl, 1993; Oliver and Quegan, 1998). Dekker (1998) used the ratio method for detecting change using repeat-pass ERS-1 data. His approach is discussed at length in section 4.2.3. Weydahl (1991; 1992) used ratio images and thresholding to detect significant changes from one imaging date to the next using ERS-1 data. The ratio of two amplitude images was then converted to dB values and a probability of false alarm was computed to determine the appropriate threshold values. ERS-2 data have also been used in the Brazilian Amazon to detect change. However, the thresholded ratio image produced after the intensity (power) images were averaged and speckle filtered was not as complete as a thresholded Landsat image acquired about the same time. Much of the difference was attributed to localized topographic effects on the strength of the radar returns (Oliver and Quegan, 1998; Grover et al., 1999). 2.1.4 Backscatter Change Detection A number of researchers have focused on parcel- or region-based change detection. Couturier et al. (1999) mapped the extent of fire damage in the eastern part of Borneo using intensity differences derived from regions identified on multi-temporal E R S and 15 J E R S images. Ruecker and Siegert (2000) did much the same using a much larger number of ERS-2 images. Two other approaches stand out. In the first, Bao (1999) applied an autoregressive model to a set of images from which speckle noise had been removed using a 3-dimensional wavelet shrinkage algorithm. He discovered that changes beyond a given range in the value of the two coefficients of a second order autoregressive model could be used to identify significant changes in backscattering intensity at the pixel level. A second approach that has been developed segments speckle-filtered images into objects of homogenous intensity using disc shaped templates fitted into regions generated by an edge filter. Two segmented images are then overlaid and areas that have the same intensity level are removed. Other areas are also removed if the difference in intensity in the overlapping area is within a defined tolerance (White, 1991). Areas left over are deemed to have experienced significant change. 2.1.5 Classification Various classification approaches have also been developed for change detection. A straightforward approach simply classifies two images, one acquired before the change of interest has occurred and the other afterwards, and maps the differences by overlaying the two classifications in a GIS environment (Dutra et al., 1999). However, more sophisticated approaches have taken advantage of the growing number of radar images available. Bruniquel et al. (1999) developed a method whereby a maximum likelihood classification map was extracted from a segmented image whose pixel values represented the widest difference in intensity among 8 ERS-1 and ERS-2 single-look complex images. Overlaying the maximum intensity difference image over a mean image derived from the source data enhanced visualization of these intensity differences. Other researchers have extended the ratio method to include classification of homogenous and near homogenous regions of intensity change. Siegert and Rucker 16 (1999) constructed single date R G B composites from two ERS-2 images acquired over eastern Borneo, one acquired before extensive fire damage occurred in the area and the other after. A Gamma MAP and 20x20 texture filtered image was assigned to the red gun, a Gamma MAP filtered image was assigned to the green gun and a Gamma MAP and 10x10 texture filtered image was assigned to the blue gun to create each composite. A ratio composite was generated to highlight the areas of significant change. A maximum likelihood classifier was then used to classify the ratio image. A more complicated approach involved the application of a Bayesian classifier to segment multitemporal radar data into regions of homogenous and similar backscatter changes (Rignot and Chellappa, 1992). Neural nets have also been applied to the whole process of change detection and, in a more restricted sense, to detecting and identifying objects that are indicative of significant change (White, 1991). 2.1.6 Coherence Repeat-pass SAR imaging has facilitated the application of interferometry to detect change. The magnitude of the amplitude and phase correlation between two or more radar images of the same area is referred to as the degree of coherence or coherence magnitude (Touzi et al., 1999). The coherence magnitude is sensitive to movements of scatterers within a resolution cell amounting to about a half a wavelength of the radar signal (Bruniquel et al., 1999). Coherence is also significantly affected by topography and by both seasonal and local weather conditions (Castel et al., 2000). In general terms, coherence is low over forested areas and high over areas where scatterers are more stable over time. S A R interferometry uses single-look complex data acquired at two or more closely separated times. The satellite tracks themselves must also be close together, typically less than 600 metres (Coulson, 1995). The required orbital separation has been achieved with the ERS-1 satellite. However, the long 35-day repeat cycle usually results in significant decorrelation of both amplitude and phase. In recent years this problem has been overcome by operating the ERS-1 and ERS-2 satellites in tandem mode. In this mode, the two satellites acquire data of the same area a day apart (Coulson, 1995). 17 In one study the accuracy achieved in separating forested from non-forested areas with one ERS-1/2 tandem pair was 94%. Change detection could be achieved with two successive tandem pairs; however, the significant decorrelation over the 35-day interval mentioned earlier contributes to a lower accuracy in separating forested from non-forested areas (Castel et ai, 2000). Repeat pass interferometry with Radarsat is not always possible because the satellite's orbital position is not accurately known or as well controlled. However, it has been determined that interferometry can be done with data acquired during a period of minimum longitudinal drift which occurs halfway through the 90 day interval between orbital boosts. Radarsat's fine beam mode is the most useful for interferometry as relatively large baselines are possible with the large chirp bandwidth and the wide range of incidence angles (Armour et al., 1997). Consequently it should be possible to separate forested areas from non-forested areas (i.e. clearcuts and withthrown areas) using one pair of Radarsat single look complex images acquired at the appropriate times. 2.1.7 Summary The techniques described above are generally successful at separating forested areas from non-forested areas. However, the search continues for new and perhaps more successful approaches to the problem of reliably detecting structural change using multi-temporal radar data. One approach that has yet to be applied to multi-temporal radar data is the use of spatial autocorrelation to detect changes in the landscape. The next section describes how spatial autocorrelation has been applied for detecting change in multi-temporal optical imagery. 18 2.2 Spatial Autocorrelation Spatial autocorrelation is defined by Goodchild (1986) as the degree to which objects at some place on the earth's surface are related to other objects located nearby. Chrisman (1997) defines it more precisely as the degree of correlation between the value of a variable and that of its neighbours. He also defines it as the tendency of spatial data to vary smoothly over distance. Global measures of spatial autocorrelation, such as Moran's I index and Geary's c index, generate a single value that expresses the strength of the spatial autocorrelation of a variable over a whole region. In the context of this study, the region would be the entire image. Local measures of spatial autocorrelation, such as the local form of Moran's I index and the Getis statistic, provide a measure of spatial autocorrelation for each location of that variable using one or more defined neighbourhoods (Anselin, 1995; Ord and Getis, 1995). The Getis statistic is of particular interest as it can be used to assess the strength of inter-pixel relationships in defined neighbourhoods and indicate whether or not spatial autocorrelation is strictly local or more extensive in nature (Wulder and Boots, 1998, 2001; LeDrewefa/ . , 2000). There are two forms of the Getis statistic, Gi (d) and G i * (d). Gi (d) excludes the value of the variable, /, while G i * (d) includes it. The "d" refers to the radius of the neighbourhood around the variable for which the Getis statistic is being computed. The Getis statistic was initially developed by Getis and Ord (1992) for use with point data. In that context, a radius of 10 metres around a particular point on a map that indicates the severity of a particular disease might, for example, enclose another four points with the same type of data. The Gi form of the Getis statistic would use the scores of those four points to compute the Getis statistic for the point at the centre of the neighbourhood whereas the G i * form would utilize all five disease scores to compute the Getis statistic for the point at the centre of the neighbourhood. A standardized version of the second form, Gi* , was later derived in the context of its application to raster data (Ord and 19 Getis, 1995). The G i * statistic (hereinafter referred to as the Getis statistic) is stated as follows (Ord and Getis, 1995): G*(d) = s{[(nS'u)-W;2]/(n-l)} j 1/2 (1) where (2) with x and s being the usual mean and standard deviation of the dataset, ^ w ^ ( J ) x y being the sum of the products of the weight assigned to the pixel times the corresponding pixel value for all pixels within distance d of /, including that of /' itself (i.e., the entire neighbourhood), W* being the sum of the weights assigned to all pixels within the neighbourhood, 5,* being the sum of squared weights assigned to all pixels within the neighbourhood, and n being the total number of pixels in the dataset. In its application here, binary weights are used where the weight is one for all pixels within distance d of /' and zero otherwise. In its application to raster data, such as satellite imagery as discussed here, the Getis statistic is normally computed for each pixel in an image using a sampling window or kernel passed over the entire image. A 3x3 kernel is used to compute the Getis statistic for the centre pixel in the kernel using the values of all of the pixels within a distance, d, of 1 pixel as measured from the centre pixel. A 5x5 kernel is used for a of of 2 pixels while 7x7, 9x9, and 11x11 kernels are used to compute the Getis statistic using all of the pixels within a d of 3, 4, and 5 pixels of the centre of the kernel. A total of 9, 25, 49, 81, and 121 pixels are therefore used to calculate the Getis statistic for the centre pixel in a 3x3, 5x5, 7x7, 9x9, and 11x11 kernel respectively. When a particular kernel size is used to compute the Getis value, a new image containing the Getis value for each pixel can be generated. High positive values of this 20 statistic indicate a clustering of high pixel values in the kernel whereas high negative values indicate a clustering of low pixel values. Clustering of the highest and lowest Getis statistics in the image reflect high local spatial autocorrelation among the associated pixels in the source data. However, a weakness of the Getis statistic is that it cannot be used to identify the strength of spatial autocorrelation among pixels with intermediate values since a mid-range Getis value could be the result of a clustering of intermediate valued pixels or a clustering of pixels that have a mix of values (Wulder and Boots, 1998). If the statistic is calculated using the five kernels mentioned above, a maximum Getis or MaxGetis value can be chosen from among the five Getis values computed for each pixel. The MaxGetis value indicates both the magnitude of the source data and the strength of local spatial autocorrelation of the data within the kernel. The size of kernel that generated the MaxGetis value indicates the extent of the local spatial autocorrelation around the pixel in the source data. The MaxGetis value is chosen by first ranking the Getis values by kernel size, converting each to an absolute value and then comparing the values in successive pairs starting with the values generated by the 3x3 and 5x5 kernels. As soon as the magnitude of the first absolute value exceeds that of the second, the comparison stops. The actual Getis value corresponding to the bigger absolute value is then written to a MaxGetis image. For example, if the ordered list of Getis values for a particular pixel are -6, - 9 , 1 , 4 , 1 1 , the Getis value generated by the 5x5 kernel (-9) would be chosen as the MaxGetis value since it is followed by a smaller absolute value. By choosing the first maximum Getis value, maximum spatial autocorrelation is determined in the smallest possible neighbourhood or spatial domain (Mike Wulder, personal communication)4. The distance corresponding to the kernel that generated the MaxGetis value is written to a maximum Getis distance or MGD image. Aggregations of the same maximum Getis distances in the MGD image indicate areas of homogeneity in the source data whereas a mixture of these distance values indicate some degree of heterogeneity (Wulder and Boots, 1998). 4 Research Scientist, Pacific Forestry Centre, Natural Resources Canada, Victoria, B.C. 21 The Getis statistic has been used in a number of remote sensing applications. It was used with LANDSAT TM data to assess the within channel and within species spatial autocorrelation trends in a managed forest area in New Brunswick (Wulder and Boots, 2001) and with passive microwave imagery to determine maximum hemispheric snow extent and degree of snow cover variability (Derksen et al., 1998). It was also used to detect the extent of flooding in China (Chen etal., 1999). In the context of change detection, it has been used to differentiate healthy coral reefs from dead coral reefs using multi-temporal S P O T imagery (LeDrew et al., 2000; Holden et al., 2000). Being heterogeneous, the maximum Getis statistic for healthy coral areas was achieved using a 3x3 kernel. Dead or damaged coral appears as homogeneous areas on S P O T imagery and consequently the maximum Getis statistic in these areas was obtained by using a 9x9 kernel. In other words, a decline in the health of the coral was revealed by a widespread increase in the size of kernel responsible for generating the maximum Getis statistic. 2.3 Developing a Change Detection Technique using the Getis Statistic It is postulated here that the standardized version of the Getis statistic can be used with multi-temporal Radarsat data to detect windthrow damage in riparian zones. The rationale for this is as follows. Cut-blocks on northern Vancouver Island are typically bright on R A D A R S A T images when the soil and overlying slash are soggy at the time of image acquisition (Figure 3 and Figure 4) (Murtha, 2000c). Under such conditions, pixels in an area of standing timber are darker than those in an adjacent clearcut area. When the standing timber is windthrown, the pixels at the location of the windthrow are noticeably brighter in an image acquired after the windthrow occurs (Figure 4). This localized change in brightness should be revealed by an increase in the maximum value of the Getis statistic computed for the affected pixels. It may or may not be accompanied by a change in the size of the kernel responsible for deriving this higher value. 22 Chapter 3 Methods The study area for this thesis was dictated by both the availability of multi-temporal radar data and the fact that monitoring the condition of riparian zones was already the subject of a long-term study funded by Grant #384 under the Canadian Space Agency Application and Development Research Opportunity (ADRO) programme (Murtha, 1997, 1998a, 1998b, 1998c, 2000a, 2000b, 2000c). The study area on Northern Vancouver Island is ideal for monitoring with radar data in particular because it is frequently obscured by cloud cover, making it impossible to regularly assess forest conditions with optical data. 3.1 Location and Climate The location of the study area is indicated in Figure 5. The cut-blocks and riparian areas examined here are situated within a triangular area with the corners located at Port McNeill, Port Hardy, and Port Alice. The study area is located within the Coastal Western Hemlock biogeoclimatic zone. This zone is dominated by western hemlock (Tsuga heterophylla (Raf.) Sarg.), western red cedar {Thuja plicata Donn), amabilis fir (Abies amabilis (Dougl.) Forbes) and Sitka spruce (Picea sitchensis (Bong.) Carr) (Klinka etal., 1991; Krajina, 1965). The climate is typically wet with relatively mild winters and somewhat drier but cool summers. The persistent cloud cover and attendant rainfall is a result of the westerly warm Pacific air flows being forced up over mountainous terrain on the west side of Vancouver Island. 23 The surficial geology consists of ground moraines, drumlinized till plains, till over bedrock, and glacio-marine and lacustrine deposits. Both the frequent precipitation and the surficial geology are considered to be responsible for cut-blocks being lighter in tone in the winter months when they are soggy and as dark or darker than the surrounding forest when dry in the summer (Murtha, 2000b). Figure 5: Map of study area showing footprint of the 4 RADARSAT images. 3.2 Data Preparation The four R A D A R S A T datasets made available for this research through A D R O 384 are calibrated Level 1 Extra Fine Resolution (SGX) products generated by the Canadian Data Processing Facility (CDPF) of the Canadian Space Agency. The 3.125 metre pixel spacing of the S G X product was kept throughout all of the data processing steps. Such a small spacing fully utilizes the spatial resolution capabilities of the S A R instrument by 24 retaining all of its input image information. This is achieved at the expense of a significantly larger dataset size than other R A D A R S A T products (RADARSAT International, 2000). All of the datasets used in this research were processed and analysed using the EASI /PACE image processing software developed by PCI Geomatics. Single pixel registration of the datasets was done using ERMapper. When the C D P F generates a 16-bit product for distribution, it applies an output scaling gain and offset to the processed data in order to more fully utilize the dynamic range of the data. The degree of scaling and offsetting varies from one image to another, thereby making it difficult to ensure consistency when performing multi-temporal analysis of images acquired under different weather conditions. Removing the output scaling gain and offset reverts the data back to its original radar brightness values that were used as the base data in this research. These reverted data accurately represent the original, signal amplitudes received at the satellite antenna (RADARSAT International, 2000). To convert the digital number data in the S G X product back to its original radar brightness values, the following formula was used in the PCI Xpace program, S A R B E T A (PCI Geomatics 2000): {3°J=[(DNJ2+A3)/A2J] where DNj is the digital number of the j t h pixel in a range line /3°j is the corresponding radar brightness value in power format A2j is the scaling gain value for the j t h pixel A3 is the fixed offset value Radar brightness data can be expressed in three formats: power (intensity), decibel, and amplitude. Of the three, amplitude format data are the most appropriate for 25 arithmetic operations. These data are characterized by a fairly wide, relaxed Rayleigh distribution curve with all values above zero. Power format data are also larger than zero; however, the data are confined to a much narrower exponential distribution that gives the corresponding image a fairly dark appearance. Decibel format data are on a logarithmic scale and consequently the distribution can include both positive and negative values. Computing certain statistics such as the mean will yield an incorrect result and therefore these data are not used for change detection unless the data are ratioed (PCI Geomatics, 2000; Oliver and Quegan, 1998). Much of the research reported in peer reviewed journals utilizes radar data that are expressed in terms of the backscatter coefficient, a 0 . These values usually account for the angle of illumination only and assume that the terrain is at the earth ellipsoid (PCI Geomatics, 2000). However, the correct approach is to take account of both the illumination angle and the local incidence angle as governed by the slope of the terrain (Raney et al., 1994). The illumination incidence angle is available from the SAR leader file that accompanies the R A D A R S A T data; however, the local incidence angle can only be derived from a digital elevation model (DEM). Although a D E M was available for the entire study area, the cut-blocks and associated riparian zones that are discussed here are all situated on relatively flat ground. In addition to this, the orbital track and the direction the satellite was travelling were identical during acquisition of the four datasets utilized in this research. All of these factors suggested that there was no need to transform the radar brightness data further to backscatter coefficient values. In fact, inconsistencies between datasets could have arisen had the data been corrected for slope and registered to local forest cover maps. Relevant information concerning the data used in this research is summarized in Table 1. All of the data were acquired in Fine 2 mode on ascending orbits between 0215 and 0230 Coordinated Universal Time (UTC). The local time of acquisition is 8 hours earlier (6:15 to 6:30 pm on the day before) (Murtha, 2000b). 26 Date of Acquisition Weather during Acquisition Mean Brightness of Forested Area December 4, 1996 Heavy rain, strong winds 0.328386 August 25, 1997 Light rain 0.342076 November 24, 1998 Rain 0.405267 March 24, 1999 Showers 0.322444 Table 1: Date and weather at time of acquisition of RADARSAT data used in this research. Although acquired at different times of the year, the weather conditions during data acquisition were similar. Of particular importance is the moisture content of materials that interact with the incident radar wave, and their density and position relative to other scatterers nearby. Noted before was the observation that soggy cut-blocks are typically brighter than the surrounding forest. At other times, much drier cut-blocks are either indistinguishable or even darker than the surrounding forest (Murtha, 2000b). The brilliance of soggy cut-blocks likely results from the incident energy being bounced between pieces of wet slash and off very wet humus and soil, and again off sodden slash, vegetation, tree stumps, rocks, and other debris. Similar responses are seen when the incident energy double-bounces off the ground and nearby tree trunks along the edge of riparian zones and road rights-of-way facing the satellite. Because the forested areas are much darker in comparison, the images that should be used for this particular kind of change detection are those in which the scene is quite wet during data acquisition (i.e., the cut-blocks are bright and the riparian zones and forested areas are relatively dark). The four datasets listed in Table 1 appear to meet this requirement. It should be noted that the mean radar brightness (and standard deviation) of a well defined even-textured forest area northeast of the Port Hardy diesel generating station was quite a bit higher during the 1998 acquisition than at other times (Table 1). This suggests that another dataset with a mean radar brightness closer to the other three 27 should have been substituted for the 1998 dataset. Unfortunately, a suitable replacement was not available. 3.3 Registering the Datasets Once the digital number data in the delivered products were converted to radar brightness values, the 1997, 1998, and 1999 datasets were registered against the 1996 dataset. It has been determined that the data processing at the C D P F is so geometrically consistent that for the area of the study all of the datasets could be registered without orthorectification using just one pixel (Peter Murtha, personal communication).5 A target common to all of the datasets is a steel tank at the Port Hardy diesel generating station. Enlargements of the target were printed on transparent plastic sheets and then fitted on a light table. The UTM coordinates of the registration pixel in the 1996 dataset were then assigned to the matching pixels in the other three datasets. Once registration had been completed, the datasets were subsetted into several smaller ones ranging in size from 256 pixels by 256 lines to 1024 pixels by 1024 lines. Cut-block and riparian zone boundaries were either sketched by hand or scanned and fitted to the imagery rather than fitting the imagery to existing mapping. This approach ensured that pixels properly overlaid one another; had orthorectification been done, there was a possibility that the process could have distorted the shape of the images, leading to problems with image differencing. Once the datasets had been subsetted, one of the subsets was subjected to available change detection methods. This established some results that were used as a benchmark to compare the results obtained from the application of the Getis statistic to change detection. 5 Professor Emeritus, Faculty of Forestry, University of British Columbia 28 Chapter 4 Evaluating Change Detection Approaches 4.1 Comparing Change Detection Approaches As indicated earlier, the main objective of this thesis is to determine how effective the Getis statistic is in detecting windthrow damage in narrow riparian zones. To do so requires a benchmark against which the results can be evaluated. Such a benchmark can be established by applying existing change detection methods to the source data. The subject of the benchmarking described in this section is cut-block S501 that makes its first appearance on the 1998 R A D A R S A T image. The cut-block is situated just north and east of the Port Hardy diesel power plant (Figure 6). Its boundaries above and below road R400 were transferred from a logging completion map produced in July 1998. The east bank of the Keogh River is clearly visible northeast of the cut-block. Although the two access roads in the cut-block north of R400 appear in the 1995 aerial photograph, harvesting did not commence until after the August 25, 1997 R A D A R S A T image was acquired. Harvesting was completed in late June 1998, well before the second R A D A R S A T image was acquired on November 24 t h of that year. 29 Figure 6: 1995 aerial photograph showing the approximate boundary of the cut-block as taken from a logging completion map. Portions of the 1997 and 1998 images that were used to evaluate the change detection methods are shown in Figure 7 and outlined in Figure 6. Both of the images in Figure 7 have a dimension of 256 pixels by 256 lines. The cursor (white cross) is resting on road R400 in the 1997 image; it and the Rupert Main Line meet at the lower right edge in both images. Although there are noticeable variations in radar brightness in some areas such as those circled in white, the cut-block is considered to be the only significant structural change to have occurred in the landscape between the times the two images were acquired. 30 Figure 7:1998 (left) and 1997 images of the area in and around cut-block S501. The 1997 image on the right is appreciably darker in the area where the cut-block is in the 1998 image. However, the 1998 image has greater contrast with the east bank of the Keogh River showing up as a ribbon of bright pixels. There is also a distinct shadow along the west boundary of the cut-block in the 1998 image. Although both images were acquired under similar weather conditions, there is considerable variability in brightness. Figure 8 shows that the radar brightness values, ranging from just above 0 to a little over 2, are Rayleigh distributed (Oliver and Quegan, 1998). The coincidence of the two histograms is remarkable, given that the 1998 image appears to contain a greater proportion of bright pixels. 31 4.2 Change Detection Using Established Techniques Detecting changes in the landscape with R A D A R S A T images acquired on separate dates can be done with a multi-temporal composite of three images, by image differencing or by ratioing two images. 4.2.1 Multi-temporal Composites A simple method used for detecting change involves the construction of a multi-temporal composite of three images acquired on different dates and the application of colour theory to interpreting colour and brightness at each pixel location. The composite image in Figure 9 was constructed by assigning the 1998 data to the red colour gun, the 1997 data to the green colour gun and the 1996 data to the blue colour gun. Areas of high backscatter in 1998, 1997, and 1996 are indicated by the bright red, bright green and bright blue pixels respectively. Areas bright in both 1998 and 1997 are yellow while areas bright in both 1997 and 1996 are cyan. Areas indicative of high backscatter in 1996 and 1998 are magenta while areas bright at all three times are white. 32 Figure 9: Multi-temporal composite image showing areas of high backscatter in 1998 (red), 1997 (green), and 1996 (blue). When interpreting multi-temporal composites of forested areas, one must be aware of the scattering mechanisms at work in the landscape being imaged. In these areas, C-band radar is characterized by canopy volume scattering, whereby the coherent electromagnetic energy emitted by the antenna on the satellite penetrates the canopy of the vegetation only a short distance before it is scattered in all directions. The proportion of the incident energy reflected back to the antenna from a particular target can vary considerably from one imaging date to the next. Some of these variations may be due to slight changes in the orbital track and altitude of the satellite and the incidence angle at which the data are acquired. Other variations may be introduced by a shift in the position of scatterers and by changes in the dielectric constant of the materials interacting with the electromagnetic energy (Villasenor et al., 1993). A combination of these factors may be responsible for the highly variable and somewhat inconsistent pattern of bright pixels in the composite shown in Figure 9. A distinct pattern of bright pixels can be associated with structural change in the scene. For example, red pixels in the composite identify targets that were bright in 1998 but not in 1996 or 1997. The majority of these pixels are associated with the highlighted cut-block that first appears in the 1998 image. When a stand is harvested, the usual canopy volume scattering of the incident energy is replaced by a double bounce of this energy 33 off pieces of slash, the ground and stumps. This double bounce effect increases the strength of the signal return. In the case of the 1998 image, the strength of the return was further enhanced by an increase in the dielectric constant of the surface materials and soil resulting from heavy rainfall before and during image acquisition (Murtha, 2000c; Villasenor et al., 1993). Despite the obvious correspondence of the majority of red pixels in the composite to the new cut-block on the 1998 image, it can be argued that it is not always possible to associate increases in backscatter with structural change. Aside from the many bright red pixels that are inside the cut-block area, there are also a significant number outside of that area. Bright green and bright blue pixels are also numerous and widespread. However, very few of the three colours of pixels overlap as evidenced by the low number of white pixels in the composite. Undoubtedly, some of these bright pixels are indicative of structural changes in the landscape. However, there are simply too many of them for all to be associated with such change. Consequently, any indication of structural change in the landscape, such as windthrow damage in riparian zones, would always have to be verified in the field or by some other data source. Without such verification, the technique can only be said to work consistently where extensive change has occurred as in the case of the cut-block discussed here. 4.2.2 Image Differencing When dealing with multi-spectral data collected by a satellite such as LANDSAT, the location and extent of change can be determined by simply subtracting a recent image from one acquired earlier. The two images must first be carefully registered to use this technique effectively. To accommodate a certain degree of variability in the data acquired at two widely separated times, a threshold is applied to the difference image in order to identify areas where significant changes in radiometric brightness have occurred. The distribution of this difference image is usually nearly normal and consequently the standard deviation of the difference data can be used to set the threshold level (Eastman, 2001). 34 Although Rignot and van Zyl (1993) determined that the difference method is not suitable for application to radar data (see Section 2.1.2), the method was applied to the radar data to examine how it behaves in relatively flat terrain. A difference image derived from the 1997 and 1998 intensity format data is shown in Figure 10. the corresponding histogram on the right. A normal curve has been overlaid on the histogram to show that the data in the image on the left is not normally distributed. The difference image shows a distinct rise in the strength of the signal returned to the sensor in the area where the cut-block appears for the first time in the 1998 image. It also shows a number of clusters of bright pixels outside the cut-block area. At first glance these appear to be located at random. However, there is in fact a distinct pattern of bright pixels along the east bank of the Keogh River. A frequency distribution of the difference values is also shown in Figure 10. The frequency of each discrete difference value has been converted to a per cent based on the total number of pixels in the image. This facilitates the overlay of a normal curve on top of the frequency distribution. Since the frequency distribution of the image in Figure 10 is not normal, the standard deviation of the difference values cannot be used as a threshold value for identifying where significant change has occurred. In any case, it appears that any thresholding applied to this image may confuse logging-induced 35 increases in backscatter in the cut-block with increases in backscatter in the undisturbed area outside the cut-block. 4 . 2 . 3 Image Ratioing As an alternative to differencing, the ratio of two images could be taken and appropriate threshold levels applied to reveal significant changes on the ratio image. Unlike image differencing, image ratioing is unaffected by calibration errors (as long as they are the same for both images) and by slope-induced effects on backscattering (providing that both images have the same imaging geometry) (Rignot and van Zyl, 1993; Oliver and Quegan, 1998). The images used here meet these requirements. Although a simple ratio image is possible, the scale of relative change is not symmetric about 1 (the no change value) because ratio image data exhibit a steep negative exponential distribution. The solution is to perform a logarithmic transformation of the ratio image thereby converting the data into dB format. For amplitude format data, the appropriate formula is f = 20 logio(r) where r is the ratio between two pixels (Rignot and van Zyl, 1993; Dekker, 1998). The logarithmically scaled ratio image derived from the 1997 and 1998 amplitude images appears in Figure 11. Figure 11: Logarithmically scaled ratio image of the 1997 and 1998 amplitude data. 36 The shape of the cut-block below the road can be inferred by the pattern of bright pixels along the road and the radar shadow along the west boundary. Bright areas above the road are relatively small and only marginally brighter than other pixels in the image. Thresholding this image could either lead to the identification of changes that are not there or to some changes being missed altogether (Dekker, 1998). It may be possible to overcome this problem by applying an adaptive speckle filter to the ratio image before thresholding it. The resulting threshold bitmaps can then be overlaid on the image shown in Figure 11. This approach is the basis of the C H D E T change detection algorithm in PCI Geomatic's P A C E Radar Analysis Module. The ratio image can either be filtered by an averaging filter or an extended Kuan Minimum Mean Square Error filter (MMSE) adapted for detection of structures (Dekker, 1998). Figure 12 shows the 1998/1997 logarithmically scaled ratio image filtered with the Kuan M M S E filter using a 7 x 7 kernel. Figure 12: Logarithmically scaled ratio image filtered with a 7 x 7 Kuan MMSE filter. The appearance of Figure 12 is unlike other speckle filtered images commonly seen in the literature. However, it resembles a ratio image filtered in the same way in Dekker's (1998) paper and therefore it was concluded that the filtering performed by the PCI Geomatics software was done correctly. 37 To set the threshold variables in the C H D E T algorithm, the variance and standard deviation of the logarithmically scaled image must be computed. Unlike a non-scaled ratio image, the variance of a logarithmically scaled ratio image is only dependent on the number of looks (Dekker, 1998). The formula for the variance is (Dekker, 1998): where L is the number of looks and £(n, m) is the Rieman zeta-function (Hoekman, 1991). For a ratio image generated from one-look images, £(2, 1) is simply — 6 (Hoekman, 1991). Consequently, the standard deviation of the logarithmically scaled ratio image derived from the 1997 and 1998 one-look data is 7.877 dB. The positive threshold bitmaps corresponding to this threshold level applied to the filtered image have been overlaid on the logarithmically scaled ratio image in Figure 13. Only those pixels that have at least 3 neighbours with a value higher than the threshold level are included in the bitmap coverage. As expected, all of the bright pixels located within the cut-block area are over the threshold level. This could lead to a conclusion that the algorithm works relatively well with data processed in this way. However, what is troubling is the considerable number of pixels beyond the cut-block area that are equally bright, indicating significant increases in backscatter in areas where structural change is 200 C(2,L) w i t h a 2 , £ ) = ^ - | V 6 k=\ In2 (10) not known to have occurred. 38 Figure 13: Bitmaps corresponding to +1 standard deviation overlaid on the logarithmically scaled ratio image. The significance of such errors associated with correctly identifying structural changes in the scene can be characterized by the probability of false alarm (PFA). The conditional probability distribution of a logarithmically scaled ratio image closely approximates that of the normal distribution, even for a ratio image derived from 1-look images (Dekker, 1998; Dekker, personal communication)6. A frequency distribution graph of the image in Figure 13 with a normal curve superimposed on it suggests that this is indeed the case (Figure 14). Appropriate values corresponding to different levels of the probability of false alarm can now be extracted from a normal table to calculate a threshold level corresponding to a specific PFA. 6 Research Scientist, T N O Physics and Electronics Laboratory, The Hague, Netherlands. 39 O N i o o i - ' - n i p N O i LO TT CO CN ^ " T - C M C O - * I I I I I Ratio (dB) Figure 14: Frequency distribution of the logarithmically scaled ratio image with the normal curve superimposed on it. The threshold level of 7.877 dB corresponds to a P F A of approximately 16%. For a smaller P F A of 5%, the thresholds would be 1.645a or ±12.958 dB, the critical value of 1.645 being obtained from a normal distribution table. Figure 15 (left) shows the bitmaps overlaid on the ratio image that correspond to a threshold of 12.958 dB. As was the case with thresholding at the lower level of 7.877 dB, the size and number of bits was thinned by setting the minimum number of neighbours to 3. On the right is the same bitmap coverage generated with the minimum number of neighbours set to 0. 40 Figure 15: Positive threshold bitmaps corresponding to a PFA of 5% overlaid on the logarithmically scaled ratio image. The bitmaps on the right are unthinned while those on the left have been thinned by setting the minimum number of neighbours to 3 in the CHDET algorithm. The bitmaps in Figure 15 are considerably less extensive than in Figure 13 and are now entirely confined to the cut-block area where structural change is known to have occurred. The apparent reduction of false alarms is a desirable outcome; however, the bitmaps now only indicate where highly significant change has occurred. Even at a PFA of 16%, the CHDET algorithm does not adequately reveal the known extent of change associated with harvesting the cut-block. At best, the algorithm reveals the location of the most significant change in the scene. The nature and extent of this change would have to be inferred from other data sources or by a field reconnaissance. Regardless of the results obtained here, the CHDET algorithm has proven to be successful for detecting change in urban areas (Dekker, 1998). Improved methods for detecting change in more challenging areas such as forested scenes are still being actively pursued (Couturier et al; 1999). 4.2.4 Summary Several existing change detection methods were used to detect the presence of significant change in a forested area due to timber harvesting in one cut-block. The three-colour composite revealed the highly inconsistent nature of the three data sets under examination (Figure 9). One indication of how different the three images are is 41 the very low proportion of pixels that are bright on all three images. Coincident pixels that are bright in all three images would appear white on the composite. Care had been taken to co-register the images using a strong point target at the nearby Port Hardy power plant. Therefore, a higher degree of coincidence of bright pixels would be expected, especially along the east bank of the Keogh River where the illumination from the left double bounces off the exposed water, ground surface and neighbouring tree trunks. That it does not makes it more difficult to sort out what is and what is not structural change in the multi-temporal composite. The apparent ambiguity in the data makes it difficult to successfully map where structural change has occurred. There was moderate success with thresholding the filtered logarithmically scaled ratio image, however, the probability of false alarm had to be lowered to 5% in order for significant change to be limited to the cut-block area where the only structural change is known to have occurred. The CHDET algorithm is therefore potentially useful for indicating where significant change has occurred; however, a secondary data source or a field visit would be required to adequately map the extent of the change that is not revealed by the algorithm. It should be noted that the results discussed above used a rather small 256 pixel by 256 line dataset. The same procedures were also applied to a 1024 pixel by 1024 line dataset of the same area. The larger dataset revealed that cut-block S501 constituted the only area where extensive physical changes in the scene appeared to have occurred. Consequently the cut-block stands out rather well on multi-temporal composites. The C H D E T algorithm produced an identical result since the threshold level is the same for both sizes of datasets. As was the case with the smaller dataset, the most readily apparent CHDET threshold bitmaps were confined to the S501 cut-block area. Bitmaps outside of this area were practically impossible to see due to their very small size. The CHDET algorithm was also applied to three larger cut-blocks located near the head of Rupert Inlet, some 10 km west of cut-block S501. Logging in cut-block M5307 42 commenced soon after the 1996 image was acquired and was completed by the time the 1998 image was taken. Logging in cut-blocks C99 and 594 first appears on the 1998 image and was completed by the time the 1999 image was acquired. In all cases, image pairs were selected for optimal detection of structural change in the cut-blocks. Application of the C H D E T algorithm to logarithmically scaled ratio images of these three cut-blocks produced results similar to those generated for cut-block S501 (Table A-1 , Appendix A). In all cases, the probability of false alarm was set to 5% and the minimum neighbour parameter was set to 3. The bitmaps indicated that the majority of pixels exceeding the threshold were confined mostly to areas disturbed by logging. However, their numbers were simply too small to accurately map the extent of the logging itself. Judging by the number of threshold bitmaps outside of the cut-blocks, one would still be unsure of which correctly identified structural change and which identified a false alarm. Although the number of probable false alarms is quite low, the appearance of so many bitmaps outside of the disturbed areas suggests that a more effective means of detecting structural change is required. The methods examined in this section provide a basis for comparing their performance against the results obtained in the next section where the Getis statistic is used as a change detection tool. 4.3 Change Detection using the Getis Statistic In Section 4.2, it was shown that existing methods for detecting change on radar images were not entirely successful when applied to the single-look, high spatial resolution R A D A R S A T images of the study area. This section will evaluate the performance of the standardized form of the Getis statistic as an alternative for detecting both large- and small-scale changes on these images. An EASI (Engineering Analysis and Scientific Interface) script (Getis EASI) was written in PCI Geomatics' Geomatica image processing software to generate Getis statistics 43 images using the data from one radar brightness image. The EASI scripting language provides a convenient method for passing the centre of an odd-size symmetric sampling window or kernel over every pixel in a radar brightness image and calculating a Getis statistic or value at that location. The equation used in the Getis EASI script to compute the Getis value is given in Section 2.2. In its application in that script, all of the pixels in the kernel are used in the computation of the Getis value and each pixel in the kernel is assigned a weight of one. Pixels at the edge of the image are replicated outwards as far as necessary to ensure the kernel is filled with data. Five kernels (3x3, 5x5, 7x7, 9x9, and 11x11) are used to generate five Getis value images from one radar brightness image. Another EASI script (MaxGetis EASI) subsequently examined the five Getis values at each pixel location and chose the highest Getis value that occurred in the smallest spatial domain. The distance value associated with the kernel that generated the MaxGetis value was written to a MaxGetis Distance image. The MaxGetis image is essentially a spatial autocorrelation map while the MaxGetis Distance image can be used to assess the extent of spatial autocorrelation associated with each pixel in the source data as well as to determine the degree of homogeneity of these data. The next few sections will examine the utility of the Getis value, MaxGetis and MaxGetis Distance images for determining the extent of structural change within the same areas examined in the previous section. 4.3.1 Getis Value Images Figures Figure 16, Figure 17, Figure 18, Figure 19, and Figure 20 show the Getis value images generated by the Getis EASI program using the data from the 1997 and 1998 images first shown in Figure 7. The magnitude of the Getis values in these images closely corresponds with the magnitude of radar brightness values in the source data. Getis values tend to be highly positive where radar brightness values within the kernel are both very high and homogeneous. Conversely, the Getis values are highly negative where radar brightness values within the kernel are both very low and homogeneous. 44 The Getis statistic is unusual in that high spatial autocorrelation is indicated only by those values that are at both ends of the Getis value distribution. Consequently a potentially serious weakness of the Getis statistic is the inability to differentiate Getis values computed from homogeneous radar brightness values of intermediate value from those that are computed from a heterogeneous mix of values within the kernel. However, this is of no concern here since the indicator of change is a brightening of clusters of radar brightness values in areas where trees have either been windthrown or removed by timber harvesting. These events replace the relatively weak volume scattering returns typical of closed canopy forest and no return in radar shadow areas with stronger radar returns generated by an array of newly created targets exposed to the radar sensor including ground vegetation, stumps, large pieces of slash, tree boles and root balls. As indicated earlier, the Getis statistic is a tool for evaluating the strength of spatial autocorrelation over localized areas. Since the pixel spacing of the R A D A R S A T data used in this study is 3.125 metres, a 3x3 kernel computes the Getis value for the centre pixel using the B° data from 9 pixels that cover an area that is just over 9 metres on a side. At the other extreme, the 11x11 kernel computes the Getis value for the centre pixel using a total of 121 pixels in a kernel that measures roughly 34 metres on the side. The magnitude of the Getis value for a particular pixel relative to others around it indicates both the magnitude of the source data and the strength of the spatial autocorrelation associated with that pixel. The calculation of the Getis value using progressively larger kernels is responsible for the increased fuzziness of features in the output images shown below. The fine detail of bright features that is obvious in the source data, such as the one circled on the 1998 images on the left, is lost as the kernel grows larger. It is even degraded to a small extent in the image generated by the 3x3 kernel. The visibility of very small clusters of bright pixels against a dark background slowly diminishes until it is lost altogether once the kernel size reaches 9x9. This stands to reason as the value of the Getis statistic 45 calculated for the centre pixel in the kernel is influenced by the radar brightness values of a progressively larger number of neighbouring pixels as the kernel increases in size. Consequently, the series of Getis value images shown below appear to be the source image subjected to ever larger averaging filters. Despite the blurring, however, it is possible to perceive the extent and location of the cut-block on all of the images on the left. Figure 17: Getis value images generated by a 5x5 kernel using the 1998 (left) and 1997 (right) radar brightness data. 46 Figure 18: Getis value images generated by a 7x7 kernel using the 1998 (left) and 1997 (right) radar brightness data. Figure 19: Getis value images generated by a 9x9 kernel using the 1998 (left) and 1997 (right) radar brightness data. 47 Figure 20: Getis value images generated by an 11x11 kernel using the 1998 (left) and 1997 (right) radar brightness data. Because of the blurring effect, kernel size should play a significant role in determining the scale or size of features that will be retained after the Getis EASI program processes the source data. If the Getis value images themselves are used for detecting small scale changes by differencing two images, then a 3x3 or 5x5 kernel would appear be the most appropriate because small scale features are better preserved than in the images generated with the larger kernels. Larger kernels, particularly the 11x11, may be more desirable for revealing larger scale changes in the landscape such as the appearance of whole cut-blocks in a set of multi-temporal R A D A R S A T images. Another consideration may be to use larger kernels to overcome the undesirable influence of speckle noise on the detection of small and medium scale changes in the landscape. Of the Getis value images shown above, those generated by the 3x3 kernel (Figure 16) most closely resemble the source images shown in Figure 7 and would be the images of choice for detecting change by visual means. However, appearances can be deceiving. Although the histograms in Figure 21 have roughly the same shape as those of the source data (Figure 8), the range of values is wider and drops further below zero. That the corresponding Getis value images exhibit the same appearance as the source images is a function of how the image processing software displays them. In both cases, the full range of data has been rescaled to a range of zero to a number that 48 corresponds with the bit depth of the graphics card driving the video screen. This rescaling, combined with linear enhancement to improve overall brightness and contrast, can lead to erroneous comparisons being made between images whose source data have a somewhat different range of values. Consequently, the information content of the Getis value images (and the MaxGetis images discussed later) will be assessed quantitatively. First, the data will be examined through profiling. Change will then be assessed by differencing Getis value data generated by the same sized kernel. This will be followed by thresholding the difference images. 350 300 250 3 200 o o _ 150 x bl 100 50 0 -4.4 -2.3 -0.3 1.8 3.8 5.9 7.9 10.0 12.9 Getis Values 1998 3x3 1997 3x3 Figure 21: Histograms of the 1997 and 1998 Getis value images generated by a 3x3 kernel. 4.3.1.1 Profiling Profiling facilitates the comparison of pixel values from two or more images in the same PCI database. A line is drawn through the area of interest and the images from which the data are to be obtained are selected. A graph is then generated to show the pixel values under that line. 49 To illustrate how this works, a yellow line has been drawn in Figure 16 through a portion of the cut-block west of road R400. Figure 22 is the resulting graph that compares the Getis values computed from the 1997 data by the 3x3 and 11x11 kernels with the corresponding source data, the 1997 radar brightness data. Figure 23 shows the corresponding profiles for the 1998 radar brightness data. Sample Profile of Selected Channelfs) Distance (m) Figure 22: Profile of pixel values on the 1997 radar brightness image (black), 3x3 Getis value image (red), and 11x11 Getis value image (blue). These profiles show the extent to which the source data have been altered in their conversion to Getis values. The black lines show the relatively low dynamic range of the radar brightness values while the red and blue lines show the highly dynamic Getis values that were generated by the 3x3 and 11x11 kernels respectively. Although the profiles indicate the values of single pixels directly under the line, adjacent pixels that are within the confines of the kernel contribute to the asynchronous fluctuations in the Getis values along the profile. Clusters of bright and dark pixels within the 3x3 kernel are responsible for the sharp peaks and valleys seen in the Getis values generated by that kernel. The differential between the peaks and the valleys in the 11x11 data is not as great because the clusters of bright and dark pixels that are highly influential in a 3x3 kernel are generally smaller than the 11x11 kernel itself. The larger kernel size blurs the 50 image considerably and consequently bright and dark areas become larger and more homogeneous in appearance thus resulting in the flattened peaks and valleys. In Figure 22, the 11x11 profile is generally below that of the 3x3 profile reflecting the source image's generally dark appearance in the cut-block area. In Figure 23, the opposite is the case since the cut-block area is noticeably brighter overall after logging. As a consequence, the Getis calculations in the 11x11 kernel are influenced by more extensive groupings of bright pixels in the disturbed area. r e V V a I u e s Sample Profile of Selected Channelfs) " « A n f\ • A / M A / W M Legend ii r U T A V \ MM 1 Channel 1 i l u W "i W/uiLA 1/rtW \\\\ I Channel 4 | Channel 12 100 200 300 Distance [m] 400 500 Figure 23: Profile of grey (pixel) values on 1998 radar brightness image (black), 3x3 Getis value image (red), and 11x11 Getis value image (blue). The profiles in Figure 22 and Figure 23 appear to suggest that the magnitude of spatial autocorrelation, as exhibited by the highest absolute Getis value, is for the most part the greatest in the 11x11 data. However, the MaxGetis data (discussed in Section 4.3.1.3) show that this is generally not the case unless the source data are from relatively homogeneous areas of the image such as water bodies and radar shadow areas. In heterogeneous areas of the image, the first peak in the absolute Getis values normally occurs in a kernel much smaller than 11x11. 51 4.3.1.2 Thresholding Thresholding has been discussed earlier in the context of differencing and ratioing radar brightness data. There it was used with limited success to differentiate areas of significant change on the 1998-97 difference image of the same cut-block and surrounding area. The Getis value images themselves suggest that differencing followed by thresholding may yield better results than differencing radar brightness images because the former exhibit greater overall contrast. This contrast stretch is reflected in the wider range of grey level values in the histograms of the Getis value images (compare Figure 8 and Figure 21). Despite the stretch, the 3x3 Getis value images are remarkably similar to the source images, right down to the brightness inconsistencies noted between the 1997 and 1998 source images (compare Figure 7 with Figure 16). Difference images were created from the 1997 and 1998 images processed by the 3x3 and 11x11 kernels (Figure 24). In both instances, the location and extent of the cut-block can be deduced after the cut-block boundary has been added to the image. In fact, the cut-block is more easily discerned in the image produced from the 3x3 Getis value data than it is in the image created by simply differencing the radar brightness data (Figure 10). It is also easily discerned as a bright, large scale feature on a much higher contrast image generated from the 11x11 Getis value data (Figure 24, right). 52 Figure 24: Getis value difference images generated from the 3x3 data (left) and the 11x11 data (right). To differentiate the extent of significant change, a threshold can be applied to difference data. In practice, thresholding identifies a population of pixels on the difference image that have a value higher than a certain critical value on the x-axis of the image's histogram. Pixels in this population reflect the greatest spread in grey value between a bright pixel on one image and a dark pixel in the same location on the second image. Since Getis values extend below zero, significant change in a Getis difference image usually results from the subtraction of a negative Getis value, corresponding to a dark pixel, from a positive Getis value corresponding to a bright pixel. An examination of the frequency distributions in Figure 25 suggests that standard deviation-based thresholding may not be appropriate as both distributions deviate somewhat from the normal. After some trials, it was determined that a 95% threshold value produced consistent results upon which comparisons could be made. This value is the point along the x-axis of the difference distribution at which the area under the curve amounts to 95% of the total. Pixels with values higher than this critical value would then be considered to be indicative of 'significant' change (Jensen, 1996). It is important to note that the 95% threshold level implies that 5% of the total area of the image would then be flagged as indicative of significant change. The distribution of the thresholded image relative to where change is known to have occurred greatly 53 influences the degree to which false alarms are generated by this method of thresholding. -19.4-12.0 -7.6 -3.2 1.2 5.6 10.0 14.4 Getis Difference Value Figure 25: Frequency distributions of the 3x3 (left) and 11x11 Getis difference images with a normal curve overlaid on them in black. The results of applying a 95% threshold to the 3x3 and 11x11 Getis value difference images are shown in Figure 26. Thresholding the 11x11 difference image at this level produced the results that were expected as the substantially larger bitmaps are almost entirely confined to the area where structural changes in the landscape are known to have occurred. The size of the larger bitmaps goes up only slightly when the threshold is lowered to 90%. However, isolated bitmaps begin to appear in areas beyond the cut-block boundary where structural change is not known to have occurred. Further fine tuning in this instance may lead to an ideal threshold value between 90 and 95%. 54 Figure 26: 3x3 (top) and 11x11 (bottom) Getis value difference images with 95% threshold bitmaps overlaid on them (left). Images on the right show the bitmap coverages after 5x5 mode filtering. Initially, the result of thresholding the 3x3 difference image was not as encouraging. First of all, the combined area of the bitmaps inside the cut-block area was considerably less than 5% of the total area of the image although a sufficient number are present to visually deduce the location and extent of the cut-block. The balance of the bitmap area, distributed among a high number of very small bitmaps, lay outside of the cut-block in areas where structural change is not known to have occurred. Consequently, many of these bitmaps could be the result of a false alarm. Later trials with speckle filtered data (Section 4.3.3) suggested that speckle noise may be partly responsible for generating many of these apparent false alarms, not only outside of the cut-block area 55 but inside as well. Consequently, it was decided that some means of removing the smallest bitmaps was needed to reduce the number of these apparent false alarms. The C H D E T algorithm explored in section 4.2.3 offers an attractive option to retain only those bitmaps that have a specified minimum number of neighbours within a kernel. After thresholding has been completed in that algorithm, the bitmap coverage is processed iteratively until the number of neighbours remains constant around each bit in the bitmap coverage. It thereby reduces the number of very small bitmaps in images that have been generated by oversampling the sensor data (PCI Geomatics, 2000). Unfortunately, this minimum number of neighbours feature is not available outside of the CHDET algorithm program. However, a mode filter applied to the bitmap coverage produced almost the same result. While the minimum neighbour algorithm preserves the original geometry of the larger bitmaps, the mode filter alters the appearance of these bitmaps somewhat by smoothing their edges and filling internal holes. These alterations can be kept to a minimum by mode filtering the bitmaps no more than once. The images on the right in Figure 26 show the effect of mode filtering the 95% threshold bitmaps with a 5x5 kernel. The images clearly emphasize where significant change has occurred in areas bigger than about 245 m 2 . Even single bitmap pixels remaining after one round of mode filtering indicate a significant area where change has occurred. In this way, mode filtering clears up some of the ambiguities associated with thresholding Getis value difference images derived from unfiltered source data. The mode filtered 95% threshold bitmaps derived from the 3x3 and 11x11 Getis value difference images appeared to indicate significant areas of change along the east bank of the Keogh River east of the cut-block. However, the nature of these changes is unknown. Thresholding the 5x5, 7x7, and 9x9 difference data produced bitmaps that were transitional both in terms of their number and size. The thresholded 7x7 image is shown in Figure 27. The location and size of the cut-block is easily discernible on the left; 56 however, there are a significant number of what might be false alarms outside of the cut-block. Their number was considerably reduced after the bitmap coverage was mode filtered with a 5x5 kernel (Figure 27, right). Figure 27: 7x7 Getis value difference image with 95% threshold bitmaps on the left. On the right are the same bitmaps mode filtered by a 5x5 kernel. As Figures Figure 26 and Figure 27 show, there is a definite trend when Getis value images generated by ever larger kernels are differenced and thresholded. As the kernel increases from 3x3 to 11x11, the proportion of the total threshold bitmap area increases in the cut-block area where the only known structural change in the scene has occurred. As this proportion changes, the number of small bitmaps drops and the size of larger bitmaps increases. This is in response to a suppression of very small bright targets whose significance is lost as they become outnumbered by an ever increasing number of darker pixels in the kernel. The opposite effect occurs where large clusters of bright pixels occur. In these areas, the outnumbered dark pixels at the centre of large kernels are converted to bright pixels in the Getis value image, thus contributing to a filling-in and "spreading" of clusters of bright pixels. Thus, as the kernel size increases, significant clusters of bright pixels are emphasized at the expense of small point targets. This suggests that the Getis value images generated by the larger kernels may not be suitable for detecting small isolated windthrow events if these are signified by small bright point targets in the source data. 57 Ratioing is another change detection technique that was explored earlier with the radar brightness data. The same procedure was attempted with Getis value images; however, this did not produce anything meaningful because the 32-bit Getis values range from negative to positive. Ratioing two negative and two positive numbers both yield a positive fraction, thereby leading to results that are confusing. 4.3.1.3 MaxGetis Images As indicated earlier, a MaxGetis EASI script was developed to select the maximum Getis value from among the five generated for each pixel in an image and to write the selected values to a MaxGetis image of the same size and shape. Although Getis values were generated for all pixels in the source images, the kernel processing algorithm used in the EASI script replicated pixel data along the edges of the images to ensure kernels were completely filled with data. For some purposes, such as image filtering, the consequent generation of artificial data from replicated data does not have a serious impact. However, it was deemed essential in this application that all Getis values considered in the selection of the MaxGetis value be derived from non-replicated data. To ensure that this would be the case, the values for the five outermost pixels in each Getis image were discarded. This action left a narrow band of grey pixels around the edge of the MaxGetis images seen in the figures below. The logic behind the selection of a MaxGetis value can be illustrated with the data in Table 2. It shows the five Getis values computed for 10 pixels under the yellow line shown in Figure 16. The table also shows the MaxGetis value derived by the MaxGetis EASI script. The MaxGetis EASI script first converts all five Getis values to an absolute value. It then conducts a series of pair-wise comparisons. If the first number in the pair is larger than the second, the Getis value generated by the kernel corresponding to the first number is written to the MaxGetis image. If the second number is bigger, then it and the next number are compared in the next pair-wise comparison. The pair-wise comparisons continue until a Getis value corresponding to either the first value in the 58 comparison or, when the comparisons have been exhausted, the 11x11 Getis value is written to the MaxGetis image. The selection of the MaxGetis value in this manner is consistent with the observation that Getis values closer to the negative and positive ends of the distribution indicate higher spatial autocorrelation than those closer to the mid-range value of zero (Wulder and Boots, 1998). To illustrate how this observation is put into practice in the MaxGetis EASI script, consider the first pixel in Table 2. The Getis value of -0.727894 generated by the 3x3 kernel represents a higher degree of spatial autocorrelation than the value generated by the 5x5 kernel, 0.431740. Although the absolute values of the data generated by the second through fifth kernel steadily increase, the value generated by the 3x3 kernel is higher than that generated by the next largest kernel. Selecting the 3x3 Getis value as the MaxGetis value therefore preserves a localized maxima in spatial autocorrelation. In the case of the second pixel, the absolute value continues to increase from the first kernel on until it reaches a maximum with the 11x11 kernel. At this point, the pair-wise comparisons are exhausted and the 11x11 Getis value is selected as the MaxGetis value. 3x3 Getis Value 5x5 Getis Value 7x7 Getis Value 9x9 Getis Value 11x11 Getis Value MaxGetis Value -0.727894 0.431740 0.913595 1.738873 j 3.636315 -0.727894 0.167225 -0.452083 -0.526995 1.106282 3.095275 3.095275 0.519464 -1.120023 -0.348780 1.221612 2.215084 -1.120023 -0.330167 -1.067662 -0.936104 0.234050 1.583690 -1.067662 0.085204 0.334647 -0.324215 -0.342095 0.458235 0.334647 1.004708 0.209331 -0.504218 -0.626444 -0.307791 1.004708 -0.291433 -1.077540 -1.980569 -1.879974 -1.361409 -1.980569 -1.922766 -2.331922 -2.868374 -2.981281 -2.240482 -2.981281 -2.645383 -3.231967 -3.735220 -4.008456 -2.554649 -4.008456 -2.277489 -3.172531 ^.173574 -3.507102 -1.850815 -4.173574 Table 2: Example showing which Getis Value is chosen as the MaxGetis Value. The MaxGetis images generated from the 1997 and 1998 Getis value images of cut-block S501 and its surrounding area are shown in Figure 28. The sharp detail that is 59 characteristic of the source images shown in Figure 7 is largely absent in both of these images. Instead, they exhibit extensive areas of very bright and very dark pixels where 24-28% of the Getis values were generated by the 11x11 kernel. The intervening areas are quite variable and blurry, reflecting the influence of Getis values generated by a diverse mix of kernels. Although more than 30% of the Getis values were generated by the 3x3 kernel which earlier showed the least propensity for producing a blurry image (see Figure 16), they were so scattered over both images that their contribution towards preserving the detail of the source images was minimal. Figure 28: MaxGetis images generated from the 1998 (left) and 1997 source data. The histograms of the two MaxGetis images are shown in Figure 29. The black vertical line rising out of the valley in the histograms indicates the position of the zero MaxGetis value on the x-axis. The double hump of these histograms is a characteristic feature of all MaxGetis images generated by the EASI script. It arises from the decision rule described earlier that is used to select the maximum Getis value from the five computed for each pixel in the source image. 60 350 300 250 | 200 o o 1 150 a. 100 50 0 ,l J I / i f \ J \_ •1998 1997 -13.2 -9.8 -6.5 -3.1 0.2 3.5 6.9 10.2 13.5 16.9 MaxGetis Value Figure 29. Histograms of the 1997 and 1998 MaxGetis images. Although not obvious, the two MaxGetis images include Getis values associated with the brightest pixels in the source data. These pixels, with a p° value greater than 0.95, are aggregated in clusters on the images in Figure 7. Approximately 1,800 pixels are associated with these clusters in each of the two images. Since the clusters are quite small, the 3x3 and 5x5 kernels generated 71% and 78% of the associated Getis values that appear on the 1998 and 1997 MaxGetis images respectively. However, only 37% and 30% of these pixels retain their comparatively high brilliance in the 1998 and 1997 MaxGetis images. The remainder lose their brilliance because they are either attached to the edges of bright clusters or are scattered throughout the source images as small, isolated clusters averaging three pixels in size. In those cases, a predominance of darker pixels in the kernel has a dampening effect on the magnitude of the Getis value. These pixels either lose their brilliance altogether or acquire an intermediate value that has the effect of blurring transitions between adjoining bright and dark clusters of pixels. In other areas, the apparent contrast of the source image is enhanced. This is a result of smoothing done by larger kernels where closely spaced bright or dark pixels 61 predominate over extensive areas. The effect of this smoothing is illustrated in Figures Figure 30 and Figure 31. Not only are the bright and dark areas more homogeneous on the MaxGetis images than in the source images on the left, they also appear to be brighter or darker in the respective areas. Figure 30:1997 p° image (left) and the MaxGetis image of the same area. Figure 31:1998 p° image (left) and the MaxGetis image of the same area. 62 Although much of the detail seen in the source images is lost during the creation of the MaxGetis image, the localized optimization of spatial autocorrelation enhances the differentiation of significant areas of both high and low backscatter. This suggests that MaxGetis images could perhaps be suitable for detecting change using the image differencing and thresholding approaches applied to the Getis value images described earlier. To determine if this might be the case, a difference image was created from the 1997 and 1998 MaxGetis images (Figure 32). This image shows the S501 cut-block with considerably less ambiguity than on the 3x3 difference image in Figure 24. Figure 32 also shows that the distribution of values in the MaxGetis difference image does not fit the normal distribution overlaid in red on the histogram. Consequently, thresholding the image to identify significant increases in backscatter must use the same procedure applied to the Getis value difference images. Figure 33 (left) illustrates the result of applying a 95% threshold to the MaxGetis difference image. It shows remarkably fewer apparent false alarms outside of the cut-block area than the thresholded 3x3 difference image (Figure 26). Mode filtering the bitmaps in Figure 33 (left) removed only the smallest of them outside of the cut-block area (Figure 33, right). Figure 32:1998-1997 MaxGetis difference image and corresponding histogram. 63 Figure 33: 95% threshold applied to the 1998-1997 MaxGetis difference image (left). On the right are the same bitmaps mode filtered with a 5x5 kernel. An analysis of the interactions between the MaxGetis values whose differences exceeded the 95% threshold level shows that almost 60% were between the 1998 MaxGetis values generated by the 11x11 kernel and the 1997 MaxGetis values generated by all five kernels. This suggested that the magnitude of the MaxGetis values and the extent of spatial autocorrelation were quite variable in the 1997 MaxGetis image in areas that exhibited high MaxGetis values and extended areas of high spatial autocorrelation in the 1998 MaxGetis image. It is important to note that once spatial autocorrelation has been optimized in the MaxGetis image, it does not play any further role during image differencing. The image differencing applied here is at the pixel level and consequently the spatial arrangement of MaxGetis values is what governs the outcome of image differencing and thresholding. The assumption made here is that differencing will detect significant increases in pixel brightness in the 1998 B° image, i.e. pixels that were bright in the 1998 image where they were previously dark in the 1997 image will make up the bulk of pixels identified by thresholding as experiencing a significant increase in the Getis value. 64 To determine if the location of the 95% bitmaps were reasonable, the mode filtered bitmaps were overlaid as polygons on multi-temporal composites of the MaxGetis and B° images in Figure 34. The 1998, 1997, and 1996 images of each were assigned to the red, green, and blue colour guns, respectively. The overlay on the B° composite indicates that generating the 1997 and 1998 MaxGetis images followed by differencing and thresholding appears to have captured the most significant increases in backscatter as signified by the clusters of bright red pixels. These clusters are more pronounced on the MaxGetis composite on the left in Figure 34. The cut-block itself is more readily discernible on the MaxGetis composite than it is on the B° composite, making the former particularly useful for visual detection of change. However, it is important to note that not all of the clusters of bright pixels visible in the B° composite were captured by the mode filtered threshold bitmaps nor were they captured by the unfiltered threshold bitmaps. It also became apparent that many of the small bitmaps removed by mode filtering and some of the larger ones that survived the filtering are not associated with clusters of bright red clusters in the B° composite. Consequently the bitmaps cannot be used to consistently identify areas where backscatter has increased by more than the threshold value. Figure 34: Multi-temporal composites of MaxGetis and radar brightness images The Getis value-based techniques described here were also applied to the three cut-blocks described in section 4.2.5. Table A-2 (Appendix A) shows that the largest 65 proportion of the cut-block captured by the unfiltered 95% threshold bitmaps was on the 11x11 Getis value difference image while the lowest proportion was on the 3x3 Getis value difference image. The proportion captured on the 11x11 Getis value difference image was the highest because a large kernel tends to brighten pixels around the edges of some clusters, thereby enlarging associated threshold bitmaps (compare Figures Figure 26, Figure 27, and Figure 33). The proportion captured by thresholding on the MaxGetis difference image is much closer to that of the 3x3 difference image; however, they are far fewer in number and much closer to the size and number of bitmaps on the 7x7 Getis value difference image. The MaxGetis difference image also appears to contain a fairly low number of apparent false alarms. Dataset size has a significant influence on the total threshold bitmap area inside disturbed areas. This is illustrated by cut-block S501 for which difference images were generated using 256 line by 256 pixel and 1024 line by 1024 pixel datasets (Table A-2). Although both sets of data were subsets of much larger B° images, the slight differences in global mean and standard deviation between a small dataset and a large dataset derived from the same parent image did have an effect on the magnitude of the Getis values. In the case of the two datasets derived from the 1998 6° image, the mean and standard deviation were somewhat lower in the larger dataset. The opposite was the case for the datasets derived from the 1997 B° image. These had a corresponding effect on the magnitude of the Getis value computed for each pixel. However, of greater influence on the total threshold bitmap area inside S501 was the threshold value itself. This value was lower in the larger dataset resulting in a much higher proportion of the cut-block being covered by the threshold bitmaps, a desirable result. However, this was also accompanied by a greater proportion of the area immediately surrounding the cut-block also being covered by the threshold bitmaps. This was not a desirable result since it merely increased the number of apparent false alarms, particularly as S501 was the only area of disturbance evident in both the small and large datasets. A s expected, the number of apparent false alarms was much lower in datasets that contained more than one recently logged cut-block. This observation, as well as the 66 increased bitmap area in S501 resulting from the use of larger datasets, demonstrates a weakness in the thresholding approach used here. Thresholding at 95% simply means that the brightest 5% of pixels in the difference dataset will be highlighted by the threshold bitmaps, regardless of whether or not any actual disturbance had been responsible for the increase in backscatter. The solution would be to raise the threshold level to reduce the number of false alarms at the expense of perhaps not detecting all of the known structural change. This was attempted with the 1024 line by 1024 pixel MaxGetis difference image containing the S501 cut-block. The result was that a threshold of 99% considerably reduced the number of apparent false alarms without seriously affecting the size of the bitmaps in the cut-block area. In fact, the bitmaps in the cut-block area were quite similar in shape and number to those generated by the 95% threshold level on the 256 line by 256 pixel MaxGetis difference image. A similar result was obtained by thresholding the 1024 line by 1024 pixel MaxGetis difference image containing the M5307 cut-block at 99%. However, the ideal threshold limit for the 1024 line by 1024 pixel dataset containing cut-block 594 appeared to be 97.5%. Since 95% thresholding followed by mode filtering seemed to be the most effective at picking up all of the significant change while at the same time reducing the number of apparent false alarms, the decision was made to conduct thresholding in riparian areas and along cut-block margins at the 95% level. An examination of the mean and standard deviation within cut-blocks found them to be always higher than in the entire difference dataset generated from either Getis value or MaxGetis difference images (Table A-2). The higher mean was due to the prominence of bright pixels within the cut-blocks whereas the higher standard deviation was attributable to the smaller number of pixels in the cut-blocks relative to the whole dataset and the wide array of very bright and very dark pixels in these features. This suggests that change could perhaps be detected by evaluating significant changes in the mean and standard deviation of the population of pixels within a feature that is tracked over time. 67 This section has demonstrated that MaxGetis difference images and multi-temporal composites are both potentially useful for identifying where significant increases in backscatter have occurred as a result of a disturbance. In addition, MaxGetis images enhance backscatter differences in multi-temporal images through local optimization of spatial autocorrelation, while at the same time reducing the prevalence of apparent false alarms. The next section evaluates the utility of the maximum Getis distance (MGD) data for detecting local changes in the extent of spatial autocorrelation and data homogeneity from one imaging date to the next. 4.3.2 Maximum Getis Distance Images The MaxGetis EASI script also creates a Maximum Getis Distance (MGD) image at the same time a MaxGetis image is produced. The MGD image indirectly indicates which kernel generated the maximum Getis value. Values in the M G D image range from 1 to 5 with each value corresponding to the radial distance over which the maximum Getis value was attained. As indicated earlier, the practical method used for computing the Getis value over a particular distance was to use a kernel. Therefore, 3x3, 5x5, 7x7, 9x9, and 11x11 kernels were used to compute the Getis value over radial distances of 1, 2, 3, 4, and 5 pixels from the centre of the kernel. The significance of MGD images is that they can be used to determine if spatial autocorrelation is localized or spatially extensive. They can also be used to determine if the extent of spatial autocorrelation has changed in a particular area over time (Wulder and Boots, 1998; Holden et al., 2000). If the MaxGetis value for a particular pixel in the source data is generated by a kernel smaller than 11x11, then the extent of spatial autocorrelation is usually confined to a local area and the region is considered to be heterogeneous. However, if it is generated by an 11x11 kernel, then the extent of spatial autocorrelation may extend much further resulting in a region that is homogeneous. When the M G D values are displayed as an image, the extent and heterogeneity of spatial autocorrelation can be visualized. 68 Figure 35 shows the M G D images generated from the 1997 and 1998 (3° images of cut-block S501. The various kernels that contributed the MaxGetis values to the 1997 and 1998 MaxGetis images are coded as yellow (3x3), magenta (5x5), blue (7x7), green (9x9), and red (11x11). Immediately noticeable on the two MGD images are the extensive areas in red. The majority of the red area corresponds to radar shadow areas in the source data while the remainder corresponds to sizable areas populated by closely spaced bright pixels. Spatial autocorrelation is spatially extensive in both of these areas, reflecting the relatively homogeneous nature of the corresponding source data. In other areas populated by a mixture of distance values, the extent of spatial autocorrelation is much more limited thus reflecting the heterogeneous nature of the source data in these areas. mi Figure 35.1998 (left) and 1997 Maximum Getis Distance images. The 1997 MGD image shows homogeneous areas in red in the middle of cut-block S501 outlined in white that later appear to be heterogeneous in the 1998 M G D image. Conversely, areas in the cut-block that are heterogeneous on the 1997 MGD image apparently became homogeneous after the cut-block was logged. However, in both cases the homogeneity of the source data was largely due to the low radar returns from areas obscured from the sensor. In the 1998 MGD image, the homogeneous area along 69 the west boundary of the cut-block was the result of a radar shadow cast by the timber left standing along that boundary. The extensive homogeneous area further east on the 1997 M G D image disappeared after logging, suggesting that this area was in a radar shadow before it was logged. Aside from these areas, it appears from the M G D data that the clear-cut area on the 1998 6° image is just as heterogeneous as it was when it appeared as a forested area on the 1997 B° image. This is borne out by the coefficient of variation images in Figure 36. The two images were generated by passing a 5x5 kernel over the two MGD images and calculating the coefficient of variation for the centre pixel. The 25 pixels in the kernel were used to compute a local standard deviation and mean. The coefficients of variation were then grouped into classes: 0-20% (red), 20-40% (green), 40-60% (blue), and 60% and over (magenta). Both of these images indicate that the extent of spatial autocorrelation is highly variable in both before and after the cut-block was logged. Based on the MGD data, the coefficient of variation of the cut-block area as a whole is roughly the same in the two images, 61% on the 1998 M G D image and 62% on the 1997 MGD image. Figure 36. 1998 (left) and 1997 coefficient of variation images derived from the MGD images. 70 Figure 37 provides a close-up view of a portion of the two M G D images. A mode filtered 95% threshold bitmap derived from the 1998-1997 MaxGetis difference image is overlaid in vector format7. p i Figure 37: MGD values over a portion of the 1998 (left) and 1997 MGD images. These two close-ups clearly show that the majority of the MaxGetis values within the bitmap area on both the 1997 and 1998 MaxGetis images were contributed by the 11x11 kernel. The bitmap itself indicates that the difference between the 1998 and 1997 MaxGetis values exceeded the 95% threshold level, i.e., the 1998 MaxGetis values were more than 10 units higher than what they were in the 1997 MaxGetis image. In some areas, the increase in MaxGetis value was accompanied by an increase in the size of kernel (area in green in Figure 38). However, in other areas the significantly higher Getis values were attained with either a smaller kernel (area in blue in Figure 38) or a kernel of the same size (area in red in Figure 38). 7 The location of this bitmap has been circled in Figure 33. 71 correlated with changes in MaxGetis values indicates that M G D images can only be used to visually assess changes in the homogeneity of the underlying source data. Increased homogeneity was observed in some areas that experienced an increase in MaxGetis values over the threshold limit; however, there were other similar areas that experienced increased heterogeneity when they were no longer obscured from the sensor. Therefore, it was concluded that MGD images cannot be reliably used for detecting structural changes in the landscape. 4.3.3 Speckle Filtered Data Several references have already been made to the influence that speckle may have on the success of change detection involving both the source data and the data that have been transformed by the Getis statistic. In many quantitative applications of radar data, including change detection, the objective is to remove the unfavourable influence of speckle in such a manner that it preserves the underlying texture and fine structures in the image. A great many speckle filters have been developed in an attempt to do just that. While a review of the subject is beyond the scope of this paper, a recent performance evaluation of more than 30 different speckle filters applied to a standard set of SAR images shows that the maximum a posteriori (MAP) filter, based on a 72 Gamma distributed model of texture and speckle and further enhanced for local detection of structures, performs the best since it preserves both texture and fine features in the scene (Fjortoft era/., 2000). With this advice, the 1997 and 1998 radar brightness data of the cut-block area was filtered using a Gamma MAP filter enhanced for detection of local structures. Structure detection was carried out by the filter software using second order statistics. The limit set for searching for structures was four pixels from the centre of the kernel (Edmond Nezry, personal communication)8. Although not shown here, the 1997 and 1998 speckle filtered images exhibit greater contrast than the corresponding (3° images. Although they are distinctly blurry, fine details such as road rights-of-way are well preserved. A better-looking product would probably have resulted had the data been multi-look rather than single-look (Edmond Nezry, personal communication)9. Figure 39 shows the MaxGetis difference image derived from the speckle-filtered radar brightness data. The difference image itself is not as sharp as the one in Figure 32; however, there is greater contrast overall. This makes it easier to identify areas where significant increases in radar brightness have occurred. It also results in much the same bitmap coverage as in Figure 33 (right) when the data have been thresholded at 95% (Figure 39, right). Very small, isolated bitmaps that were thought to represent false alarms in the unfiltered MaxGetis difference image are almost completely absent. Mode filtering therefore removed only a few of the smallest bitmaps and trimming some of the larger bitmaps (Figure 40). Senior Scientist, ParBleu Technologies Inc., Montreal 'ibid. 73 Figure 39: MaxGetis difference image derived from speckle filtered data (left) with 95% threshold bitmaps overlaid on the right. Figure 40: MaxGetis threshold bitmaps from Figure 39 mode filtered with a 5x5 kernel. The difference image generated from the 11x11 data appears in Figure 41. This image is much smoother in appearance than the equivalent image generated from non-filtered data. Gone are the peculiar vertical and horizontal banding seen in the right-hand image in Figure 26. The 95% threshold bitmaps on the right in Figure 41 are remarkably similar to the mode filtered bitmaps overlaid on the right in Figure 26. They are also quite similar in appearance to the mode filtered threshold bitmaps in Figure 40. 74 Figure 41: Difference image (left) generated from 1997 and 1998 11x11 Getis value images and overlaid with the 95% threshold bitmaps on the right. The 11x11 Getis value images were generated from speckle filtered data. Although MaxGetis and 11x11 difference images derived from speckle filtered data produced much cleaner looking threshold bitmaps than those derived from unfiltered data, it was noted earlier that the Getis statistic already appears to behave as a filter when it transforms radar brightness data to Getis values. Speckle filtering the data beforehand could erase signs of small scale change associated with windthrow in riparian zones and along cut-block margins. Since the desire here is to detect such changes wherever possible, the decision was made not to undermine the performance of the Getis statistic by speckle filtering the data. Instead, mode filtering would be used as a post-processing method to reduce the number of apparent false alarms associated with very small bitmaps in MaxGetis difference images derived from unfiltered data. 4.3.4 Change in the Next Period Logging in the S501 cut-block was the only event between late August, 1997 and late November, 1998 that effected significant structural change in the area examined thus far. This activity substantially increased backscatter on such a large area that the majority of the 95% MaxGetis difference threshold bitmaps were confined to the cut-block area. However, no further anthropogenic changes are known to have occurred in the cut-block or surrounding area after logging was completed in the cut-block. This 75 raised the possibility that incremental changes of a less extensive nature, perhaps brought about by windthrow, could be detected in a subsequent image. With this in mind, the 1999 6° image of the same area was converted to a MaxGetis image. Part of the east boundary of the S501 cut-block runs quite close to the Keogh River. However, it is set far enough west from the river's edge to accommodate the required riparian zone. The exposure of the trees along this boundary and the fact that they are on soft ground suggests there may be a significant risk that windthrow damage could be sustained during winter months when strong south-westerly winds cross the open area of the cut-block. Again image differencing and thresholding at the 95% level, followed by 5x5 mode filtering, was used to assess increases in backscatter between late November, 1998 and late March, 1999 when the two R A D A R S A T images were acquired. The challenge was to correctly identify significant structural changes in the landscape in places where one would expect to find them. Of particular interest was any indication of increased backscatter in the riparian zone along the east boundary of the cut-block that might be attributable to windthrow. Figure 42 shows the 1998 and 1999 radar brightness images. Careful examination of these two images reveals significant differences in the detail, composition, shape and brightness of certain clusters of pixels in areas where no change in backscatter had been expected. One example is highlighted with a yellow circle in the two images. The expectation had been that this area would be somewhat similar in appearance in both images considering that both R A D A R S A T images were acquired in winter, well before any vegetation of any significant size could have become established after logging. It was also expected that the backscatter from areas outside of the cut-block boundary which have not been disturbed would have been similar as weather conditions at the time the two images were acquired were essentially the same. 76 Figure 42:1999 (left) and 1998 p° image of the cut block area. The 1998 and 1999 MaxGetis images were somewhat better for perceiving differences between the two B° images (Figure 43). Because the MaxGetis images have enhanced differences between bright and dark areas on the (3° images, they are better for spotting some of the areas where a significant increase in radar brightness has occurred. Some of the brighter areas have been circled in yellow on the 1999 image. Figure 43: 1999 (left) and 1998 MaxGetis images of the cut-block area. Both of these MaxGetis images were combined in a multi-temporal composite in Figure 44. The 1999 image has been assigned to the red colour gun while the 1998 image has 77 been assigned to both the green and blue colour guns. The composite shows that areas brighter on the 1999 image are red while areas brighter on the 1998 MaxGetis image are cyan. Areas bright at the time both images were acquired are white. The latter signature is particularly evident along the east bank of the Keogh River. While the multi-temporal composite is helpful in identifying areas where significant increases in brightness have occurred, its effectiveness for change detection is compromised by the large number of distinctly red pixel groups scattered everywhere on the image. A more effective approach is to generate a MaxGetis difference image. The 1999-98 MaxGetis difference image is shown on the right in Figure 44. It is overlaid in yellow by the 95% threshold bitmap vectors from the 1998-97 MaxGetis difference image. The 95% threshold bitmaps generated from the 1999-98 difference image are overlaid in red. All of the overlays have previously been mode filtered with a 5x5 kernel. Figure 44: Multi-temporal composite of the 1999 and 1998 MaxGetis images of the cut-block area (left) and the 1999-98 difference image on the right. Both the red bitmaps and the yellow polygons reveal that the backscattering characteristics of particular areas in the landscape were quite different at the time the 1998 and 1999 R A D A R S A T images were acquired. The red bitmaps indicate areas 78 where the backscatter was significantly higher in the 1999 data than it was in the 1998 data. Many of the larger red bitmaps can be easily correlated with patches of bright red pixels on the MaxGetis multi-temporal composite on the left in Figure 44. Raising the threshold to 99% reduced the size and number of bitmaps; however, bitmaps still existed wherever the 95% bitmaps were fairly large. Since there was no further disturbance in the cut-block after the 1998 image was acquired, it was a surprise to see a significant increase in backscatter from that area. Equally surprising was the significant decline in 1999 MaxGetis values within the yellow polygons, the same areas in which significant increases in backscatter had been found on the 1998 data. Had the areas inside the yellow polygons been bright at the time the 1999 image was acquired, the expected colour would have been grey, i.e. similar to that of the area inside the west boundary of the cut-block, and not black as seen in the MaxGetis difference image. The cause of many of these anomalies may be localized differences in soil and surface moisture at the time the two images were acquired (Leckie, 1998). However, it is plausible that in areas of standing timber, localized increases in backscatter could arise as a result of windthrow exposing stems of both standing and windthrown trees to the incident radar wave. In this situation, the increased backscatter could be the result of standing and fallen tree stems acting as corner reflectors. It has also been observed that increased HH-polarized backscatter does occur from horizontally oriented tree trunks (Leckie, 1998). It is therefore possible that the bitmap circled in black in Figure 44 signifies increased backscatter due to windthrow since it occurs on the exposed side of the riparian zone established along the Keogh River. 79 Chapter 5 Testing the MaxGetis Differencing Technique The previous chapter demonstrated that the Getis statistic could perhaps provide a viable means of detecting both large- and small-scale structural change in the landscape. The best results were obtained by differencing two co-registered MaxGetis images, thresholding the difference image at 95% and then mode filtering the change bitmaps to eliminate the smallest ones that were considered to be largely false alarms. The results of applying this approach to four areas where windthrow is known to have occurred will be described in this chapter. The first two areas, Strip 31a and Strip 16, are riparian zones that were established to protect streams running through the middle of their respective cut-blocks. The other two cases involve sizable areas where windthrow occurred after the cut-block had been logged. Since one of these areas has been salvage logged while the other has not, an opportunity was provided to compare the strength and consistency of backscatter from areas that have experienced different levels of disturbance. These two cases also provided an opportunity to model damage detection in riparian zones that are located along cut-block boundaries. The results of applying the MaxGetis differencing technique to six other riparian zones where windthrow damage has been reported will also be discussed briefly. 5.1 Strip 31a Strip 31a, first illustrated in Figure 3, is a riparian zone that cuts through the middle of cut-block M5304. The red-green-blue multi-temporal composites (MTC) of the B° and MaxGetis images in Figure 45 clearly show a dark riparian zone cutting through a 80 comparatively bright cut-block. The cut-block itself is predominantly blue which is consistent with its highest reflectivity in the 1996 image, the first image acquired after logging had been completed in 1995. Clusters of red, yellow, and green pixels within and adjacent to the riparian zone are in areas that were dark on the 1996 image but were brighter on images acquired later. Bright clusters of red and yellow pixels within the zone are apparently indicative of windthrow damage that occurred between the time the 1996 and 1999 images were acquired (Peter Murtha, personal communication) 1 0. Figure 45: Multi-temporal composite of the 1999, 1997, and 1996 p° images (left) and the equivalent as a MaxGetis multi-temporal composite. The 1999-96 MaxGetis difference image with its 95% threshold bitmaps overlaid in red is shown on the left in Figure 46. Without the bitmaps, the riparian zone appears as a bright ribbon running through a much darker cut-block. Mode filtering these bitmaps removed many of the small ones while at the same time preserving the larger ones (Figure 46, right). The cyan, magenta, and yellow polygons were created from the mode-filtered 95% threshold bitmaps generated from the 1997-96, 1998-97, and 1999-98 MaxGetis difference images. The high degree of overlap between the 1999-96 bitmaps in red and the 1999-98 bitmaps shown as yellow polygons suggests that much of the significantly increased backscatter in the riparian zone occurred between 1998 1 0 Professor Emeritus, Faculty of Forestry, University of British Columbia 81 and 1999. On their own, the red bitmaps appear to correlate reasonably well with the bright red and yellow clusters in the riparian zone in Figure 45, left. Although this could be construed as being consistent with the location of the reported windthrow, the presence of significant bitmaps in the adjacent cut-block suggests that other factors may be responsible for the increased backscatter in the riparian zone. Figure 46: 1999-96 MaxGetis difference image (left) overlaid with 95% threshold bitmaps. On the right, the same difference image has been overlaid with the mode filtered 1999-96 bitmaps in red and the mode filtered bitmaps from the 1997-96, 1998-97, and 1999-98 difference images as cyan, magenta, and yellow polygons. 82 5.2 Strip 16 Strip 16, shown on the 1995 aerial photograph in Figure 47, is located in cut-block S505 close to the Port Hardy Airport. Figure 47: Aerial photo showing location and condition of Strip 16 in 1995. Although the photograph clearly shows that Strip 16 existed after logging had been completed in 1995, it was reported that the riparian zone had been decimated by the time the December 1996 image was acquired (Murtha, 2000b). An examination of the four R A D A R S A T images acquired between December 1996 and March 1999 appears to confirm this. There is no definite band of dark pixels where the riparian zone should 83 be as was the case with Strip 31a. However, a 25-metre resolution LANDSAT image, acquired in August 1997 and shown in Figure 48, suggests that there were perhaps some trees remaining at the eastern end of Strip 16. This is more strongly supported by a 15-metre resolution A S T E R image of the same area acquired in November 2000 (Figure 48). In addition, a somewhat darker pixel at the western end of Strip 16 on the LANDSAT image indicates the possibility a few trees may have been still standing in that area in 1997. A careful examination of the 1996 and 1997 R A D A R S A T images shows a rather dark patch in that area which could be a radar shadow cast by some trees. The same area is somewhat brighter on the 1998 image and very bright on the 1999 image. The significant brightening in 1999 could perhaps be due to the disappearance of the radar shadow as the trees were blown over. The pixels at the location of this supposed windthrow event show up quite clearly in red on the multi-temporal composite in Figure 49, consistent with the substantially higher reflectivity of this area on the 1999 image. The A S T E R image acquired a year later confirms that there are no standing trees at the west end of the riparian zone. Consequently, it was assumed that the changes in reflectivity recorded in that area on the RADARSAT images were due to a windthrow event that occurred after the 1997 R A D A R S A T image was acquired. 84 Figure 49: Multi-temporal composite of the 1999,1997, and 1996 p° images. The multi-temporal composite in Figure 49 also shows two other bright red patches of pixels at the eastern end of Strip 16. There are also some bright red patches of pixels along the east boundary of the cut-block, both north and south of Strip 16. These may also be indicative of windthrow damage in Strip 16 and the riparian zones situated along the cut-block boundary. Of the several MaxGetis difference images thresholded at 95%, the 1999-96 and 1998-97 difference images provided a bitmap that matched the location of the increased brightening at the western and eastern ends of Strip 16. Figure 50 shows the 1999-96 mode filtered bitmaps overlaid on the difference image and on the aerial photograph. The particularly bright patches of pixels seen at both ends of the riparian zone in Figure 49 are indeed highlighted by fairly large red bitmaps. However, very few bitmaps match up with the bright red patches of pixels seen along the east boundary on the multi-temporal composite. 85 Figure 50: 1999-96 MaxGetis difference image with mode filtered 95% threshold bitmaps overlaid in red. The bitmaps have been overlaid on the aerial photo on the right. Also evident in Figure 50 are a large number of sizable bitmaps scattered throughout the cut-block as well as in adjacent areas that have well established advanced regeneration. The associated pixels on the 1999 B° image are indeed much brighter than the corresponding ones on the 1996 B° image. This is in spite of the fact that the cut-block is generally darker on the 1999 B° image due to increased scattering from vegetation that has become established in the intervening period. Raising the threshold to 99% reduced the number of these bitmaps without removing the significant ones in Strip 16; however, a considerable number of large bitmaps remained in the general area. Surprisingly, the bitmaps generated from the 1998-97 difference image were considerably larger and more prevalent in the cut-block. The presence of so many large bitmaps in both the 1999-96 and 1998-97 difference images raised doubts about whether or not the increased backscatter at both ends of Strip 16 was in fact due to windthrow or to some other factor that was also responsible for increased backscatter in the adjacent areas. 86 5.3 Edge Blowdown Further south of Strip 16 is a cut-block that suffered serious windthrow damage in early 1997 along the eastern half of its northern boundary (Peter Murtha, personal communication) 1 1. The attendant change in spectral signature of the normally green forest in that area is outlined in Figure 51 on a LANDSAT image acquired in August 1997. Also visible but less obvious is a right-of-way for a road that was under construction at the time to salvage the timber. Figure 51: 1997 LANDSAT image showing location of new cut-block and road right-of-way leading to it. A multi-temporal composite of the 1999 (red), 1997 (green), and 1996 (blue) R A D A R S A T images shows the boundary of the cut-block in red and the extent of the blowdown within the green boundary (Figure 52). 1 1 Professor Emeritus, Faculty of Forestry, University of British Columbia 87 Figure 52: Multi-temporal composite of the RADARSAT p° images showing the same blowdown damage in the northeast corner of the cut-block. The preponderance of red pixels in the blowdown area in Figure 52 indicates that the majority of the increased backscatter from that area occurred after the 1997 radar data were collected. This was unexpected since the 1997 R A D A R S A T image was acquired in the same month as the LANDSAT image and therefore it was expected that the blowdown area in the multi-temporal composite would be dominated by bright yellow pixels in response to elevated backscatter levels in both 1997 and 1999. However, when the radar images are viewed separately, the blowdown area is generally quite dark on both the 1996 and 1997 RADARSAT images. Therefore it is reasonable for that part of the composite to be dominated by red pixels. The 1997-96 MaxGetis difference image in Figure 53 suggests that there was a substantial increase in backscatter in 1997 in the blowdown area in a very small area. Since such increases have occurred in presumably undisturbed areas nearby, as indicated by the bitmaps in Figure 54, it is not certain that the increased backscatter from the blowdown area can in fact be attributed to windthrow. 88 Figure 53. 1997-96 MaxGetis difference image of the blowdown area overlaid with mode filtered 95% threshold bitmaps. As Figure 51 shows, the LANDSAT image revealed the extent of the blowdown very well because dead trees exhibit a distinct spectral signature that allows them to be separated from living trees. Radar sensors operate on an entirely different principle in that they record the strength of backscatter from features in the landscape. Therefore, to explain the rather weak backscatter from the area, it is speculated that the radar signal interacted with the windthrown trees in much the same manner as it does with standing trees that also return low levels of backscatter. It is not until later when the trees were removed that suitable conditions were created for much higher levels of backscatter to be returned to the sensor (i.e., when there is a clear shot to the ground and there are sufficient targets to reflect the radar signal back to the sensor). From this example, one could conclude that windthrown trees are not easily detectable since radar returns from these targets appear to be quite weak. However, significant increases in backscatter due to the salvage logging are easily detected as illustrated by the large bitmaps on the 1999-96 and 1998-97 MaxGetis difference images in Figure 54. The bitmaps on the 1998-97 MaxGetis difference image indicate that the salvage logging appeared to have been completed before the 1998 radar image was acquired. 89 Figure 54:1999-96 MaxGetis difference image (left) and 1998-97 MaxGetis difference image with mode filtered 95% threshold bitmaps. The bitmaps in Figure 54 were generated by thresholding the difference image at 95%. Raising the threshold to 99% slightly reduced the size of the bitmaps in the blowdown area while considerably reducing the number of minor bitmaps elsewhere. However, neither image shows the extent of the disturbance as well as the LANDSAT image. Nevertheless, its significance is clearly flagged by the bitmaps, although the disturbance detected was logging rather than windthrow. 5.4 Edge Blowdown at M5307 The last area examined in detail was the northeast corner of cut-block M5307. Logged just before the 1997 R A D A R S A T image was acquired, it is located to the south of cut-block M5304 where Strip 31a is situated. A multi-temporal composite of the 1999, 1997, and 1996 MaxGetis images, assigned to the red, green, and blue colour guns respectively, shows the cut-block boundary in red and a second boundary in green demarcating the area where the blowdown is known to have occurred (Figure 55, left) (Peter Murtha, personal communication) 1 2. Professor Emeritus, Faculty of Forestry, University of British Columbia 90 The cut-block area is dominated by patches of bright red, green, and yellow pixels while the blowdown area features distinct areas of red, green, and blue pixels. The red, green, and blue pixels indicate that the highest MaxGetis values were on the 1999, 1997, and 1996 MaxGetis images respectively. Although blurry, the colouration in the blowdown area correlates quite well with the distribution of similarly coloured pixels on the B° composite in Figure 55, right. Figure 55: Multi-temporal composites of the 1999, 1997, and 1996 MaxGetis (left) and B° (right) images of the northeast corner of cut-block M5307. Figure 56, left, illustrates the extent of bitmaps generated from three MaxGetis difference images that used the 1999 image as a reference. The cyan polygons correspond with the extent of the 95% threshold bitmaps generated from the 1999-98 MaxGetis difference image while the white and black polygons correspond with the extent of the bitmaps generated at the same threshold level from the 1999-97 and 1999-96 MaxGetis difference images. 91 Figure 56: The 1999-1997-1996 MaxGetis multi-temporal composite overlaid with threshold bitmap vectors derived from various MaxGetis difference images to highlight areas of significant change. On the right are the same bitmap vectors after 5x5 mode filtering. As expected, the 1999-96 MaxGetis difference image bitmaps (outlined in black) are largely confined to the cut-block which in 1999 still exhibited much higher backscatter than that captured by the 1996 RADARSAT image. However, the extent of the bitmaps in the cut-block area is significantly less than that obtained by thresholding the 1997-96 difference image. The reduced extent of high backscatter may be due to a rise in canopy volume scattering as vegetation progressively obscured slash, stumps, and other logging debris which had earlier been responsible for larger areas exhibiting strong radar returns on the 1997 image. The bitmaps within the blowdown area indicate that the 1999-97 MaxGetis difference image generated the most extensive 95% threshold bitmaps (outlined in white in Figure 56, left). Thresholding the 1999-98 and 1999-96 MaxGetis difference images also generated bitmaps in the same general area (1999-98 in cyan, 1999-96 in black). However, when checked against two-colour composites of the respective (3° images, they did not correspond with any of the high backscatter pixels in the 1999 image. Figure 57 shows that 95% threshold bitmaps generated from the 1999-98 and 1998-97 MaxGetis difference images were not consistent with those derived from the 1999-97 92 MaxGetis difference image. All three MaxGetis images show that the northeast corner of the blowdown area was bright only on the 1998 image; on the other two images, that corner is rather dark, even in 1999. The high backscatter in 1998 is also quite evident in the cut-block where a substantial number of bitmaps were generated from the 1998-97 MaxGetis difference data. However, the cut-block area is noticeably darker in the 1999 image. This led to the conclusion that some factor other than windthrow is probably responsible for the high backscatter from the blowdown area in 1998 since the elevated level of backscatter is not present in the 1999 data. Figure 57: The 1999-1998-1997 MaxGetis multi-temporal composite overlaid with the 1999-97 (white), 1999-98 (cyan), and 1998-97 (black) threshold bitmap polygons to highlight areas of significant change on the 1999 image. On the right are the same bitmap polygons after 5x5 mode filtering. Figure 57 also shows that the left half of the blowdown area exhibited the highest reflectivity in 1996 as indicated by the dark blue pixels in that area. The same area is quite dark on the 1997, 1998, and 1999 images. Taking a cue from the low reflectivity of the area immediately above, the explanation may be that once the trees had fallen, fewer scattering elements were visible to the satellite because of steep terrain obscuring the area. However, this explanation has not been verified in the field. 93 In summary, the blowdown in the area adjacent to the northeast corner of cut-block M307 could not be identified with certainty because higher levels of backscatter from both the eastern half of the blowdown area and adjacent areas recorded in the 1998 data were followed by noticeably lower levels in the 1999 data. There was also insufficient information available to confirm that reduced levels of backscatter in the western half of the blowdown area were the result of windthrow reducing the number of scatterers visible to the sensor. 5.5 Other Areas Six other riparian areas were examined to determine if there were significant MaxGetis value differences that could perhaps be associated with windthrow damage. Strips 17 and 18 are located along the east boundary of cut-block S505, running north and south off the eastern end of Strip 16. The other four are located a short distance away to the southeast. Strips 13 and 25 are situated along the southwest and eastern boundaries of a cut-block that was logged in 1995. Strip 26 runs through the middle of another cut-block nearby that was also logged in the same year. Strip 28 is situated along the eastern boundary of a cut-block that was logged several years before. The 1999-96 difference image was chosen for the sake of consistency since it was used for Strip 16. Except for Strip 17, less than 10% of the area of each of the six riparian zones showed significant increases in MaxGetis values between the 1996 and 1999 image (Table A-3, Appendix A ) 1 3 . Although windthrow to varying degrees has been reported in all of these riparian areas (Peter Murtha, personal communication) 1 4, it is not known how well the location of the pixels exceeding the 95% threshold level agree with the specific location of windthrow or if they indeed indicate any disturbance at all. Although some of the bitmap area did correspond with areas that exhibited high backscatter in the 1999 image, there was a significant area that did not, an issue that is discussed further in the next chapter. The proportions shown in the table would be lower had the bitmaps been mode filtered. Professor Emeritus, Faculty of Forestry, University of British Columbia 94 Table A-3 also shows that the mean MaxGetis difference value within some of the features examined here is generally higher than for the image as a whole. However, in Strips 16, 25, and 26 the mean is very close to that of the image itself. In these cases, the areas were not easily distinguished from the surrounding area. This highlighted the need to overlay riparian zone and cut-block boundaries first before any meaningful interpretations could be made of high backscatter from any particular area. The only area where there was a sizable clustering of bright pixels that could be associated with a known disturbance with certainty was in the edge blowdown area described in Section 5.3. In that case, the data suggest that it was the salvage logging that was responsible for the elevated backscatter and not the windthrow itself. The few examples discussed in this chapter indicate that the radar data used here could not be used to reliably detect windthrow in areas where it might occur because it was uncertain that factors responsible for elevated backscatter in adjacent undisturbed areas were not also responsible for elevated backscatter levels from areas where windthrow may or may not have occurred. In one area where windthrow had occurred, a LANDSAT image detected the spectral signature of dead timber in the windthrow area; however, a radar image acquired at about the same time did not record any significant increase in backscatter from that area. It was only later when salvage logging had been completed that the affected area returned elevated levels of backscatter to the radar sensor. However, this one example cannot be used to rule out windthrow causing elevated levels of backscatter from areas being monitored for such a disturbance. Furthermore, it was found that the MaxGetis differencing and thresholding method did not reliably identify all areas that were bright on one image and correspondingly dark on an image acquired under similar conditions earlier. Had it been more consistent, it could have been used to identify areas which experienced significant increases in backscatter. These bitmaps could have acted as flags indicating the possibility of windthrow in areas where it is expected it could occur. The next chapter discusses 95 some of the issues that may have influenced the lack of success in identifying windthrow damage with the aid of the Getis statistic. 96 Chapter 6 Discussion This thesis set out to evaluate the utility of the Getis statistic for detecting and mapping the extent of windthrow damage in riparian management areas using multi-temporal radar images acquired under wet conditions. Intact riparian zones and forested areas are typically much darker on these images than adjacent freshly logged areas because of considerable canopy volume scattering losses and significant areas of radar-induced shadow. On an image acquired after the event, windthrow has been characterized in some instances by a significant brightening of pixels in the area where it has occurred (Murtha, 1998a, 1998b, 1998c, 2000a, 2000b). It was postulated that the Getis statistic, designed to detect both an increase in the magnitude of values and degree of spatial autocorrelation associated with them, could be used to detect the event and its extent. A predicted increase in spatial autocorrelation associated with windthrow was based on the premise that a group of pixels associated with standing timber exhibit heterogeneous values, and therefore low spatial autocorrelation, while brighter areas associated with windthrow damage and logging activity exhibit comparatively homogeneous pixel values and therefore higher spatial autocorrelation. It was suggested that the increased brightness and higher spatial autocorrelation could be used to detect and map the extent of windthrow damage. Based on trials first conducted with cut-block S501, it was found that the Getis statistic could indeed be used for detecting large-scale structural changes in the landscape caused by timber harvesting. In its application here, the Getis statistic behaved somewhat like an averaging filter with the enhancement of both bright and dark areas being an added effect. Highly positive Getis values were associated with bright areas in the image whereas highly negative Getis values were associated with dark areas. Areas 97 of increased brightness in the cut-block were effectively revealed by differencing two MaxGetis value images. This was found to be sufficient to perceive the location and extent of the freshly logged area. The difference image could also be thresholded at 95% to reveal the area where the most significant increases in MaxGetis value had occurred in the cut-block between one imaging date and the next. The results obtained from this approach were considered more meaningful for detecting change than the next best alternative, the CHDET algorithm, because no data were lost to speckle filtering. The absence of such filtering allowed thresholding to capture most of the significant increase in scene brightness that had occurred between two imaging dates. However, thresholding did not adequately identify the extent of the disturbance that could be perceived by eye because not all areas exhibited significantly higher MaxGetis values. Unexpectedly, a corresponding increase in the extent of high spatial autocorrelation did not occur in the cut-block, something that was presumed to occur and would be useful for mapping the extent of structural change. For many pixels, it was found that the size of kernel which generated the MaxGetis value could just as easily decrease as increase after the trees had been removed. In fact, the proportion of MaxGetis values contributed by the various kernels was about the same in the cut-block area before and after logging with the majority being split between the 3x3 and 11x11 kernels. It also revealed that freshly logged areas could be as heterogeneous as forested areas for this type of radar data. MaxGetis differencing and thresholding of the S501 cut-block appeared to offer promising results for detecting a significant proportion of structural change known to have occurred in the scene and therefore the approach was adopted in the work that followed. However, this initial success was largely due to the fact that almost all of the change detected by differencing and thresholding such a small dataset (256 line by 256 pixel) was confined to the cut-block area. An assessment of significant change in the following period, as well as with changes associated with timber harvesting and windthrow on all of the other images that were 1024 pixels by 1024 lines in size, 98 demonstrated that thresholding MaxGetis difference images produced an unacceptably high number of what can only be described as false alarms since no structural change could be associated with the majority of them. The number of apparent false alarms varied from one set of images to another depending on the extent of structural change in the scene, i.e., a set of images with one recently logged cut-block had a far higher number of apparent false alarms than another set that had several recently logged cut-blocks. The most obvious source of the problem may be the threshold level itself. A 95% threshold level means that 5% of the pixels on the difference image will be identified as representing significant change irrespective of the reason for the significant increase in MaxGetis value. Raising the threshold level would be an obvious solution, however, it was found in all cases that raising it to as high as 99% did not entirely eliminate false alarms. There were still a considerable number of instances where there was a significant increase in MaxGetis value with no known structural change associated with them. Part of the problem can be attributed to MaxGetis differencing itself. MaxGetis values range from positive to negative with 90% of the values typically distributed between -6 and +8, although the range can be wider and more asymmetrical in some datasets. Depending on the spatial arrangement of MaxGetis values, it is possible to obtain significant differences between MaxGetis values at the 95% threshold level of 7.5 for pixels that appear to be grey on one MaxGetis image and black on another. The same could arise between pixels that are very bright on one image and not so bright on the other. The MaxGetis values relate back to the corresponding values on the 3° image from which the MaxGetis values were generated. In some cases, the MaxGetis value is higher or lower than what would be suggested by the (3° value simply because the MaxGetis value for that pixel is influenced by the B° values of the adjacent pixels in the kernel. 99 Another contributor may have been speckle noise, something that is inherent in all radar data, particularly in datasets where higher spatial resolution is provided at the cost of lower radiometric resolution (Raney, 1998). It was decided early on that the data should not be speckle filtered, lest it alter data necessary for detecting small scale changes that would commonly be associated with isolated windthrow events. The way in which the Getis value is computed for each pixel in the image certainly suggests that it would be suitable for application to unfiltered radar data, since the kernel method of processing the data would reveal the nature of the data in the absence of speckle 1 5 . This was supported with the observation that the slightly larger bitmaps generated from speckle filtered data are quite similar in number and appearance to those derived from unfiltered data. The relatively minor reduction of small bitmaps in the speckle filtered difference data suggests that speckle noise did not seriously interfere with the consistent detection of significant increases in MaxGetis value. However, it did indicate that speckle noise may be responsible for some of the apparent false alarms. Some of the false alarms may have also been generated as a result of the data being single look. Due to Gaussian scattering within a particular resolution cell or pixel, the true reflectivity of a target can only be reliably determined from a large number of measurements or samples (Raney, 1998). For the 6° amplitude format data used here, these measurements will have a Rayleigh distribution, the same distribution that characterizes the entire image (Oliver and Quegan, 1998). Since the reflectivity of a pixel is derived from just a single sample (or look), its observed brightness can vary considerably from one image to the next due to Gaussian scattering of the incident radar wave (Raney, 1998). This variability is further compounded by the pulse repetition frequency (PRF) of the radar instrument. For the R A D A R S A T data used here, the P R F is around 1254 Hz. This parameter influences the perspective from which the satellite 'sees' the target each time the backscatter data are acquired. Consequently, the target In fact, the 5x5 Getis value images appeared to be identical in appearance to the Gamma M A P filtered images. 100 could conceivably be dark on one image and bright on the next, leading to a potential false alarm in the difference image. Despite the fact that all four parent images were acquired under similar conditions, weather, and other factors may have had more of an influence on backscatter than initially thought. Seasonal and environmental factors, particularly phenological state and moisture content of the vegetation, as well as variability in soil and surface moisture conditions, are known to greatly influence backscatter signal strength. They can also influence the characteristics of speckle noise in the data (Leckie, 1998). It is certainly possible that local surface wetness was sufficiently different at each of the times the data were acquired. This could certainly vary across the area during a particular scene acquisition, leading to localized brightening and darkening which may in itself generate false alarms. Perhaps a more fundamental issue that influenced the appearance of backscatter in multi-temporal B° composites was one particular characteristic of the radar data used in this research. All of the B° images are quite grainy with strings and blocks of bright and not-so-bright pixels separated by valleys filled with dark pixels. This pattern of graininess, quite unlike the more homogeneous optical data collected by satellites such as LANDSAT, differs significantly among images as can be seen in Figure 4 on page 11. It is particularly obvious when rapidly flipping back and forth between two B° images on a computer screen. When three images are put together into a multi-temporal composite, the graininess is preserved in the arrangement of the pixels (see Figure 55, right). The brightness and colour of each pixel in the composite is very much dependent upon the strength of the backscatter recorded in all three images at that particular pixel location. Although converting B° to Getis value images reduces this effect, the inherent graininess of the data could significantly influence the outcome of MaxGetis image differencing. The stability of the RADARSAT sensor is itself not in question. The calibration and stability of the sensor is verified continuously using RADARSAT-1 Precision 101 Transponders and multi-temporal analysis of data collected over a large swath of Amazonian rainforest (Srivastava et al., 2003). A recent analysis of 780 datasets acquired over that area showed that the R A D A R S A T sensor is remarkably consistent in all beam modes, the radiometric error being no more than about ±0.5 dB. A noticeable annual cyclic variability of lesser magnitude was thought to be due to inherent variations in backscatter from the site or variability in instrument gain (Luscombe, 2003). In the end, the number of apparent false alarms was arbitrarily reduced so that the analysis could more correctly identify where extensive increases in backscatter had occurred. The method used here was to post-process the bitmap coverage so that only bitmaps over a certain size would be retained. The reasoning behind this was that the larger ones were always correlated with structural change while the smaller ones removed in the post processing were not in the vast majority of cases. The preferred method would have been to use something similar to the MINNEIGH feature in the CHDET algorithm. However, this functionality was not available outside of that algorithm. Consequently, two alternative approaches were evaluated. The first was SIEVE, a tool commonly used to post-process bitmaps generated by classification algorithms. This tool amalgamates small bitmaps under a specified size into adjacent bitmaps of another class. It preserves the geometry of bitmaps over the specified size, an attractive feature; however, narrow bitmaps that are only one pixel wide are also retained if they contain a sufficient number of adjoining pixels. Retaining such narrow bitmaps (one pixel wide is equivalent to just over 3 metres) was considered to be undesirable in this application and therefore SIEVE was dropped in favour of the second approach. The second approach makes use of the mode filter. This filter closely mimics the MINNEIGH functionality in that it takes into account the number of neighbours that a particular pixel has before removing the pixel from the coverage. In its implementation here, the bitmap segment was converted to an image plane. In this state, the bitmap pixel had a value of 1 and the pixels in between the bitmaps had a value of 0. Although 102 the mode filter can use a kernel from 1x3 to 7x7 in size, trials suggested that the most useful size was 5x5. Using a 5x5 kernel centred upon a particular pixel, the mode filtering program counted the number of neighbouring pixels within the kernel that had a value of 1. If the number was 12 or more, the centre pixel was retained. Unlike MINNEIGH, which is run iteratively until changes in the bitmap coverage no longer occur, the mode filter was run only once in its application here. The 5x5 kernel removed a substantial number of small bitmaps while smoothing the edges and filling in holes in the larger bitmaps. This effect on larger bitmaps was considered inconsequential as the bitmaps were to have been used to indicate where backscatter had increased significantly from one imaging date to the next. However, the number of apparent false alarms remained far too great to be acceptable for reliable detection of structural change in the landscape. Therefore, it was concluded that differencing Getis and MaxGetis images and thresholding the resulting difference images was inappropriate for the type of radar data used in this research. Since areas that exhibited elevated levels of backscatter did not exhibit a greater degree of homogeneity than forested areas as originally thought, it appeared that the Getis statistic had no value whatsoever for change detection. However, structural features, areas of high backscatter and contrast are all enhanced on multi-temporal MaxGetis composites and consequently they are much easier to interpret than multi-temporal B° composites. Therefore, these could potentially be quite useful for visual detection of change in areas where it is expected to occur. 103 Chapter 7 Conclusion 7.1 Suggestions for Further Research By comparing the spectral signature of dead trees in a blowdown area on a LANDSAT image against the extent and magnitude of backscatter captured from the same area at about the same time on a R A D A R S A T image, it was determined that windthrow damage generally did not exhibit noticeably higher backscatter than adjacent standing trees in the data used here. This of course does not rule out the possibility that such damage can be detected more effectively by other types of radar data. Fine 2 Mode, single look R A D A R S A T data was used here because it was the highest spatial resolution data available and therefore offered the greatest opportunity to detect small, isolated windthrow events in riparian zones. Other data at approximately the same spatial resolution could be obtained by the C-band Advanced Synthetic Aperture Radar (ASAR) instrument on board the ENVISAT satellite (European Space Agency, 2002). These data can be collected in W , V H , and HV polarization in addition to the HH polarization of the R A D A R S A T Fine 2 Mode data used here. Similar data at 3 and 10 metres resolution could also be obtained from RADARSAT-2 , due to be launched in 2005 (RADARSAT International, n.d.). Data collected at a polarization other than HH may be more sensitive for detection of increased backscatter associated with windthrow damage and therefore its potential for detecting windthrow damage should be evaluated. Multi-look data at a higher spatial resolution than previously available could also be acquired by RADARSAT-2 . Multi-look data reduce uncertainty over the true magnitude of backscatter from a target by averaging several samples obtained 104 from different vantage points and therefore they have a potential for providing more consistent measurements of backscatter from areas being monitored on a continuous basis. Previously available multi-look data have had a maximum spatial resolution of 30 metres which is more appropriate for differentiating change at the cut-block scale. However, the forthcoming availability of 10 metre resolution "multi-look fine" imagery from RADARSAT-2 could perhaps be usable for monitoring small scale windthrow events, provided these areas exhibit significantly higher backscatter than adjacent forest areas. If exploratory studies suggest that other types of data are able to detect significantly higher backscatter resulting from windthrow, then a rigorous field-based approach should be undertaken over several years to correlate elevated levels of backscatter in riparian zones and along cut-block boundaries with windthrow damage and other factors. A selection of riparian zones identified on the local licensee's cutting permit development plan maps should be monitored from before the associated cut-block is logged to several years after. Since surface and soil moisture have a significant influence on the characteristics and level of backscatter returned to the satellite, the radar data used in the study should all have been acquired under similar weather conditions. Data acquired during the winter appear to be the best for the North Vancouver Island area since the soil is typically saturated and the vegetation is frequently quite wet thus ensuring that cut-blocks are highly visible. A number of weather stations should be established in the immediate area to record weather conditions at the time the images are acquired. Soil and humus moisture content should also be measured in the riparian zones and adjacent cut-blocks at the time of image acquisition. Site conditions can conveniently be captured by camera. Vertical and oblique aerial photos taken with a hand-held digital camera taken from a helicopter or light aircraft should be taken just prior to the radar data being acquired. To facilitate registration of both aerial photos and radar images, registration markers 105 and trihedral corner reflectors should be set up and their positions determined with a G P S receiver. Once a radar image has been received from the agency responsible for acquiring it, the forest cover map showing the location and extent of the riparian zones being monitored should be registered to the image and the location of areas exhibiting high backscatter should be digitized. The registration of the map to the image can then be reversed so that the correct ground coordinates of the high backscatter areas can be obtained. A G P S receiver can then be used to visit the areas on the ground so that factors that may be responsible for the elevated backscatter, such as soil and surface moisture conditions; density, phenological state and morphological characteristics of vegetation; configuration and characteristics of scatterers; and microtopography, can be determined. Although MaxGetis multi-temporal composites are useful for visualizing areas exhibiting high backscatter, the characteristics of radar data are such that the Getis statistic does not provide any additional benefits that would allow it to be used for detecting and mapping the extent of structural changes in the landscape. Consequently, some other approach will be required. 7.2 Summary of Results The research conducted for this thesis was principally concerned with determining the effectiveness of the Getis statistic for detecting structural change in the landscape associated with timber harvesting and windthrow. Both of these types of disturbance were presumed to be detectable by an increase in backscatter in the area in which they occurred. Elevated levels of backscatter can be visually perceived in multi-temporal composites or by thresholding difference images. One established method for elucidating significant increases in backscatter applied a threshold to a speckle filtered, logarithmically scaled ratio image. However, it was found that this approach did not capture the full 106 extent of timber harvesting in a cut-block chosen for testing change detection methods. Based on results obtained from the same test cut-block, Getis value difference images proved to be more helpful than logarithmically scaled ratio images for determining the extent of the increased backscatter caused by timber harvesting. Difference images generated from the 11x11 Getis value images appeared to reveal the greatest extent of change in the test cut-block. However, such a large kernel also tended to obscure small scale changes. Difference images generated from MaxGetis value images appeared to be better as they enabled optimal detection of both small and large scale changes in backscatter. Thresholding these images to highlight areas where backscatter had increased significantly also appeared to provide more promising results, particularly after the 95% threshold bitmaps had been mode filtered to reduce the number of very small bitmaps associated with apparent false alarms. However, changes in the distance over which the highest local spatial autocorrelation associated with the maximum Getis value for a particular pixel were inconsistent and resulted in the finding that the pattern of backscatter from freshly logged areas was not homogeneous as originally presumed but that it was just as heterogeneous as it was from forested areas. The lack of homogeneity therefore made it impossible to map the extent of the disturbance. Reducing speckle noise through filtering was also determined to be inappropriate because it diminished the magnitude of differences between pixels needed for detecting small scale changes. Despite the initial success in using MaxGetis images to reveal the location and extent of significant increases in backscatter caused by logging in the test cut-block, applying the technique to other images to detect structural change resulted in a significant number of false alarms, even after the thresholds had been raised to 99%. The prevalence of these false alarms made it difficult to be sure that the significant increases in backscatter that emanated from riparian zones were in fact due to windthrow, especially after it was discovered that blowdown evident 107 on a LANDSAT image was not particularly evident on a R A D A R S A T image acquired at about the same time. Therefore, it was concluded that thresholding MaxGetis difference images was not a reliable means of detecting significant increases in backscatter due to windthrow. However, it was found that multi-temporal composites of MaxGetis images could be particularly useful for interpretation of features and visualizing macro scale changes in the landscape attributable to logging and road construction due to their enhancement of contrast and structure in the scene. 108 References Anselin, L. 1995. Local indicators of spatial association - LISA. Geographical Analysis, 27, 93-115. Armour, B., Farris-Manning, P., Vachon, P.W. 1997. The suitability of RADARSAT 's orbit for repeat pass S A R interferometry. Proceedings Geomatics in the Era of RADARSAT (GER '97), 25-30 May 1997, Ottawa, Canada. C D - R O M . Bao, M. 1999. Backscattering change detection in SAR images using wavelet techniques. Proceedings IGARSS '99, 28 June - 2 July, 1999, Hamburg, Germany. C D -ROM. B C M O F and BC Environment. 1995. Riparian management area guidebook. British Columbia Ministry of Forests and BC Environment, Victoria, B.C. http://www.for.gov.bc.ca/tasb/legsregs/fpc/fpcguide/riparian/Rip-toc.htm Bruniquel, J . , Sassier, H., Peraudeau, S., Rognant, L., Goze, S., Chust, G. , Ducrot, D., Fellah, K., Laur, H. 1999. Assessment of multi-temporal products for multi-thematic applications with E R S S A R data. Proceedings CEOS SAR Workshop (CEOS '99), 26-29 October, 1999, Toulouse, France. http://www.estec.esa.nl/conferences/99b02/ Castel, T., Martinez, J-M. , Beaudoin, A., Wegmuller, U., Strozzi, T. 2000. E R S INSAR data for remote sensing hilly forested areas. Remote Sensing of Environment, 73, 73-86. Chen, Y., Zheng, F-H., Zhang, N-X., Jin, Y -Q. 1999. Monitoring and evaluation of flooding extent by using the Getis statistic of the SSM/I data. Proceedings IGARSS '99, 28 June - 2 July, 1999, Hamburg, Germany. C D - R O M . Chrisman, N. 1997. Exploring Geographic Information Systems. John Wiley & Sons Inc. New York. 298 pp. Coulson, S.N. 1995. SAR Interferometry with E R S . http://earth.esa.int/printer_friendly/rootcollection/sysutil/inter_steve.html 109 Couturier, S., Liew, S . C , Nakayama, M., Lim, H. 1999. Monitoring vegetation regeneration in fire-affected tropical forests using E R S / J E R S synthetic aperture radar. Proceedings IGARSS '99, 28 June - 2 July, 1999, Hamburg, Germany. C D - R O M . Dekker, R. J . 1998. Speckle filtering in satellite S A R change detection imagery. Int. J. Remote Sensing, 19, 1133-1146. Derksen, C , Wulder, M., LeDrew, E., Goodison, B. 1998. Application of the Getis Statistic to hemispheric and regional scale passive microwave derived snow water equivalent imagery, Proceedings IGARSS '98, 6-10 July 1998, Seattle, USA. C D - R O M . Dierking, W., Schou, J . , Skriver, H. 2000. Change detection of small objects and linear features in multi-temporal polarimetric images. Proceedings IGARSS '00, 24-28 July 2000, Honolulu, USA. C D - R O M . Dutra, L. V., Hernandez F., P., Mazzocato, M.E., de Souza, R.C.M. , Oliver, C. 1999. Land cover classification based on multi-date JERS-1 imagery as a basis for deforestation detection. Proceedings IGARSS '99, 28 June - 2 July, 1999, Hamburg, Germany. C D - R O M . Eastman, J . R. 2001. Idrisi 32 Release 2. Guide to GIS and image processing, Vol. 2. Clark University, Worcester, MA. 151 pp. European Space Agency. 2002. The A S A R user guide, Issue 1.1. http://envisat.esa.int/dataproducts/asar/toc.htm Fjortoft, R., Lopes, A., Adragna, F. 2000. Radiometric and spatial aspects of speckle filtering. Proceedings IGARSS '00, 24-28 July 2000, Honolulu, USA. C D - R O M Getis, A., Ord, J.K. 1992. The analysis of spatial association by use of distance statistics. Geographical Analysis, 24, 189-206. Goodchild, M. 1986. Spatial autocorrelation. Concepts and Techniques in Modem Geography, 47, Geo Books, Norwich. Gong, P., Marceau, D.J., Howarth, P .J . 1992. A comparison of spatial feature extraction algorithms for land-use classification with S P O T HRV data. Remote Sensing of Environment, 40, 137-151. Grover, K., Quegan, S., da Costa Freitas, C. 1999. Quantitative estimation of tropical forest cover by SAR. IEEE Trans. Geoscience and Remote Sensing, 37, 479-490. Haralick, R.M., Shanmugan, K., Dinstein, I. 1973. Texture features for image classification. IEEE Trans. Syst. Man. Cybernet, 3, 610-621. 110 Henderson, F. M., Lewis, A . J . , Eds. 1998. Principles & applications of imaging radar. Manual of Remote Sensing, Third Edition, Volume 2. John Wiley & Sons, New York. 866 pp. Hoekman, D. H. 1991. Speckle ensemble statistics of| logarithmically scaled data. IEEE Trans. Geoscience and Remote Sensing, 29, 180-182. Holden, H., Derksen, C , LeDrew, E. 2000. Coral reef ecosystem change detection based on spatial autocorrelation of multispectral satellite data. Proceedings ACRS 2000, 4 - 8 December, 2000, Taipei, Taiwan. http://www.gisdevelopment.net/aars/acrs/2000/ts3/cost006.shtml Jensen, J . R. 1996. Introductory digital image processing. 2 n d Edition. Prentice Hall, Upper Saddle River, N.J. 318 pp. Khorram, S., Ed. 1999. Accuracy assessment of remote sensing-derived change detection. American Society of Photogrammetry and Remote Sensing, Bethesda, Maryland. 64 pp. Kurvonen, L., Hallikainen, M.T. 1999. Textural information of multitemporal ERS-1 and JERS-1 S A R images with applications to land and forest type classification in boreal zone. IEEE Trans. Geoscience and Remote Sensing, 37, 680-689. Klinka, K., Pojar, J . , Meidinger, D.V. 1991. Revision of biogeoclimatic units of coastal British Columbia. Northwest Science, 65, 32-47. Krajina, V. 1965. Biogeoclimatic zones in British Columbia. Ecology of Western North America, 1, 1-17. Kux, H.J.H., dos Santos, J.O., Lacruz, M.S.P., Ahern, F.J., Pietsch, R.W. 1998. Capabilities of R A D A R S A T data for land use/land cover analysis in S W Amazonia (test site Acre, Brazil). In F.J. Ahern (Ed.), RADARSAT for Amazonia: Results of ProRADAR Investigations (pp. 35-44), Canada Centre for Remote Sensing, Ottawa. Leckie, D.G. 1997. The effect of environmental influences on R A D A R S A T images of forest land. Proceedings Geomatics in the Era of RADARSAT (GER '97), 25-30 May 1997, Ottawa, Canada. C D - R O M . Leckie, D.G., Hill, D.A., Yatabe, S.M., Copis, P.L., D'eon, S.P. , Robinson, C.F., Banner, A., Landry, R. 1998. Temporal dynamics of R A D A R S A T imagery over a forest site. Proceedings ADRO Final Symposium, 13-15 October, 1998, Montreal, Canada. C D - R O M . Leckie, D. G. 1998. Forestry applications using imaging radar. In F.M. Henderson & A . J . Lewis (Ed.), Principles and Applications of Imaging Radar (pp. 435-509), John Wiley & Sons, New York. 111 LeDrew, E.F., Wulder, M. and Holden, H. 2000. Change detection of satellite imagery for reconnaissance of stressed tropical corals. Proceedings IGARSS '00, 24-28 July 2000, Honolulu, USA. C D - R O M . Lopes, A., Nezry, E., Touzi, R., Laur, H. 1993. Structure detection and statistical adaptive speckle filtering in S A R images. Int. J. Remote Sensing, 14, 1735-1758. Luscombe, A. P. 2003. Radiometric stability of RADARSAT-1 images of the Amazon calibration site. Proceedings ASAR 2003 Workshop, 25-27 June 2003, St. Hubert, Canada. C D - R O M . v Mitchell, S . J . 1995. A synopsis of windthrow in British Columbia: occurrence, implications, assessment and management. In M. Coutts and J . Grace (Ed.) Wind and Trees, (pp. 448-459), Cambridge University Press. Mitchell, S .J . , Hailemariam, T., Kulis, Y. 2001. Empirical modelling of cutblock edge windthrow risk on Vancouver Island, Canada, using stand level information. Forest Ecology and Management, 154, 117-130. Murtha, P.A. 1997. RAdar Imaging Natural Systems (RAINS): stumps and dead trees on RADARSAT. Proceedings Geomatics in the Era of RADARSAT (GER '97), 25-30 May 1997, Ottawa, Canada. C D - R O M Murtha, P.A. 1998a. Monitoring riparian leave strips in forest clearcuts with multi-temporal R A D A R S A T Fine 2 Mode image data. Proceedings ADRO Final Symposium, 13-15 October, 1998, Montreal, Canada. C D - R O M Murtha, P.A. 1998b. Weather or not: monitoring riparian strips in clearcuts with multi-temporal Fine 2 Mode R A D A R S A T data. Proceedings 20th Canadian Symposium on Remote Sensing, 11-14 May, 1998, Calgary, Canada. Canadian Aeronautics and Space Institute, Ottawa, 269 pp. Murtha, P.A. 1998c. Interesting images from the RAINS study, northern Vancouver Island. Proceedings ADRO Final Symposium, 13-15 October, 1998, Montreal, Canada. C D - R O M Murtha, P.A., Mitchell, S . J . 1998. Riparian leave strip monitoring with multi-temporal Fine 2 Mode R A D A R S A T data. Proceedings ASPRS Annual Meeting, 1-3 April, 1998, Bethesda, USA. C D - R O M Murtha, P.A. 2000a. Monitoring riparian leave strips with multi-temporal R A D A R S A T C-band satellite data. In L. Darling (Ed.), At Risk: Proceedings of a conference on the biology and management of species and habitats at risk (pp. 553-559). B.C. Ministry of Environment Lands and Parks, Victoria, Canada. 112 Murtha, P.A. 2000b. Monitoring cut-blocks, riparian strips and windthrow on northern Vancouver Island with R A D A R S A T F2 data. Proceedings BC RADARSAT Workshop, 9 May 2000, Victoria, Canada, (http://www.space.qc.ca/asc/pdf/bc14.pdf) Murtha, P.A. 2000c. Surficial geology and climatic effects on forest clearcut tone in R A D A R S A T images of northern Vancouver Island. Canadian Journal of Remote Sensing, 26, 253-262. Nezry, E., Leysen, M., De Grandi, G. 1995. Speckle and scene spatial estimators for SAR image filtering and texture analysis: some applications to agriculture, forestry and point target detection. Proceedings of SPIE, 2584, 110-120. Oliver, C , Quegan, S. 1998. Understanding synthetic aperture radar images. Artech House, 479 pp. Ord, J.K., Getis, A. 1995. Local spatial autocorrelation statistics: distributional issues and an application. Geographical Analysis, 27, 286-306. PCI Geomatics. 2000. P A C E radar analysis, version 7.0. PCI Geomatics Group, Richmond Hill, Ontario. 170 pp. R A D A R S A T International, n.d. Innovation in earth observation: RADARSAT-2 overview. http://www. radarsat2. info/rs2_satellite/overview. asp R A D A R S A T International. 1995. R A D A R S A T illuminated: your guide to products and services. R A D A R S A T International, Richmond, Canada. ftp://ftp.rsi.ca/USER/USGUIDE.PDF R A D A R S A T International. 2000. R A D A R S A T data products specifications. RSI -GS-026, Vers. 3/0, May 8, 2000, R A D A R S A T International, Richmond, Canada. 126 pp. ftp ://ftp. rsi. ca/ceos_dfn/d4_3-0 .zip Raney, R. K. 1998. Radar fundamentals: technical perspective. In F.M. Henderson & A . J . Lewis (Ed.), Principles and Applications of Imaging Radar (pp. 435-509), John Wiley & Sons, New York. Raney, R.K., Freeman, T., Hawkins, R.W., Bamler, R. 1994. A plea for radar brightness. Proceedings IGARSS '94, 8- 12 August, 1994, Pasadena, USA, 2, 1090-1092. Rignot, E. J . M., van Zyl, J . J . 1993. Change detection techniques for ERS-1 S A R data. IEEE Trans. Geoscience and Remote Sensing, 31, 896-906. Rignot, E., Chellappa, R. 1992. A Bayes classifier for change detection in synthetic aperture radar imagery. Proceedings Int. Conf. on Acoustics, Speech and Signal Processing, 23 - 26 March 1992, San Francisco, USA, 3, 25-28. 113 Ruecker, G. , Siegert, F. 2000. Burn scar mapping and fire damage assessment using ERS-2 S A R images in East Kalimantan, Indonesia. Proceedings IAPRS 2000, 17-21 August 2000, Amsterdam, The Netherlands. http://www.rssgmbh.de/pdf/ISPRS-2000.pdf Siegert, F., Rucker, G . 1999. Evaluation of the 1998 forest fires in East-Kalimantan (Indonesia) using multitemporal ERS-2 S A R images. Earth Observation Quarterly, 61, 7-12. Smits, P.C., Annoni, A. 2000. Toward specification-driven change detection. IEEE Trans. Geoscience and Remote Sensing, 38, 1484-1488. Srivastava, S.K., Le Dantec, P., Gray, R., Hawkins, R.K., Murnaghan, K. 2003. RADARSAT-1 calibration performance maintained in extended mission. Proceedings ASAR 2003 Workshop, 25-27 June 2003, St. Hubert, Canada. C D - R O M . Touzi, R., Lopez, A., Bruniquel, J . , Vachon, P.W. 1999. Coherence estimation for S A R imagery. IEEE Trans. Geoscience and Remote Sensing, 37, 135-149. Villasenor, J . D., Flatland, D.R., Hinzman, L.D. 1993. Change detection on Alaska's North Slope using repeat-pass ERS-1 S A R Images. IEEE Trans. Geoscience and Remote Sensing, 31, 229-236. Walpole, R. E. 1974. Introduction to statistics. 2 n d Edition. MacMillan Publishing, New York. 340 pp. Wells, M.L., Getis, A. 1999. Spatial characteristics of stand structure in Pinus torreyana. Plant Ecology, 143, 153-170. Weydahl, D. J . 1991. Change detection in S A R images. Proceedings IGARSS '91, 3-6 June 1991, Helsinki, Finland, 3, 1421-1424. Weydahl, D. J . 1992. Temporal change detection in ERS-1 S A R images. Proceedings IGARSS '92, 26-29 May 1992, Houston, USA, 2, 1346-1348. White, R.G. 1991. Change detection in S A R imagery. Int. J. Remote Sensing, 12, 339-360. Wulder, M., Boots, B. 1998. Local spatial autocorrelation characteristics of remotely sensed imagery assessed with the Getis statistic, Int. J. Remote Sensing, 19, 2223-2231. Wulder, M., Boots, B. 2001. Local spatial autocorrelation characteristics of Landsat TM imagery of a managed forest area. Canadian Journal of Remote Sensing, 27, 67-75. 114 Appendix A 115 Table A-1. Feature-based threshold results for datasets processed by the C H D E T algorithm. Feature Dataset Feature Imaging Data Arithmetic Thresholding Proportion of Proportion of Size Size Dates for Format Operation Applied Pixels in Feature Pixels in Feature Cutblock (pixels) (pixels) Arithmetic After Arithmetic Exceeding Positive Exceeding Positive Operations Operation Threshold (%) and Negative Threshold (%) Cutblock S501 (Small dataset) 65,536 18,930 1998, 1997 dB Filtered log ratio 1.65o (PFA 5%) 3.0 4.0 Cutblock S501 (Large dataset) 1,048,576 18,930 1998, 1997 dB Filtered log ratio 1.65a (PFA 5%) 3.0 4.0 Cutblock 594 1,048,576 46,971 1999, 1996 dB Filtered log ratio 1.65a (PFA 5%) 2.7 3.9 Cutblock M5307 1,048,576 40,698 1998, 1996 dB Filtered log ratio 1.65a (PFA 5%) 5.3 5.7 Cutblock C99 1,048,576 32,889 1999, 1997 dB Filtered log ratio 1.65a (PFA 5%) 3.0 5.1 116 Table A-2. Feature-based threshold results and summary statistics for Getis format data. Feature Dataset Feature Imaging Data Arithmetic Thresholding Proportion of Following Mean Standard Resulting Size Size Dates for Format Operation Applied Pixels in Feature Statistics Deviation (a) Dataset Cutblock (pixels) (pixels) Arithmetic After Arithmetic Exceeding Positive Apply To Normally Operations Operation Threshold (%) Distributed? Cutblock S501 (Small dataset) 65,536 18,930 1998, 1997 3x3 Getis Difference 95% 11.1 Whole Dataset 0.003 0.991 yes Feature 0.341 1.110 yes 11x11 Getis Difference 95% 15.8 Whole Dataset 0.013 0.480 no Feature 0.325 0.600 no MaxGetis Difference 95% 12.1 Whole Dataset 0.007 1.031 no Feature 0.423 1.163 no Cutblock S501 (Large dataset) 1,048,576 18,930 1998, 1997 3x3 Getis Difference 95% 16.7 Whole Dataset 0.000 0.913 yes Feature 0.442 1.151 almost 11x11 Getis Difference 95% 42.1 Whole Dataset 0.000 0.348 yes Feature 0.425 0.623 no MaxGetis Difference 95% 20.4 Whole Dataset -0.001 0.936 no Feature 0.548 1.202 no Cutblock 594 1,048,576 46,971 1999, 1996 3x3 Getis Difference 95% 11.3 Whole Dataset -0.0011 0.896 no Feature 0.251 1.050 almost 11x11 Getis Difference 95% 24.0 Whole Dataset -0.003 0.391 no Feature 0.248 0.603 no MaxGetis Difference 95% 13.3 Whole Dataset -0.007 0.909 no Feature 0.326 1.107 no Cutblock M5307 1,048,576 40,698 1998, 1996 3x3 Getis Difference 95% 19.0 Whole Dataset 0.000J 0.911 no Feature 0.541 1.162 almost 11x11 Getis Difference 95% 44.6 Whole Dataset 0.001 0.380 yes Feature 0.531 0.601 no MaxGetis Difference 95% 23.3 Whole Dataset 0.003 0.917 no Feature 0.679 1.197 no Cutblock C99 1,048,576 32,889 1999, 1997 3x3 Getis Difference 95% 14.9 Whole Dataset 0.000 0.874 no Feature 0.420 1.050 yes 11x11 Getis Difference 95% 37.6 Whole Dataset 0.001 0.360 almost Feature 0.411 0.566 no MaxGetis Difference 95% 18.2 Whole Dataset 0.000 0.881 no Feature 0.524 1.092 no 117 Table A-3. Proportion of pixels in riparian zones and blowdown areas indicating significant change on MaxGetis difference images. Feature Riparian Zone or B l o w d o w n Area Dataset Size (Pixels) Feature Size (Pixels) Imaging Dates for MaxGetis Difference Image T h r e s h o l d i n g A p p l i e d T o Difference Image Proportion of Pixels in Feature Exceeding Positive Threshold (%) Fol lowing Statistics Apply T o Mean Standard Deviation (a) Strip 31a 1,048,576 3160 1999, 1996 95% 12.1 Whole Image -0.001 0.909 Feature 0.365 1.032 Strip 16(S505) 1,048,576 922 1999, 1996 95% 8.9 Whole Image -0.001 0.943 Feature 0.051 1.077 Edge Blowdown 1,048,576 1731 1999, 1996 95% 23.6 Whole Image -0.001 0.943 Feature 0.827 1.023 Edge Blowdown (M5307) 1,048,576 1649 1999, 1997 95% 3.3 Whole Image 0.000 0.881 Feature 0.100 0.763 Strip 17 (S505) 1,048,576 986 1999, 1996 95% 14.1 Whole Image -0.001 0.943 Feature 0.324 1.078 Strip 18 (S505) 1,048,576 888 1999, 1996 95% 7.5 Whole Image -0.001 0.943 Feature 0.367 0.920 Strip 13 1,048,576 1106 1999, 1996 95% 8.7 Whole Image Feature -0.001 0.278 0.953 0.950 Strip 25 1,048,576 4515 1999, 1996 95% 6.2 Whole Image -0.001 0.953 Feature -0.006 0.998 Strip 26 1,048,576 1029 1999, 1996 95% 7.8 Whole Image -0.001 0.953 Feature -0.049 1.260 Strip 28 1,048,576 2733 1999, 1996 95% 5.7 Whole Image -0.001 0.953 Feature 0.139 0.918 1 18 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0075005/manifest

Comment

Related Items