Implementación de algoritmos para la monitorización de calidad de aguas y de espacios costeros mediante la utilización de imágenes de teledetección satelital de muy alta resolución

   UniversidaddeLasPalmasdeGranCanaria DepartamentodeSeñalesyComunicaciones  ProgramadeDoctorado CibernéticayTelecomunicación  TesisDoctoral Implementación de algoritmos para la monitorización de calidad de aguas y de espacios costeros mediante la utilización de imágenes de teledetección satelital de muy alta resolución AUTOR: Javier Martín Abasolo DIRECTORES: Dr. Francisco Eugenio González Dr. Javier Marcello Ruiz El Director El Codirector El Doctorando Las Palmas de Gran Canaria, Marzo de 2016 Dedicatoria  A mi familia y en especial a mi pareja Anabella, por todo su apoyo incondicional. Agradecimientos Este agradecimiento está dirigido a todas aquellas personas que han hecho posible la realización de esta Tesis Doctoral. En primer lugar a mis tutores Francisco Eugenio González y Javier Marcello Ruiz, gracias por todo vuestro apoyo, conocimientos y esfuerzo derrochados durante este tiempo. Muchas gracias. También a todos los compañeros con los que he tenido el honor de trabajar durante este tiempo en el Grupo de Procesado de Imágenes y Teledetección (GPIT). Toda mi gratitud a los responsables de los diferentes proyectos de investigación en que se basa esta tesis doctoral, por todos los recursos materiales y conocimientos proporcionados. Especial agradecimiento al Observatorio Ambiental de Granadilla (OAG), y a su director Dr Antonio Machado Carrillo, por ser el precursor de esta investigación, al apostar por la teledetección de alta resolución en la monitorización de los espacios costeros. Este trabajo ha sido apoyado por los siguientes proyectos: 1. OAG. Elaboración de algoritmos específicos para interpretar la turbidez y clorofila en imágenes del satélite Worldwiew2 en el litoral de Granadilla, en Tenerife (CN-47/11142001). 2. VULCANO. Volcanic erUption at El Hierro IsLand. Sensitivity and ReCovery of the mAriNe EcOsystem (CTM2012-36317). 3. TELECAN. Programa para el Desarrollo de Redes Tecnológicas y de Aplicación de Datos de Teledetección en África Occidental (MAC/3/C181). 4. TECHMARAT. Tecnologías de Vegetales Marinos para la Región Atlántica (0111_TECHMARAT_2_A). 5. ARTeMISat. Desarrollo de técnicas avanzadas de procesado de imágenes de satélite de alta resolución para la gestión sostenible de los recursos naturales marinos y terrestres (CGL2013-46674-R).  VI Figura 63. (a) Imagen del volcán submarino. (b) Localización de las diferentes tomas in-situ realizadas durante la erupción submarina. (Fuente: IEO, proyecto Vulcano). .............................. 117 Figura 64. (a) Composición RGB de las imágenes MODIS, y (b) MERIS del mismo día (9 noviembre 2011). (c) Producto kd(490) para el sensor MODIS, y (d) MERIS. ............................... 117 Figura 65. Composición RGB de la imagen WV2 del 27 de octubre de 2011 (a), concentración de kd(490) obtenida mediante el algoritmo operativo ecuación (96). .................................................. 118 Figura 66. Resultados del algoritmo de calidad de aguas para la imagen WV2 del 27 de octubre de 2011. (a) Atenuación debida a la concentración de clorofila, (b) atenuación debida a la materia disuelta, (c) difusión debida a la materia suspendida, y (d) atenuación difusa kd(490). ................ 119 Figura 67. (a) Transecto costero utilizado en el ajuste. (b) Ajuste del ratio algorithm mediante regresión con la batimetría in-situ: Ajuste lineal (línea azul) y ajuste cuadrático (curva roja). ...... 121 Figura 68. Resultado del ratio algorithm: (a) ajuste lineal, (b) ajuste cuadrático, (c) batimetría sonar. .............................................................................................................................................. 121 Figura 69. Resultados de batimetría para la costa de Granadilla. (a) Imagen RGB de la zona de estudio, (b) gráfica de dispersión de datos satélite-sonar, (c) mapa de batimetría sonar, y (d) mapa de batimetría satélite. ..................................................................................................................... 123 Figura 70. Resultado del algoritmo RTM de batimetría para la costa de Maspalomas: (a) Batimetría sonar de la zona de estudio del año 2003, (b) composición RGB del 20 de noviembre de 2011, (c) resultado del algoritmo de batimetría para el 20 de noviembre de 2011, (d) composición RGB del área de estudio para el 17 de enero de 2013, (e) resultado del algoritmo de batimetría para el 17 de enero de 2013. ........................................................................................ 124 Figura 71 Resultados de batimetría para la costa de Corralejo: (a) Imagen RGB de la zona de estudio, (b) gráfica de dispersión de los datos satélite-sonar, (c) mapa de batimetría sonar, y (d) mapa de batimetría satélite. ........................................................................................................... 126 Figura 72. Reflectividad normalizada de las clases puras más usuales en los fondos costeros (arena, algas, y sedimentos) [103]. ................................................................................................ 126 Figura 73. Resultados del albedo submarino para la costa de Granadilla. (a) Composición RGB del albedo submarino, (b) abundancia de la arena, (c) abundancia de las algas, (d) abundancia del sedimento. ...................................................................................................................................... 128 Figura 74. Resultados del albedo costero de Corralejo-Lobos: (a) Composición RGB del albedo costero, (b) abundancia de la arena, (c) abundancia de las algas, (d) abundancia del sedimentorocas. .............................................................................................................................................. 130 Figura 75. (a) Clasificación supervisada de comunidades bénticas de alta resolución mediante datos de albedo de fondo de imágenes WV2. (b) Clasificación de comunidades bénticas CIMA 2008. ............................................................................................................................................... 138 VII ÍNDICE DE TABLAS Tabla 1. Bandas de paso del satélite WorlView-2 en micrómetros. ................................................ 24 Tabla 2. Ancho de banda efectivo de las bandas del WorldView-2. ............................................... 25 Tabla 3. Irradiancia solar espectral promedio para las bandas del WorldView-2. .......................... 26 Tabla 4. Principales métodos de eliminación del brillo solar. .......................................................... 55 Tabla 5. Resumen de los resultados de la regresión lineal de las bandas WV2. ........................... 65 Tabla 6. Resultados de la irradiancia normalizada directa y difusa del modelo 6S. ....................... 73 Tabla 7. Coeficientes de absorción del fitoplancton. (Fuente: Lee et al. 1998). ............................. 87 Tabla 8. Datos in-situ de turbidez obtenidos para la imagen de granadilla del 1/12/2011 ............ 110 9 Capítulo 1 1. Introducción 1.1. Antecedentes La necesidad de realizar una monitorización de los entornos costeros es una labor de gran interés medioambiental, que requiere de un esfuerzo continuo, debido al impacto que supone la actividad humana en estas zonas, donde se concentra la mayor diversidad marina. El estudio de los entornos marinos mediante imágenes multiespectrales de teledetección espacial se ha enfocado históricamente en entornos oceánicos a una escala global y en zonas de aguas profundas, utilizándose sensores de muy baja resolución espacial, pero con una gran resolución temporal, como son, por ejemplo, el CZCS, AVHRR y SeaWiFS. El estudio de estas aguas, denominadas Tipo I, requiere de un procedimiento bien conocido gracias a décadas de experiencia en la parametrización de sus componentes, dominados por el fitoplancton, detectable a través de su concentración de clorofila (chl-a). Con el desarrollo tecnológico y el lanzamiento de sensores con mayor resolución espectral, como, por ejemplo MODIS y MERIS, ha sido posible el estudio de entornos más cercanos a la costa, denominadas aguas de Tipo II, y en ellas se concentran distintos tipos de materias en suspensión y disolución. En estos entornos más complejos, se ha requerido aplicar importantes mejoras en los algoritmos que permitan modelar la presencia de estos componentes, logrando obtener parámetros de calidad del agua y su eutrofización. Sin embargo, la baja resolución espacial que proporcionan estos satélites, de cientos de metros o entorno al kilómetro, limitan la monitorización de las aguas más cercanas a la costa, así como la vigilancia de aguas interiores, como lagos y presas. Introducción 10 Gracias al avance de los satélites de media resolución como, por ejemplo, LANDSAT, con resoluciones cercanas a 30 metros, se ha iniciado un proceso de adaptación de los algoritmos desarrollados para sensores como el MODIS, con objeto de lograr una monitorización de costas a una escala más regional. Sin embargo, la resolución espacial aún es insuficiente para la monitorización precisa de zonas litorales. En esta última década, la teledetección de alta resolución espacial se ha consolidado como herramienta indispensable para la monitorización en múltiples aplicaciones terrestres y marinas a escala local. En el contexto de la monitorización costera, la alta resolución ha abierto un novedoso campo de aplicación que permite estudiar con gran detalle las primeras decenas de metros de la costa, en donde la biodiversidad alcanza su mayor cota. Además, con la incorporación de bandas multiespectrales de alta penetración en el agua, estos sensores permiten, por primera vez, la monitorización del fondo costero hasta profundidades cercanas a 30 metros. El modelado conjunto de los parámetros inherentes del agua y del albedo del fondo costero conllevan una gran complejidad, situando esta temática en el más puntero estado del arte. Así, se requiere del desarrollo e implementación de nuevos y sofisticados algoritmos para su resolución. A la gran complejidad del uso de datos de alta resolución, se une el hecho de que dichos satélites suelen disponer de pocos canales espectrales con una baja relación señal a ruido, debido a la gran absorción de la luz producida por el agua. Por este motivo, estas imágenes han de ser adquiridas en condiciones idóneas, desde el punto de vista de calidad atmosférica, marina y de oleaje, debiéndose aplicar complejos procesados para posibilitar su uso. El lanzamiento del satélite WorldView-2 (en lo sucesivo WV2), a finales del 2009, aportó un nuevo hito en el estado del arte de los satélites de muy alta resolución, proporcionando una resolución espacial de 0.5 m en la banda pancromática (PAN) y 2 m en sus ocho canales multiespectrales (MS). WV2 introduce una cantidad de canales inusualmente alta, para este tipo de satélites, entre ellos un azul de alta penetración que permite mejorar las capacidades de este satélite en la monitorización de aguas costeras. En esta misma línea, DigitalGlobe ha apostado por dar continuidad a estos servicios al lanzar el nuevo satélite WorldView-3 (WV3) en agosto de 2014. WV3 mantiene el mismo número de bandas en el rango óptico-NIR mejorando la resolución espacial (0.31 m PAN y 1.24 m en los canales MS) y añadiendo nuevas bandas en el IR cercano. Los entornos insulares como las Islas Canarias, contienen una gran riqueza y diversidad marina en sus costas, donde reside la mayor parte de su población, siendo además un reclamo para su importantísima industria turística. Por este motivo, se requiere de una monitorización exhaustiva de sus zonas litorales, en donde la teledetección de alta resolución es una herramienta fundamental. Estas necesidades de monitorización incluyen la obtención de parámetros de calidad de las aguas en entornos naturales y playas, la monitorización de la biodiversidad bentónica del lecho costero, así como la determinación de la variación de la batimetría en contextos como obras de infraestructuras portuarias. Asimismo, debido a la naturaleza volcánica de las islas, la utilización de imágenes de alta resolución tuvo gran impacto en la monitorización de la erupción submarina acaecida en la Isla de El Hierro en octubre de 2011. El interés del Grupo de Procesado de Imágenes y Teledetección (GPIT) de la Universidad de Las Palmas de Gran Canaria (ULPGC) y de organismos públicos de vigilancia y gestión costera ha quedado patente durante este tiempo con su participación en el desarrollo de diferentes proyectos y contratos de I+D+i europeos (TELECAN1, TECHMARAT2), nacionales (ARTEMISAT3, 1 Programa para el Desarrollo de Redes Tecnológicas y de Aplicación de Datos de Teledetección en África Occidental (MAC/3/C181). 2 Tecnologías de Vegetales Marinos para la Región Atlántica (0111_TECHMARAT_2_A). CAPÍTULO 1 11 VULCANO4) y autonómicos (Observatorio Ambiental Granadilla, OAG5) estableciendo las bases para la realización de esta Tesis Doctoral. Para poner en antecedente, se procederá a describir con detalle las actuaciones que se han ido llevando a cabo en materia de monitorización costera en diferentes áreas de Canarias, con gran interés medioambiental. La costa de Granadilla, situada en el sur-este de la isla de Tenerife, es un entorno rico en un tipo de fanerógama marino denominadas sebas (Cymodocea nodosa) que crean colonias llamadas sebadales. Este entorno de alto valor ecológico se vio comprometido por causa de la construcción del puerto de Granadilla, cuyas obras se iniciaron en el 2009. Para garantizar la conservación medioambiental se creó la Fundación Pública Observatorio Ambiental Granadilla (OAG), la cual se encarga de la monitorización de este entorno no solo mediante los métodos in-situ tradicionales, sino también mediante el procesado de imágenes WV2 de muy alta resolución. Ésta última tarea fue encargada al Grupo de Procesado de Imágenes y Teledetección (GPIT) del Instituto universitario de Oceanografía y Cambio Global (IOCAG), en el marco del Programa Europeo de Monitorización Ambiental de Granadilla, establecido en 2010, con el fin de garantizar la calidad ambiental adecuada, dentro y fuera del puerto, durante su construcción. Por otro lado, la erupción del volcán submarino cercano a la costa sur de la Isla de El Hierro proporcionó una oportunidad única para los grupos de investigación de las universidades Canarias, relacionadas con las ciencias marinas, de estudiar las variaciones de las propiedades físico-químicas, biológicas y geológicas a consecuencia de dicha erupción. Para ello, el Instituto Español de Oceanografía (IEO) coordinó un proyecto multidisciplinar, VULCANO, donde diversos grupos de investigación pudieron participar. En particular, el GPIT, llevó a cabo la tarea de monitorizar la variación de la mancha eruptiva, mediante imágenes de teledetección espacial procedentes de sensores de baja resolución, principalmente MODIS y MERIS. Asimismo, se utilizaron imágenes de muy alta resolución WV2 para generar con gran detalle espacial datos del contenido físico-químico en el epicentro de la mancha. La reserva natural de Maspalomas y Playa del inglés, ubicada en la costa sur de la Isla de Gran Canaria, alberga una laguna costera conocida como “la charca de Maspalomas” situada junto a un sistema de dunas móviles de gran interés medioambiental. Este entorno protegido mantiene un complejo equilibrio entre el agua marina y el agua aportada por la desembocadura del barranco, lo que hace que esta agua sea salobre. Una de las principales características de este ecosistema es su estacionalidad, lo cual se traduce en importantes variaciones en la salinidad y volumen de agua, así como en sus niveles de fitoplancton, materias disueltas y en suspensión. Por este motivo, desde el GPIT se promovió un proyecto de cooperación transnacional MAC llamado TELECAN. A partir del cual se realizó un estudio de calidad de aguas en este tipo de entorno costero, mediante el uso de imágenes de alta resolución WV2. Finalmente, la costa noreste de la Isla de Fuerteventura, concretamente el área Corralejo-Lobos, declarada en su conjunto ‘Reserva de la Biosfera’ por la UNESCO (mayo 2009) es otro entorno natural de gran valor ambiental debido a sus fondos marinos de baja profundidad. A través de la colaboración entre el GPIT y el Banco Español de Algas se logró la adquisición de imágenes de esta zona de aguas someras, logrando obtener el albedo en grandes áreas del fondo marino. 3 Desarrollo de técnicas avanzadas de procesado de imágenes de satélite de alta resolución para la gestión sostenible de los recursos naturales marinos y terrestres (CGL2013-46674-R). 4 Volcanic erUption at El Hierro IsLand. Sensitivity and ReCovery of the mAriNe EcOsystem (CTM2012-36317). 5 Elaboración de algoritmos específicos para interpretar la turbidez y clorofila en imágenes del satélite Worldwiew2 en el litoral de Granadilla, en Tenerife (CN-47/11-142001). Introducción 12 Recientemente y en el marco del proyecto del Plan Nacional ARTEMISAT (2014-2016), coordinado por el GPIT, se está continuando con el desarrollo y la optimización de los diversos algoritmos implementados. Estas áreas estratégicas, ver Figura 1, de elevado interés para la conservación costera del Archipiélago Canario y descritas en las líneas de actuación anteriores, han proporcionado el marco para el desarrollado esta Tesis Doctoral. (a) (b) (c) (d) (e) Figura 1. (a) Localización de la zona de estudio (Islas Canarias). (b) Área de Granadilla. (c) Área de Maspalomas. (d) Área del norte de Fuerteventura (Corralejo). (e) Área del sur de El Hierro (La restinga). CAPÍTULO 1 13 En el proceso de desarrollo de esta tesis, el primer paso de esta investigación ha sido la identificación del sensor y plataforma apropiados para realizar este tipo de estudio. Debido a la complejidad del modelado físico-radiativo del medio marino se hace necesario, no sólo una gran resolución espacial, sino también espectral y radiométrica. El satélite WV2 ha demostrado tener gran potencial en la monitorización de aguas costeras e interiores, con ocho bandas multiespectrales de alta resolución, codificadas en 11 bits de cuantificación, siendo el único satélite, junto al reciente WV3, que dispone de estas características. El segundo paso ha requerido el estudio de las correcciones atmosféricas y del brillo solar, definido en inglés con el término glinting, de la superficie del agua. Este tipo de correcciones son imprescindibles, debido a la gran aportación tanto de la atmósfera como del glinting a la reflectividad obtenida por el satélite, respecto a la baja reflectividad del agua. En esta fase se han analizado diferentes algoritmos y modelos atmosféricos, algunos de ellos más adecuados para entornos marinos. Así, una importante parte de esta investigación se ha focalizado en la obtención de la mejor corrección atmosférica y del brillo solar, que permita la obtención de la reflectividad marina de una forma fidedigna. El tercer paso ha consistido en el modelado radiativo del agua, centrándose en los dos principales efectos en la luz: la absorción y la retro-difusión, conocido en inglés como back-scattering. El efecto del agua y los diferentes componentes en suspensión y disolución producen un efecto de absorción y back-scattering para cada banda del satélite que puede ser modelado según su profundidad mediante sistemas de ecuaciones no lineales. De esta forma, el modelado mediante las ecuaciones de transferencia radiativa (RTE, Radiative Transfer Equations) permite la obtención de la concentración de los principales parámetros de la calidad del agua, Clorofila (chl-a), materia en disolución (CDOM) y materia suspendida (TSS). Además de estos parámetros, la profundidad (batimetría) y el albedo del fondo costero son variables modelables en el sistema de ecuaciones, lo que nos permite obtener a su vez mapas batimétricos y de albedo del fondo costero de alta resolución. El cuarto y último paso ha consistido en la obtención del albedo del fondo costero a partir de las bandas de mayor penetración. Para este cometido se ha hecho uso de técnicas avanzadas de desmezclado lineal de clases bentónicas puras (conocido en inglés como endmembers). Finalmente, indicar que en la realización de esta Tesis Doctoral ha sido necesaria la investigación de múltiples campos heterogéneos del conocimiento, como son: (i) el estudio del funcionamiento de los nuevos satélites y sensores ópticos de muy alta resolución, tanto los aspectos físicos y geométricos, que permiten entender el funcionamiento y adquisición de las imágenes por el sensor, como las limitaciones de este tipo de plataforma; (ii) conocer los fenómenos físicos relacionados con la transferencia radiativa de la atmósfera y de los entornos marinos; (iii) analizar los fenómenos químicos y biológicos relacionados con la calidad de aguas y las especies bentónicas que coexisten en los fondos costeros; (iv) procesado de imágenes orientado a la eliminación de ruidos presentes en las imágenes de teledetección los cuales son más evidentes tras la aplicación de los algoritmos de eliminación del brillo solar, y (v) aplicar métodos numéricos avanzados para la resolución de modelos matemáticos complejos. 1.2. Objetivos Como se ha mencionado, el uso de imágenes satelitales de muy alta resolución para el modelado radiativo de las aguas costeras y sus componentes, así como el modelado del albedo del fondo costero y su batimetría resulta ser un problema muy complejo y en fase incipiente de investigación. La comprensión detallada de esta temática puede llegar a ser un ejercicio muy complejo para los investigadores que no están intensamente especializados en ella. Por lo tanto, resulta más útil Introducción 14 dividir un problema de esta naturaleza en una serie de tareas más pequeñas y manejables. De esta forma, esta metodología nos permitirá presentar una descripción clara del proceso de modelado realizado en esta investigación. El propósito de esta sección es proporcionar una descripción del problema a resolver para seguidamente presentar los objetivos de esta Tesis Doctoral y las tareas necesarias para la consecución de los objetivos planteados. 1.2.1. Planteamiento del problema El propósito de esta investigación es el uso de imágenes de satélite de muy alta resolución espacial WoldView-2 para la monitorización de entornos acuáticos costeros mediante la implantación del modelado de transferencia radiativa. Dicho modelado físico permitirá obtener parámetros de calidad de aguas, batimetría y albedo del fondo costero con una muy alta resolución espacial. Aunque nuestro objetivo es el modelado radiativo del agua, debido a las grandes perturbaciones (atenuación y retro-difusión) que origina la atmósfera respecto a la baja cantidad de energía aportada por el entorno acuático, se ha requerido del modelado conjunto de estos dos medios en pos de la eliminación del componente atmosférico de la reflectividad obtenida por el satélite a lo alto de la atmósfera (ToA: Top of Atmosphere). Junto a estos fenómenos es necesario modelar la reflectividad especular de la superficie marina con el objetivo de eliminar dicha componente. Además de lograr modelar radiatívamente el comportamiento de la atmósfera y de los elementos contenidos en el agua, es de gran interés traducir estos datos a una estimación de concentraciones de clorofila, materia suspendida y materia disuelta. Para ello se ha realizado un esfuerzo para la estimación del comportamiento radiativo de estos elementos. El cálculo de la profundidad de las áreas costeras, conocido como batimetría, es un proceso costoso que típicamente es realizado normalmente mediante barcos equipados de sonares de barrido. Evidentemente, la obtención de mapas de batimetría de alta resolución mediante imágenes de teledetección tiene un gran interés gracias a que minimiza los costes. Sin embargo, la complejidad requerida para la obtención de buenos mapas batimétricos, con errores aceptables, requiere de modelos complejos y de condiciones ideales del posicionamiento del satélite, de la calidad atmosférica, oleaje y transparencia del agua. Asimismo, la obtención del albedo del fondo costero, con aplicación directa a la detección y clasificación de los diferentes tipos de elementos bentónicos que se encuentran en la superficie del lecho marino, proporciona a los biólogos una importante herramienta para la gestión medioambiental de estos entornos. Para ello se han de modelar los diferentes tipos de fondos mediante el método de desmezclado lineal de firmas espectrales puras de las principales clases bentónicas presentes en cada uno de los fondos costeros. En la Figura 2, se presenta el diagrama de flujo de la arquitectura propuesta en esta Tesis Doctoral. CAPÍTULO 1 15 Figura 2. Diagrama de flujo del sistema del modelado radiativo de imágenes de alta resolución WV2 para la obtención de parámetros de calidad de agua, batimetría y albedo del fondo. Las imágenes adquiridas del satélite WV2 son calibradas radiométricamente para obtener los valores físicos de radiancia a partir de los valores digitales de la imagen. A continuación, mediante un modelado atmosférico, estos datos son convertidos a valores de reflectividad superficial, en donde se normalizan los valores de radiancia según las condiciones de iluminación y se corrigen los efectos de absorción y back-scattering de la atmósfera, lo que permite obtener la reflectividad de la superficie (ToC: Top of Canopy). El siguiente módulo hace uso de la reflectividad superficial para corregir los efectos de brillo especular de la superficie del agua, el cual no aporta información sobre los fenómenos de absorción y retro-difusión producidos por debajo de la superficie. A continuación, el siguiente módulo hace uso de los canales corregidos del WV2 para realizar el cálculo de las ecuaciones de transferencia radiativa (RTE). De dicho cálculo se obtienen las propiedades inherentes de las substancias contenidas en el agua (IOPs: Inherent Optical Properties), la profundidad de la columna de agua (batimetría) y una estimación del albedo del fondo costero. Finalmente, gracias a un post-procesado final, en el que se infiere el comportamiento de las principales substancias presentes en las aguas costeras, se pueden obtener mapas de concentraciones de clorofila (chl-a), la concentración total de sólidos en suspensión (TSS) y materia disuelta (CDOM). A su vez, a partir de la información de reflectividad de los canales de mayor penetración y a la información de batimetría, se pueden obtener mapas del albedo del fondo marino y las abundancias de las clases puras modeladas en el desmezclado lineal. Imágenes Multiespectrales de Alta Resolución Espacial 22 1.85 m y una banda pancromática de 0.46 m de resolución espacial, proporcionando un ancho de imagen de 16.4 km. Las imágenes se suministran comercialmente con resoluciones de 0.5 y 2 m, respectivamente. WV2 está dotado de un sensor tipo push-broom, el cual funciona como un escáner que obtiene imágenes de una dimensión de forma continua para cada instante de su órbita. WV2 proporciona imágenes de 11 bits de resolución radiométrica para sus nueve bandas (pancromática, costal-blue, blue, green, yellow, red, red-edge, NIR1, NIR2). La introducción de cuatro nuevas bandas espectrales respecto a las usuales (R, G, B e IR cercano), de otros satélites de muy alta resolución, es el rasgo más diferenciador del satélite, proporcionando un gran potencial en el procesado espectral de las imágenes. Las principales ventajas frente a otros satélites de teledetección de alta resolución son las siguientes: 1. Mayor resolución espectral. WV2 es el primer satélite comercial en proporcionar alta resolución espacial e imágenes multiespectrales de 8 bandas. Así, además de las cuatro bandas multiespectrales estándares (R, G, B e IR cercano), incluye 4 bandas adicionales para mejorar el análisis espectral y permitir nuevas aplicaciones. WV2 permite la adquisición por separado de cuatro u ocho bandas gracias a la utilización de dos subsistemas independientes. 2. Mayor agilidad. La serie de satélites WorldView son las primeras plataformas comerciales con capacidad para controlar los momentos de fuerza generados por los giróscopos. Esta tecnología de alto rendimiento proporciona una aceleración hasta 10 veces superior que la de otros actuadores de control de actitud y mejora la maniobrabilidad y la capacidad de orientación. El tiempo de giro se reduce de más de 60 segundos a sólo 9 para cubrir 300 km, permitiendo obtener imágenes de diferentes zonas en un pase orbital único. 3. Mayor capacidad de revisita. Gracias a su mayor agilidad, WV2 puede obtener imágenes multiespectrales de diferentes áreas en un solo pase. WV2 tiene una capacidad de exploración de hasta 975.000 km2 por día. La combinación de una mayor agilidad y la altura orbital le permite realizar una revisita de cualquier zona en 1.1 días. 4. Mejor precisión. La tecnología avanzada de posicionamiento de WV2 es lo que permite mejoras significativas en su precisión de geolocalización. La especificación de precisión se ha mejorado hasta 6.5 m CE90 sin ningún tratamiento adicional, uso de modelo de elevación, ni puntos de control en tierra, pudiendo llegar hasta 2 m con datos adicionales. 2.1.1. Características espectrales WV2 lleva un instrumento que genera una imagen pancromática de 0.46 m de resolución espacial, con una banda de paso que abarca desde el azul hasta el infrarrojo cercano, y ocho bandas espectrales estrechas (entre 40 y 60 nm), para este tipo de satélites, de 1.85 m de resolución espacial. Las ocho bandas multiespectrales son capaces de proporcionar una precisión de color excelente, permitiendo el desarrollo de nuevas aplicaciones. Las nuevas bandas proporcionadas por WV2, comparativamente con los satélites usuales, se distribuyen de la siguiente manera: una centrada en longitudes de onda más corta que el azul, aproximadamente en 427 nm; una banda amarilla, a 608 nm; una banda en el borde del rojo, centrada estratégicamente en, aproximadamente, 724 nm al ser el inicio de la parte de alta reflectividad de la respuesta de la vegetación, y una adicional en el infrarrojo cercano, pero a mayor longitud de onda, centrada aproximadamente en 949 nm, que es sensible al vapor de agua atmosférico. CAPÍTULO 2 23 Las principales características de las bandas multiespectrales del WV2 son las siguientes [1]:  Azul costero: Ayuda a la realización de análisis vegetativo debido a su alta absorción. Mayor penetración en el agua, muy útil en los estudios batimétricos. Tiene el potencial para mejorar las técnicas de corrección atmosférica debido a la mayor absorción del ozono y por su alto nivel de difusión de Rayleigh.  Azul: Es equivalente a la banda azul de otros satélites de alta resolución como el QuickBird. Al igual que el azul costero es afectado fuertemente por la absorción de la clorofila y tiene una alta penetración en el agua.  Verde: Esta banda es más estrecha que la proporcionada por QuickBird. Permite la detección de la vegetación gracias al pico de reflectividad dentro del rango visible. Combinado con la banda amarilla permite discriminar ente tipos de vegetaciones terrestres. Esta banda permite una amplia penetración en el agua, incluso de mayor profundidad en aguas costeras con altos niveles de materia disuelta y/o suspendida.  Amarillo: Muy importante para la clasificación, dado que detecta la "amarillez" particular de la vegetación. Para el caso del agua, permite una penetración moderada debido al incremento de la absorción del agua en esta longitud de onda.  Rojo: Esta banda es más estrecha que la proporcionada por QuickBird. Permite discriminar la vegetación terrestre, junto al rojo borde, gracias a su baja reflectividad. La penetración de esta banda en el agua se reduce a unos pocos metros.  Rojo borde: Muy valiosa para medir la salud de plantas y ayudar en la clasificación de la vegetación. La penetración de esta banda en el agua es muy reducida.  Infrarrojo cercano (NIR1): Proporciona una alta separación con la banda rojo borde, proporcionando un alto valor de reflectividad en la vegetación terrestre siendo muy útil en la clasificación de cubiertas terrestres. La penetración en el agua es casi nula siendo muy importante en los algoritmos de eliminación de brillo solar.  Infrarrojo cercano (NIR2): Tiene un alto valor de reflectividad en la vegetación terrestre siendo muy útil en la clasificación de cubiertas terrestres. Se ve menos afectada por la influencia de la atmósfera. Permite el análisis de la vegetación y estudios de la biomasa. La penetración en el agua es casi nula siendo muy importante en los algoritmos de eliminación de brillo solar. Las imágenes Worldview-2 se pueden adquirir en 3 niveles de procesamiento [2]:  Basic, que incluyen únicamente la corrección radiométrica.  Standard/Ortho-ready, que incluyen la corrección radiométrica y geométrica.  Stereo, que incluyen 2 escenas superpuestas con ángulo de visión estéreo. En la realización de esta Tesis Doctoral se han utilizado imágenes Standard/Ortho-ready. 2.2. Calibración radiométrica de los datos WorldView-2 Las correcciones radiométricas son aquellas técnicas que tienen por objeto modificar los niveles digitales de las imágenes procedentes de los sensores de observación de la Tierra, con el fin de corregir los problemas derivados del funcionamiento de los mismos. La respuesta de la radiancia espectral relativa se define como el rango del número de fotoelectrones medidos por el sistema, convertidos en radiancia espectral a una longitud de onda concreta, presente a la entrada de la apertura del telescopio. Esto no sólo incluye la eficiencia cuántica del detector, sino también las pérdidas de transmisión debidas a la óptica del telescopio y a los filtros ópticos multiespectrales. A continuación se describirán las características radiométricas de las imágenes procedentes del satélite WV2 [3], cuya respuesta espectral se muestra en la Figura 3. Las bandas de paso del sistema, determinadas a partir de las respuestas espectrales, se proporcionan en la Tabla 1. Imágenes Multiespectrales de Alta Resolución Espacial 24 Figura 3. Respuesta espectral de las bandas del satélite WorldView-2. (Fuente: DigitalGlobe). Tabla 1. Bandas de paso del satélite WorlView-2 en micrómetros. Banda espectral Longitud de onda central 50% Banda de paso 5% Banda de paso Pancromático 0.632 0.464 – 0.801 0.447 – 0.808 Costera 0.427 0.401 – 0.453 0.396 - 0.458 Azul 0.478 0.448 – 0.508 0.442 – 0.515 Verde 0.546 0.511 – 0.581 0.506 – 0.586 Amarilla 0.608 0.589 – 0.627 0.584 - 0.632 Rojo 0.659 0.629 – 0.689 0.624 – 0.694 Rojo borde 0.724 0.704 – 0.744 0.699 - 0.749 Infrarrojo cercano 1 0.831 0.772 – 0.890 0.765 – 0.901 Infrarrojo cercano 2 0.908 0.862 – 0.954 0.856 - 1.043 El ancho de banda efectivo para cada banda se define como: ∆′󰇛󰇜∙ ∞  (1) donde ∆ es el ancho de banda efectivo medido en  para una banda, y ′󰇛󰇜 es la respuesta de radiancia espectral relativa. El ancho de banda efectivo debería ser utilizado en la conversión a radiancia espectral en el nivel de la atmósfera para cada banda. Este dato está incluido en el fichero de metadatos (.IMD) que acompaña a cada producto y se proporciona en la Tabla 2. CAPÍTULO 2 25 Tabla 2. Ancho de banda efectivo de las bandas del WorldView-2. Banda espectral Ancho de banda efectivo [μm] Pancromática 0.2846 Costera 0.0473 Azul 0.0543 Verde 0.0630 Amarillo 0.0374 Rojo 0.0574 Rojo borde 0.0393 Infrarrojo cercano 1 0.0989 Infrarrojo cercano 2 0.0996  Irradiancia solar El instrumento WV2 es sensible a las longitudes de onda de la luz en el visible hasta áreas del infrarrojo cercano del espectro electromagnético. En esta región, la medida de radiancia en el nivel de la atmósfera, medida por el WV2, está dominada por la radiación solar reflejada, donde la irradiación de la superficie y atmósfera, como cuerpo negro, es despreciable. La irradiancia espectral está definida como la energía por unidad de área que incide sobre la superficie en función de la longitud de onda. Como el Sol actúa como un radiador de cuerpo negro, la irradiancia solar espectral puede ser aproximada con las curvas de cuerpo negro de Planck a 5900 K, corregido por el área del disco solar y la distancia entre la Tierra y el Sol. Sin embargo, un modelo de irradiancia solar [4] fue creado por el World Radiation Center (WRC) a partir de una serie de medidas solares y es el utilizado para la conversión a reflectancia, como se muestra en la Figura 4. Figura 4. Curva de irradiancia espectral solar estándar WRC. (Fuente: NASA). Imágenes Multiespectrales de Alta Resolución Espacial 26 Como se muestra en la figura, la curva de irradiancia espectral solar WRC, alcanza su máximo en 450 nm (azul-costero y azul) y disminuye ligeramente para las demás longitudes de onda. En general, la irradiancia solar espectral promedio está definida como la media ponderada de los valores de máxima irradiancia efectiva normalizada sobre la banda de paso del detector, definida por:  󰇛󰇜∙′󰇛󰇜∙ ∞ ′󰇛󰇜∙ ∞   (2) donde  es la irradiancia solar espectral promedio en [Wm-2µm-1] para una banda dada, 󰇛󰇜 es la curva de irradiancia espectral solar WRC [Wm-2µm-1], mostrada en la Figura 4, y 󰆒󰇛󰇜 es la respuesta de radiancia espectral relativa para una banda dada. En concreto, para el WV2 los valores de irradiancia solar espectral promedio para una distancia Tierra-Sol de 1 Unidad Astronómica, normal a la superficie iluminada, se proporcionan en la Tabla 3. Tabla 3. Irradiancia solar espectral promedio para las bandas del WorldView-2. Banda espectral Irradiancia solar [W*m-2*µm-1 ] Pancromática 1580.8140 Costera 1758.2229 Azul 1974.2416 Verde 1856.4104 Amarillo 1738.4791 Rojo 1559.4555 Rojo borde 1342.0695 Infrarrojo cercano 1 1069.7302 Infrarrojo cercano 2 861.2866  Conversión de valores digitales a radiancia Asumiendo que los detectores tienen una respuesta lineal en función de la radiancia de entrada, los niveles digitales de la imagen vendrá dada por, ∙  (3) donde  es la radiancia de la banda [Wm-2sr-1µm-1 ],  es la ganancia absoluta, y  es el offset del instrumento. De esta manera se determina una única ganancia para cada banda, y luego cada detector es escalado respecto a los otros detectores en la misma banda. Separando estas ganancias absoluta y relativa se llega a la siguiente expresión:  ∙∙(4) donde  es la ganancia absoluta y  es la ganancia relativa del detector. Por definición la media relativa de las ganancias tiene que ser igual a uno. Para normalizar la nomenclatura, redefinimos  como , ∙ como  y como . Así, podemos expresar la ecuación (4) como, CAPÍTULO 2 27  ∙   (5) donde  son los datos crudos,  son los datos del detector radiométricamente corregidos, los cuales han sido linealmente escalados respecto a los valores de radiancia,  es la ganacia relativa del detector, y  es el offset. Los valores de configuración de la ganancia para la imagen pancromática (PAN) y las bandas multiespectrales (MS) del WV2 dependen de varios parámetros como la banda, el tiempo-retardo de integración (TDI, time-delayed-integration), la resolución radiométrica del producto, etc. Los valores de ganancia apropiados, teniendo en cuenta la combinación de estos parámetros, son proporcionados en el archivo de metadatos (.IMD) que acompaña a cada producto.  Corrección radiométrica de los productos WorldView-2 La calibración y corrección radiométrica relativa es necesaria debido a que una escena uniforme no crea una imagen cruda uniforme en términos de niveles digitales. La mayoría de las causas de esta no uniformidad son la variabilidad en la respuesta de los detectores, variabilidad en la ganancia y offset, caída de las lentes y partículas contaminantes sobre el plano focal. Estos fallos producen un efecto de rayas verticales y bandeadas en la imagen. La Figura 5 muestra un ejemplo del efecto de bandeado de la imagen producida, en este caso, por la diferencia entre las ganancias y offset de los registros de entrada. La corrección radiométrica relativa es realizada sobre los datos crudos provenientes de todos los detectores en todas las bandas, durante los momentos iniciales de generación de los productos del satélite. Esta corrección incluye la eliminación del offset y una corrección no uniforme. A partir de la ecuación (5), los datos del detector, radiométricamente corregidos, son:     (6) Los resultados de las correcciones radiométricas, aplicadas a la imagen WV2 de la Figura 5 (a), se muestran en la Figura 5 (b). Se puede observar la desaparición virtual de los defectos de la imagen sin corregir radiométricamente. (a) (b) Figura 5. (a) Datos crudos de una imagen WorldView-2 sin corrección radiométrica. (b) Imagen WorldView-2 tras la corrección radiométrica. (Fuente: DigitalGlobe). Imágenes Multiespectrales de Alta Resolución Espacial 28  Radiancia en la parte superior de la atmosfera (Top of Atmosphere: ToA) La radiancia ToA está definida como la radiancia reflejada por la superficie de la tierra y la columna vertical de la atmosfera que entra por la apertura del telescopio a la altura del satélite, 770 km para el caso del WV2. La conversión desde los datos del satélite, radiométricamente corregidos, a valores de radiancia se realiza mediante,  ∙, ∆ (7) donde  representa la imagen de radiancia,  es el factor de calibración radiométrica para una banda, ,es la imagen radiométricamente corregida y ∆ es el ancho de banda efectivo para cada banda específica. La eliminación del offset no es necesaria en este punto, dado que ya ha sido realizado en el paso de corrección radiométrica durante la generación del producto. La conversión a radiancia es un proceso simple que implica la realización de dos pasos: 1. Multiplicar los valores de píxel de la imagen corregida radiométricamente por el factor de calibración K. 2. Dividir el resultado por el ancho de banda efectivo apropiado para cada banda. Los valores del factor de calibración y el ancho de banda efectivo para cada banda son proporcionados con cada producto del WV2, y se encuentran localizados en el archivo de metadatos (.IMD). Las variaciones en la irradiancia solar están dominadas por la geometría solar durante la adquisición de una imagen específica. El Sol puede ser aproximado como un punto dado en que la distancia entre la Tierra y el Sol es mucho mayor que el diámetro del Sol. La irradiancia de un punto es proporcional a la inversa del cuadrado de la distancia. Así, la irradiancia de un punto a la distancia deseada, puede ser calculada dando la irradiancia de la fuente a una distancia específica, mediante, ∙  (8) donde  es la irradiancia buscada a la distancia deseada , y  es la irradiancia conocida de la fuente a la distancia específica . La distancia media entre la Tierra y el Sol es una unidad astronómica (UA). Así la ecuación queda como,   (9) dado que el dato de irradiancia solar a la distancia de 1 UA ha sido anteriormente definido como , se puede reescribir la ecuación anterior como,    (10) donde  es la irradiancia solar media a una distancia dada entre la Tierra y el Sol,  es la irradiancia solar dada en la Tabla 3 y   es la distancia entre el Sol y la Tierra, dada en unidades astronómicas, en el momento de adquisición de la imagen. CAPÍTULO 2 29 La irradiancia solar definida se refiere a la normal a la superficie que está siendo iluminada. A medida que el ángulo cenital solar se mueve fuera de la normal, aparece el efecto de área proyectada y, consecuentemente, el mismo haz de luz ilumina un área más grande. Este efecto es una función del coseno del ángulo de iluminación y viene dado por, ∙󰇛󰇜 (11) donde  es la irradiancia solar media para un ángulo cenital del Sol,  es la irradiancia solar media normal a la superficie que está siendo iluminada, que coinciden con los datos de la Tabla 3 y  es el ángulo cenital del Sol. Las ecuaciones (10) y (11) pueden ser combinadas para obtener la irradiancia solar media, teniendo en cuenta la geometría del Sol, en el momento de adquisición de la imagen:   ∙󰇛󰇜 (12) Para calcular la distancia entre la Tierra y el Sol, para un producto concreto, es necesario calcular el Día Juliano (JD) a partir de los datos de la hora de adquisición de la imagen [5]. La hora de adquisición se encuentra en el archivo de metadatos y representa la hora UTC. Del formato UTC, se extrae el año, el mes, el día y se calcula la hora universal (UT, Universal Time),  60.0  3600.0 (13) El Día Juliano, se puede calcular a partir de la siguiente ecuación:  󰇟365.25∙󰇛ñ4716󰇜󰇠󰇟30.6001∙󰇛1󰇜󰇠í 24.01524.5 (14) donde 󰇡ñ 󰇢, 2󰇡 󰇢. int hace referencia al truncado del valor decimal, quedándose sólo con la parte entera del número en cuestión. Si la imagen fue adquirida en Enero o Febrero, antes de aplicar la ecuación del día juliano, el año y el mes debe ser modificado como: (año = año - 1) y (mes = mes +12). Una vez se ha calculado el JD, se puede obtener la distancia entre la Tierra y el Sol a partir de la siguiente ecuación [6]: 1.000140.01671∙󰇛󰇜0.00014∗󰇛2󰇜 󰇟.  .󰇠 (15) donde, 357.5290.98560028∗󰇟󰇠, y 2451545.0. El valor medio del ángulo cenital del Sol de una imagen de muy alta resolución es una aproximación válida para toda la imagen debido a que la variación del ángulo es despreciable dentro de la escena. La media del ángulo de elevación (), expresado en grados, para un producto dado, es calculada para el centro de la escena y puede ser obtenida a partir del archivo de metadatos. El ángulo cenital del Sol se puede obtener, a partir de este dato, aplicando la siguiente ecuación: Imágenes Multiespectrales de Alta Resolución Espacial 30 90.0 (16) La forma de las curvas de radiancia espectral a la altura de la atmósfera, como una función de las longitudes de onda del WV2, están dominadas por la forma de la curva solar. Esta radiancia puede ser modelada con la contribución de tres radiaciones principales:  (17) donde  es la radiancia total a lo alto de la atmósfera, es la radiancia reflejada por la superficie,  es la radiación dispersada por la atmósfera y reflejada en la superficie y es la radiancia dispersada por las moléculas de la atmósfera que depende de la longitud de onda. Expandiendo la radiación reflejada por la superficie, y asumiendo una superficie Lambertiana, se puede expresar como: ∙∙,∙󰇛󰇜  ∙ (18) donde  es la reflectancia espectral difusa,  es la transmisividad de la atmósfera en la dirección del sensor,  es la transmisividad de la atmósfera que atraviesa la radiación solar dirección descendente,  es la irradiancia solar,  es el ángulo cenital del Sol, y   es la distancia entre la Tierra y el Sol. Sustituyendo en la ecuación (17), obtenemos lo siguiente: ∙∙,∙󰇛󰇜  ∙  (19) Sin tener en cuenta los efectos de la atmosfera se puede obtener la siguiente ecuación: ,∙󰇛󰇜  ∙ (20) Finalmente, reorganizando los términos, para despejar el valor de reflectancia, obtenemos la ecuación que permite calcular la reflectividad difusa como: ∙ ∙ ,∙󰇛󰇜 (21) 2.3. Resumen En este capítulo se han introducido las principales características del satélite WV2, indicándose las mejoras proporcionadas en la resolución espacial de las imágenes, en la agilidad de la adquisición de las escenas, mejorándose a su vez el tiempo de revisita y la precisión geométrica de las imágenes. A continuación, se han descrito las características espectrales del sensor, donde se han presentado las ocho bandas multiespectrales, así como las principales aplicaciones de dichas bandas, siendo el incremento en la resolución espectral en ocho bandas multiespectrales más estrechas, el factor más determinante en la selección del satélite WV2. Finalmente, se han descrito CAPÍTULO 2 31 los pre-procesados necesarios para la corrección radiométrica de las imágenes de alta resolución, detallándose, a su vez, los pasos necesarios para la obtención de los parámetros físicos de radiancia y de reflectividad difusa a lo alto de la atmósfera a partir de los valores digitales de las imágenes. En el capítulo 3, una vez que hemos obtenido la reflectividad en la parte superior de la atmósfera, se va a proceder a describir la corrección atmosférica de las imágenes de teledetección necesaria para la eliminación de las perturbaciones que dicho medio introduce en la radiación recibida por el sensor. De esta manera es posible corregir estos fenómenos lográndose determinar la reflectividad superficial del agua marina que se pretende monitorizar. Modelado Atmosférico de las Imágenes Multiespectrales del Satélite WorldView-2 38 Junto a la propuesta anterior, también existe otra simplificación del método COST, que consiste en aproximar el efecto multiplicativo de la transmisividad, que afecta al rayo incidente 󰇛󰇜, por el coseno del ángulo cenital del Sol y la transmisividad de la atmósfera para el flujo ascendente 󰇛󰇜, por el coseno del ángulo de observación. Finalmente, podemos expresar la ecuación de la reflectividad superficial mediante: 󰇛 󰇜∙∙ cos∙󰇛,∙cos∙cos󰇜 (32) 3.3. Algoritmo de corrección atmosférica basado en parámetros físicos: Modelo 6S El Second Simulation of a Satellite Signal in the Solar Spectrum (6S) es un modelo avanzado de transferencia radiativa diseñado para simular la reflexión de la radiación solar en condiciones de una atmósfera libre de nubes, según condiciones específicas geométricas y espectrales. Este modelo tiene en cuenta los principales parámetros atmosféricos para modelar la dispersión y la absorción que produce la atmósfera en la longitud de onda del canal del satélite. El algoritmo 6S es utilizado para generar las LUTs (Look-Up Tables) en los algoritmos de corrección atmosférica del sensor MODIS de la NASA [11] [12]. El código fuente de este algoritmo es libre bajo licencia GNU programado en el lenguaje de programación Fortran, frente a otros modelos comerciales como el FLAASH (Fast Line-of-Sight Atmospheric Analysis), basado en el código de transferencia radiativa MODTRAN encapsulado en el software ENVI [13], y el ATCOR integrado en el paquete comercial de procesado ERDAS [14]. El código de transferencia radiativa 6S estima la reflectividad aparente ToA teniendo en cuenta los efectos de absorción de los gases, la difusión de las moléculas y aerosoles presentes en la atmosfera, y la falta de homogeneidad de la reflectividad de la superficie terrestre [15]. 6S define  como la reflectividad superficial del objetivo, rodeado de un entorno homogéneo de reflectividad , mientras que la reflectividad ToA se define como , la cual es definida de la siguiente manera. 󰇛,,Δ󰇜󰇛,󰇜󰇛,,Δ󰇜󰇣 󰇛󰇜󰇤 󰇛󰇜 1 (33) La referencia a la longitud de onda ha sido eliminada para una mayor claridad de la ecuación. El significado de los términos de la ecuación se explica a continuación:  cos󰇛󰇜,cos󰇛󰇜.  Δ representa la diferencia entre el azimut solar y del satélite.   representa la transmisibilidad total de los gases (en la trayectoria de bajada y subida), teniendo en cuenta la absorción de los diferentes gases de la atmosfera.   representa a la reflectividad atmosférica, el cual depende de las propiedades moleculares y de los aerosoles presentes en la atmosfera.   representa al espesor atmosférico (Atmospheric Optical Depth, AOD)  󰇛󰇜,󰇛󰇜 representa a la transmitancia difusa de la atmósfera.   representa el albedo esférico de la atmosfera. CAPÍTULO 3 39 El término 󰇛1󰇜 tiene en cuenta las difusiones múltiples entre la superficie y la atmósfera. Como se puede observar, en la ecuación 6S se abordan separadamente los procesos de absorción y difusión. Las difusiones producidas por las moléculas y por los aerosoles son, análogamente, diferenciados. La reflectividad total de la atmósfera se obtiene mediante la introducción de coeficientes procedentes del scattering de Rayleigh y de los aerosoles. Estos coeficientes son obtenidos mediante aproximaciones de primer orden. 6S es un modelo de una única capa, en donde no se tiene en cuenta las variaciones de los parámetros en la columna vertical de la atmósfera. La ecuación (33) es una expresión monocromática. Para obtener el resultado de una banda multiespectral, 6S calcula la ecuación para todo el rango de longitudes de onda con un paso de 5 nm, integrando todos estos resultados en el resultado final de la banda, según la respuesta espectral del sensor para dicha longitud de onda. 3.3.1. Configuración del modelo de corrección 6S La configuración de entrada del modelo se divide en 5 partes principales: condiciones geométricas, modelado atmosférico y de aerosoles, alturas del área de estudio y del sensor, condiciones espectrales, y reflectancia del suelo. A continuación, procederemos a su descripción.  Condiciones geométricas El modelo 6S define diferentes modos de insertar las condiciones geométricas según el tipo de satélite o plataforma aerotransportada utilizada. Define condiciones para satélites geoestacionarios como el METEOSAT, así mismo parametriza de forma sencilla las condiciones de los satélites polares de observación de la Tierra como el AVHRR (NOAA), el HRV (SPOT) y el TM (LANDSAT). Para definir las condiciones geométricas de otros satélites polares de alta resolución como el WV2 es necesario utilizar el tipo genérico denominado User’s. Los parámetros introducidos en el modelo para satélites genéricos son los siguientes: 1. Mes 2. Día 3. Ángulo cenital del Sol 4. Ángulo acimutal del Sol 5. Ángulo cenital del sensor 6. Ángulo acimutal del sensor  Modelado atmosférico Para modelar la atenuación de la atmósfera en el momento de adquisición de los datos, el modelo 6S requiere de la definición del tipo atmosférico-climático de la zona. Estos tipos son los siguientes: 1. Sin absorción gaseosa 2. Tropical 3. Latitud media de verano 4. Latitud media de invierno 5. Sub-ártico de verano 6. Sub-ártico de invierno 7. Estándar 62 US Modelado Atmosférico de las Imágenes Multiespectrales del Satélite WorldView-2 40 8. Definición de la concentración de vapor de agua y de ozono 9. Perfil de usuario por base de datos a. Altura [km] b. Presión [mb] c. Temperatura [K] d. Densidad de H2O [g/m3] e. Densidad de O3 [g/m3]  Modelo de aerosoles Para modelar el fenómeno del scattering de la atmósfera se hace necesaria la selección del modelo de tipos y concentración de aerosoles que se va a utilizar de entre los siguientes: 1. No aerosoles 2. Modelo continental 3. Modelo marítimo 4. Modelo urbano 5. Aerosoles zonas desérticas 6. Quema de biomasa 7. Modelo estratosférico  Para los modelos del 2 al 7 se ha de definir la profundidad óptica (AOD) a 550 nm o su equivalente en visibilidad en km 8. Definición de los valores del modelo  C(1) = volumen % de polvo  C(2) = volumen % de agua suspendida  C(3) = volumen % de agua oceánica suspendida  C(4) = volumen % de hollín Los valores en un escenario típico, por ejemplo, para las Islas Canarias, serán: latitud media y modelo marítimo. La utilización de valores obtenidos por radiosondas permite modelar de una manera más precisa la atmósfera. Al disponer de bases de datos de medidas diarias, con radiosondas, permite mejorar los resultados de la corrección atmosférica [16] [17]. El parámetro de profundidad óptica (AOD: Aerosol Optical Depth) necesaria para los modelos de aerosoles puede ser obtenida en los productos atmosféricos generados por la NASA a partir del sensor MODIS [18].  Altura del área de estudio y del satélite En este paso se proporciona al modelo los datos de altura media del terreno en estudio y de la altura de la plataforma: 1. Altura del área en estudio en [km] 2. Altura del sensor en [km]  Condiciones espectrales Para la definición de las condiciones espectrales es necesario conocer el ancho de banda y la forma de la banda de paso del sensor, existiendo perfiles de los sensores más utilizados de media y baja resolución como el SEVIRI del MSG y el TM del LANSAT. Para otros sensores de alta resolución como el WV2 es necesario obtener la respuesta del filtro paso banda de los canales CAPÍTULO 3 41 para introducirla en el modelo. Si no conocemos la forma del canal se puede optar por simular un canal filtro paso banda perfecto. A continuación, se muestran las diferentes opciones: 1. Utilizar una banda de paso igual a 1 2. Canal monocromático 3. Definición de un filtro por el usuario a. Introducción del inicio y el fin de la banda de paso b. Valores de la banda de paso en paso de 2.5 nm 4. Selección de perfiles de sensores predefinidos por el modelo  Reflectividad superficial del entorno El modelo 6S permite describir la reflectividad del entorno del objetivo, para ello se puede definir un modelo de suelo homogéneo o no homogéneo. En nuestro contexto, para el modelado de la superficie marina, seleccionamos la primera opción. Con este modelo se puede seleccionar perfiles de suelo predefinidos entre los que se encuentran los siguientes: 1. Introducir firma espectral de forma manual con paso de 2.5 nm 2. Valor medio de vegetación 3. Valor medio de agua marina de baja turbidez 4. Valor medio de agua de lago de baja turbidez 5. Valor medio de la arena Una vez configurados los archivos necesarios, para la corrección de cada uno de los canales del satélite, se ejecuta el modelo, obteniéndose los valores de los principales parámetros atmosféricos que intervienen en los fenómenos de absorción y scattering producidos por la atmosfera. Resaltar que la expresión genérica para la corrección atmosférica de las imágenes de radiancia del modelo 6S viene dada por: , 1󰇛∗󰇜, 󰇛∗󰇜 (34) donde,  , es la reflectancia superficial corregida atmosféricamente (ToC)   es la radiancia medida por el satélite en [w/m2/sr]   es la inversa de la transmitancia atmosférica   es el scattering de la atmósfera   es el albedo atmosférico para la luz isotrópica De esta forma se realiza una conversión de la radiancia medida por el satélite a reflectividad de la superficie terrestre, ya corregida atmosféricamente, mediante el uso de tres variables generadas por el modelo 6S, para las condiciones geométricas, espectrales y atmosféricas indicadas en los archivos de configuración. Modelado Atmosférico de las Imágenes Multiespectrales del Satélite WorldView-2 42 3.4. Evaluación de los algoritmos de corrección atmosférica de imágenes WV2 En este apartado se presentan los resultados obtenidos en la corrección atmosférica de imágenes multiespectrales WV2 mediante la evaluación de los algoritmos DOS y COST y el modelo atmosférico basado en parámetros físicos 6S. Primeramente, además de la localización del área de estudio, se realizará un análisis comparativo entre los valores de reflectividad obtenidos en la parte superior de la atmósfera (ToA) y los correspondientes obtenidos tanto por los algoritmos DOS y COST, basados en propiedades de la imagen, como por el modelo atmosférico 6S. Seguidamente, se presentarán los resultados obtenidos respecto a la reflectividad espectral corregida mediante el modelo atmosférico 6S en comparación de la reflectividad espectral obtenida in-situ mediante un espectro-radiómetro de campo y, finalmente, se expondrán las principales conclusiones del modelado atmosférico de las imágenes WV2. 3.4.1. Comparativa entre los métodos DOS y COST respecto al modelo 6S Para realizar la comparativa entre los diferentes métodos de corrección atmosférica se ha hecho uso de una imagen WV2 del área del puerto de Granadilla (Isla de Tenerife), específicamente, se ha utilizado la imagen del pase del 29 de octubre del 2011, como se muestra en la Figura 8. Figura 8. Imagen WV2 del puerto de Granadilla (Tenerife), 29 de octubre del 2011. CAPÍTULO 3 43 Para esta imagen se han seleccionado 6 puntos de interés, de zonas homogéneas y suficientemente amplias, que se resaltan en la figura, dos situados en tierra y otros cuatro en el mar, intentando obtener valores de agua con diferentes grados de turbidez. Los puntos de interés son los siguientes:  Punto 1: Tejado de edificio  Punto 2: Zona arenosa  Punto 3: Agua de mar de baja turbidez (Sin Turb)  Punto 4: Agua de mar con turbidez (Turb 1)  Punto 5: Agua de mar con turbidez (Turb 2)  Punto 6: Agua de mar con turbidez (Turb 3) En la Figura 9 se muestran los resultados de reflectividad obtenidos en la parte superior de la atmósfera (ToA) y los correspondientes obtenidos, tanto por los algoritmos DOS y COST, basados en propiedades de la imagen, como por el modelo atmosférico 6S, para los diferentes puntos de interés y para cada banda multiespectral del WV2. La configuración del modelo 6S para la adquisición de la imagen WV2 del día 29 de octubre de 2011 es la siguiente:  Condiciones geométricas: o Ángulo cenital solar = 41.8º o Ángulo acimutal solar =170.5º o Ángulo de visión del sensor = 20.4º o Ángulo acimutal del sensor = 283.8º  Modelo atmosférico = Latitud media de invierno  Modelo aerosoles = Marítimo  Espesor óptico de la atmósfera (AOD) = 0.22  Reflectividad superficial del entorno = superficie homogénea de agua de baja turbidez Modelado Atmosférico de las Imágenes Multiespectrales del Satélite WorldView-2 44 Figura 9. Correcciones atmosféricas para los puntos de interés en el área del puerto de Granadilla (29 de octubre del 2011). Podemos observar como las mayores diferencias obtenidas se encuentran en el canal azul costa, debidas al mayor nivel de scattering y absorción. Adicionalmente, se deduce como la reflectividad difusa en lo alto de la atmósfera (ToA) tiene valores superiores a los obtenidos con las correcciones atmosféricas para reflectividades bajas, como en el caso del agua marina. Por el contrario en áreas de alta reflectividad las correcciones atmosféricas 6S y COST corrigen al alza los valores de reflectividad, mientras que se puede observar como el algoritmo DOS no tiene en cuenta el fenómeno de absorción de la atmósfera, corrigiendo siempre a la baja las reflectividades. Como se muestra en la Figura 10, los diferentes niveles de turbidez del agua varían la reflectividad de los ocho canales del WV2, produciendo mayores alteraciones en el visible, en especial en los canales azul y verde. Por otro lado, las diferencias en el rojo e infrarrojo son bajas. Figura 10. Variación de la reflectividad corregida 6S del agua en función de la turbidez. CAPÍTULO 3 45 3.4.2. Resultados de la corrección atmosférica 6S respecto a datos in-situ Para verificar el correcto funcionamiento de la corrección atmosférica 6S, se procedió a realizar una campaña para la obtención de datos in-situ en la zona de adquisición de las imágenes, con una proximidad temporal alta lo que permitió obtener estas muestras en condiciones similares de la superficie a medir, iluminación, y condiciones atmosféricas. Para garantizar una correcta obtención de datos in-situ se eligieron áreas homogéneas suficientemente grandes que permitiese una correcta identificación de los puntos en las imágenes. Para las muestras se hizo uso de un espectro-radiómetro Vis/NIR ASD FieldSpec 3, proporcionado por el grupo GOTA (Grupo de Observación de la Tierra y la Atmósfera) de la Universidad de La Laguna. Las principales características del espectro-radiómetro son las siguientes: rango espectral entre 350-2500 nm; una resolución espectral de 3.5 nm a partir de la longitud de onda de 700 nm, y de 10 nm a partir de los 1500 hasta los 2100 nm; para ello hace uso de anchos de banda de 1.4 nm desde los 350 hasta los 1050 nm y de 2 nm desde los 1050 hasta los 2100 nm. Una de las características fundamentales del espectro-radiómetro es su bajo peso 5.2 Kg y su alta portabilidad, lo que permite la realización de campañas de muestreo in-situ de forma más sencilla. En la Figura 11, se muestra el modelo de espectro-radiómetro utilizado en la campaña. Figura 11. Espectro-radiómetro Vis/NIR ASD FieldSpec 3.  Simulación de reflectividad multiespectral a partir de datos hiperespectrales Las bandas multiespectrales del WV2 proporcionan reflectividades promediadas para todo su ancho de banda, mientras que la reflectividad del radiómetro corresponde a anchos de banda cercanos a 1 nanómetro (sensores hiperespectrales). Por este motivo, los resultados de reflectividad corregida de las bandas WV2 y los resultados obtenidos por el radiómetro no resultan ser cuantitativamente comparables. Para obtener un resultado cuantitativo de los valores de reflectividad corregida del WV2 y los datos del radiómetro se hace uso de la simulación de bandas multiespectrales a partir de datos hiperespectrales (SMS, Simulated Multi-Spectral). Dicha simulación consiste en obtener la banda de reflectividad superficial simulada, pseudo WV2, (,  ) a partir de multiplicar la función normalizada de respuesta del filtro paso banda (NMRF, Normalized Multispectral Response Modelado Atmosférico de las Imágenes Multiespectrales del Satélite WorldView-2 46 Function) de la banda (Figura 3) por el valor de reflectividad monocromáticas del radiómetro. Así, se obtiene una integración de la reflectividad del radiómetro para la banda Multiespectral, desde  hasta  [19]. ,   ∗   (35)  Resultados de la corrección atmosférica 6S respecto a la campaña realizada en Granadilla La campaña para la obtención de la reflectividad superficial, mediante el espectro-radiómetro fue realizada el día 7 de febrero del 2012, de forma que existiese la menor diferencia temporal con la adquisición de la imagen WV2. En la Figura 12 se muestra la localización de las muestras in-situ respecto a la imagen WV2 obtenida el día 18 de febrero. Se puede observar la presencia de cuatro muestras de diferentes tipos de arenas (Arena1, Arena2, Arena3, Arena4), así como la obtención de una muestra en un parking, y en una rotonda. Finalmente se obtuvieron dos muestras más en el mar, la primera en una playa y la segunda en el puerto. La configuración del modelo 6S para la adquisición del día 18 de febrero de 2012 es la siguiente:  Condiciones geométricas: o Ángulo cenital solar = 43.5º o Ángulo acimutal solar = 153.1º o Ángulo de visión del sensor = 17.5º o Ángulo acimutal del sensor = 104.1º  Modelo atmosférico = Latitud media de invierno  Modelo aerosoles = Marítimo  Espesor óptico de la atmósfera (AOD) = 0.1  Reflectividad superficial del entorno = superficie homogénea de agua de baja turbidez Figura 12. Emplazamiento de las muestras in-situ obtenidas con el radiómetro en la costa de Granadilla. CAPÍTULO 3 47 En la Figura 13 (a) se muestran las firmas espectrales, reflectividades a lo largo de la longitud de onda que define a los elementos presentes en la superficie terrestre, de los datos in-situ obtenidos con el radiómetro. En la Figura 13 (b) se muestran las reflectividades simuladas de los ocho canales multiespectrales pWV2 obtenidos a partir de los valores de reflectividad del radiómetro. En la Figura 13 (c) se muestran los valores de reflectividad corregidos de los ocho canales WV2 para los píxeles asociados a las posiciones geográficas de los datos in-situ. (a) (b) (c) Figura 13. (a) Firmas espectrales de las muestras in-situ. (b) Pseudo bandas WV2 generadas a partir del radiómetro para las muestras in-situ. (c) Reflectividad de las bandas WV2 corregidas atmosféricamente para los píxeles asociados a las ubicaciones in-situ. Corrección del Reflejo Solar en Imágenes de Alta Resolución 54 Podemos observar como los valores de reflectividad son muy reducidos para ángulos de incidencia menores a 30º, valores inferiores a 0.025. Continúan siendo muy bajos hasta los 60º, donde los valores son inferiores a 0.06. A partir de este ángulo, los valores de reflectividad crecen exponencialmente, llegando a 0.14 a los 70º, 0.35 a los 80º y 1.00 a los 90º. Estos valores nos indican que el brillo solar adquiere importancia sólo en ángulos de incidencia elevados, aunque debido a que la reflectividad del agua es muy baja siempre se ha de tener en cuenta. El brillo especular se produce dependiendo de las condiciones de iluminación (hora del día, estación y latitud de la zona de estudio), y en condiciones de fuerte oleaje. Normalmente, el oleaje es el mayor problema en la generación de brillo especular en la superficie marina. El cálculo riguroso de la reflectividad especular mediante el uso de la ecuación de Fresnel resulta inviable debido a que no es posible conocer la orientación de la superficie de las olas, parámetro completamente necesario para su cálculo. Además, dicha ecuación se basa en la existencia de una superficie plana, por lo que sería necesario trabajar con una resolución espacial muy elevada para poder aproximar la superficie del píxel a una superficie plana. Tampoco se puede despreciar el tiempo de exposición en la adquisición del píxel, el cual ronda las milésimas de segundo, por lo que la superficie puede variar y con ello el valor de reflectividad especular. Por ese motivo, la superficie especular de un píxel puede ser modelada, estadísticamente, mediante una función de distribución de densidad de probabilidad (Probability Density Function, PDF), como la probabilidad de que la pendiente de la ola esté orientada hacia el satélite, teniendo en cuenta que la fuente de luz no es puntual sino una elipse de 0.53º de diámetro angular [21] [23] [24]. En la Figura 17 se representa la reflectividad especular promediada para un píxel de 2x2 m2 con un nivel de oleaje concreto. Se puede observar como los haces de luz inciden con un ángulo similar (fuente no puntual) sobre una superficie ondulada que representa el oleaje del mar. Debido a las diferentes pendientes de la superficie marina, los ángulos de los haces de luz reflejada se orientan en múltiples direcciones. Por lo tanto, sólo una proporción de la luz incidente es reflejada de forma especular al satélite (flechas naranjas). La intensidad de la reflectividad especular para dicho haces de luz puede ser calculada mediante la ecuación (40) dado que se conocen los ángulos de incidencia del Sol, los ángulos de visión del satélite y el ángulo de tolerancia de la elipse Solar [23]. Figura 17. Reflectividad especular promedio de la superficie de un píxel. CAPÍTULO 4 55 De esta forma, se puede observar como el parámetro más importante en el cálculo de la reflectividad especular resulta ser la forma de las olas. Analizando detalladamente la representación simplificada de la figura anterior, podemos concluir que una opción para el cálculo del brillo solar pasa por modelar estadísticamente el comportamiento del oleaje y la probabilidad de que se generen ángulos especulares, teniendo en cuenta parámetros físicos como la velocidad y dirección del viento, los cuales son parámetros suficientes para modelar de forma estadística la rugosidad del oleaje. En este contexto, si en el interior de un píxel existe suficiente superficie marina con oleaje para modelar estadísticamente sus pendientes, se podrá estimar la reflectividad del brillo solar mediante la información de los ángulos de incidencia del Sol y visión del satélite, velocidad del viento y su dirección. Esta aproximación es válida para satélites de baja resolución como el SeaWiFS, MERIS o el MODIS, en escenarios de aguas oceánicas abiertas. 4.2. Estrategias de corrección del reflejo solar para longitudes de onda del visible e infrarrojo La mejor manera para abordar el problema del brillo solar es evitar su presencia mediante la selección apropiada del lugar y el tiempo de adquisición. Satélites como el SeaWiFS permiten variar el ángulo de visión respecto el nadir (tilt) hasta 20º para eludir ángulos problemáticos de visión del satélite según el ángulo de incidencia solar. Sin embargo, hay otros muchos satélites que no permiten esta opción como el MODIS o el MERIS, mientras que satélites de alta resolución como WV2, que tienen agilidad angular, no contemplan esta opción en las adquisiciones debido a que son utilizados principalmente en aplicaciones terrestres. Según lo descrito en [25] [26], los principales métodos de eliminación del brillo solar publicados en artículos científicos se pueden resumir en la Tabla 4. Existen otros métodos más complejos de deglinting que requieren de mucha más información espectral, la cual solo puede ser proporcionada por sensores hiperespectrales [27] [28] [29] [30]. Tabla 4. Principales métodos de eliminación del brillo solar. Método Autores Sensor Metodología Hipótesis Aguas abiertas Wang & Bailey [31] SeaWiFS El brillo solar es predicho mediante la velocidad del viento (ECMWF data), y sustraído de la reflectividad cuando este valor se encuentra entre dos niveles predefinidos. No se tiene en cuenta la dirección del viento ni la posibilidad de la dispersión múltiple de la atmósfera. Aguas abiertas Montagner, Billat & Belanger [24] MERIS El brillo solar es predicho mediante la velocidad y dirección del viento (ECMWF data), y sustraído de la reflectividad cuando este valor se encuentra entre dos niveles predefinidos. No se tiene en cuenta la posibilidad de la dispersión múltiple de la atmósfera. Aguas abiertas Fukushima et al. [32] GLI Similar al SeaWiFS, pero obteniendo la velocidad del viento a partir de datos del medidor de scatterer microondas (ADEOS-II). No se tiene en cuenta la dirección del viento ni la posibilidad de la dispersión múltiple de la atmósfera. Aguas abiertas Ottaviani et al. [33] SeaWiFS Obtiene la solución completa de las ecuaciones de transferencia radiativa, incluyendo efectos de la dispersión múltiple de la atmósfera, múltiples reflexiones y sobras. No se tiene en cuenta la dirección del viento. Corrección del Reflejo Solar en Imágenes de Alta Resolución 56 Aguas costeras Hochberg, Andrefouet & Tyler [34] IKONOS La banda NIR es usada para determinar la variación espacial del brillo en la imagen, escalando los resultados según la pendiente obtenida del uso de tan solo dos puntos del NIR. Los índices de refracción son independientes de la longitud de onda, no se tiene en cuenta la reflectividad del agua de la banda NIR. Aguas costeras Hedley, Harborne and Mumby [35] IKONOS Mejora el método anterior utilizando un área de valores de NIR obteniendo la pendiente mediante regresión en lugar de mediante dos puntos. Para ello usa una zona de aguas profundas. Los índices de refracción son independientes de la longitud de onda, no se tiene en cuenta la reflectividad del agua de la banda NIR. Aguas costeras Lyzenga, Malinas and Tanis [36] IKONOS El factor de corrección se basa en la medida de covariancia en la banda del NIR en un área de aguas profundas. Presupone que la radiancia de las aguas profundas es nulas. 4.2.1. Métodos basados en modelos estadísticos del estado de la superficie marina Los primeros métodos utilizados para la eliminación del brillo solar en la superficie marina fueron implementados para sensores oceanográficos de baja resolución, por ejemplo el SeaWiFS, estando basados en modelos estadísticos que predicen el estado de la superficie marina. Una vez predicha la cuantía de la reflectividad del brillo solar, ésta es sustraída de la reflectividad captada por el sensor, o bien si el valor es muy elevado es enmascarada en la imagen. A continuación, se va a presentar, resumidamente, el funcionamiento del algoritmo propuesto por Cox y Munk [23], que es la base de los demás algoritmos avanzados para sensores de baja resolución (100-1000 m), como el MODIS y el MERIS. Cox y Munk, desarrollaron la función de densidad de probabilidad para el brillo solar, basado en la velocidad del viento en la superficie marina, haciendo uso de 29 fotografías aéreas en un periodo de 20 días de una misma zona de estudio. En su trabajo, obtenían la probabilidad de pendiente especular para predecir la cantidad de brillo solar mediante la siguiente expresión, ρ       ∗󰇛,,󰇜 4cos (42) donde ρ es la reflectividad especular solar,  es la radiancia del brillo solar con dirección al sensor,  es la irradiancia solar que incide sobre la superficie,   es la reflectividad de Fresnel para el ángulo respecto a la normal , 󰇛,,󰇜 es la probabilidad de pendiente especular, teniendo en cuenta la geometría de adquisición y la tolerancia de la fuente solar en forma de elipse. Esta probabilidad depende tanto de la geometría de iluminación del Sol como de la de visión del sensor. El valor del PDF es aproximado a una función gaussiana, y su expresión hace uso de una expansión de Gram-Charlier de cuarto grado [37], dada por: 󰇛,󰇜1 2exp󰇧󰇛󰇜 2󰇨112󰇛1󰇜16󰇛3󰇜 1 24󰇛36󰇜14󰇛1󰇜󰇛1󰇜1 24󰇛36󰇜 (43) CAPÍTULO 4 57 donde  y  son pendientes normalizadas de las olas que dependen de la geometría de iluminación-visión.  representa las pendientes en la dirección del viento y  representa las pendientes en la dirección perpendicular al viento. En la ecuación el termino unitario “1” modela el comportamiento gaussiano, mientras que los múltiples términos cxx modifican su comportamiento. Dichos coeficientes son constantes o funciones que dependen de la velocidad del viento: c12 y c30 indican la asimetría, mientras que c40, c22 y c04 representan el grado de apuntamiento de la función. Gracias a que existen una relación lineal entre las pendientes normalizadas y la velocidad del viento es posible obtener el valor estadístico de la pendiente especular [38]. El método de Cox y Munk ha demostrado ser suficientemente robusto para imágenes de baja resolución siendo ampliamente utilizado, sin grandes modificaciones, hasta la actualidad. Así, sigue siendo el algoritmo base de referencia para la corrección de brillo solar de los sensores oceanográficos de baja resolución espacial, como el SeaWiFS y el MODIS, en donde se modela tanto la magnitud como la dirección del viento. Es importante resaltar que los métodos estadísticos de corrección del brillo solar no producen buenos resultados para imágenes de alta resolución (resoluciones espaciales menores a 10 metros) debido a que no hay suficiente población estadística de las diferentes pendientes del oleaje dentro de un píxel [26]. 4.2.2. Métodos para imágenes de alta resolución para aguas costeras poco profundas Los métodos de corrección para imágenes de alta resolución se basan en explotar la elevada absorción del agua en la banda NIR, lo que nos permite hacer uso de la premisa relacionada con que la reflectividad del agua a esas longitudes de onda es despreciable y, por lo tanto, la reflectividad tras la corrección atmosférica es debida únicamente al brillo especular de la superficie. Dado que la reflectividad de Fresnel puede ser considerada constante para todo el rango de longitudes de onda del visible, se puede buscar una relación numérica que vincule la proporción de brillo especular de la banda NIR con el brillo especular de cada una de las bandas del óptico. Relacionado con este concepto, Hochberg et al. [34], obtuvieron unos coeficientes lineales basados en la obtención de un pixel brillante y otro oscuro para determinar la ecuación de la recta que proporcione la pendiente utilizada en la corrección. La utilización de sólo dos píxeles hace que este método sea vulnerable a errores y ruidos al adquirir píxeles con nubes, sombras o contaminados con espuma de mar (White-cups). Por ese motivo, como se analizará a continuación, Hedley [35] y Lyzenga [36] mejoraron la obtención de estas pendientes mediante la utilización de regiones con oleaje para obtener, mediante regresión lineal y covarianza, el valor de la pendiente, permitiendo así eliminar los píxeles erróneos o contaminados.  Método de Lyzenga et al. El método de Lyzenga [36] calcula la pendiente que relaciona el brillo especular de la banda NIR con la óptica mediante el cálculo del valor de covarianza obtenido entre la banda NIR y la banda del visible para una región de interés seleccionada con los mismos requisitos que en el método de Hedley, como se analizará posteriormente. Así, 󰇛,󰇜1 󰇛󰇜󰇛󰇜1 󰇛󰇜1 󰇛󰇜       (44) Corrección del Reflejo Solar en Imágenes de Alta Resolución 58 donde  es la reflectividad superficial del mar, VIS es la banda visible, NIR es la banda del infrarrojo cercano, y N es el número de píxeles del área de interés. De esta forma, se obtiene el coeficiente como sigue, ,󰇛,󰇜 󰇛󰇜 (45) donde 󰇛󰇜󰇛,󰇜 es la varianza de los valores del NIR. Una vez obtenidos los coeficientes lineales se puede definir la expresión para la corrección del brillo solar mediante: 󰇛󰇜󰆒󰇛󰇜,󰇛󰇜󰇛󰇜            (46) donde 󰇛󰇜           es el valor medio de la reflectividad de la banda NIR para la región de interés. La pendiente ,, calculada mediante la covarianza, tiene la misma función en el algoritmo de corrección que la pendiente obtenida por regresión lineal del método de Hedley, siendo métodos equivalentes. La única diferencia es que en el método de Hedley se utiliza, como veremos seguidamente, el valor mínimo de reflectividad NIR, mientras que en el método de Lyzenga se hace uso del valor medio de la reflectividad del canal NIR.  Método de Hedley et al. En esta aproximación [35] se hace uso de una o más regiones con oleaje para obtener la escala lineal del brillo solar en la banda NIR y ópticas de la imagen a corregir. Las regiones seleccionadas han de tener un valor mínimo de reflectividad en el canal NIR, típicamente en áreas de alta profundidad y sin turbidez, y a su vez deben contener zonas de oleaje constante y continuo. Así, se obtiene una pendiente que relaciona los valores de brillo solar de los píxeles de la banda NIR con los píxeles de cada una de las bandas del óptico visible. Estos valores de pendiente sólo son válidos para la imagen procesada, dado que este valor depende de las condiciones atmosféricas en el momento de la adquisición. A su vez se asume que las condiciones atmosféricas son constantes para toda la imagen. En la Figura 18 se muestra, gráficamente, un ejemplo de curva de regresión para una región de oleaje ideal. Así, se puede observar una nube de puntos azules, los cuales representan los valores de reflectividad de la banda azul (eje de ordenadas) respecto a la reflectividad de la banda NIR (eje de abscisas). CAPÍTULO 4 59 Figura 18. Representación de la corrección del brillo solar mediante el uso de una recta de regresión. La nube de puntos tiene una forma aproximada a la de una función lineal (con una cierta variabilidad), en donde la pendiente (bi) de la recta, obtenida mediante el algoritmo de regresión lineal, es representada mediante la línea roja gruesa. El punto rojo representa un píxel con brillo especular a ser corregido Ri(VIS), mientras que el punto verde representa al pixel corregido una vez eliminado el brillo especular Ri(VIS). Podemos ver como a estos puntos les corresponde los valores de reflectividad del eje de abscisas del infrarrojo R(NIR) y Rmin(NIR). En este procedimiento, cada uno de los píxeles de la banda visible es corregido asumiendo que el valor de reflectividad en la banda NIR libre de brillo solar tiene un valor igual a min󰇛󰇛󰇜󰇜. De esta forma se reduce el valor de reflectividad en el visible mediante el valor de pendiente bi obtenido mediante regresión lineal. Asimismo, se define la siguiente expresión para la eliminación del brillo solar de la banda visible i, 󰇛󰇜󰆒󰇛󰇜󰇟󰇛󰇜min󰇛󰇛󰇜󰇜󰇠 (47) Podemos concluir que la corrección basada en regresión lineal es más robusta frente a los efectos de contaminación de los píxeles por elementos presentes en el agua como los white-cups. Aunque la eliminación del brillo solar se realiza después de la corrección atmosférica, pequeños errores en el modelado de los aerosoles pueden introducir ciertos niveles residuales de reflectividad en las bandas. Este hecho tendrá efecto en el valor mínimo de la banda infrarrojo min󰇛󰇛󰇜󰇜, pero gracias a que este valor residual es similar en las dos bandas y gracias al uso del parámetro min󰇛󰇛󰇜󰇜 en la ecuación, el valor de la pendiente permanecerá estable. En el siguiente apartado se va a describir el método de corrección de brillo solar implementado y adaptado para las imágenes WV2, basado en el método descrito de Hedley et al. Corrección del Reflejo Solar en Imágenes de Alta Resolución 60 4.2.3. Método de corrección del brillo solar para imágenes WV2 basado en el método Hedley El método de Hedley et al. es el más ampliamente utilizado en la corrección del brillo solar en imágenes de alta resolución, al hacer uso del método de regresión que permite un mejor ajuste de la recta, permitiendo visualizar y eliminar valores atípicos de píxeles ruidosos o contaminados. Por ese motivo se ha seleccionado este método en la implementación del algoritmo de deglinting para imágenes WV2. El paso más importante y crítico para la adaptación del algoritmo a las bandas del WV2 es comprobar que existe un buen ajuste entre las bandas NIR y las bandas ópticas. Una diferencia reseñable que tiene el WV2 es que dispone de ocho bandas, dos de las cuales están en el rango del infrarrojo cercano (banda 7 NIR1, y banda 8 NIR2). Por lo tanto, se va a proceder a realizar el ajuste con las dos bandas NIR debido a que el instrumento dispone de 2 grupos de sensores independientes y no sincronizados espacio-temporalmente. Para evaluar esta adaptación del modelo se ha seleccionado una imagen en donde existe un oleaje elevado con un brillo solar sobre las olas que resulta evidente (área de Granadilla, Isla de Tenerife, 22 de agosto de 2011). Se ha seleccionado una región profunda sin turbidez aparente y baja presencia de espuma de mar (white-cups) que contamine la zona de interés. El tamaño de la región de estudio es de 500 x 500 píxeles, o lo que es lo mismo 1000 x 1000 metros, por lo que se tiene suficiente variabilidad en el brillo solar debido a las diferentes pendientes del oleaje para calcular la regresión lineal en las bandas. En la Figura 19 se muestran tanto la imagen seleccionada así como el área de interés. Una vez obtenida el área de interés de la imagen se ha procedido a calcular las pendientes (bi) a partir del cálculo de regresión lineal. El método de regresión lineal es ampliamente conocido y utilizado por múltiples programas y librerías matemáticas. Para esta prueba se ha optado por utilizar la toolbox cftool de Matlab [39] que permite una fácil visualización de los resultados obtenidos. Figura 19. Imagen de Granadilla (sur de Tenerife, 22-08-2011) contaminada con brillo solar debido al elevado oleaje, y área de interés para el cálculo de la pendiente de ajuste entre banda NIR y bandas ópticas. CAPÍTULO 4 61 En las siguientes figuras se presentan los resultados obtenidos en las regresiones lineales entre las seis bandas ópticas con respecto a las dos bandas NIR del satélite WV2. Siendo b y c la pendiente y la constante de la recta, respectivamente, R2 es el nivel de correlación (cuadrado de la correlación de Pearson [40]) de la nube de puntos con la recta obtenida, y RMSE es el error cuadrático medio obtenido entre la recta de regresión y los puntos de la nube. R1 vs NIR1 R1 vs NIR2 b 0.4247 b 0.7654 c 0.05386 c 0.03432 R 2 0.3815 R 2 0.8513 RMSE 0.009928 RMSE 0.004868 Figura 20. Ajuste lineal de las bandas: azul costa R1 con la banda NIR1 R7 (izquierda), y azul costa R1 con la banda NIR2 R8 (derecha). Valores del ajuste lineal (tabla inferior). R2 vs NIR1 R1 vs NIR2 b 0.8743 b 0.6341 c 0.03918 c 0.05772 R 2 0.936 R 2 0.4323 RMSE 0.004083 RMSE 0.01216 Figura 21. Ajuste lineal de las bandas: azul R2 con la banda NIR1 R7 (izquierda), y azul R2 con la banda NIR2 R8 (derecha). Valores obtenidos del ajuste lineal (tabla inferior). Corrección del Reflejo Solar en Imágenes de Alta Resolución 62 R3 vs NIR1 R3 vs NIR2 b 0.9673 b 0.6848 c 0.01505 c 0.03742 R 2 0.9573 R 2 0.422 RMSE 0.003645 RMSE 0.0134 Figura 22. Ajuste lineal de las bandas: verde R3 con la banda NIR1 R7 (izquierda), y verde R3 con la banda NIR2 R8 (derecha). Valores del ajuste lineal (tabla inferior). R4 vs NIR1 R4 vs NIR2 b 0.6074 b 0.9848 c 0.03497 c 0.01034 R 2 0.3747 R 2 0.8036 RMSE 0.01419 RMSE 0.008153 Figura 23. Ajuste lineal de las bandas: amarillo R4 con la banda NIR1 R7 (izquierda), y amarillo R4 con la banda NIR2 R8 (derecha). Valores del ajuste lineal (tabla inferior). CAPÍTULO 4 63 R5 vs NIR1 R5 vs NIR2 b 1.007 b 0.6913 c 0.004162 c 0.0277 R 2 0.958 R 2 0.3966 RMSE 0.00374 RMSE 0.01418 Figura 24. Ajuste lineal de las bandas: rojo R5 con la banda NIR1 R7 (izquierda), y rojo R5 con la banda NIR2 R8 (derecha). Valores del ajuste lineal (tabla inferior). R6 vs NIR1 R6 vs NIR2 b 0.601 b 1.032 c 0.03016 c 0.00362 R 2 0.3665 R 2 0.722 RMSE 0.01409 RMSE 0.009332 Figura 25. Ajuste lineal de las bandas: rojo borde R6 con la banda NIR1 R7 (izquierda), y rojo borde R6 con la banda NIR2 R8 (derecha). Valores del ajuste lineal (tabla inferior). Corrección del Reflejo Solar en Imágenes de Alta Resolución 70 Figura 30. Post-procesado Histogram Matching de una imagen con alto nivel de brillo solar que ha sido corregida mediante el algoritmo de deglinting (22/08/2011) y haciendo uso de los valores estadísticos de reflectividad de una imagen sin brillo solar aparente (01/12/2011). En la Figura 31, se presenta el resultado final de la corrección del brillo solar y de los ruidos asociados al procesado de la imagen de ejemplo (Figura 27). Se puede observar como se ha logrado eliminar casi la totalidad del ruido con forma de cuadrícula, así como buena parte de los white-cups presentes en el área de interés. (a) (b) (c) (d) (e) (f) Figura 31. Seis primeras bandas del WV2, correspondientes con (a b c d e f), tras la corrección del brillo solar y la eliminación del ruido residual procedente del deglinting. CAPÍTULO 4 71 4.3. Implementación de un nuevo algoritmo automático para la eliminación del reflejo solar en imágenes multiespectrales de alta resolución Worldview-2 Como se ha analizado previamente, el método adaptado de Hedley et al., que hace uso de la relación lineal existente entre el brillo solar especular de las bandas del visible y las bandas NIR, proporciona buenos resultados en imágenes de alta resolución WV2. Sin embargo presenta importantes inconvenientes: (i) la obtención de la pendiente de ajuste que relaciona el brillo en la banda NIR y la banda óptica a corregir es una parte complicada y tediosa durante la ejecución del algoritmo; (ii) aunque esta pendiente puede ser calculada mediante algoritmos de regresión lineal, no siempre las imágenes proporcionan las condiciones necesarias para obtener estos coeficientes mediante regresión. Así, en zonas costeras de muy baja profundidad y aguas interiores continentales no se pueden encontrar zonas de alta profundidad y baja turbidez para calcular la regresión lineal y, (iii) la no existencia de oleaje genera un brillo especular casi constante para toda la imagen, que es óptimo para la corrección, pero no proporciona información suficiente de la pendiente de la recta, generando importantes errores en la obtención de la misma. Por lo anteriormente reseñado y con el objetivo de desarrollar un algoritmo automático para la eliminación del reflejo solar en las imágenes WV2, seguidamente se analiza y determina el significado físico de las pendientes de escala utilizadas en la corrección de brillo solar especular en imágenes de alta resolución. El objetivo es reemplazar la etapa empírica del algoritmo de regresión lineal por una expresión física, permitiendo así no depender de las condiciones cambiantes del oleaje y proporcionando, al mismo tiempo, un método novedoso y completamente automático. En este contexto, en la Figura 32 se muestra un modelo simplificado con las diferentes rutas, y cómo la luz procedente del Sol, que penetra en la atmósfera hacia la interface aire-agua, alcanza al sensor del satélite óptico. Concretamente:  La ruta A representa el rayo solar con dirección descendente hacia la tierra, el cual es afectado por fenómenos de absorción y difusión de la atmósfera.  La ruta B representa la parte de la luz que ha sido difundida por la atmósfera, que incide de forma especular en la superficie marina y es rebotada hacia el sensor (sky-glint).  La ruta C representa el porcentaje de luz que procedente del rayo directo incide sobre la superficie, penetrando en el agua y mediante el proceso de retro-difusión del agua parte de esta luz es enviada hacia el sensor (reflectividad del agua).  La ruta D representa al rayo solar directo que incide sobre la superficie marina y es rebotada especularmente hacia el sensor (sun-glint).  La ruta E representa la parte de la luz que ha sido difundida por la atmosfera hacia la dirección del sensor sin alcanzar la superficie del agua (debido al scattering de Rayleigh y a la presencia de aerosoles). Tras la corrección atmosférica, los fenómenos de absorción y de retro-difusión (ruta E) son compensados. Siguiendo con la hipótesis de que en la banda NIR la reflectividad es debida a la contribución especular del brillo solar, existen dos contribuciones distintas en la reflectividad especular, el sun-glint (ruta D) y el sky-glint (ruta B). Estas contribuciones proceden de la luz directa y de la luz difusa, principalmente causada por la difusión de Rayleigh. Corrección del Reflejo Solar en Imágenes de Alta Resolución 72 Figura 32. Esquema simplificado de las rutas por las cuales la luz procedente del Sol, que atraviesa la atmósfera hacia el mar, alcanza al sensor óptico del WV2. Como es conocido, la difusión de Rayleigh tiene un mayor impacto sobre las bandas del azul (de menor longitud de onda), mientras que en las bandas de rojo y el infrarrojo cercano este fenómeno es casi despreciable. Sin embargo, la reflectividad especular del brillo solar es mayoritariamente causada por la componente de la luz directa debido a que la luz difusa es dispersada isotrópicamente, debido a la función de fase del scattering de Rayleigh [47], iluminando difusamente a la superficie marina en múltiples ángulos. Por lo tanto, el brillo difuso del cielo es muy poco dependiente del estado del mar, permitiendo modelar este valor como un parámetro que depende del ángulo cenital solar y la velocidad del viento [48], mediante. ρ  ρ  󰇛 θ ,󰇜∗   (48) donde ρ es la reflectividad del reflejo solar debido a la luz difusa, ρ󰇛θ,󰇜 es el coeficiente de corrección del brillo solar difuso que depende del ángulo cenital solar y de la velocidad del viento () y   es la radiancia difusa descendente producida por la atmósfera. De esta forma, teniendo en cuenta que el brillo solar se compone de una aportación de luz directa, ecuación (42), y otra difusa, ecuación (48), y haciendo uso de los valores de reflectividad normalizada directa y difusa obtenidos del modelo 6S, podemos obtener la siguiente expresión: ρ  , ,∗ρ, ,∗ρ  , ρ, ,∗ρ , ,∗ρ  (49) CAPÍTULO 4 73 donde , , representa a la irradiancia normalizada directa, , , representa a la irradiancia normalizada difusa y , representa la reflectividad inherente del agua (water-leaving reflectance). De esta manera, si utilizamos la ecuación de la pendiente de una recta a partir de dos puntos, e introducimos las expresiones anteriores obtenemos la siguiente ecuación para la pendiente b en relación de la longitud de onda. ∗ρ  󰇛󰇜  ∗ρ  ρ  󰇛∗ρ󰇛󰇜∗ρρ 󰇜  ∗ρ 󰇛󰇜 ∗ρ  ∗ρ 󰇛󰇜 ∗ρ    ∗ρ󰇛󰇜ρ󰇛󰇜 󰇛ρ 󰇛󰇜ρ 󰇛󰇜󰇜   (50) donde ρ󰇛󰇜 es la reflectividad mínima del píxel generada por un menor número de superficies planas con pendientes orientadas especularmente y ρ󰇛󰇜 es la reflectividad máxima del píxel generada por un mayor número de superficies planas especulares. Nótese como ρρ  debido a que la reflectividad de Fresnel, de los diferentes canales ópticos, varían mínimamente debido a que el índice refractivo es muy similar para todas las bandas del visible. De esta forma podemos ver que mediante esta aproximación es posible obtener la pendiente de ajuste de cada una de las bandas, a partir de sus valores de irradiancia normalizada directa. Gracias a la utilización de un modelo atmosférico avanzado como el 6S, podemos acceder a estos valores al calcular y proporcionar este modelo la irradiancia normalizada en la superficie. En la Tabla 6 se muestran los resultados de la irradiancia normalizada directa y difusa obtenida por el modelo para la corrección atmosférica de la imagen de ejemplo (Granadilla 22/08/2011). Tabla 6. Resultados de la irradiancia normalizada directa y difusa del modelo 6S. R1 R2 R3 R4 R5 R6 R7 R8 , 0.774 0.846 0.888 0.911 0.927 0.938 0.943 0.951 , 0.226 0.154 0.112 0.093 0.073 0.062 0.057 0.049 ,  , 0.813 0.897 0.941 0.958 0.983 0.986 - - Podemos observar como las bandas con valores de irradiancia difusa corresponden a los canales de menor longitud de onda, donde el scattering de Rayleigh tiene un mayor impacto. A su vez, el cociente entre las bandas ópticas respecto a su correspondiente banda NIR es un valor creciente hacia la unidad como sucede en los parámetros de pendiente de la regresión lineal (ver Tabla 5). En la Figura 33 se muestran los resultados obtenidos para las pendientes de escala mediante el uso del método de regresión lineal y mediante el nuevo método físico del cociente entre la irradiancia normalizada directa. Corrección del Reflejo Solar en Imágenes de Alta Resolución 74 Figura 33. Comparativa entre el cálculo de la pendiente mediante regresión lineal y la pendiente mediante el cociente de irradiancias normalizadas directas para cada banda visible del WV2. Se puede observar como los resultados obtenidos mediante el método físico del cociente de las irradiancias proporcionan valores muy similares a los obtenidos mediante el método de regresión lineal. La diferencia cuadrática media entre dichos resultados se sitúa en torno a 0.01. La desviación más elevada se encuentra en el canal azul costa, en donde por un lado, es el canal más afectado por el scattering de Rayleigh y, por lo tanto, el modelo puede cometer mayores errores y, por otro lado, es una banda del sub-sensor MS2, el cual proporciona un menor ajuste en el cálculo de regresión por lo que puede introducir cierto error en su cálculo. El siguiente error más elevado puede observarse en la banda rojo borde, en donde el método de regresión lineal proporciona un resultado ligeramente superior a uno (1.032), mientras que en el canal rojo con una longitud de onda muy similar la pendiente de regresión es prácticamente 1. Teniendo en cuenta que en la banda rojo borde (MS2) se obtiene el peor resultado de correlación es probable que éste sea el causante del incremento en la diferencia entre los dos métodos. Por lo tanto, haciendo uso del cociente entre irradiancias normalizadas directas se presenta la siguiente expresión que permite la corregir el brillo solar especular de las imágenes mediante una aproximación física: 󰇛󰇜󰆒󰇛󰇜, , , ,∗󰇛󰇜 (51) En la ecuación se ha eliminado el valor de reflectividad mínima de la banda NIR dado que no ha sido utilizada en ningún momento para el cálculo de la pendiente y debido a que este valor residual de reflectividad en la banda NIR es mayormente causada por errores en la configuración del modelo atmosférico 6S, al utilizarse valores de AOD procedentes de imágenes de baja resolución espacial y debido a ciertas aproximaciones en el modelado de los aerosoles. De esta forma el nuevo método para la eliminación del brillo solar permite una ejecución completamente automática del procesado de la imagen. CAPÍTULO 4 75 4.3.1. Resultados obtenidos en la corrección del brillo solar especular en imágenes WV2 Como se ha podido apreciar en la Figura 33, los resultados obtenidos en la pendiente (bi) son muy similares tanto para la adaptación del método de Hedley et Al. para imágenes WV2, ecuación (47), como para el nuevo algoritmo automático de corrección de brillo solar propuesto en esta Tesis, ecuación (51). A continuación se presentan los resultados obtenidos en la corrección del brillo solar especular de las imágenes WV2 generadas por el nuevo algoritmo automático de corrección de brillo solar, en donde se ha hecho uso de las funciones de post-procesado previamente descritos. Para la evaluación del funcionamiento del algoritmo desarrollado para la eliminación del reflejo solar (glinting) se han procesado imágenes con alta contaminación por brillo solar, en las diferentes áreas bajo estudio, con distintas condiciones de iluminación, según la estación del año.  Área de Granadilla (Tenerife) En la Figura 34 se muestra el resultado del algoritmo de deglinting para la imagen WV2 del 22 de agosto 2011 del área bajo estudio. Para ello se muestra, en la parte superior, una composición RGB de la imagen de reflectividad corregida atmosféricamente y con alta contaminación de brillo solar y, en la parte inferior, la composición RGB de la imagen de reflectividad una vez eliminado el brillo solar. Obsérvese que los valores de contraste en el visualizador son idénticos en las dos imágenes, sin embargo, la reflectividad del agua en la imagen (a) es muy superior a la imagen corregida (b). Se puede observar que el color del agua en las zonas profundas tiene una tonalidad mucho más azulada en la imagen corregida. Así, en la imagen corregida los detalles de la reflectividad del fondo marino, en aguas poco profundas, así como plumas de turbidez, cerca de la costa, presentan mayor contraste. Estas características a analizar posteriormente, fondo marino y calidad del agua, son difíciles de distinguir, aunque sea visualmente, en la imagen no corregida. Finalmente, resaltar la máscara de tierra en color gris, en la imagen corregida, que nos permite eliminar las áreas terrestres de la imagen. Este enmascaramiento se realiza en el paso de eliminación de valores atípicos en el proceso de eliminación de ruidos residuales de la imagen. (a) Corrección del Reflejo Solar en Imágenes de Alta Resolución 76 (b) Figura 34. Imagen de Granadilla 22 de agosto de 2011: (a) reflectividad superficial con alta contaminación de brillo solar, y (b) reflectividad superficial sin brillo solar tras su corrección.  Área de Maspalomas (Gran Canaria) En la Figura 35 se muestra el resultado del algoritmo de deglinting para la imagen WV2 del 9 de mayo 2012 de la zona sureste de Gran Canaria (área de Maspalomas). En la Figura 35 (a) se proporciona una composición RGB de la imagen de reflectividad corregida atmosféricamente y con alta contaminación de brillo solar y, en la Figura 35 (b), la composición RGB de la imagen de reflectividad una vez eliminado el brillo solar. Al igual que el caso de estudio previo, se puede observar que el color del agua en las zonas profundas tiene una tonalidad mucho más azulada en la imagen corregida, evidenciando los detalles de la reflectividad del fondo costero. Finalmente, son observables diferentes estructuras generadas por la variación del viento, añadidas a las formas del oleaje en la Figura 35 (a), que son parcialmente eliminadas en la imagen corregida, pero aún son apreciables en la Figura 35 (b). Nótese que esta imagen ha sido tomada con condiciones de viento muy elevados apareciendo conjuntamente al oleaje una gran cantidad de espuma marina. CAPÍTULO 4 77 (a) (b) Figura 35. Imagen de Maspalomas 9 de mayo 2012: (a) reflectividad superficial con alta contaminación de brillo solar, y (b) reflectividad superficial sin brillo especular solar.  Área de la Restinga (El Hierro) Finalmente, en las Figura 36 (a) y (b) se muestra el resultado del algoritmo de deglinting para una imagen WV2 de la zona de La Restinga en la isla de El Hierro (imagen del 2 de marzo 2012), adquirida en la finalización del proceso eruptivo submarino acontecido en octubre de 2011. Se puede observar como en la imagen corregida sin brillo solar, Figura 36 (b), los detalles del resto de la pluma de turbidez (zona inferior derecha de la imagen) procedente de la erupción del volcán submarino comienzan a ser más evidentes respecto a la imagen contaminada de brillo solar, Figura 36 (a), en la que no se puede apreciar la ligera mancha verde en el mar. Se puede observar como en la imagen corregida se ha enmascarado tanto el área terrestre como las nubes de la imagen. Corrección del Reflejo Solar en Imágenes de Alta Resolución 78 (a) (b) Figura 36. Imagen de La Restinga 2 de marzo de 2012: (a) reflectividad superficial con alta contaminación de brillo solar, y (b) reflectividad superficial sin brillo solar tras su corrección. CAPÍTULO 4 79 4.4. Resumen En el presente capítulo se ha tratado el fenómeno físico de la reflexión especular del Sol en la superficie marina, describiéndose los principales métodos de eliminación del brillo solar en imágenes de teledetección. Por un lado, los algoritmos para imágenes de baja resolución utilizan una aproximación estadística para cuantificar el brillo solar, según las condiciones de iluminaciónvisión en la adquisición de la imagen, así como para las condiciones de oleaje, gracias a la información de la magnitud y dirección del viento. Por otro lado, los algoritmos para imágenes de alta resolución, se basan en la hipótesis de que la reflectividad propia del agua en las bandas NIR es despreciable, por lo que la reflectividad presente en dicha banda es debida al brillo solar. Esta aproximación utiliza este valor de reflectividad y un parámetro de escala que relaciona el brillo solar en la banda NIR con la banda visible para eliminar la aportación del brillo solar en la banda visible, siendo necesario el cálculo de las pendientes de ajuste, las cuales pueden ser calculadas mediante regresión lineal (Hedley et al.) o mediante el cálculo de la covarianza (Lyzenga et al.). A continuación, se ha descrito la metodología para la corrección de imágenes del WV2 mediante la adaptación del método de Hedley a las bandas multiespectrales del satélite, incluyéndose una etapa de post-procesado para eliminar el ruido residual entre diferentes bandas de la imagen. En dicho post procesado se han realizado ajustes de correlación de las bandas mediante técnicas de ventana deslizante, se han eliminado valores de píxeles atípicos mediante el uso de un umbral, se ha realizado técnicas de inpainting para el rellenado de los píxeles eliminados y se ha introducido un paso de filtro gaussiano bilateral para el suavizado del ruido de altas frecuencias respetando la información de la imagen de bajas frecuencias. Indicar que la adaptación del algoritmo Hedley ha requerido del conocimiento del funcionamiento de las bandas del WV2, dado que las ocho bandas trabajan agrupadas en dos sub-sensores (MS1 y MS2), los cuales no están correlacionados temporalmente entre sí. Esta metodología ha sido publicada, además de en diferentes congresos nacionales e internacionales, en el artículo: “High-resolution maps of bathymetry and benthic habitats in shallowwater environments using multispectral remote sensing imagery,” IEEE Transactions on Geoscience and Remote Sensing 2015. Posteriormente, se ha presentado un nuevo algoritmo de deglinting que permite calcular la pendiente de ajuste mediante una aproximación física sobre la aportación de la luz directa de cada canal respecto al brillo especular. Así, se ha obtenido la pendiente como cociente entre la irradiancia normalizada directa del canal visible y la irradiancia normalizada directa del canal NIR. Gracias al uso del modelo atmosférico 6S es posible obtener esta información generada en el proceso de corrección atmosférica. El nuevo algoritmo propuesto, basado en una aproximación física, permite la obtención de la pendiente de ajuste y es una alternativa muy interesante a los métodos anteriores dado que, aunque proporciona resultados muy similares, permite ser utilizado de una manera completamente automática y en todo tipo de zonas de estudio, como aguas continentales, en donde no se puede obtener una zona óptima para el cálculo de regresión lineal. Además, el cálculo de regresión lineal requiere de un oleaje mínimo que permita proporcionar información suficiente de la pendiente de la recta. Si no existe este oleaje los datos se concentrarán en un solo punto generando errores importantes en el cálculo de la pendiente y una correlación muy baja. El nuevo método es idóneo para este tipo de situaciones dado que no requiere de esta información. A su vez, el hecho de proporcionar un método completamente automático permite automatizar completamente la cadena de procesamiento sin necesidad de intervención humana, siendo idóneo en entornos en donde se han de procesar un gran número de imágenes. Finalmente, se han presentado los resultados del procedimiento de eliminación del brillo solar para imágenes WV2 obtenidos en diferentes zonas de estudio con diferentes condiciones de Determinación de Parámetros de Calidad del Agua, Batimetría y Albedo del Fondo Costero 86 Figura 40. Atenuación de los restos (detritus) dependiendo de . (Fuente: Lee-IOPs_Lecture-1). Figura 41. Espectro de absorción del plancton según su tamaño. (Fuente: Ciotti et al. 2002). Figura 42. Absorción del fitoplancton según su especie o grupo funcional. (Fuente: Carr et al. 2006). CAPÍTULO 5 87 Figura 43. Espectro de absorción de los diferentes pigmentos del fitoplancton. (Fuente: Bricaud et al. 2004). En la Figura 43 se muestra el espectro de absorción de los diferentes pigmentos fotosintéticos presentes en el fitoplancton [58], pudiéndose observar una multitud de pigmentos fotosintéticos como la chl-a (el más predominante y conocido), chl-b y chl-c. Es posible modelar un espectro de absorción promediado según el tamaño medio, especies y contenido de pigmentos promedios del fitoplancton. Para ello existen varios métodos que utilizan uno, dos, o múltiples parámetros, siendo el método más conocido y utilizado el propuesto por Lee et al. [59], mediante el uso de dos parámetros, dado por: 󰇛󰇛󰇜󰇛󰇜ln󰇛󰇜󰇜∗ (58) donde 0 󰇛  󰇜 y 1 󰇛  󰇜 son los dos parámetros dependientes de la longitud de onda, mientras que  es el valor de la absorción del fitoplancton a los 440 nm (󰇛440󰇜). Por lo tanto, se puede modelar la absorción del fitoplancton mediante un solo valor de referencia. En la Tabla 7 se muestran los coeficientes de absorción del fitoplancton propuestos por Lee et al., para el modelado del espectro de absorción del fitoplancton. Finalmente, en la Figura 44 se muestran los espectros de absorción del fitoplancton según la simulación del método propuesto por Lee et al. para diferentes concentraciones del chl-a. Tabla 7. Coeficientes de absorción del fitoplancton. (Fuente: Lee et al. 1998). Determinación de Parámetros de Calidad del Agua, Batimetría y Albedo del Fondo Costero 88 Figura 44. Simulación del espectro de absorción del fitoplancton según diferentes concentraciones de chl-a. (Fuente: Lee et al. 1998).  Atenuación de la materia disuelta en el agua La absorción por la materia disuelta en el agua (gelbstoff) ha sido modelada gracias a múltiples pruebas experimentales, constatándose que la dependencia espectral de la atenuación puede ser descrita mediante una función exponencial dada por [60], 󰇛󰇜󰇛󰇜∗󰇛  󰇜 (59) donde el parámetro  describe el grado de decrecimiento de la función exponencial según la longitud de onda y puede tener valores que oscilan entre 0.01 y 0.03 nm-1. El parámetro 󰇛󰇜 indica la magnitud de la atenuación a una longitud de onda de referencia. En la Figura 45 se muestra los diferentes espectros de absorción producidos por la presencia de materia disuelta en el agua para diferentes localizaciones, según las concentraciones de materia disuelta [61]. Figura 45. Absorciones producidas por la materia amarilla en diferentes lugares de test. (Fuente: Kirk 1994). CAPÍTULO 5 89 Podemos observar que el comportamiento espectral es idéntico a la absorción debida a los detritos, en donde solo varía el rango del parámetro Sg el cual alcanza un valor algo más elevado. Debido a que los espectros de absorción de los dos parámetros son muy similares se suele unificar dichos parámetros utilizando un valor de exponente que es un compromiso entre estos dos parámetros (típicamente 0.015). El valor de  utilizando como valor de referencia se fija a 440 nm y se define el parámetro G como el valor de 󰇛440󰇜 mediante, 󰇛󰇜∗.󰇛󰇜 (60)  Comparativa entre los diferentes valores de absorción del agua En la Figura 46 se pueden observar tanto las diferentes contribuciones de atenuación en el espectro del visible, para un entorno oceánico con baja concentración de IOPs, como para un entorno costero con mayores contribuciones de materia suspendida, disuelta y mayores niveles de clorofila. Donde a_tot representa la suma de las diferentes absorciones. Se puede observar como la atenuación del agua pura son idénticas en las dos figuras, siendo el factor de atenuación dominante en aguas oceánicas limpias. Sin embargo, en aguas costeras, con altos contenidos de materia disuelta y suspendida, la absorción en el rango visible cercano al azul y el verde se eleva de una manera evidente. (a) (b) Figura 46. Contribución a la absorción de los diferentes IOPs, según su concentración: (a) entorno oceánico y (b) entorno costero. (Fuente: Lee-IOPs_Lecture-1). 5.1.1.2. Propiedades de difusión de los IOPs Al igual que la absorción, las propiedades de difusión del agua pueden ser calculadas como la suma de las difusiones de los diferentes elementos del agua. Existen innumerables elementos que producen difusión en el agua [51]. Sin embargo, estos parámetros pueden ser clasificados de una forma más funcional y práctica como la suma de la difusión del agua pura (), debido al scattering de Rayleigh, y la debida a partículas suspendidas ().  (61) Determinación de Parámetros de Calidad del Agua, Batimetría y Albedo del Fondo Costero 90 A su vez, la difusión de las partículas se compone por partículas inorgánicas  (Particulate Inorganic Matter) y por partículas orgánicas  (Particulate Organic Matter) con una respuesta espectral a la difusión muy similar [62].  (62) En la Figura 47 se muestra un diagrama con los diferentes elementos existentes en el agua. Se puede observar como hasta las 0.2 micras se identifican a los elementos como disueltos, mientras que si superan este tamaño son considerados como materia suspendida, la cual comienza a tener propiedades de difusión. Las bacterias están en el límite de las 0.2 micras, mientras que el fitoplancton se encuentra dentro de la materia suspendida. Se puede observar como existe una gran variedad de detritus, partículas inorgánicas y orgánicas, con una gran variedad de tamaños. Hay que recalcar que el fenómeno de difusión es fuertemente dependiente del tamaño y la forma de las partículas, generando mayor difusión las partículas de mayor tamaño, mientras que la forma de las partículas incide en su espectro de difusión. Figura 47. Diagrama esquemático de los componentes presentes en el agua marina y su tamaño. (Fuente: Stramski 2004).  Difusión del agua marina pura La difusión del agua marina pura se debe al scattering molecular tipo Rayleigh el cual tiene una función volumétrica de dispersión VSF constante, desde el punto de vista angular. Morel et al. (1974) y Shifrin et al. (1988) [63] modelaron el comportamiento de este parámetro mediante: 450 . (63) donde  es el valor de referencia a una longitud de onda de 450 nm. La retro-difusión del agua () es el 50% del valor de dispersión (), pudiéndose modelar espectralmente con la ecuación (65), propuesto por Morel et al. (1974) y Zhang et al. (2009) [64]. CAPÍTULO 5 91 0.5∗ (64) 󰇛󰇜0.002450 . (65)  Difusión de las partículas presentes en el agua El scattering generado por las partículas suspendidas en el agua es de una naturaleza muy heterogénea y dependiente de su tamaño y forma. El grado de difusión hacia delante respecto a la retro-difusión es muy elevado en contraposición al caso del agua pura donde este valor es idéntico. De esta forma, se puede definir   como la proporción de retro-difusión generado en el fenómeno de difusión,     (66) Para el caso del agua pura el cociente   es igual a 0.5, mientras que para las partículas disueltas en el agua el valor   varía típicamente entre 0.005 y 0.05. Múltiples autores han intentado modelar la función volumétrica de dispersión de las partículas suspendidas en el agua marina [65] [66], aunque estos modelos son muy complejos y dependen de un gran número de parámetros que no suelen estar disponibles para este tipo de aplicaciones. Sin embargo, es posible obtener un valor promediado de back-scattering generado por las partículas en suspensión mediante la siguiente expresión [67]: 󰇛󰇜  󰇛󰇜400 . (67) donde  es el valor de referencia de back-scattering para la longitud de onda , e  es el valor del exponente que modela la función. Dicho valor varía típicamente entre 0 y 2. Un valor muy utilizado para la referencia de los 400 nm es  = 1.7, pudiéndose definir el parámetro  como 󰇛400󰇜.  es el encargado de modelar el espectro de retro-difusión, siendo afectado por el tamaño y la forma de las partículas. Normalmente, cuando las partículas son muy pequeñas,  alcanza valores elevados, mientras que si las partículas son de gran tamaño  tiende a valores cercanos a 0.  Comparativa entre las diferentes valores de back-scattering del agua En la Figura 48 se muestran las diferentes contribuciones de back-scattering en el espectro del visible, para un entorno oceánico con baja concentración de partículas, y para un entorno más costero con mayor contribución de materia suspendida. Como se puede observar, en aguas oceánicas la mayor contribución de back-scattering es de origen molecular del agua, con un valor más elevado en las bandas del azul, confiriendo el color azulado típico del agua. Mientras que en entornos costeros el mayor contribuyente es la retro-difusión de las partículas que tiene un espectro más plano. Determinación de Parámetros de Calidad del Agua, Batimetría y Albedo del Fondo Costero 92 Figura 48. Back-scattering generado por el agua pura (azul), y por la materia suspendida en el agua (roja aguas oceánicas) y (verde aguas costeras). (Fuente: Lee-IOPs_Lecture-1). 5.1.1.3. Reflectividad inherente del agua calculada a partir de los IOPs y las condiciones de contorno Las condiciones de contorno hacen referencia a cuáles son las condiciones de iluminación del Sol y de visión del satélite en el instante de adquisición de la imagen. Estas condiciones de contorno son, conjuntamente con la contribución de absorción y retro-difusión de los IOPs, las responsables de la reflectividad del agua. La reflectividad del agua calculada mediante los IOPs se basan en el cálculo numérico de la transferencia radiativa mediante el parámetro  designado por Gordon [68], y dado por,    (68) Podemos observar como el parámetro  depende directamente de la cantidad de back-scattering generado en el agua, e inversamente de la cantidad de atenuación difusa generada, como la suma de la atenuación y de la retro-difusión. Haciendo uso del parámetro  múltiples autores han modelado la reflectividad del agua profunda, cuando no existe interacción con el fondo marino. El primero fue Gordon en 1988 [69], correspondiente con la ecuación (69), siendo actualizada por Lee et al. [67] mediante la ecuación (70). Debido a que estas ecuaciones no tienen en cuenta las dependencias de los ángulos de iluminación del Sol y visión del satélite, actualmente [70] [71], se han implementado unas nuevas expresiones que modelan la geometría de adquisición. Estas ecuaciones (71) y (72) vienen dadas por: 󰇛0.09490.0794󰇜 (69) 󰇛0.08400.017󰇜 (70) 0.0512󰇛14.66597.83875.4571󰇜∗10.1098  ∗10.4021   (71) 0.287410.28211.0190.4561∗󰇭10.4021  1.4021 󰇮 (72) CAPÍTULO 5 93 donde  es la reflectividad por debajo de la superficie del agua medida, con respecto a la unidad del ángulo sólido sr-1, obtenida mediante el sensor de teledetección (remote sensing reflectance).  y  representan a los cosenos de los ángulos de incidencia solar y de visión del satélite, justo por debajo de la superficie del agua, modificados por la refracción según la ley de Shell, dada por la ecuación (73). 󰇧󰇛󰇜 1.34󰇨 (73) donde la constante 1.34 representa el índice de refracción del agua marina. A su vez, se define  como la reflectividad en la parte superior de la superficie del agua medida con respecto a la unidad del ángulo sólido sr-1 obtenida mediante el sensor de teledetección (remote sensing Reflectance). La conversión entre ambas reflectividades, parte inferior y superior de la superficie marina, requiere tener en cuenta tanto el fenómeno de reflexión-refracción como la reflexión total interna de la luz retro-difundida con dirección a la superficie que alcanza la interfaz agua-aire con un ángulo superior al ángulo límite (típicamente 48.6º para el agua marina). En la Figura 49 se muestra el fenómeno de reflexión total interna al pasar del medio de mayor índice de refracción (agua) a otro medio con menor índice (aire). Varios autores han modelado la reflectividad  a partir del valor de reflectividad  mediante aproximaciones numéricas. Inicialmente, Mobley et al. [72] calculó un modelo simplificado para condiciones de iluminación especulares, siendo posteriormente actualizado por Lee et al. [73]. A continuación, Loisel [74] implementó un modelo más complejo, considerando el ángulo de incidencia de la luz solar. Las diferentes expresiones citadas son: Figura 49. Reflexión total interna del agua. Determinación de Parámetros de Calidad del Agua, Batimetría y Albedo del Fondo Costero 94 0.54∗ (74) 0.52∗ 11.17∗ (75) 󰇛󰇜∗ 1󰇛󰇜∗ (76) donde 󰇛󰇜󰇟0°0.5236,30°0.5169,0°0.4933󰇠 y 󰇛󰇜󰇟0°2.1941,30°2.3001, 0°2.6796󰇠 son coeficientes que depende del ángulo de incidencia solar. Dichos resultados están basados en el simulador comercial hydrolight [75], el cual realiza el cálculo de inversión de IOPs y el uso de los datos del IOCCG informe 5 [76]. Finalmente, la relación que existe entre la reflectividad superficial obtenida mediante teledetección  y el parámetro de reflectividad del agua superficial  viene dada por la ecuación (77), presuponiendo un comportamiento lambertiano de la superficie (Figura 50), donde la radiancia incidente en una superficie plana es dispersada por igual por todo el casquete esférico que representa el ángulo sólido.   (77) Gracias al modelado del comportamiento de los IOPs y de las condiciones de contorno de iluminación y visión, como se ha analizado previamente, se puede modelar la reflectividad superficial del agua debido a sus propiedades inherentes, o al contrario, se pueden inferir los valores de los principales parámetros de los IOPs que corresponden con la reflectividad obtenida por un sensor remoto. Esto último es la base en la que se apoya el modelado de transferencia radiativa para invertir las ecuaciones que permiten obtener la reflectividad, a partir de los IOPs, en pos de calcular los parámetros intrínsecos asociados al agua marina. Como se ha descrito anteriormente, es posible obtener la reflectividad inherente al agua marina en condiciones normales, además de los valores constantes de atenuación y retro difusión del agua pura, mediante tres parámetros: la atenuación debida a la materia disuelta y a los detritos suspendidos en el agua G; la atenuación del fitoplancton P, y la retro difusión generada por las partículas suspendidas en el agua X. Figura 50. Suposición de superficie lambertiana de la reflectividad inherente del agua en la superficie. CAPÍTULO 5 95 En la Figura 51 se muestran diferentes tipos de reflectividades del agua según el modelado de los tres principales parámetros IOPs que son las responsables del color inherente del agua, para unas condiciones de iluminación y visión cenitales. Figura 51. Diferentes colores (RGB) del agua marina modelados según los parámetros GPX de los IOPs. (Fuente: http://www.exo.net/~pauld/colorofwater/ColorofWater.html). 5.1.1.4. Métodos de inversión del RTM basado en IOPs Los modelos analíticos o semi-analíticos que estudian la composición de los elementos presentes en el agua, mediante el modelado del comportamiento óptico de los elementos bioquímicos, requieren de la implementación de algoritmos complejos para el cálculo mediante inversión numérica de estos parámetros a partir de los valores de observación proporcionados por las bandas del sensor remoto. Las dos principales estrategias para la implementación de estos algoritmos son la de Bottom Up Strategy (BUS), y la Top Down Strategy (TDS). Para la implementación de la estrategia BUS, ampliamente utilizada, se asume que conocemos perfectamente el comportamiento espectral de los parámetros envueltos en los IOPs, teniéndose que calcular todos los parámetros simultáneamente. Mientras que para la estrategia TDS, el modelado espectral de los parámetros sólo es utilizado cuando es necesario en el modelo y se hace uso, normalmente, de expresiones empíricas obtenidas mediante relaciones entre bandas para obtener ciertos parámetros, permitiendo un cálculo menos complejo del problema al no necesitar resolver todos los parámetros de forma simultánea. En la Figura 52 se muestra el funcionamiento de las dos estrategias respecto al cálculo de los parámetros necesarios en la estrategia BUS y TDS. Figura 52. Estrategias BUS y TDS para el cálculo de los parámetros inherentes del agua IOPs. Determinación de Parámetros de Calidad del Agua, Batimetría y Albedo del Fondo Costero 102 Teniendo en cuenta que la atenuación de la luz en el agua es exponencial, este algoritmo realiza una linealización mediante la utilización del cociente del logaritmo de dos bandas y la profundidad. Teniendo en cuenta que la atenuación debida al agua es mucho mayor que la del albedo del fondo y las propiedades del cociente del logaritmo de las bandas, el autor afirma que es de esperar una variación mínima en el ajuste lineal debido a la variación del albedo marino. La ecuación utilizada para el cálculo de batimetría, mediante el ratio algorithm, es la siguiente: ln󰇛󰇜 ln󰇛󰇜 (90) donde  es el parámetro utilizado para ajustar el cociente a la variación en la batimetría,  es la constante utilizada para ajustar la profundidad inicial a 0 m y  es un valor fijo que permite asegurar que el logaritmo va a resultar positivo para cualquier valor de reflectividad (típicamente 1000). De esta manera, mediante la utilización de dos bandas de alta penetración, normalmente azul y verde, se pueden obtener los parámetros  y  mediante regresión lineal entre el valor del cociente del logaritmo de las bandas y un mapa de batimetría conocido. En la Figura 55 se muestra el ajuste obtenido por Stumpf et al. para las bandas azul y verde, realizado para el satélite IKONOS. Como se puede observar, este algoritmo no tiene en cuenta el coeficiente de atenuación debido a los IOPs, el cual puede variar mucho en aguas costeras modificando en gran medida los resultados de las constantes de ajuste. Tampoco tiene en cuenta los ángulos de incidencia y visión que influyen en la distancia recorrida por la luz en el agua. Sin embargo, dicho algoritmo proporciona resultados aceptables con errores reducidos para aguas poco profundas y en condiciones similares a los del ajuste. Figura 55. Ajuste de los parámetros lineales del algoritmo del cociente. (Fuente: Stumpf [99]). CAPÍTULO 5 103 5.1.3.3. Inversión del modelo de transferencia radiativa en aguas costeras Como se ha descrito anteriormente, las ecuaciones radiativas para aguas costeras modelan la reflectividad inherente del agua mediante los IOPs, así como la reflectividad procedente del albedo del fondo, siendo la profundidad el parámetro de ajuste entre una y otra fuente de reflectividad. De esta forma se puede reescribir el sistema de ecuaciones presentado en la ecuación (78), el cual relaciona las variables del modelado de transferencia radiativa, de la siguiente manera.         󰇛󰇜󰇛󰇜,󰇛󰇜,󰇛󰇜,󰇛󰇜,󰇛󰇜,,󰇛󰇜󰇛,,,,󰇜  󰇛󰇜󰇡󰇛󰇜,󰇛󰇜,󰇛󰇜,󰇛󰇜,󰇛󰇜,,󰇛󰇜󰇢󰇛,,,,󰇜 ...  󰇛󰇜󰇡󰇛󰇜,󰇛󰇜,󰇛󰇜,󰇛󰇜,󰇛󰇜,,󰇛󰇜󰇢󰇛,,,,󰇜 (91) donde la reflectividad del  puede ser modelado mediante un único valor de reflectividad normalizada a 555 nm, o bien puede ser modelado mediante una mezcla lineal de dos o tres EndMembers. La resolución de los valores de abundancia de la mezcla lineal ha de ser resuelto mediante optimización por mínimos cuadrados en el sistema de ecuaciones, al igual que las demás incógnitas. Por lo que el nuevo número de incógnitas pasa a ser de 5 o 6, según el modelado del albedo costero realizado. 5.2. Modelo de transferencia radiativo mejorado para aguas costeras y adaptado a los canales multiespectrales WV2 Teniendo en cuenta las limitaciones descritas y el número de bandas multiespectrales disponibles, dentro del rango óptico, será necesario el desarrollo e implementación de un modelo robusto que permita la utilización de las ecuaciones de transferencia radiativa para imágenes multiespectrales de alta resolución espacial. Es importante resaltar, como científicamente se constata, que el modelado de transferencia radiativa fue ideado, inicialmente, para sensores hiperespectrales con una multitud de bandas monocromáticas [59] [67] [100] [101] [102] [103]. Por este motivo, el mayor desafío en la implementación del RTM para imágenes WV2 es la adaptación del problema de la resolución de las ecuaciones de transferencia radiativa a las necesidades multiespectrales y la configuración de los parámetros a modelar, teniendo en cuenta la limitación del número de ecuaciones disponibles para la resolución del sistema de ecuaciones de transferencia radiativa. De esta manera, se comenzó con la implementación del diagrama clásico de bloques aislados interconectados entre sí [104] [105] [106]. Estos bloques representan (i) la corrección atmosférica, (ii) corrección de brillo solar y (iii) el modelado de transferencia radiativo, como se muestra en la Figura 56. Sin embargo, el hecho de hacer uso de imágenes de alta resolución espacial, en donde se tiene acceso a áreas de muy baja profundidad y alta concentración de IOPs, ha hecho necesaria la introducción de ciertas modificaciones a esta cadena de procesamiento estándar, intentando así mitigar desajustes generados en los diferentes bloques de la cadena de procesamiento. A continuación, se procede a describir detalladamente la metodología de adaptación y resolución de las ecuaciones de transferencia radiativa para la obtención de los parámetros de calidad del agua, batimetría y albedo del fondo con imágenes multiespectrales WV2, para diferentes entornos costeros. Determinación de Parámetros de Calidad del Agua, Batimetría y Albedo del Fondo Costero 104 Figura 56. Cadena estándar de procesado con módulos separados. 5.2.1. Adaptación del modelo de transferencia radiativa a las bandas multiespectrales WV2 Como se ha descrito en el capítulo 2, las bandas multiespectrales del satélite WV2 integran la respuesta espectral de la radiancia recibida para unos anchos de banda en torno a 50 nm. Consecuentemente, los procesos no lineales a los que es sometida la luz, tanto en la atmósfera como en el agua, son integrados mediante el filtro paso banda del sensor. El comportamiento no lineal hace necesario modelar la reflectividad de cada longitud de onda y, posteriormente, integrarlo mediante su producto con la función normalizada de respuesta del filtro paso banda. Este modelado requeriría el cálculo de unas 50 ecuaciones de transferencia radiativa para cada banda, lo que ralentizaría en exceso el cómputo del modelo. Por ese motivo, al igual que para el modelo 6S, se va a proceder a calcular el modelo para el rango de longitudes de onda con un paso de 5 nm, aproximación correcta debido a la variación suave de la respuesta espectral, permitiéndonos acelerar el cálculo del modelado radiativo. La integración de los resultados de transferencia radiativa monocromáticas, de anchos de banda de 1 nm en pasos de 5 nm, para la respuesta de las bandas multiespectrales WV2, vendrá dada por:  󰇛󰇜  󰇛󰇜∗5     (92) donde  representa el número de banda (1 a 8) dentro de las bandas del WV2,  󰇛󰇜 constituye el resultado de la ecuación radiativa para la banda ,   representa la longitud de onda inicial del ancho de banda para la banda ,   representa la longitud de onda final del ancho de banda para la banda  y 5 representa la función normalizada de respuesta del filtro paso banda, paso 5 nm, para la longitud de onda  determinada. CAPÍTULO 5 105 De esta forma se puede reformular la ecuación (91), para las ocho ecuaciones multiespectrales del WV2, como se muestra en la ecuación.                       󰇛󰇜󰇛󰇜,󰇛󰇜,󰇛󰇜,󰇛󰇜,󰇛󰇜,,󰇛󰇜 󰇛,,,,󰇜  󰇛󰇜󰇡󰇛󰇜,󰇛󰇜,󰇛󰇜,󰇛󰇜,󰇛󰇜,,󰇛󰇜󰇢 󰇛,,,,󰇜  󰇛󰇜󰇡󰇛󰇜,󰇛󰇜,󰇛󰇜,󰇛󰇜,󰇛󰇜,,󰇛󰇜󰇢 󰇛,,,,󰇜  󰇛󰇜󰇡󰇛󰇜,󰇛󰇜,󰇛󰇜,󰇛󰇜,󰇛󰇜,,󰇛󰇜󰇢 󰇛,,,,󰇜  󰇛󰇜󰇡󰇛󰇜,󰇛󰇜,󰇛󰇜,󰇛󰇜,󰇛󰇜,,󰇛󰇜󰇢 󰇛,,,,󰇜  󰇛󰇜󰇡󰇛󰇜,󰇛󰇜,󰇛󰇜,󰇛󰇜,󰇛󰇜,,󰇛󰇜󰇢 󰇛,,,,󰇜  󰇛󰇜󰇡󰇛󰇜,󰇛󰇜,󰇛󰇜,󰇛󰇜,󰇛󰇜,,󰇛󰇜󰇢 󰇛,,,,󰇜  󰇛󰇜󰇡󰇛󰇜,󰇛󰇜,󰇛󰇜,󰇛󰇜,󰇛󰇜,,󰇛󰇜󰇢 󰇛,,,,󰇜 (93) donde ∑󰇛󰇜  representa la integración del ancho de banda multiespectral de las bandas WV2 para el rango de longitudes de onda del canal . Como se puede observar, podemos modelar la reflectividad sobre la superficie marina de las ocho bandas multiespectrales del WV2, sin embargo, no todas las bandas proporcionan la misma cantidad ni tipo de información: i. Las primeras dos bandas, azul costa y azul (1, 2), son las que proporcionan mayor cantidad de señal, siendo muy sensibles a la atenuación producida por el fitoplancton y materia disuelta, así como al albedo del fondo marino. ii. La banda verde (3) es menos sensible a la absorción del fitoplancton y la materia disuelta, pero aún tiene una alta penetración en el agua por lo que es sensible al albedo costero. iii. La banda amarilla (4) es mucho menos sensible a la atenuación debida al fitoplancton y la materia disuelta, permitiendo ser detectadas variaciones de la retro-difusión debida a la materia suspendida. La influencia del albedo del fondo comienza a ser más tenue. iv. La banda roja (5) proporciona información sobre la absorción del fitoplancton debido a que la clorofila tiene su segundo pico de absorción dentro de esta banda. Esta banda es sensible a las variaciones de la retro-difusión debida a la materia suspendida, teniendo el albedo del fondo una influencia baja en la reflectividad de la banda. v. La banda rojo borde (6) es prácticamente inmune al fitoplancton y a la materia disuelta, respondiendo levemente a la variación de la retro-difusión debida a la materia suspendida. La influencia del albedo del fondo es mínima. vi. Las bandas NIR (7, 8), debido a su altísima atenuación, son prácticamente inmunes a las variaciones de los IOPs y al albedo del fondo. Solamente en situaciones de turbidez muy elevadas o profundidades inferiores a uno o dos metros se pueden apreciar variaciones sensibles en estas bandas. De esta manera, y con el objetivo de resolver el sistema de ecuaciones del modelo de transferencia radiativa, se puede hacer uso de las seis primeras bandas del WV2, considerando Determinación de Parámetros de Calidad del Agua, Batimetría y Albedo del Fondo Costero 106 que cada uno de los canales es afectado de manera muy diferente por los IOPs y el albedo de fondo, siendo viable su resolución mediante optimización para las variables definidas. Por otro lado, las ecuaciones de transferencia radiativa modeladas de las bandas NIR no pueden ser utilizadas en la resolución de las RTE debido a su baja reflectividad y a la imposibilidad de eliminar el brillo solar en dichas bandas NIR. Sin embargo, pueden ser de gran ayuda en el algoritmo de eliminación del brillo solar. Considerando los datos proporcionados en la Tabla 4, según lo descrito por [26] [34] [35], la mayor debilidad de los métodos basados en la sustracción del brillo solar a partir de la banda NIR reside en que no se tiene en cuenta la reflectividad intrínseca del agua de las bandas NIR. De esta forma, en los casos en que esta reflectividad no sea despreciable se producirán errores importantes en el algoritmo de eliminación del brillo solar, al eliminarse en todas las bandas del óptico una mayor cantidad de reflectividad de la que le correspondería, influyendo negativamente en el modelado radiativo. Debido a que disponemos de imágenes de muy alta resolución espacial, la hipótesis de tratar píxeles de muy baja profundidad y/o alta turbidez en áreas muy próximas a la costa comienza a ser muy plausible, por lo que resulta necesario modelar la reflectividad intrínseca del agua para las bandas NIR. Teniendo en cuenta el análisis previo, se propone una nueva expresión, basada en el modelo físico de corrección del brillo solar (ecuación (51)), donde el valor de reflectividad inherente modelada de la banda del infrarrojo cercano  󰇛󰇜 es introducida en la ecuación, en pos de eliminar la proporción de reflectividad de la banda NIR que no es debida al brillo solar superficial. Esta expresión vendrá dada por: 󰇛󰇜󰆒󰇛󰇜, , , ,∗󰇛󰇜 󰇛󰇜 (94) Es muy importante resaltar que la introducción del valor modelado de reflectividad NIR implica que la corrección del brillo solar será realizada conjuntamente en el modelo de transferencia radiativa en la resolución iterativa del problema. Como resultado, las bandas modeladas son comparadas con los valores de reflectividad obtenidas en los canales WV2, eliminado el brillo solar, haciendo uso de la reflectividad modelada de los canales NIR. En la Figura 57 se muestra la nueva metodología de procesado propuesta para la obtención de los mapas de alta resolución de parámetros de calidad de agua, a partir de los IOPs, de batimetría y de albedo del fondo, para los ecosistemas litorales. En este nuevo procedimiento, con las tres etapas de procesado interrelacionadas, resaltar: 1. El modelo de corrección atmosférica 6S (i) proporciona información de la irradiancia directa utilizada en el algoritmo de eliminación de brillo solar. 2. El módulo de eliminación del brillo solar (ii) es introducido junto a las ecuaciones de transferencia radiativa en el algoritmo de optimización para el modelado de transferencia radiativa (iii). 3. Los resultados del modelado radiativo de las bandas NIR son utilizados en el módulo de eliminación del brillo solar (ii) para eliminar el exceso de reflectividad del NIR. Si consideramos la cadena de procesado como un sistema global e interdependiente, y teniendo en cuenta la teoría de las limitaciones, en donde la suma de los óptimos locales no tiene por qué proporcionar el óptimo global, el procesado conjunto de estos módulos ha de mejorar el resultado de todo el sistema, dado que los desajustes generados en los diferentes módulos pueden ser compensados en otros módulos en el cálculo global del problema. CAPÍTULO 5 107 Figura 57. Esquema propuesto de funcionamiento de la nueva cadena de procesado de imágenes multiespectrales WV2 de alta resolución espacial. 5.2.2. Resolución de las ecuaciones de transferencia radiativa mediante optimización numérica El modelado de las ecuaciones de transferencia radiativa de los canales del WV2 nos ha permitido obtener ochos ecuaciones, de las cuales las seis primeras proporcionan información suficiente para calcular una solución óptima de los parámetros a modelar. Gracias a la utilización de la optimización espectral (ver ecuación (79)), como función de coste en algoritmos de minimización por mínimos cuadrados, se puede tratar de calcular los valores óptimos de estos parámetros que minimizan el error espectral. Para ello existen múltiples algoritmos de optimización numérica, los cuales proceden de técnicas de análisis numéricos como el método de Newton, interpolación polinómica de Lagrange, el método de eliminación de Gauss, o el método de Euler [107]. Sin embargo, para la resolución de sistemas de ecuaciones no lineales mediante el método de mínimos cuadrados optimizado a la computación por ordenador, uno de los métodos más utilizados es el del Levenberg-Marquardt (LMA) [108]. LMA es un método iterativo que calcula el mínimo local de la suma de los cuadrados de la función de coste. Simplificadamente, es una combinación entre el método de descenso por gradiente (steepest descent), que realiza una minimización a través de la dirección del gradiente, y el método de Gauss-Newton, mediante la utilización del modelo cuadrático para acelerar el proceso iterativo de búsqueda de la solución. El algoritmo se comporta como un método de máxima pendiente, lento pero que garantiza la correcta convergencia, cuando la solución obtenida en el paso iterativo está lejos de la óptima, y se convierte en método de Gauss-Newton cuando la solución actual está cerca de la óptima, obteniendo una mayor velocidad de convergencia. Los métodos basados en gradientes necesitan calcular la matriz de los diferenciales de las variables respecto a la función, conocida como Jacobiano, por lo que se requiere obtener las derivadas parciales de primer orden de las ecuaciones de transferencia radiativa. Para ello existen dos posibilidades: (i) obtener la expresión analítica de las derivadas parciales de la función, lo cual es solo aconsejable en funciones poco complejas y con derivadas conocidas o, (ii) utilizar Determinación de Parámetros de Calidad del Agua, Batimetría y Albedo del Fondo Costero 108 diferenciación numérica, en donde se obtiene el diferencial mediante la resta del valor de la función para la variable x y el resultado de la función para la variable x+delta. En nuestro caso, aunque existe una derivada conocida de las ecuaciones de transferencia radiativa (ecuación exponencial), la complejidad de dicha ecuación, respecto a las variables implicadas en ella, hace que el resultado de la derivada parcial sea tan complejo que requeriría de centenares de líneas de código para representarlas, siendo a su vez una limitación en la velocidad de cómputo del algoritmo de optimización. Por ese motivo, y debido a que los métodos de diferenciación numérica proporcionan buenos resultados en este tipo de ecuaciones, se decidió utilizar el algoritmo LMA con diferenciación numérica en la resolución de las ecuaciones de transferencia radiativa. Como se ha mencionado, el algoritmo LMA sólo asegura la obtención de un mínimo de la función, el cual puede ser local o global. Teniendo en cuenta que estamos trabajando con un sistema de ecuaciones no lineal mal condicionado, es de esperar que existan más de un mínimo posible en el cálculo del sistema de ecuaciones. Para limitar este problema, LMA proporciona herramientas que permiten restringir y orientar el resultado mediante cuatro tipos fundamentales de parámetros: I. Valor inicial de las variables: inicializar las variables cerca del mínimo global permite al algoritmo converger al resultado correcto, siendo la profundidad el parámetro más sensible dado que funciona como parámetro regulador, entre la reflectividad generada por los IOPs y la generada por el albedo del fondo. II. Condiciones de contorno: fijan un rango en donde los valores de las variables pueden fluctuar, por ejemplo, limitar que la variable de profundidad no adquiera valores negativos pues físicamente no tiene sentido. III. Valor de escala de las variables: estos parámetros permiten normalizar los valores de gradientes del LMA, por ejemplo, la variabilidad de la profundidad con valores que se mueven típicamente entre 0 y 20 m no puede ser comparable con la variable de materia suspendida que varía entre 0.001 y 0.100. Trabajar con gradientes no normalizados hace que el optimizador tienda a minimizar más el error de las variables de mayor variabilidad, en nuestro ejemplo la profundidad. IV. Salto máximo de las variables: Teniendo en cuenta que trabajamos con una ecuación exponencial, si no se fija un valor máximo de salto de las variables, las variables que se encuentra en el exponente, por ejemplo, la profundidad tendería a obtener variaciones mucho más rápidas que las demás variables. Este hecho produce una rápida optimización del error de esta variable impidiendo a las restantes minimizar su error, alcanzando un mínimo local en donde sólo se corrige el error de profundidad. 5.3. Resultado del modelo de transferencia radiativa propuesto en entornos costeros Como se ha detallado previamente, sólo se dispone de seis canales con información suficiente para ser utilizados en el modelado de transferencia radiativa. Si bien se puede utilizar una configuración genérica, que intente calcular la mayoría de los parámetros, se puede reducir el conjunto de variables para tener un mayor número relativo de ecuaciones, introduciendo mayor robustez en el cálculo de las variables restantes. En nuestro contexto, se configurarán los parámetros prioritarios para cada estudio específico, dentro de la monitorización de áreas litorales, concretamente: análisis de calidad del agua, batimetría y albedo del fondo marino, a partir de las CAPÍTULO 5 109 imágenes multiespectrales de alta resolución del satélite WV2, previamente pre-procesadas radiométrica y atmosféricamente. A continuación, para todos los ecosistemas costeros seleccionados del Archipiélago Canario se proporcionan los resultados obtenidos, con los nuevos algoritmos implementados, en la monitorización de los parámetros vinculados con la calidad del agua, batimetría, y albedo de fondo costero, en áreas litorales con imágenes multiespectrales WV2, previamente adquiridas y procesadas. 5.3.1. Monitorización de la calidad de agua en zonas costeras de Canarias Para el cálculo de calidad de aguas vamos a seleccionar los tres parámetros asociados PGX, así como la profundidad . A su vez, se va a utilizar únicamente el albedo normalizado de la arena a los 555 nm para modelar la reflectividad del fondo marino. La elección de una sola firma espectral, al igual que en Lee et al. [67], es un intento de modelar el fondo costero con una única variable, lo que genera pequeños desajustes en el modelado, los cuales son insignificantes en el cálculo de los IOPs debido a que la aportación de reflectividad de las propiedades inherentes del agua es muy superior a la del albedo costero. Mientras que los beneficios de calcular un sistema de ecuaciones con menos incógnitas que ecuaciones proporciona una robustez adicional a los resultados. De esta forma se va a generar un sistema de seis ecuaciones y cinco incógnitas.  󰇛󰇜,,,, (95) En el desarrollo de esta investigación se ha trabajado en la estimación de los parámetros vinculados a la calidad del agua en cuatro entornos litorales canarios, con características específicas, concretamente: (i) Entorno del nuevo puerto de Granadilla, en el marco del ‘Programa Europeo de Monitorización Medioambiental’, Isla de Tenerife; (ii) Espacio natural protegido de Maspalomas, Isla de Gran Canaria; (iii) Reserva de la Biosfera de Corralejo e Isla de Lobos, Isla de Fuerteventura y, (iv) Zona de la Restinga, al sureste de la isla de El Hierro, durante la erupción del volcán submarino. A continuación, se proporcionan los resultados obtenidos para las diferentes áreas de estudio y parámetros específicos, indicadores de la calidad del agua marina, en zonas litorales: concentración de materia suspendida, disuelta y de clorofila  Costa de Granadilla El objeto principal era el estudio, análisis y cuantificación, con imágenes de satélite de alta resolución, de los parámetros oceanográficos vinculados a la determinación de la calidad del agua en la zona costera de obras del puerto de Granadilla y el cumplimiento de la evaluación del estado de conservación de las especies y los hábitats recogidos en los anexos de la Directiva Hábitat europea (directiva 2008/105/CE del Parlamento Europeo y del Consejo de 16 de diciembre de 2008 relativa a las normas de calidad ambiental en el ámbito de la política de aguas). En este contexto, comprende la evaluación de la gestión de las ZEC marinas de la red Natura 2000 en Canarias y la propuesta de criterios e indicadores de seguimiento y de medidas de protección. El interés de la investigación residía en la monitorización de la construcción de un puerto comercial, en el marco de un Programa de Vigilancia y Monitorización en esa zona industrial. Las imágenes WV2 de la zona han permitido una monitorización de la calidad del agua, en donde miles de toneladas de tierra y escombros fueron arrojadas a la costa. La imagen que mejor ilustra Determinación de Parámetros de Calidad del Agua, Batimetría y Albedo del Fondo Costero 110 los resultados de monitorización de la calidad del agua durante su construcción fue adquirida el día 1 de diciembre del 2011 (Figura 58 (a)), cuando las obras de expansión del puerto se encontraban en plena actividad. Dicha imagen, con una gran pluma de turbidez, ha sido de gran ayuda para la validación del algoritmo, gracias a que en ese mismo día se tomaron muestras insitu de los niveles de turbidez del agua (ver Figura 58 (b)). En la Tabla 8 se muestran los valores de turbidez óptica (Nephelometric Turbidity Unit, NTU) y su equivalente en concentración de materia suspendida [109]. Tabla 8. Datos in-situ de turbidez obtenidos para la imagen de granadilla del 1/12/2011 Turbidez (NTU) TSM (g/m3) P1 0.4 1.36 P2 5.6 19.16 P3 2.7 9.23 P4 3.5 11.97 En la Figura 58 se muestran los resultados obtenidos en la estimación, mediante el modelo de transferencia radiativa propuesto, de los parámetros de calidad del agua para la imagen de la costa de Granadilla del 1 de diciembre de 2011. Concretamente, en la Figura 58 (a) se puede observar las obras del puerto y la gran pluma de turbidez; en la Figura 58 (b) se muestra el mapa obtenido de concentración de materia suspendida y la ubicación exacta de las muestras in-situ contemporáneas con el paso del satélite; en la Figura 58 (c) se muestran las concentraciones estimadas de clorofila y, finalmente, en la Figura 58 (d) se muestran las concentraciones de materia disuelta presentes en el área de estudio. (a) (b) CAPÍTULO 5 111 (c) (d) Figura 58. Resultados parámetros de calidad del agua. (a) Área de estudio (Granadilla 1/12/2011), (b) mapa de concentración materia suspendida, (c) concentración de clorofila, (d) concentración de materia disuelta. Se puede observar una importante correlación entre los datos de concentración de materia suspendida in-situ y los puntos de la imagen: P1 tiene un valor de turbidez normal de aguas abiertas, P2 es tomado en el punto exacto de vertido de materiales, con lo que el valor de turbidez es muy elevado (color rojo cercano a 20 g/m3) y P3-P4 se han tomado dentro del remolino generado en la pluma, con valores cercanos a 10 g/m3. De la misma manera se puede observar que la concentración obtenida de clorofila es nula, lo cual era de esperar en estas condiciones, mejorando los resultados obtenidos mediante algoritmos empíricos como el OC3, los cuales confunden el valor verdoso del agua con presencia de clorofila. Resaltar que, como se observa en la Figura 58 (d), una parte de la arena vertida al mar es lo suficientemente fina para disolverse en el agua, produciendo cierta concentración de materia disuelta asociada a la pluma de turbidez. Seguidamente, en la Figura 59 se muestran los resultados del algoritmo de calidad de agua desarrollado cuando se aplica a la imagen completa, donde se aprecia, claramente, como el vertido de materiales genera valores de turbidez que se extienden por toda la costa de Granadilla. Se aprecia una pequeña interferencia entre el valor del albedo de alta reflectividad del fondo calcáreo situado a unos 25 m de profundidad, en la parte central de la imagen, con los resultados de materia suspendida, que adquiere un incremento de entrono a 0.5 g/m3 respecto a los valores cercanos. Los valores de clorofila obtenidos son de carácter residual, inferiores a 0.1 mg/m3, consistentes con los obtenidos a baja resolución espacial mediante los satélites MERIS-MODIS para dicho día (pueden ser consultados en la Web de Ocean Color [110]). Determinación de Parámetros de Calidad del Agua, Batimetría y Albedo del Fondo Costero 118 El ajuste mediante regresión (función exponencial) del parámetro kd(490) de la imagen WV2 se realizó con cuatro bandas (R2, R3, R4, y R5) [112]. Se introdujo una función de peso () en la ecuación, encargada de balancear el uso de las bandas de mayor longitud de onda, que permite un mejor ajuste cuando existen altos contenidos de turbidez, pero en cambio son más ruidosas en ausencia de reflectividad. La función de ajuste implementada, para la obtención del parámetro kd(490), viene dada por: 󰇛490󰇜10󰇛󰇜󰇧.󰇡      󰇢 . .󰇨󰇧.󰇡      󰇢 . .󰇨 2∗ 0.4,  0.4 0.9 (96) En la Figura 65, se muestra el resultado del algoritmo operacional para el cálculo de kd(490), para la imagen WV2 adquirida el 27 de octubre de 2011 (Figura 65 (a)), durante el proceso eruptivo submarino. En la Figura 65 (b) podemos observar que las concentraciones de kd(490) alcanzan valores de 1.5 m-1 en el centro (color marrón) de la imagen. En la Figura 66 se muestra el resultado del algoritmo de calidad del agua que permite obtener el valor de kd(490) (Figura 66 (d)). Se puede observar como se alcanzan valores cercanos a 1.5 m-1, con una elevada correlación respecto al algoritmo de ajuste. A su vez, gracias al modelado radiativo podemos distinguir entre los diferentes parámetros de calidad del agua: en la Figura 66 (a) se puede observar como la atenuación de la clorofila, en la mancha eruptiva, es despreciable, como se había constatado en la campaña, pero que producían grandes errores en los algoritmos de clorofila para aguas abiertas; un elevado contenido de sustancias disueltas en el agua es detectado en la Figura 66 (b), siendo el parámetro dominante con valores de atenuación cercanos a 1 m-1 y, finalmente, en la Figura 66 (c) se pueden apreciar los valores de atenuación difusa debida a la materia suspendida, la cual alcanza valores elevados en el centro eruptivo, (~ 0.3 m-1). (a) (b) Figura 65. Composición RGB de la imagen WV2 del 27 de octubre de 2011 (a), concentración de k d (490) obtenida mediante el algoritmo operativo ecuación (96). CAPÍTULO 5 119 (a) (b) (c) (d) Figura 66. Resultados del algoritmo de calidad de aguas para la imagen WV2 del 27 de octubre de 2011. (a) Atenuación debida a la concentración de clorofila, (b) atenuación debida a la materia disuelta, (c) difusión debida a la materia suspendida, y (d) atenuación difusa k d (490). 5.3.2. Batimetría en zonas litorales de Canarias En este caso, con el objeto de determinar la batimetría vamos a seleccionar la profundidad z, parámetro que va a ser inicializado con el resultado obtenido mediante el ratio algorithm, así como los tres parámetros asociados a la calidad del agua PGX, aunque el parámetro de absorción de la clorofila P podría ser obviado debido al bajo contenido de clorofila en las aguas costera en estudio. Determinación de Parámetros de Calidad del Agua, Batimetría y Albedo del Fondo Costero 120 A su vez, se va a utilizar únicamente el albedo normalizado de la arena a los 555 nm para modelar la reflectividad del fondo marino. De esta forma se va a generar un sistema de ecuaciones de seis ecuaciones y cinco incógnitas.  󰇛í󰇜,,,, (97) En el desarrollo de esta investigación se ha trabajado en el cálculo de batimetría en tres entornos litorales canarios, con características específicas: (i) costa de Granadilla, (ii) la costa de Maspalomas, y (iii) costa de Corralejo. A continuación, se procede a mostrar los resultados del modelado de transferencia radiativa para el cálculo de batimetría en los diferentes escenarios.  Implementación del ratio algorithm para el cálculo de batimetría en imágenes WV2 En el marco de la investigación realizada para la monitorización del puerto de Granadilla en el 2011, y debido a la necesidad urgente de monitorización de la profundidad en el área de construcción del puerto, fue necesario la implementación de un algoritmo de cálculo de batimetría basado en el método ratio algorithm. Para implementar el ratio algorithm se realizó un estudio detallado de la variación de la reflectividad de las bandas del WV2 respecto a la profundidad del píxel, mediante curvas de regresión con los datos de batimetría de la costa de Granadilla para una adquisición de la imagen con un ángulo de visión e iluminación medio. La zona de estudio tiene un albedo de fondo que es representativo de las costas canarias con arenales de reflectividad media, áreas de sebadales y alta reflectividad en áreas submarinas más profundas. Para la implementación del algoritmo de dos bandas se llegó a la conclusión que las mejores bandas son el verde y el azul. Los principales motivos fueron:  Las bandas del verde y el azul tienen la mayor penetración en aguas costeras, cercano a los 25 m. La banda azul costa y el amarillo, en este tipo de aguas de mayor turbidez, tienen una menor penetración, en especial el amarillo que no supera los 15 metros.  Las bandas del verde y el azul son obtenidas por el mismo sub-sensor, por lo que no se produce ruido ligado a la no correlación temporal de las bandas obtenidas en los dos dispositivos sensores del WV2, como se analizó en el Capítulo 3. En la Figura 67 (a) se muestra el transecto costero utilizado en el ajuste por regresión, y en la Figura 67 (b) se muestra el ajuste de la relación numérica 󰇡󰇛∗󰇜 󰇛∗󰇜󰇢 representado en el eje X y los valores de batimetría en el eje Y. La línea azul representa el ajuste lineal, mostrando tal comportamiento entre 1 y 20 metros, mientras que se puede observar cierta saturación sobre todo en aguas más profundas, a partir de los 20 metros, proporcionando un ajuste con un R2 = 0.917. En la Figura 55 se observaba como los albedos de alta reflectividad, a partir de los 18 metros, generaba ese tipo de comportamiento no lineal a esas profundidades. Este fenómeno es muy común en las costas canarias, donde se pueden encontrar albedos de fondo costero muy elevados (angileras), a partir de los 20 metros de profundidad. Para lograr una mayor adaptación en estas condiciones, proponemos un ajuste cuadrático (Figura 67 (b) en rojo) que se ajusta mejor a estos tipos de fondos, proporcionando un ajuste con un R2 = 0.979. CAPÍTULO 5 121 (a) (b) Figura 67. (a) Transecto costero utilizado en el ajuste. (b) Ajuste del ratio algorithm mediante regresión con la batimetría in-situ: Ajuste lineal (línea azul) y ajuste cuadrático (curva roja). A continuación, en las ecuaciones (98) y (99), se muestra el ajuste lineal y cuadrático del ratio algorithm realizados para las bandas del WV2. 103.2∗ln󰇛1000∗󰇜 ln󰇛1000∗󰇜106.8 (98) 507.5∗ln󰇛1000∗󰇜 ln󰇛1000∗󰇜1017∗ln󰇛1000∗󰇜 ln󰇛1000∗󰇜510.4 (99) Finalmente, en la Figura 68 (a) y (b) se muestran los resultados obtenidos mediante el ratio algorithm tanto para el ajuste lineal como para el cuadrático, en el área bajo estudio, en comparación con la batimetría obtenida mediante sonar, que se muestra en la Figura 68 (c). Se puede observar como el resultado del ajuste cuadrático se adapta mejor en el rango de profundidades superiores a los 20 metros. (a) (b) (c) Figura 68. Resultado del ratio algorithm: (a) ajuste lineal, (b) ajuste cuadrático, (c) batimetría sonar. Determinación de Parámetros de Calidad del Agua, Batimetría y Albedo del Fondo Costero 122  Costa de Granadilla La costa de Granadilla situada en la isla de Tenerife, con una plataforma insular poco desarrollada, adquiere rápidamente profundidades elevadas a pocas decenas de metros del litoral. Sin embargo, existen zonas de profundidades más bajas, como por ejemplo las playas. Esto nos ha permitido comparar los resultados del algoritmo de batimetría en amplias zonas de la imagen. La imagen de estudio tiene una calidad media, debido a unas condiciones de oleaje moderado, lo que ha producido una imagen con cierto contenido de brillo solar especular. La presencia moderada de brillo especular hace que su corrección genere un cierto grado de ruido debido al cálculo entre bandas del algoritmo de corrección. En la Figura 69 se muestran los resultados obtenidos en el modelado de transferencia radiativa para el cálculo de batimetría: La Figura 69 (a) muestra una composición RGB de una zona de baja profundidad al sur de la costa de Granadilla; en la Figura 69 (b) se muestra la comparativa entre los valores de batimetría satélite y la batimetría sonar, para más de 300 puntos de test, mientras que en la Figura 69 (c) se muestra la imagen de batimetría sonar de alta resolución de la zona de estudio. Finalmente, en la Figura 69 (d) se muestra el resultado de la batimetría del algoritmo RTM obtenida para la zona de estudio. Se puede observar como la gráfica de dispersión genera una función lineal, con una pendiente cercana a la unidad, con una constante (sesgo) de 1.94 metros asociado con el nivel de marea. El valor de ajuste R2 = 0.94 es muy alto, así como el error cuadrático medio es de 1.94 metros, siendo un error bastante reducido. El parecido entre dichos resultados se pueden observar en la similitud de los dos mapas de batimetría de las imágenes Figura 69 (c) y (d). A su vez, en la Figura 69 (d) se puede apreciar un cierto nivel de ruido asociado con el generado en la eliminación del brillo solar debido al oleaje y a la baja reflectividad del fondo costero. Se puede observar que la obtención de la información del fondo marino es más dependiente del ruido de la imagen, lo que es lógico teniendo en cuenta la baja reflectividad del píxel de agua costera, y a que la aportación del albedo marino es mínimo en comparación a la reflectividad del agua marina generada por el back-scattering. (a) (b) CAPÍTULO 5 123 (c) (d) Figura 69. Resultados de batimetría para la costa de Granadilla. (a) Imagen RGB de la zona de estudio, (b) gráfica de dispersión de datos satélite-sonar, (c) mapa de batimetría sonar, y (d) mapa de batimetría satélite.  Costa de Maspalomas La costa de Maspalomas y Playa del Inglés tiene gran interés para el estudio de batimetría debido a que, de la misma forma que existe una acumulación y pérdida de arena en el campo de dunas, dicha variabilidad se extiende en el fondo marino, generando un entorno dinámico en donde la arena es desplazada por las corrientes costeras y el arrastre del oleaje. Este hecho es notable dado que el área de estudio se sitúa en el extremo sur de la isla de Gran Canaria, en donde se genera un vórtice debido a la corriente superficial con dirección norte-sur producido por los vientos Alisios. El hecho de disponer de un fondo dinámico nos va a permitir detectar posibles cambios en el fondo arenoso, una vez que hemos probado que el algoritmo produce resultados robustos en fondos más estáticos. En la Figura 70 se muestran los resultados obtenidos del algoritmo de batimetría para dos fechas diferentes, 20 de noviembre de 2011 y 17 de enero de 2013. Dichos resultados son comparados con los resultados de batimetría de alta resolución obtenida en el año 2003 (Figura 70 (a)). Se puede observar un cierto parecido entre los resultados de la batimetría sonar con los resultados obtenidos en el 2011 y el 2013. Sin embargo, existen variaciones significativas en las batimetrías satelitales, pudiéndose apreciar una mayor acumulación de arena en las aguas someras situadas por debajo de la costa para la imagen del 20 de noviembre de 2011. Como consecuencia de este dinamismo, la línea de costa en ambas imágenes tiene ciertas diferencias. Dichas variaciones en la línea de costa en la esquina sureste de Playa del Inglés es algo recurrente en todas las imágenes satelitales adquiridas pues, como se ha indicado, es una zona conocida por su acumulación y pérdida estacional de arena. Determinación de Parámetros de Calidad del Agua, Batimetría y Albedo del Fondo Costero 124 Figura 70. Resultado del algoritmo RTM de batimetría para la costa de Maspalomas: (a) Batimetría sonar de la zona de estudio del año 2003, (b) composición RGB del 20 de noviembre de 2011, (c) resultado del algoritmo de batimetría para el 20 de noviembre de 2011, (d) composición RGB del área de estudio para el 17 de enero de 2013, (e) resultado del algoritmo de batimetría para el 17 de enero de 2013. CAPÍTULO 5 125  Costa de Corralejo-Isla de Lobos La costa de Corralejo tiene gran interés para el estudio de batimetría debido a que se encuentra situado al lado del canal interinsular que separa Fuerteventura y Lanzarote, cerca de la isla de Lobos, con amplias áreas de aguas someras con profundidades inferiores a los 25 metros. Esto nos ha permitido comparar los resultados del algoritmo de batimetría en amplias franjas de la imagen. La imagen de estudio tiene una gran calidad, debido a las condiciones de oleaje casi inexistente, lo que ha proporcionado una imagen libre de todo tipo de brillo solar especular. La no existencia de brillo solar hace que la corrección sea prácticamente superflua, lo que se traduce en un mínimo grado de ruido debido al cálculo entre bandas. Otro elemento que contribuye a la calidad de la imagen son los fondos arenosos de alta reflectividad, los cuales pueden ser apreciables en la imagen a simple vista y junto al bajo nivel de ruido permite el correcto modelado de la profundidad del fondo costero. En la Figura 71 se muestran los resultados obtenidos en el modelado de transferencia radiativa para el cálculo de la batimetría: En la Figura 71 (a) muestra una composición RGB de una zona de baja profundidad entre la costa de Corralejo y la isla de Lobos; la Figura 71 (b) muestra la comparativa entre los valores de batimetría satélite y la batimetría sonar para más de 300 puntos de test mientras que en la Figura 71 (c) se muestra la imagen de batimetría sonar de alta resolución de la zona de estudio. Finalmente, en la Figura 71 (d) muestra el resultado de la batimetría del algoritmo RTM obtenida para la zona de estudio. Se puede observar como la gráfica de dispersión genera una función lineal muy cercana a y=x, con un valor de pendiente muy cercano a 1. El valor del parámetro de ajuste R2 = 0.93 es muy alto, de la misma manera el error cuadrático medio es de 1.2 metros, lo que supone un error bastante pequeño, siendo inferior al tamaño del píxel de la imagen. La correlación entre dichos resultados se pueden observar en el parecido entre los dos mapas de batimetría de las imágenes Figura 71 (c) y (d). (a) (b) Determinación de Parámetros de Calidad del Agua, Batimetría y Albedo del Fondo Costero 126 (c) (d) Figura 71 Resultados de batimetría para la costa de Corralejo: (a) Imagen RGB de la zona de estudio, (b) gráfica de dispersión de los datos satélite-sonar, (c) mapa de batimetría sonar, y (d) mapa de batimetría satélite. 5.3.3. Estimación del albedo del fondo costero en zonas litorales de Canarias En el cálculo del albedo marino, usando el modelo de transferencia radiativa, se han utilizado siete parámetros. Como en las configuraciones anteriores, inicialmente, se ha hecho uso de los parámetros de calidad de agua PGX y el parámetro de profundidad z. Sin embargo, para el modelado del albedo marino, adicionalmente, se hace uso del desmezclado lineal de las tres clases puras más comunes en fondos costeros. Dichas clases son la arena de alta reflectividad, la reflectividad de las algas, utilizándose un valor promedio de las clases más comunes y, finalmente, la clase sedimento con una baja reflectividad en todos los canales (ver Figura 72). Figura 72. Reflectividad normalizada de las clases puras más usuales en los fondos costeros (arena, algas, y sedimentos) [103]. CAPÍTULO 5 127 Por lo tanto, la reflectividad final será el resultado de la mezcla lineal de estas tres clases puras mayoritarias en los fondos costeros:  󰇛󰇜,,,,,, (100) A continuación, y específicamente para los ecosistemas costeros de Granadilla y Corralejo-Lobos, dada su singularidad para este estudio, se presentan los resultados del modelado de transferencia radiativa para el cálculo del albedo de fondo para estos escenarios.  Costa de Granadilla En la costa de Granadilla se puede encontrar zonas con fondos de baja profundidad, permitiéndonos obtener el albedo del fondo costero, como se muestra en la Figura 73. En la Figura 73 (a) se muestra una composición RGB del albedo obtenido mediante la mezcla lineal de las tres clases puras y sus abundancias. Se puede observar como la reflectividad del fondo costero tiene una tonalidad más apagada, correspondiente a la presencia de un fondo sedimentario y de arena volcánico de menor reflectividad. A su vez, se puede apreciar la presencia de un tono verdoso asociado a la presencia de sebadales. En la Figura 73 (b) se puede observar la abundancia de la arena de alta reflectividad, siendo esta clase minoritaria en este fondo costero. En los afloramientos rocosos la abundancia de la arena desaparece por completo y, al igual que la zona donde se encuentran los sebadales, existe una mezcla de reflectividad de la arena que llega casi al 40 %, lo cual junto a que la profundidad media de estas zonas nos hace pensar que es congruente con la existencia de sebadales. Se puede observar como en las zonas más profundas es donde la abundancia aumenta por encima del 75 %, esto es debido a que se tratan de fondos calcáreos de muy alta reflectividad. En la Figura 73 (c) se muestra la abundancia de las algas, en donde se observa unas abundancias elevadas superiores al 40 % en profundidades medias, en torno a 10-20 metros de profundidad. Esto junto a la mezcla de fondo arenoso es congruente con la existencia de sebadales de media y alta densidad, lo cual puede ser refrendado en los mapas bentónicos de la zona de estudio [113]. Finalmente, en la Figura 73 (d) se muestra la abundancia de los sedimentos, en este caso modela los afloramientos de rocas y parte de la reflectividad de la arena que en esta costa es de menor reflectividad. Conclusiones 134 Estos productos satelitales pasan desde la obtención de la calidad de las aguas en entornos naturales y playas, a la monitorización de la biodiversidad bentónica del lecho costero y al cálculo batimétrico de sus costas. Asimismo, debido a la naturaleza volcánica de las islas, la utilización y procesado de imágenes de alta resolución, en el contexto de este trabajo, ha permitido la monitorización de eventos extraordinarios, tales como la erupción submarina acaecida en la Isla de El Hierro en octubre de 2011. Especial énfasis se ha realizado, dentro de los diferentes niveles de procesado, en los procedimientos relacionados con la corrección atmosférica, la corrección del brillo solar y el modelado de la ecuación de transferencia radiativa del medio marino en las zonas litorales. En el contexto de la corrección atmosférica y del reflejo del brillo solar, la determinación precisa de la reflectividad superficial es crítica para estudios cuantitativos de los ecosistemas y monitorización del entorno. Con este objetivo, se ha modificado y adaptado el modelo de transferencia radiativa 6S y se ha implementado un nuevo modelo de deglinting. Para ambos procesos se han analizado, tanto los aspectos relacionados con el proceso sistemático de compilación del conjunto de datos de comparación entre medidas in-situ y observaciones coincidentes del satélite, como aquellos otros vinculados con los términos de corrección y optimización del algoritmo 6S y del algoritmo de corrección del reflejo solar automático propuesto. En relación con la implementación de un nuevo modelo de transferencia radiativa del agua, éste nos ha permitido modelar los tres parámetros principales de calidad del agua (clorofila, materia disuelta y materia suspendida), así como el modelado de la batimetría y el albedo del fondo costero. El trabajo realizado en esta Tesis, convenientemente modularizado y documentado sobre interfaces gráficas, podría permitir disponer de herramientas de procesado de imágenes del satélite WorldView-2 que mejoran las precisiones del análisis y de la interpretación de datos y facilitarían su adaptación a futuros sensores de observación de la Tierra u otras aplicaciones de la tecnología de la teledetección, como pueden ser la clasificación multitemporal, la fusión de datos multisensoriales, la detección de movimiento y el reconocimiento de estructuras. A continuación, se analizarán con más detalle tanto las principales contribuciones de esta Tesis Doctoral, como las líneas de investigación que han quedado abiertas. Finalmente, se detallarán las publicaciones indexadas, y las contribuciones a Congresos Nacionales e Internacionales vinculadas directamente a esta Tesis. 6.1. Principales contribuciones El trabajo realizado en la Tesis puede ser dividido en tres bloques principales: corrección atmosférica, corrección del brillo solar, y modelado de transferencia radiativa del medio acuático, que se corresponden con los capítulos 3, 4, y 5. A continuación se van a presentar las principales contribuciones científicas realizadas en cada uno de estos tres bloques temáticos.  Corrección atmosférica En el contexto del modelado atmosférico de las imágenes del satélite WV2, la principal aportación realizada ha sido tanto la adaptación y configuración del modelo de corrección atmosférica 6S a este nuevo sensor de alta resolución, como la validación de los resultados obtenidos mediante datos in-situ. Si bien existen múltiples métodos de corrección atmosférica, debido a la baja reflectividad del agua se requería la utilización de un modelo atmosférico avanzado que permitiese una adecuada corrección. Aunque existen varios modelos, tras una revisión del estado del arte se procedió a la elección del modelo 6S, tanto por su naturaleza de código abierto como por ser un modelo fiable y CAPÍTULO 6 135 testado en la corrección de imágenes MODIS en aplicaciones marinas (Código proporcionado por la NASA). El modelo 6S, si bien permite una fácil configuración para sensores de baja resolución, como el MODIS, no disponía de una configuración predeterminada para satélites de muy alta resolución como el WoldView-2. Por este motivo, la correcta configuración de la geometría de adquisición del satélite y la definición de las 8 bandas de paso del satélite tuvieron cierta complejidad en la configuración del modelo. Es de resaltar, en este contexto, que en el marco de un contrato de investigación con la Fundación Estatal OAG, se comenzó el procesado de imágenes WV2 en el año 2011, apenas un año después de su lanzamiento, siendo pioneros en la corrección atmosférica de imágenes WV2 con el modelo 6S. Otra de las principales aportaciones, en el marco del modelado atmosférico (detallada en el capítulo 3), ha sido la validación in-situ de la reflectividad superficial. Para ello se llevó a cabo una campaña de medidas reales de reflectividad superficial, mediante la utilización de un radiómetro de campo, en donde se tomaron muestras de zonas terrestres y en aguas costeras. Los resultados obtenidos de dicha validación nos mostraron errores reducidos, inferiores al 8 % de la reflectividad in-situ para las muestras terrestres, y errores ligeramente superiores en las aguas costeras, debido a que la adquisición de la imagen y la toma in-situ no fueron llevadas a cabo exactamente en el mismo día ni con las mismas condiciones marinas ni de oleaje. Por este motivo se considera que los resultados obtenidos en la corrección atmosférica han sido satisfactorios.  Corrección del brillo solar especular Como se ha podido demostrar en esta Tesis, el brillo solar especular, sobre la superficie marina, es un problema muy complejo de corregir y en permanente investigación. De esta manera, tras una revisión exhaustiva del estado del arte, la primera aportación realizada, en el marco de este trabajo, fue la adaptación del algoritmo de Hedley a la corrección de imágenes de alta resolución espacial. Este tipo de algoritmo se basa en la hipótesis de que la reflectividad en la banda NIR sólo es debida al brillo solar, dada la gran absorción del agua en esas longitudes de onda. Haciendo uso de una relación lineal es posible sustraer el brillo solar de las bandas del rango visible mediante el valor de reflectividad del canal NIR. Para lograr esta adaptación, se tuvo en cuenta la falta de sincronización temporal entre los dos grupos de bandas del sensor (MS1: Azul, verde, rojo y NIR1 y MS2: Azul costa, amarillo, rojo borde, NIR2). De esta forma, la corrección de cada banda del grupo multiespectral se realizó con su correspondiente banda NIR. Finalmente, en este contexto, se aplicaron diferentes técnicas de procesado de imágenes para la eliminación de píxeles ruidosos o atípicos y para el mitigado del ruido producido por las operaciones entre bandas debido, fundamentalmente, a pequeños desajustes en el corregistro de las bandas asociados a errores en el remuestreo de los datos en el pre-procesado de las imágenes. La aportación principal en la corrección del brillo solar especular fue la implementación de un nuevo método físico, basado en el método empírico de Hedley, que permite el cálculo del parámetro de la pendiente de la relación lineal entre las bandas ópticas y el NIR de una manera analítica, mediante la utilización de la información procedente del modelo 6S a partir de la irradiancia directa normalizada en la superficie. Este algoritmo analítico nos ha permitido dejar de depender del método de regresión que requiere de la selección manual de zonas costeras de oleaje uniforme de alta profundidad y sin variaciones en la turbidez. Así, dicho método empírico no puede ser utilizado en imágenes de aguas interiores, introduciendo a su vez cierta incertidumbre en la regresión debido a la existencia de valores atípicos (outliers), como la espuma del mar o desajustes espacio-temporales entre bandas. En consecuencia, la eliminación del paso de regresión lineal sobre una ventana de la imagen permite el cómputo del algoritmo de una forma Conclusiones 136 completamente automática, lo que es de gran interés en sistemas operacionales de procesado de imágenes de satélite.  Modelado de transferencia radiativa del agua para imágenes WV2 Como se refleja en el análisis del estado del arte, el modelado de transferencia radiativa en aguas costeras fue concebido para imágenes hiperespectrales obtenidas, típicamente, mediante sensores aerotransportados. El hecho de que las imágenes hiperespectrales de alta resolución solo estén disponibles mediante sistemas aerotransportados hace que dicha tecnología implique, lógicamente, unos costes muy elevados para las aplicaciones consideradas. Por otra parte, el lanzamiento del satélite WorldView-2, con sus ocho canales multiespectrales de alta resolución, abrió la puerta a la utilización de algoritmos pensados para datos hiperespectrales a precios muy inferiores gracias al uso de dicha plataforma espacial. Por ese motivo, la principal y novedosa aportación científica de esta Tesis Doctoral ha sido la adaptación del modelado radiativo costero a las necesidades y limitaciones impuestas por los canales multiespectrales de alta resolución del satélite WV2. Como se ha mostrado en el estado del arte del capítulo 5, existe una gran cantidad de literatura sobre el modelado radiativo, sin embargo, las restricciones en el número de bandas y en el ancho de banda multiespectral hizo necesario seleccionar los parámetros más apropiados para el modelado y los métodos aproximativos más adecuados a la escasa información multiespectral disponible. Gracias a la correcta adaptación del modelo radiativo ha sido posible estimar los tres parámetros principales de calidad de agua (clorofila, materia disuelta, y materia suspendida), así como modelar la batimetría y el albedo del fondo costero. Otra importante aportación ha sido la integración del método de eliminación del brillo solar dentro del modelo de transferencia radiativa. Como se ha descrito en el capítulo 4 y 5, una de las mayores debilidades del método de eliminación del brillo solar basado en la banda NIR, es suponer que la reflectividad del agua costera en la banda del infrarrojo es despreciable. Debido a la alta resolución del satélite se tiene acceso a áreas costeras de profundidades muy reducidas y de altos valores de turbidez. Este hecho provoca que la reflectividad del agua costera en la banda NIR en estas condiciones no sea despreciable, lo que produce importantes desajustes en el algoritmo de eliminación del brillo solar y, consecuentemente, en el posterior modelado radiativo. La integración del método de eliminación del brillo solar dentro del modelo ha permitido calcular la reflectividad intrínseca del agua en los canales NIR, y de esta forma eliminar esta aportación extra de reflectividad para obtener una mejor estimación de brillo solar. Esta integración permite una interconexión entre los tres módulos clásicos: corrección atmosférica, deglinting y modelado de transferencia radiativa, que originalmente siempre han estado inconexos. Este hecho permite compensar, en parte, los errores generados en los diferentes módulos al estar interconectados. Otra importante aportación, en este contexto, ha sido el cálculo de los parámetros de calidad de agua, debido a que los algoritmos empíricos existentes desarrollados para aguas abiertas han demostrado que generan resultados completamente erróneos en entornos costeros, con niveles altos de turbidez y con profundidades reducidas, donde el albedo del fondo influye en la reflectividad obtenida. Gracias al modelado radiativo se ha logrado obtener mapas de concentración de clorofila, materia disuelta y suspendida para diferentes zonas costeras Canarias. La aportación del modelado radiativo de alta resolución tiene también una gran importancia en el cálculo de mapas de batimetría de alta resolución. Gracias a los buenos resultados obtenidos en los mapas batimétricos, esta nueva metodología puede ser tenida en cuenta como una alternativa real y de calidad al cálculo de batimetría mediante sónares de alta resolución, a bordo de buques de investigación oceanográfica, disminuyendo su coste de forma importante. CAPÍTULO 6 137 La integración del algoritmo de desmezclado lineal de clases bentónicas puras en el modelo de transferencia radiativa ha permitido modelar diferentes tipos de fondos costeros, siendo una importante aportación a la clasificación bentónica. Debido a la poca información multiespectral solo ha sido posible modelar clases bentónicas primarias, presentes en el fondo marino, como son los arenales, algas y sedimentos o rocas. Sin embargo, los mapas de abundancia de clases y los mapas de albedo costero generados tienen un gran interés gracias a su alta resolución, pudiendo ser utilizados en la generación de mapas de especies bentónicas mediante algoritmos de clasificación supervisados correlado con la información experta de biólogos marinos. Finalmente, es de especial relevancia, en el conjunto de aportaciones científicas de esta Tesis Doctoral, el correcto modelado de las fuentes de ruido de la imagen, como valores excesivos en el espesor óptico del aire en el modelo atmosférico y el brillo solar superficial elevado en condiciones de oleaje prominente. Este correcto modelado nos ha permitido hacer uso de un mayor rango de imágenes, dado que originalmente las imágenes debían de tener unas condiciones atmosféricas y de oleaje casi perfectas para su correcta utilización. Este hecho tiene gran importancia dado que originalmente se descartaban más del cincuenta por ciento de las imágenes por sus condiciones atmosféricas o de brillo solar. En el desarrollo de esta investigación se ha trabajado en la estimación de los parámetros vinculados a la calidad del agua, batimetría y albedo del fondo, en cuatro ecosistemas costeros canarios, con características específicas, concretamente: entorno del nuevo puerto de Granadilla; espacio natural protegido de Maspalomas; la reserva de la Biosfera de Corralejo e Isla de Lobos y la zona de la Restinga, al sureste de la isla de El Hierro. Los resultados obtenidos para las diferentes áreas de estudio y parámetros específicos, comparativamente con medidas in-situ, mapas de batimetría obtenidos con sonar y mapas bionómicos, proporcionados por el Gobierno de Canarias, muestran una adecuada correlación con todos los indicadores de la calidad del agua marina y mapas satelitales de batimetría y albedo del fondo marino. 6.2. Líneas futuras de investigación Las aplicaciones asociadas a la monitorización costera mediante el modelado de transferencia radiativa de alta resolución espacial de esta Tesis Doctoral son múltiples, algunas de las cuáles constituyen, sin lugar a dudas, líneas abiertas de investigación, como se describe seguidamente.  Cartografía bentónica del fondo en zonas litorales con imágenes de satélite de alta resolución La clasificación bentónica de alta resolución es, científicamente, de gran interés, permitiendo obtener mapas del lecho costero de alta resolución espacial de las principales especies presentes en las áreas costeras. Para ello sería necesario el conocimiento previo de las diferentes clases existentes en la costa y su disposición según la profundidad, pudiendo obtenerse mapas de clases bentónicas mediante el cálculo del desmezclado lineal de las especies y la profundidad del fondo, haciendo uso de métodos de clasificación supervisado como, por ejemplo, el Support Vector Machine (SVM), o de clasificación orientada a objeto. En la Figura 75 se muestra un ejemplo de clasificación supervisada de comunidades bénticas mediante el uso del albedo costero de las imágenes WV2 para la zona de Granadilla. Para este estudio supervisado ad-hoc para el área de Granadilla se hizo uso de datos in-situ mediante inmersiones y transectos, los cuales fueron utilizados como entrada del clasificador. Conclusiones 138 (a) (b) Figura 75. (a) Clasificación supervisada de comunidades bénticas de alta resolución mediante datos de albedo de fondo de imágenes WV2. (b) Clasificación de comunidades bénticas CIMA 2008. Finalmente, los resultados de las clases bentónicas obtenidas fueron comparados con el estudio CIMA 2008 [114] realizado mediante muestreo de datos in-situ, que genera como se puede observar mapas de baja resolución espacial. La comparativa visual entre ambos mapas de comunidades bentónicas nos muestra una alta similitud, sin embargo, se puede apreciar un mayor nivel de detalle en la clasificación generada mediante imágenes satelitales de alta resolución.  Monitorización de calidad de aguas interiores Otra línea de investigación de interés sería el estudio de calidad de aguas interiores, como presas, pantanos y humedales protegidos. Debido al área reducida de las imágenes de muy alta resolución respecto a las imágenes de media y baja resolución, las aplicaciones orientadas a espacios específicos de alto interés como las costas, puertos, lagos, presas y humedales son las que más se adecúan a las necesidades de esta plataforma. A su vez, el cálculo de calidad de aguas en presas y pantanos, los cuales son los responsables de suministrar agua de consumo humano, puede adquirir gran interés. De la misma forma podría ser posible utilizar el cálculo de la batimetría para calcular el volumen de agua almacenada, en escenarios de baja turbidez de las aguas.  Integración en un modelo único de transferencia radiativa atmosférico-marino Si bien en la presente Tesis Doctoral se ha avanzado hacia la integración de los tres módulos necesarios para el modelado de transferencia radiativa costera, aún existe una desconexión entre los dos grandes modelos radiativos, atmosférico y marino. La integración de ambos sistemas en un modelo único atmosférico-marino que integre el nuevo método físico de eliminación del brillo solar nos permitiría estimar parámetros atmosféricos fundamentales como el espesor óptico de la atmósfera o la absorción de los gases atmosféricos, los cuales son obtenidos de forma aproximada mediante datos satelitales de muy baja resolución o mediante perfiles atmosféricos típicos. Dicho modelado conjunto permitiría reducir los errores cometidos en el modelo atmosférico, siendo posible lograr mejores adaptaciones en la parametrización de la retro difusión y absorción de la atmósfera. Punto de gran importancia en el área de Canarias debido a que los episodios de calima no están correctamente caracterizados en los modelos atmosféricos. CAPÍTULO 6 139 6.3. Producción científica A continuación se presenta la producción científica derivada de las investigaciones realizadas durante la Tesis Doctoral. Revistas de impacto 1. Martin, J.; Eugenio, F.; Marcello, J.; Medina, A. Automatic Sun Glint Removal of Multispectral High-Resolution Worldview-2 Imagery for Retrieving Coastal Shallow Water Parameters. Remote Sensing, 2016, vol. 8, no 1, p. 37. [Índice de impacto: 3.180] 2. Eugenio, F.; Marcello, J.; Martin, J. High-Resolution Maps of Bathymetry and Benthic Habitats in Shallow-Water Environments Using Multispectral Remote Sensing Imagery. Geoscience and Remote Sensing, IEEE Transactions on, 2015, vol. 53, no 7, p. 3539-3549. [Índice de impacto: 3.514] 3. Eugenio, F.; Martin, J.; Marcello, J.; Fraile-Nuez, E. Environmental monitoring of El Hierro Island submarine volcano, by combining low and high resolution satellite imagery. International Journal of Applied Earth Observation and Geoinformation, 2014, vol. 29, p. 53-66. [Índice de impacto: 3.470] Congresos internacionales 1. Marcello, J.; Eugenio, F.; Marques, F.; Martín, J. Precise classification of coastal benthic habitats using high resolution Worldview-2 imagery. IEEE International Geoscience and Remote Sensing Symposium (IGARSS), 2015. p. 2307-2310. 2. Eugenio, F.; Marcello, J.; Martín, J. Submesoscale structures monitoring and detection by satellite imagery: El Hierro island submarine volcano. IEEE International Geoscience and Remote Sensing Symposium (IGARSS), 2014. p. 2371-2374. 3. Eugenio, F.; Martin, J.; Marcello, J.; Bermejo, J. Worldview-2 high resolution remote sensing image processing for the monitoring of coastal areas. Signal Processing Conference (EUSIPCO), 2013 Proceedings of the 21st European. p. 1-5. 4. Martin, J.; Eugenio, F.; Marcello, J.; Medina, A.; Bermejo, J.; Arbelo, M. Atmospheric correction models for high resolution WorldView-2 multispectral imagery: a case study in Canary Islands, Spain. SPIE Remote Sensing. International Society for Optics and Photonics, 2012. p. 85340O85340O-10. 5. Eugenio, F.; Marcello, J.; Martín, J. Monitoring El Hierro submarine volcano with low and high resolution satellite images. SPIE Remote Sensing. International Society for Optics and Photonics, 2012. p. 853816-853816-11. Congresos nacionales 1. Martín, J.; Eugenio, F.; Marcello, J. Estimación de las propiedades ópticas del agua, batimetría y albedo del fondo en ecosistemas litorales mediante imágenes multiespectrales de satélite de alta resolución. Congreso de la Asociación Española de Teledetección 2015, Sevilla; 08/2015. 2. Hernández Cordero, A.; Martín-Abasolo, J.; Leite, L.; Hernández Hernández, N.; Arístegui Ruiz, J.; Álvarez Vázquez, R.; Hernández-Calvento, L. Aplicación de TIG en la generación de indicadores de calidad ambiental de sistemas playa-dunas. 2014. 141 Capítulo 7 7. Listado de símbolos Símbolo Unidades Descripción g/m3 Pigmento de clorofila tipo A (tipo mayoritario) m-1 Materia disuelta en el agua (absorción por unidad métrica) g/m3 Materia suspendida en el agua ∆  μm Ancho de banda efectivo 󰆒󰇛  󰇜- Radiancia espectral relativa Wm-2µm-1 Irradiancia espectral solar promedio de la banda 󰇛  󰇜Wm-2µm-1 Irradiancia espectral solar - Niveles digitales de la imagen para cada banda - Niveles digitales de la imagen para cada banda Wm-2µm-1 Radiancia espectral (WV2) - Ganancia absoluta de la banda Lista de símbolos 142  - Offset del instrumento  - Offset del instrumento (WV2)  - Ganancia relativa de la banda  Wm -2 µm -1 Radiancia espectral multiplicada por ganancia absoluta (WV2)  Wm -2 µm -1 Radiancia a lo alto de la atmosfera ToA   - Factor de calibración radiométrico para cada banda WV2 , - Imagen radiométricamente corregida UA - Unidad astronómica (distancia promedio de la Tierra al Sol)  Wm -2 µm -1 Irradiancia solar media entre una distancia dada Tierra-Sol  Wm -2 µm -1 Irradiancia solar media para un ángulo cenital solar    Wm -2 µm -1 Irradiancia solar para una longitud de onda dada ,   Wm -2 µm -1 Irradiancia solar a lo alto de la atmosfera ,   Wm -2 µm -1 Irradiancia solar difusa que llega a la superficie ,   Wm -2 µm -1 Irradiancia solar que llega a la superficie  º Ángulo cenital solar  º Ángulo de elevación respecto al Sol  º Ángulo de visión del satélite  º Ángulo de azimut del Sol  º Ángulo de azimut del satélite  º Diferencia entre los ángulos de azimut solar-visión   - Coseno del ángulo cenital solar   - Coseno del ángulo de visión del satélite JD - Día Juliano   UA Distancia entre la Tierra y el Sol    Wm -2 µm -1 Radiancia total a lo alto de la atmósfera    Wm -2 µm -1 Radiancia reflejada por la superficie    Wm -2 µm -1 Radiancia dispersada por la atmósfera y reflejada por la superficie CAPÍTULO 7 143 ,   Wm -2 µm -1 Radiancia dispersada por la atmosfera    Wm -2 µm -1 Radiancia recibida por el sensor    Wm -2 µm -1 Radiancia mínima de la banda de la imagen   ,% Wm -2 µm -1 Radiancia de un objeto oscuro libre de sombra    - Reflectividad espectral difusa (ToA)  - Reflectividad atmosférica  - Reflectividad superficial (ToC)  - Reflectividad superficial del entorno homogéneo ,   - Reflectividad superficial simulada de las bandas WV2    - Reflectividad monocromática del radiómetro    - Función normalizada de respuesta del filtro de paso de las bandas WV2    - Transmisividad de la atmósfera, dirección ascendente    - Transmisividad de la atmósfera, dirección descendente  - Profundidad óptica de la atmósfera ,   - Transmisividad del ozono ,   - Transmisividad de los aerosoles ,   - Transmisividad de las moléculas (scattering de Rayleigh) 󰇛󰇜 - Transmisividad difusa de la atmósfera (descendente) 󰇛󰇜 - Transmisividad difusa de la atmósfera (ascendente)  - Transmisividad total de descendente y ascendente  - Espesor atmosférico (Atmospheric Optical Deep, AOD)  - Albedo efectivo de la atmósfera  - Inversa de la transmitancia atmosférica  - Scattering atmosférico  - Albedo atmosférico para la luz isotrópica  º Ángulo de reflexión respecto a la normal 󰆒 º Ángulo de refracción de la luz en la interfaz aire-agua