AUTOMATIC CALCULATION OF OIL SLICK AREA FROM MULTIPLE SAR ACQUISITIONS FOR DEEPWATER HORIZON OIL SPILL

The Deepwater Horizon oil spill occurred in the Gulf of Mexic o in April 2010 and became the largest accidental marine oil s pill in history. Oil leaked continuously between April 20th and Jul y 15th of 2010, releasing about 780, 000m of crude oil into the Gulf of Mexico. The oil spill caused extensive economical and ecolo gical damage to the areas it reached, affecting the marine an d wildlife habitats along with fishing and tourism industries. For oil spill mitigation efforts, it is important to determi ne the areal extent, and most recent position of the contamin ated area. Satellitebased oil pollution monitoring systems are being used for mo nitoring and in hazard response efforts. Due to their high ac curacy, frequent acquisitions, large area coverage and day-and-night opera tion Synthetic Aperture Radar (SAR) satellites are a major c ontributer of monitoring marine environments for oil spill detection. We developed a new algorithm for determining the extent of th e oil spill from multiple SAR images, that are acquired with s hort temporal intervals using different sensors. Combining the multi-polarization data from Radarsat-2 (C-band), Envisa t ASAR (C-band) and Alos-PALSAR (L-band) sensors, we calculate the extent o f the oil spill with higher accuracy than what is possible fro m nly one image. Short temporal interval between acquisitions (hour s to days) allow us to eliminate artifacts and increase accur y. Our algorithm works automatically without any human interv ention to deliver products in a timely manner in time critica l operations. Acquisitions using different SAR sensors are radiometrica lly calibrated and processed individually to obtain oil spi ll area extent. Furthermore the algorithm provides probability maps of the ar as that are classified as oil slick. This probability info rmation is then combined with other acquisitions to estimate the combined p robability map for the spill.


INTRODUCTION
Marine oil spills are a common threat to all sea bordering countries. The environmental and related economical losses due to an oil spill can be large, and many methods have been developed to monitor oceans in operational conditions. In most cases it is critical to respond to spills in a timely manner, and therefore it is important for such operational systems to provide results in near real time. Once a spill is detected, its behavior can be modeled based on physical models. Synthetic Aperture Radar (SAR) systems provide a viable option of oil slick monitoring. SAR intensity images are sensitive to surface roughness which is altered in the case of an oil spill. The scattering of oil-free ocean surface is dominated by Bragg scattering (Bragg, 1913). A thin oil sheen covering the ocean will reduce the ocean-atmosphere interaction, and alter the smoothness of the surface. Therefore, oil slicks appear slightly darker in moderate wind conditions in SAR imagery. SAR systems have an ideal range of wind speed and direction where they are most sensitive, and do not perform well for oil spill detection under very windy or very calm conditions.
There are of course other methods to monitor oceans and detect oil spills, such as optical and infra-red remote sensing, and hyperspectral imaging (Brekke and Solberg, 2005). No matter how complex or advanced a method is single-handedly all have different levels of uncertainties. Therefore it is important to combine many observations complementing each other's weaknesses providing the best possible information. Results from different remote sensing sensors, as well as field observations can be combined, and statistically provide an outcome better than any individual resource.
In this paper we develop a fully automatic oil spill monitoring systems that is capable of combining data from multiple SAR sources. The focus of this paper is on radar imagery, however the algorithm can be then be expanded to optical imagery and ground measurements as discussed later. A brief background on the technical aspects is provided in Section 2. Results obtained from SAR data over the Deepwater Horizon oil spill are discussed in Section 3.

BACKGROUND
Synthetic Aperture Radar imagery is sensitive to surface roughness which is altered in case of an oil spill (Alpers and Hühnerfuss, 1988). Oil slicks change the smoothness of ocean surface and appear darker compared to surrounding oil free ocean, however the amount of damping is affected with wind and wave conditions. Furthermore the speckle effect in SAR imagery limits the reliability of point measurements in the image, causing spurious results (Brekke and Solberg, 2005). In our approach we apply a multiple step processing to limit the adverse effects of speckle, instead of filtering the data with a speckle filter.
In this paper, our focus is not on developing a method that performs better than any one of the methods mentioned earlier, but instead to develop a practical framework where results from different methods and sources can be combined to provide a joint solution, which is likely to have less uncertainty. The method can be easily extended to include any geospatial observation (e.g. optical remote sensing, ground measurements), however in this paper our focus is on SAR imagery. The algorithm is composed of three calculation steps: (1) Pixel (point) probability, (2) spatial probability, and (3) spatio-temporal probability. Point probability is calculated based on normalized radar cross section, where darker pixels get higher probabilities for oil contamination (Barni et al., 1995). Spatial probability is based on the damping ratio given the current wind conditions and imaging parameters (Gade et al., 1998, Kim et al., 2010. Point and spatial analysis results are then combined to provide a joint probability for oil slick at each acquisition. The spatio-temporal probability is estimated from multiple SAR acquisitions separated shortly in time, resulting in a time varying probability of oil slick over target area. All analysis in this paper is done over intensity calibrated, geocoded SAR imagery. Images are calibrated to normalized radar cross section (NRCS). The imagery is resampled to a common geometry using a sinc interpolator, and land-masked using Global, Selfconsistent, Hierarchical, High-Resolution Shoreline data (GSHHS) (Wessel and Smith, 1996).
First and second steps of the algorithm are iterative, and are performed at the same time. The first step of the algorithm is a simple dark-object selection routine based on intensity thresholding. At this step, dark areas of the image are assigned a higher probability for an oil spill. The initial threshold for the first step is the noise equivalent sigma zero (NESZ), which is the sensors noise floor. A probability value for each pixel is assigned based on it's intensity such that: where P (W |σ0) is probability of oil-free water given the NRCS, T is the threshold, and P (O|σ0) is the probability of oil given the NRCS is the complement of P (W |σ0). The P (W |σ0) is modified by bringing all larger values to 1, constraining the probability values between 0 and 1. Furthermore, iterations start using a multilooked imagery, to reduce the effect of speckle noise. It is worth noting that the multilooking operation changes the dynamic range of the radar imagery. In order to keep the thresholds equivalent at each iteration, multilooked images are further calibrated to have the same dynamic range as the full resolution image.
The second step of the algorithm starts with calculation of damping factor after Gade et al. (Gade et al., 1998, Kim et al., 2010. The estimated damping factors at different wind speeds as a function of Bragg wavenumbers are shown in Figure 1. Figure 1 also shows the maximum theoretical observation range for Radarsat-2, Envisat-ASAR, and Alos-PALSAR SAR systems. The Bragg wavenumber is defined as (Gade et al., 1998): where kB is the Bragg wavenumber, k0 is the radar wavenum-ber, and φ is the incidence angle. Figure 1 shows the maximum expected damping amount which would be observed when the angle between the wind, and radar wave is zero. ALOS-PALSAR is not expected to be effective in mild wind conditions, as shown in the figure. The SAR sensors on-board Radarsat-2 and Envisat are both C-Band, and therefore cover similar regions in the wavenumber domain. The slightly larger coverage of Radarsat-2 is due to it's larger range of incidence angles. Wind is an important factor to estimate the damping factor, which is also dependent on the relative angle between radar wave and wind direction. Figure 2 shows a plot of expected damping factors at speeds between 2.5 m/s and 12.5 m/s. The relative wind direction is plotted at counter-clockwise increasing angles, starting from zero at the horizontal axis. In this study we use wind speed, direction and ocean wave group velocity data from National Data Buoy Center station 42040, located at 29.122N, 88.207W, about 40km NNW of the Deepwater Horizon oil rig. The first and second steps of the algorithm are run iteratively at five different spatial resolutions. The analysis is performed in a pyramid structure with five different stages. At each stage the image is multilooked to the power of two, such that at the first stage International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XXXIX-B7, 2012 XXII ISPRS Congress, 25 August -01 September 2012, Melbourne, Australia the image is multilooked to 2 5 , and to 2 4 at the second stage. Both the intensity thresholding and damping factor methods provide results at each stage which are combined using: where JP (O) is the joint probability for oil, P (O|n), is the probability of oil given the NRCS, P (O|d) is the probability of oil given the damping factor analysis, P (W |n) is the probability of clear water given the NRCS, and P (W |d) is the probability of clear water given the damping factor analysis. The main reason for combining the two probability functions is that the damping factor analysis will return low probabilities for the center of large slicks. This is due to the fact that the damping factor is calculated against a moving window average. For oil slicks larger than the size of the moving window, the average and pixel values will be very close, returning no dampening. Of course a larger window size can be used to eliminate this problem, however, since that would require human intervention. Contrary to the damping factor analysis, the NRCS based thresholding algorithm will return high probabilities for the center of the slick.
The study area is shown in Figure 3, and is about 350km×350km, centered around the Deepwater Horizon oil spill. The locations of the ALOS-PALSAR and Radarsat-2 imagery are shown in the figure. It should be noted that even though the SAR imagery is calibrated to NRCS, there is still a gradual change in intensity along the range direction ( Figure 3). The Gulf of Mexico oil spill provides a great test case for the new algorithm, because it is localized and continuous over time. Furthermore there are many published research and ground observations available for validating the method. Figure 3: The black-box shows our designated study area. GSHHS shoreline is shown in light green. The SAR intensity images are from Radarsat-2, and the one on top is acquired on April 24th 2010. Light green boxes show the footprints for the ALOS PALSAR imagery. The legend shows distance in degrees.

RESULTS AND DISCUSSION
SAR images acquired from Radarsat-2 and Alos were processed using the proposed method. Some of the imaging parameters and environmental conditions are summarized in Table 1. Table  shows the observed wind-speed, relative angle between the wind and radar wave, and the calculated damping factor (D.F.). The damping factor calculation also takes into account the wave group velocity, and radar incidence angle which are not listed in the table. The damping factors listed for ALOS-PALSAR are rather low, however they are still above the NESZ of the instrument, which is about −29dB for the fine beam single (FBS) imaging mode. The results of oil spill detection algorithm is shown in Figure 4. The probability maps calculated for the five iterative steps using two different methods and their combination are presented.
The results are shown in ascending order of resolution, where the 2 5 multilooked image is located at the left hand side. The final solution for the processed imagery is shown in the bottom right corner. The joint probability (JP (O)) is calculated using the NRCS based oil probability (P (O|n)) and damping factor based oil probability (P (O|d)) as shown in (4). Oil Slick Probability 0 0.5 1 Figure 4: Results for the first (point probability) and second (spatial probability) processing steps, for the Radarsat-2 data acquired on April 27th, 2010.
The complementing behavior of the two methods can be seen in Figure 4. While the P (O|n) has very little noise at high multilooking, the opposite is true for P (O|d), which becomes less and less noisy with decreasing multilooking. Furthermore, the void in the center of the oil spill is visible in P (O|d) results for level five.
Final results of all the images processed in this study are shown in Figure 5. The analysis using Radarsat-2 imagery obtained better results compared to the Alos imagery. This is very likely due to the small damping factors that are obtained at L-band, as shown in Figure 1. The current algorithm does not employ any weighting to the data, therefore the combined probability of all observations are inconclusive. This can be improved however, by implementing a more complex filter to the final stage, such as a Kalman filter, or by simply adding more C-band data to dominate the results. We recently acquired additional imagery from Envisat to test our hypothesis. Utilizing a Kalman filter at the final step of the algorithm will allow for utilizing a larger spectrum of methods and data sets, which may only be useful under certain conditions. It is also worth noting that the Alos-PALSAR imagery, acquired on May 1st, 2010 at 04:10 UTC shows almost no brightness variation.

CONCLUSIONS AND FUTURE WORK
More and more remote sensing data is becoming available every day, some of which are applicable to the the detection of oil spills. In this paper, we presented a method to combine SAR data from multiple sensors, acquired at different times from different geometries, to obtain a combined oil slick probability map for the Deepwater Horizon oil spill. Furthermore, it is also possible to combine results from different algorithms applied to the same data set, to increase accuracy as presented in this paper.
Our fully automatic algorithm can be expanded to include other sensors (e.g. optical, infra-red, hyperspectral), and can be customized to the operational needs. The algorithm can operate with any sensor and any algorithm that can provide a probability map. The intuitive probability maps assist operators and ground personnel in determining which areas to prioritize. Expansion of the algorithm to utilize MODIS imagery, and a neural-network based SAR oil spill algorithm for a future study is currently under consideration.