ES3013922T3 - Method for processing insar images to extract ground deformation signals - Google Patents
Method for processing insar images to extract ground deformation signals Download PDFInfo
- Publication number
- ES3013922T3 ES3013922T3 ES20157709T ES20157709T ES3013922T3 ES 3013922 T3 ES3013922 T3 ES 3013922T3 ES 20157709 T ES20157709 T ES 20157709T ES 20157709 T ES20157709 T ES 20157709T ES 3013922 T3 ES3013922 T3 ES 3013922T3
- Authority
- ES
- Spain
- Prior art keywords
- time series
- image
- images
- pixel
- filtered
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Classifications
-
- G—PHYSICS
- G06—COMPUTING OR CALCULATING; COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/70—Denoising; Smoothing
-
- G—PHYSICS
- G06—COMPUTING OR CALCULATING; COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/50—Image enhancement or restoration using two or more images, e.g. averaging or subtraction
-
- G—PHYSICS
- G06—COMPUTING OR CALCULATING; COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/60—Image enhancement or restoration using machine learning, e.g. neural networks
-
- G—PHYSICS
- G06—COMPUTING OR CALCULATING; COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10032—Satellite or aerial image; Remote sensing
- G06T2207/10044—Radar image
-
- G—PHYSICS
- G06—COMPUTING OR CALCULATING; COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/20—Special algorithmic details
- G06T2207/20024—Filtering details
-
- G—PHYSICS
- G06—COMPUTING OR CALCULATING; COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/20—Special algorithmic details
- G06T2207/20084—Artificial neural networks [ANN]
-
- G—PHYSICS
- G06—COMPUTING OR CALCULATING; COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/20—Special algorithmic details
- G06T2207/20172—Image enhancement details
- G06T2207/20182—Noise reduction or smoothing in the temporal domain; Spatio-temporal filtering
-
- G—PHYSICS
- G06—COMPUTING OR CALCULATING; COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/20—Special algorithmic details
- G06T2207/20212—Image combination
- G06T2207/20221—Image fusion; Image merging
-
- G—PHYSICS
- G06—COMPUTING OR CALCULATING; COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/30—Subject of image; Context of image processing
- G06T2207/30181—Earth observation
Landscapes
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Image Processing (AREA)
Abstract
La invención se refiere a un método para procesar series temporales de imágenes ruidosas de una misma área, comprendiendo el método: generar un conjunto de series temporales de imágenes a partir de una serie temporal de imágenes de entrada (RI[t]) combinando mediante primeras combinaciones lineales cada píxel (PX[i,j]) de cada imagen de la serie temporal de imágenes de entrada con píxeles vecinos seleccionados en la imagen y en una imagen adyacente de la serie temporal de imágenes de entrada; aplicar operaciones de filtrado (FPs, s=1-3, 5, 7-11) en cascada al conjunto, combinando cada operación de filtrado cada píxel de cada imagen de cada serie temporal del conjunto mediante segundas combinaciones lineales con píxeles vecinos seleccionados en la imagen y en una imagen adyacente en cada serie temporal del conjunto; realizar una operación de combinación de imágenes (SLs, s=4, 6) para reducir cada serie temporal del conjunto a una única imagen; introducir una imagen modelo (EM) del área como imagen filtrada en el conjunto; y combinar cada imagen del conjunto en una imagen de salida (PI), mediante terceras combinaciones lineales. (Traducción automática con Google Translate, sin valor legal)
Description
DESCRIPCIÓN
Método para procesar imágenes de InSAR para extraer señales de deformación del suelo
Campo técnico
La presente invención se refiere a un método y dispositivo para procesar imágenes tales como imágenes de radar de apertura sintética (SAR) procedentes de satélites de observación terrestre y, más particularmente, imágenes de SAR interferométrico (InSAR). La presente invención es aplicable a la extracción a partir de tales imágenes de información útil en relación, en especial, con deformaciones del suelo causadas por fenómenos naturales, incluyendo procesos tectónicos y actividad volcánica, o causadas por actividades humanas, incluyendo agotamiento de acuíferos, inyección de agua (incluyendo inyección de aguas residuales e inyección geotérmica) y actividades de extracción (incluyendo petróleo, gas y mineral). Más particularmente, la invención es aplicable a la detección de sucesos de deslizamiento lento a lo largo de fallas, a la detección de deslizamientos precursores a lo largo de fallas activas.
Antecedentes
Se da cabida al movimiento relativo entre placas tectónicas por deslizamiento sobre fallas activas que puede medirse, en principio, mediante InSAR. El paradigma clásico es que la mayoría de las fallas están bloqueadas la mayor parte del tiempo. El esfuerzo elástico acumulado puede liberarse durante terremotos grandes, mientras que el deslizamiento en algunas fallas o porciones de fallas puede ser predominantemente asísmico y, por lo tanto, suponer poco riesgo sísmico. Los estudios realizados a lo largo de las últimas dos décadas que revelan las apariciones de sucesos de deslizamiento inesperados han mostrado que este panorama es simplista en aspectos importantes. De hecho, recientemente se han identificado diferentes clases de fenómenos de deslizamiento, incluyendo terremotos de frecuencia baja, sucesos de deslizamiento lento y temblor tectónico asociado, así como deslizamiento asísmico continuo.
Los terremotos suponen una amenaza significativa para las regiones pobladas densamente y, hasta la fecha, es imposible proporcionar un pronóstico a corto o largo plazo para la incidencia de terremotos individuales. Se ha mostrado teóricamente que los terremotos deberían ir precedidos de una fase de nucleación, durante la cual una falla tectónica comienza a deslizarse lentamente y se acelera antes de romperse rápidamente. La detección de tal fase de nucleación permitiría pronosticar a corto plazo terremotos grandes. Sin embargo, algunas fallas se deslizan lentamente, de forma episódica, sin generar terremotos grandes y, sin un conjunto de observaciones lo bastante denso, es difícil diferenciar de forma eficiente un suceso de deslizamiento lento e inofensivo de una fase de nucleación precursora de un suceso devastador. Además, el deslizamiento asísmico y el sísmico están estrechamente relacionados: se cree que la fluencia de las fallas precede y, a menudo, sigue a los terremotos, lo que subraya la necesidad de un avance importante en la descripción mecánica y física del deslizamiento de fallas en general.
La interferometría de radar de apertura sintética (InSAR), desarrollada a principios de la década de 1990, se usa en la actualidad de forma rutinaria para medir deformaciones del suelo debidas a procesos hidrológicos, volcánicos y tectónicos. Los satélites en órbita proporcionan en la actualidad imágenes diarias de la superficie de la Tierra y, usando el concepto de interferometría con un sistema de radar, pueden obtenerse mapas de la evolución de la distancia aparente entre el suelo y el sistema de satélites. Estos mapas pueden ensamblarse para proporcionar una monitorización continua de tal distancia. Esta distancia aparente está relacionada directamente con el movimiento del suelo a lo largo del tiempo, por lo tanto hoy en día es posible detectar cualquier suceso de deslizamiento lento, precursor o no de un terremoto grande, en cualquier parte de la Tierra. La constelación de satélites de InSAR Sentinel 1 tiene el potencial de posibilitar una cartografía sistemática de todas las regiones que presentan deformación a nivel mundial, revelando mecanismos subyacentes a todos los modos de deslizamiento de fallas tectónicas, con la condición de que pueda desarrollarse un método robusto y automático para separar las señales atmosféricas de la deformación de fallas. Antes de la constelación Sentinel 1, la extensión típica de una imagen de SAR era de 100 x 100 km para los satélites clásicos de banda C y L (ERS, Envisat, ALOS...), con un tamaño de píxel del orden de decenas de metros (dependiendo de la longitud de onda del radar y la estructuración de antena). Con la constelación Sentinel, el tamaño de las imágenes es, como promedio, 3 veces más grande y las adquisiciones son, como mínimo, 5 veces más frecuentes. El lanzamiento de la constelación de satélites Sentinel 1 es transformador en el sentido de que proporciona una cartografía de radar sistemática de todas las regiones que presentan deformación activa en el mundo con un período de retorno de 6 días. Los lanzamientos de satélites venideros (NISAR, etc.) proporcionarán datos de InSAR de una resolución y una frecuencia aún más altas, ahondando en la necesidad de herramientas de detección automática.
La técnica de InSAR se ha aplicado con éxito para monitorizar desplazamientos grandes debidos a terremotos, el movimiento de la capa de hielo, desplazamientos más pequeños en relación con el agotamiento de acuíferos, la deformación intersísmica, los corrimientos de tierra de movimiento lento y los sucesos de deslizamiento lento. Las señales de deformación rápidas y de gran amplitud, tales como campos de desplazamiento cosísmico o episodios volcánico-tectónicos, también se analizan de forma rutinaria mediante InSAR. De forma similar, la acumulación lenta pero ininterrumpida de deformación a lo largo de períodos de tiempo relativamente largos también se detecta fácilmente usando InSAR, debido a que la deformación resultante también es lo suficientemente grande. Sin embargo, estas mediciones generalmente solo han tenido éxito a través de una exploración manual minuciosa de los conjuntos de datos de deformación por expertos, evitando de este modo los estudios a gran escala.
La combinación de dos imágenes de SAR en un interferograma (InSAR) permite medir el movimiento del suelo a lo largo de un área amplia. La fase del interferograma varía como una función de cualquier fenómeno físico que afecte al tiempo de desplazamiento bidireccional de la onda de radar entre el soporte de radar y el suelo. Estos efectos incluyen cambios de órbita, así como cambios en la atmósfera y la deformación del suelo entre adquisiciones. Debido a que las órbitas son bien conocidas, la fuente más importante de ruido que obstaculiza las mediciones de deformación es la variación espacial y temporal de la temperatura, la presión y el vapor de agua en la atmósfera. Los retardos troposféricos e ionosféricos generan errores grandes en el procesamiento de InSAR, equivaliendo a menudo a errores de 10-20 cm, que pueden ser considerablemente más grandes que las señales de interés. Los retardos en las señales de radar se deben principalmente al vapor de agua presente en la troposfera. Estos retardos tienen dos causas identificadas principales: 1) turbulencias en la troposfera y 2) un componente altamente relacionado con la topografía local. En la bibliografía se han desarrollado dos enfoques primarios para aliviar esta cuestión y reducir el ruido. El primer enfoque consiste en alisar o filtrar en el tiempo, con la suposición de que el ruido atmosférico es de naturaleza estocástica y puede promediarse. Estos enfoques incluyen el apilamiento de análisis de series temporales de interferogramas basándose en métodos de filtrado y correlación. Sin embargo, el ruido atmosférico es, de hecho, no estocástico tanto en el tiempo como en el espacio, y el alisado conduce a una disminución de la resolución de las señales procedentes del suelo. En otras palabras, se ha mostrado que estos métodos enmascaran parcialmente las señales de interés.
El segundo enfoque se basa en conjuntos de datos adicionales para estimar retardos y mitigar el ruido en estimaciones de deformación del suelo, incluyendo datos meteorológicos, mapas de vapor de agua a partir de datos de espectrorradiómetros y espectrómetros, datos de posición geográfica (a partir de un sistema mundial de navegación por satélite) y estimaciones a partir de modelos y simulaciones meteorológicos. Todos estos diferentes conjuntos de datos tienen sus virtudes e inconvenientes. Los modelos atmosféricos pueden adolecer de errores sistemáticos o localizados. Los mapas de vapor de agua tienen una resolución espacial alta pero una resolución temporal relativamente deficiente, y son limitados cuando hay nubes presentes. En cambio, los datos de posición geográfica proporcionan las estimaciones con la resolución temporal más alta de retardos atmosféricos, pero requieren estaciones en la superficie de la Tierra, y los datos están limitados espacialmente a las estaciones disponibles (la interpolación entre estaciones no es obvia y es difícil evaluar su calidad). En la actualidad, no existe un enfoque fiable y usado universalmente para mitigar el impacto de los errores inducidos por la atmósfera en las estimaciones del movimiento del suelo de InSAR.
Por lo tanto, la detección de campos de deformación de amplitud baja y longitud de onda larga, tales como los debidos a acumulaciones de deformaciones intersísmicas o movimientos post-sísmicos, sigue siendo un desafío debido a la descorrelación interferométrica, las órbitas imprecisas de los satélites y los retardos de propagación atmosférica: si la atmósfera ralentiza la señal de satélite, parece que la distancia entre el suelo y el satélite aumenta, como si se hubiera movido el suelo. Para lograr este objetivo, deben hacerse avances significativos en la medición de InSAR de la deformación lenta.
Por lo tanto, aún no es posible la detección sistemática de estos sucesos de deslizamiento a lo largo de fallas tectónicas, aunque los datos están disponibles en las imágenes de InSAR. De hecho, las soluciones existentes se basan en la detección por seres humanos a través de inspección visual de series temporales de InSAR de imágenes de movimiento del suelo, donde la señal es una mezcla de ruidos atmosféricos y señales tectónicas.
La publicación “CNN-based InSAR Denoising and Coherence Metric” , de Suhhayan Mukherjee y col., Biblioteca de la Universidad de Cornell, 201 Biblioteca Olin, Universidad de Cornell, Ithaca, NY 14853, 20 de enero de 2020, describe un método para procesar una única imagen de InSAR usando una red neuronal convolucional para la reducción de ruido y artefactos.
La solicitud de patente US 2019/0005603 describe un método para procesar una única imagen mediante una red neuronal convolucional para aplicar un operador de procesamiento de imagen objetivo a la imagen.
La solicitud de patente EP 3229 038 describe un método para procesar una única imagen de InSAR realizando un filtrado de fase en el dominio de ondículas y una estimación de frecuencia local.
La solicitud de patente US 2008/0231504 describe un método para procesar imágenes de InSAR para reducir el ruido y preservar los bordes formando un interferograma a partir de dos imágenes de InSAR de la misma escena, estando desfasadas las dos imágenes. Después, el interferograma se somete a un complejo algoritmo de difusión anisotrópica y a un filtro de choque y a un desenvolvimiento de fase variacional bidimensional.
La solicitud de patente WO 2011/068783 describe un método para procesar una serie temporal de imágenes de IRM (formación de imágenes por resonancia magnética) para mejorar la relación señal/ruido de las imágenes,
La solicitud de patente CN 109541596 describe un método de procesamiento de una única imagen de InSAR basándose en un algoritmo de aprendizaje profundo para obtener una imagen de InSAR de alta precisión con ruido reducido.
En consecuencia, existe la necesidad de procesar automáticamente series temporales masivas de conjuntos de imágenes tales como series temporales de imágenes de InSAR para eliminar o atenuar las perturbaciones atmosféricas, con vistas a resaltar las señales de deformación más pequeñas posibles de todas las clases. También existe la necesidad de una detección y una caracterización sistemáticas y globales de las deformaciones del suelo de amplitud baja y/o longitud de onda larga a diferentes escalas, con vistas a diferenciar con seguridad entre sucesos de deslizamiento lento inofensivos y sucesos pequeños en relación con fases precursoras de terremotos, también puede ser deseable identificar la deformación de fallas sin requerir la interpretación de un experto.
Resumen
Se describe un método para procesar series temporales de imágenes de una misma área, sometidas a ruido. El método puede incluir: ejecutar por una red neuronal convolucional una pluralidad de etapas de filtrado sucesivas que comprenden una primera etapa de filtrado que recibe una serie temporal de imágenes de entrada, una última etapa de filtrado que proporciona una imagen de salida y etapas de filtrado intermedias, transformando, cada una de la primera etapa de filtrado y las etapas de filtrado intermedias, un conjunto de series temporales de imágenes filtradas que incluyen inicialmente la serie temporal de imágenes de entrada, transformándose el conjunto generando un número de series temporales de imágenes filtradas combinando mediante combinaciones lineales cada píxel de cada imagen de cada serie temporal del conjunto con píxeles vecinos seleccionados en la imagen y en una imagen adyacente de las series temporales de imágenes del conjunto, usando cada combinación lineal un conjunto respectivo de coeficientes de ponderación y dando como resultado un píxel de una de las imágenes filtradas generadas; y ejecutar por la red neuronal convolucional una o más operaciones de combinación, realizándose cada operación de combinación entre dos sucesivas de las etapas de filtrado intermedias, para reducir un número de imágenes en cada una de las series temporales de imágenes filtradas del conjunto, combinando imágenes de subconjuntos de imágenes adyacentes en cada una de las series temporales de imágenes filtradas del conjunto, reduciendo una última de las operaciones de combinación cada serie temporal de imágenes filtradas del conjunto a una única imagen filtrada.
Según una realización, el método comprende además, después de la última operación de combinación, introducir una imagen de modelo del área en el conjunto como una imagen filtrada adicional, yendo seguida la introducción de la imagen de modelo por una o más de las etapas de filtrado.
Según una realización, cada una de las combinaciones lineales incluye un coeficiente de sesgo que se añade a un resultado de la combinación lineal.
Según una realización: cada uno de los coeficientes de ponderación tiene un valor que depende de un signo, positivo o negativo, de un valor de píxel por el que se multiplica este, o el resultado de cada una de las combinaciones lineales se transforma mediante una función unitaria lineal rectificada.
Según una realización, el método comprende además: generar la primera serie temporal de imágenes a partir de un modelo de deformación del suelo usando parámetros seleccionados aleatoriamente; generar series temporales de imágenes de entrenamiento, a partir de la primera serie temporal de imágenes, usando modelos de diferentes señales de ruido; y usar la serie temporal de entrenamiento generada para ajustar los valores de los coeficientes de ponderación.
Según una realización, los valores de los coeficientes de ponderación se ajustan iterativamente aplicando un método de minimización de descenso basado en gradiente iterativo usando las series temporales de imágenes de entrenamiento y una función de coste de modelo, para minimizar el resultado de la función de coste de modelo.
Según una realización, la generación del conjunto de series temporales de imágenes filtradas a partir de la serie temporal de imágenes de entrada se realiza usando la siguiente ecuación:
en donde PX[i,j,t,f] es un píxel de una de las imágenes filtradas de la serie temporal f de imágenes filtradas, PX[i,j,t] es un píxel de una imagen de la serie temporal de imágenes de entrada, W[l,m,u,f] es uno de los primeros coeficientes para la serie temporal f de imágenes filtradas, B[f] es un coeficiente de sesgo para la serie temporal f, y LR() es una función unitaria lineal rectificada.
Según una realización, cada combinación lineal de primeras operaciones de filtrado de las operaciones de filtrado aplica la siguiente ecuación:
en donde PX[i,j,t,f'] es un píxel de una de las imágenes filtradas de la serie temporal de imágenes filtradas f', W[s,l,m,u,f,f'] es uno de los segundos coeficientes de ponderación para la serie temporal f de imágenes filtradas para la operación de filtrado s, B[s,f'] es un coeficiente de sesgo para la serie temporal f y la operación de filtrado s, y LR() es una función unitaria lineal rectificada con fugas.
Según una realización, cada combinación lineal de segundas operaciones de filtrado de las operaciones de filtrado aplica la siguiente ecuación:
en donde PX[i,j,f] es un píxel de una de las imágenes filtradas f', W[s,l,m,f,f’] es uno de los segundos coeficientes de ponderación para las imágenes filtradas f para la operación de filtrado s, B[s,f'] es un coeficiente de sesgo para las imágenes filtradas f y la operación de filtrado s, y l R() es una función unitaria lineal rectificada con fugas.
Según una realización, cada píxel de la imagen de salida (PI) se calcula mediante la siguiente ecuación:
en donde PX[i,j] es un píxel de la imagen de salida, PX[i,j,f] es un píxel de la imagen filtrada f, W[l,m,f] es uno de los terceros coeficientes, B es un coeficiente de sesgo y LR() es una función unitaria lineal rectificada con fugas.
Según una realización, la función unitaria lineal rectificada con fugas LR es de tal modo que LR(x) = x si x > 0 y LR(x) = 0,5x si x < 0.
Según una realización, cada una de las operaciones de combinación de imágenes aplica la siguiente ecuación:
en donde PX[i,j,t,f] es un píxel de una de las imágenes filtradas de la serie temporal de imágenes filtradas f, y MÁX (PX[i,j,t+u,f]) es una función que proporciona el valor máximo de entre los valores de píxel PX[i,j,t+u,f] con u = -1, 0 y 1.
Las realizaciones también pueden referirse a un ordenador que comprende un procesador, una memoria, estando configurado el procesador para llevar a cabo el método descrito anteriormente.
Según una realización, el ordenador comprende además un procesador gráfico controlado por el procesador, estando configurado el procesador para configurar el procesador gráfico para llevar a cabo algunas de las operaciones del método descrito anteriormente.
Las realizaciones también pueden referirse a un producto de programa informático cargable en una memoria informática y que comprende porciones de código que, cuando son llevadas a cabo por un ordenador, configuran el ordenador para llevar a cabo el método descrito anteriormente.
Breve descripción de los dibujos
El método y/o el dispositivo pueden entenderse mejor con referencia a los siguientes dibujos y descripción. Se describen descripciones no limitantes y no exhaustivas con los siguientes dibujos. En las figuras, símbolos de referencia semejantes pueden referirse a partes semejantes en la totalidad de las diferentes figuras, salvo que se especifique lo contrario.
La figura 1 es un diagrama de bloques de etapas de procesamiento de un ordenador que implementa un método para procesar imágenes, según una realización;
la figura 2 es un diagrama de bloques de una primera etapa de procesamiento del procesador, según una realización; la figura 3 es un diagrama de bloques que ilustra la operación realizada por la primera etapa de procesamiento para calcular un valor de píxel, según una realización;
la figura 4 es un diagrama de bloques de una segunda etapa de procesamiento del procesador, según una realización; la figura 5 es un diagrama de bloques que ilustra la operación realizada por la segunda etapa de procesamiento para calcular un valor de píxel, según una realización;
la figura 6 es un diagrama de bloques de una tercera etapa de procesamiento del procesador, según una realización; la figura 7 es un diagrama de bloques que ilustra la operación realizada por la tercera etapa de procesamiento para calcular un valor de píxel, según una realización;
la figura 8 es un diagrama de bloques de una cuarta etapa de procesamiento del procesador, según una realización; la figura 9 es un diagrama de bloques que ilustra la operación realizada por la cuarta etapa de procesamiento para calcular un valor de píxel, según una realización;
la figura 10 es un diagrama de bloques de una quinta etapa de procesamiento del procesador, según una realización; la figura 11 es un diagrama de bloques que ilustra la operación realizada por la quinta etapa de procesamiento para calcular un valor de píxel, según una realización;
la figura 12 es un diagrama de bloques de un ordenador personal clásico que puede configurarse para implementar el método de procesamiento de imagen.
Descripción detallada
Se describe un método para procesar imágenes tales como imágenes de InSAR. Este método comprende:
1) reconstrucción de imágenes coherentes en términos de fase;
2) corrección de las imágenes para tener en errores de posiciones de satélite;
3) eliminación de ruido atmosférico; y
4) reconstrucción de series temporales de deformación.
Una realización de software en un ordenador se describe a continuación como un conjunto de componentes de programa legibles por ordenador que cooperan para controlar el rendimiento de operaciones de procesamiento de datos cuando se cargan y se ejecutan en el ordenador. Será evidente para un experto en la técnica que las etapas individuales de métodos de la presente invención pueden implementarse en código de programa informático y que puede usarse una diversidad de lenguajes de programación e implementaciones de codificación para implementar los métodos descritos en la presente memoria. Además, los programas informáticos incluidos en el software no pretenden limitarse a los flujos de control específicos descritos en la presente memoria, y una o más de las etapas de los programas informáticos pueden realizarse en paralelo o secuencialmente. Una o más de las operaciones descritas en el contexto de una implementación controlada por programa informático podrían implementarse alternativamente como un componente electrónico de hardware.
La figura 1 ilustra operaciones de un método realizado por un ordenador según una realización, para eliminar el ruido atmosférico de las imágenes con los errores de posiciones de satélite corregidos (la etapa 3). La figura 1 muestra las<etapas de procesamiento FP0, FP1, FP2, FP3, SL4, FP5, SL6, FP7, FP8, FP9, FP10,>F<p>1<1, FP12 conectadas en>serie, de una red neuronal convolucional CNN configurada para procesar sucesivamente una serie temporal RI de imágenes para eliminar el ruido y, especialmente, el ruido atmosférico. Cada imagen RI[t] de la serie temporal RI comprende un conjunto de píxeles PX[i,j] donde i = 1, ... X y j = 1, ...Y. La serie temporal de imágenes RI comprende T imágenes RI[1], ... RI[T] de una misma área de escena (por ejemplo, la misma región de suelo) tomada en diferentes momentos o adquirida en intervalos de tiempo dados. La serie temporal de imágenes RI se procesa mediante una primera etapa de procesamiento, proporcionándose las imágenes calculadas por cada una de las etapas a la siguiente etapa, hasta una última etapa de procesamiento o salida FP12 que proporciona una imagen de salida PI en la que el ruido atmosférico está ausente o atenuado fuertemente, pero en la que siguen estando presentes señales de movimiento del suelo presentes en la serie temporal de entrada RI de imágenes.
La figura 2 ilustra la primera etapa de procesamiento o entrada FP0, según una realización. La etapa de procesamiento FP0 comprende F componentes de filtrado FT[0,f], con f igual a 0 a F, procesando cada componente de filtrado FT[0,f] la serie temporal de entrada RI. “ 0” en la notación FT[0,f] se refiere a un número asignado a la primera etapa de procesamiento FP0. Cada componente de filtrado FT[0,f] genera una serie temporal FI[0,f] de T imágenes filtradas FI[0,f,t]. Por lo tanto, la etapa de procesamiento FP0 genera F series temporales de imágenes filtradas FI[0,f,t].
La figura 3 ilustra la operación realizada por cada componente de filtrado FT[0,f] de la primera etapa de procesamiento FP0, según una realización. Cada píxel PX[i,j,t,f] de cada imagen filtrada FI[0,t,f] es calculado por un componente de filtrado FT[0,f] correspondiente realizando una combinación lineal de píxeles vecinos seleccionados (píxeles de color negro en las imágenes RI[t] y RI[t+1] en la figura 3) en las proximidades del píxel PX[i,j,t] en la tabla de píxeles tridimensional PX que incluye todos los píxeles de la serie temporal de entrada R<i>. Puede añadirse un coeficiente de sesgo al resultado de la combinación lineal y el resultado de la suma puede transformarse mediante una función unitaria lineal rectificada (ReLU), que puede ser del tipo “ con fugas” . La función ReLU es una función de activación responsable de transformar la entrada ponderada sumada a partir de la etapa en la activación de la etapa o salida para esa entrada. En general, la función ReLU es una función lineal por tramos que emitirá la entrada directamente si es positiva, de lo contrario, emitirá cero. El fin es eliminar un problema conocido como “ problema de desvanecimiento de gradiente” , que puede evitar que las redes profundas aprendan de forma eficaz. Podrían usarse otras funciones de activación. Una función ReLU con fugas mantiene los valores negativos con un coeficiente de atenuación lineal.
Según una realización, cada píxel PX[i,j,t,f] (i = 1, ...X, j = 1, ...Y, t = 1, ...T, f = 1, ...F) de cada imagen filtrada FI[0,t,f] se calcula mediante el componente de filtrado Ft[0,f] correspondiente, según la siguiente ecuación:
donde LR es la función ReLU con fugas, B[0] es una tabla unidimensional de coeficientes de sesgo B[0,f], representando cada uno de los mismos un valor constante usado por el componente de filtrado FT[0,f], W[0] es una tabla tetradimensional de coeficientes de ponderación W[0,l,m,u,f] que tiene las siguientes dimensiones (1..L, 1..M, 1..U, 1..F), teniendo los coeficientes de ponderación W[0,l,m,u,f] valores reales, y a, b, c son coeficientes que pueden ajustarse con los parámetros L, M y U para seleccionar en la tabla PX los píxeles vecinos que están implicados en las combinaciones lineales. Los píxeles fuera de las imágenes RI[t], es decir, cuando i+a.l > X, y/o j+b.m > Y, y/o t+c.u > T se establecen a 0.
Según un ejemplo, LR(x) = x si x > 0 y LR(x) = 0,5x si x < 0. En el ejemplo de la figura 3, L y M se establecen a 3, U se establece a 2, a y b se establecen a 3 y c se establece a 1. Por lo tanto, cada píxel PX[i,j,t,f] de cada imagen filtrada FI[0,t,f] se calcula a partir de 9 píxeles vecinos en cada una de las imágenes RI[t] y RI[t+1] de la serie temporal de imágenes RI[1..T].
La figura 4 ilustra las etapas de procesamiento FP1, FP2, FP3 y FP5, según una realización. Cada una de las etapas de procesamiento FP1-FP3, FP5 comprende F componentes de filtrado FT[s,f], con s igual a 1, 2, 3 o 5, y f igual a 1 a F, procesando cada componente de filtrado FT[s,f] la serie temporal de imágenes filtradas FI[s-1 ,t]. Cada componente de filtrado FT[s,f] genera series temporales de T imágenes filtradas FI[s,t,f]. Por lo tanto, cada etapa de procesamiento FPs (s: 1, 2, 3 o 5) genera F series temporales de T imágenes filtradas.
La figura 5 ilustra la operación realizada por cada componente de filtrado FT[s,f] de la etapa de procesamiento FP[s] (s = 1, 2, 3 o 5), según una realización. Cada píxel PX[i,j,t,f] de cada imagen filtrada FI[s,t,f] es calculado por un componente de filtrado FT[s,f] correspondiente realizando una combinación lineal de píxeles vecinos seleccionados (píxeles de color negro en las imágenes FI[s-1 ,t,f] y FI[s-1 ,t+1,f] en la figura 5) en las proximidades del píxel PX[i,j,t,f] en la tabla de píxeles tridimensional PX[i,j,t,f] en todas las F series temporales de imágenes filtradas FI[s,t,f]. De nuevo, puede añadirse un coeficiente de sesgo al resultado de la combinación lineal y el resultado de la suma puede transformarse mediante la misma función unitaria lineal rectificada (ReLU).
Según una realización, cada píxel PX[i,j,t,f] de cada imagen filtrada FI[s,t,f] (f igual a 1 a F) se calcula mediante el componente de filtrado Ft[s,f] correspondiente, según la siguiente ecuación:
donde B[s,f'] es una tabla bidimensional de coeficientes de sesgo B[s,f], cada uno con un valor constante atribuido al componente de filtrado FT[s,f], W[s] (s = 1, 2, 3 o 5) es una tabla pentadimensional de coeficientes de ponderación W[s,l,m,u,f,f] que tiene las siguientes dimensiones (1..L, 1..M, 1..U, 1..F, 1..F), teniendo los coeficientes de ponderación W[s,l,m,u,f,f] valores reales, y a, b, c son coeficientes que pueden ajustarse con los parámetros L, M y U para seleccionar en la tabla PX los píxeles vecinos que están implicados en las combinaciones lineales.
Según un ejemplo, LR(x) = x si x > 0 y LR(x) = 0,5x si x < 0. En el ejemplo de la figura 5, L y M se establecen a 3, y U se establece a 2, a y b se establecen a 3 y c se establece a 1. Por lo tanto, cada píxel PX[i,j,t,f] de cada imagen filtrada FI[s,t,f] se calcula a partir de 9 píxeles vecinos (desde el píxel PX[i,j]) en todos los conjuntos de imágenes FI[s-1 ,t,f] y Fl[s-1 ,t+1 ,f], con f = 1, ... F, es decir, cada píxel PX[i,j,t,f] de cada imagen filtrada FI[s,t,f] se calcula a partir de 3 x 3 x 2 x F píxeles.
Las figuras 6 y 7 ilustran las etapas de procesamiento SL4 y SL6 del procesador, según una realización. Cada una de las etapas de procesamiento<s>L4 y SL6 comprende F que componentes de combinación CB[s,f] (con s igual a 4 o 6 y f igual a 1 a F), procesando cada componente de combinación CB[s,f] la serie temporal de imágenes filtradas FI[s-1 ,t] y proporcionando una serie temporal de T/n imágenes filtradas FI[s,t,f], siendo T un múltiplo de n. Las operaciones de combinación realizadas por el componente de combinación CB[s,f] combina todos los grupos de píxeles PX[i,j,t], PX[i,j,t+1], ... [i,j,t+k] en cada una de las F series temporales FI[s-1,f] de imágenes filtradas, para dar píxeles resultantes individuales respectivos, considerándose cada píxel de cada una de las imágenes filtradas FI[s-1 ,t,f] solo una vez para calcular la serie temporal FI[s,f], y representando k el número de imágenes filtradas FI[s,t,f] combinadas por el componente de combinación CB[s,f].
Según una realización, cada píxel PX[i,j,t,f] de cada imagen filtrada FI[s,t,f] (f igual a 1 a F) se calcula a partir de tres imágenes filtradas FI[s-1 ,t], Fl[s-1,t+1] y FI[s-1 ,t+2] mediante el componente de combinación correspondiente CB[s,f], según la siguiente ecuación:
PX[¡,j,t,f] = MAX (PX[i,j,3(t-1 )+u,f]>
donde MÁX es una función que proporciona el valor mayor en cada grupo de tres píxeles PX[i,j,3(t-1)+1 ], PX[i,j,3(t-1)+2] y PX[i,j,3(t-1 )+3], siendo t igual a 1, 2, ... T/3. Como resultado, cada una de las series temporales de imágenes filtradas FI(s,f) proporcionadas por las etapas de procesamiento SL4 y SL6 comprende T/3 imágenes, siendo T un múltiplo de 3.
Después de la etapa de procesamiento SL6, solo queda una imagen FI(s,f) en cada una de las series temporales de imágenes filtradas FI (s,f,t).
La figura 8 ilustra la etapa de procesamiento FP7, FP8, FP9, FP10 del procesador, según una realización. Cada una de las etapas de procesamiento FP7-FP10 comprende F componentes de filtrado FT[s,f], con s igual a 7, 8, 9 o 10, y f igual a 1 a F, procesando cada componente de filtrado FT[s,f] la imagen filtrada FI[s-1,f] y generando una imagen filtrada FI[s,f].
La figura 9 ilustra la operación realizada por cada componente de filtrado FT[s,f] de la etapa de procesamiento FP[s] (s = 7, 8, 9, 10 u 11), según una realización. Cada píxel PX[i,j,f] de cada imagen filtrada FI[s,f] es calculado por un componente de filtrado FT[s,f] correspondiente realizando una combinación lineal de píxeles vecinos seleccionados (píxeles de color negro en la imagen Fl[s-1 ,f] en la figura 9) en las proximidades del píxel PX[i,j,f] en la tabla de píxeles tridimensional PX en todas las F imágenes filtradas Fl[s,f]. De nuevo, puede añadirse un coeficiente de sesgo al resultado de la combinación lineal y el resultado de la suma puede transformarse mediante la misma función unitaria lineal rectificada (ReLU).
Según una realización, cada píxel PX[i,j,f'] de cada imagen filtrada Fl[s,f'] (f igual a 1 a F) se calcula mediante el componente de filtrado Ft[s,f] correspondiente, según la siguiente ecuación:
donde B[s] es una tabla unidimensional de coeficientes de sesgo B[s,f], cada uno con un valor constante atribuido al componente de filtrado FT[s,f'], W[s] son tablas tetradimensionales de coeficientes de ponderación W[s,l,m,f,f'], con s = 7, 8, 9, 10 u 11, que tienen las siguientes dimensiones (1..L, 1..M, 1..F, 1..F), teniendo los coeficientes de ponderación W[s,l,m,f,f'] valores reales, y a, b, c son coeficientes que pueden ajustarse con los parámetros L, M y U para seleccionar en la tabla PX los píxeles vecinos que están implicados en las combinaciones lineales. En el ejemplo de la figura 9, L y M se establecen a 3, y a y b se establecen a 3. Por lo tanto, cada píxel PX[i,j] de cada una de las imágenes filtradas FI[s,f] se calcula a partir de 9 píxeles vecinos (desde la posición[i,j]) en todas las imágenes FI[s-1 ,f] (f = 1 a F).
La imagen EM del modelo de elevación del área del suelo sometida a formación de imágenes por la serie temporal de imágenes RI[t] se introduce en la etapa de procesamiento FP7 como una imagen filtrada que tiene el índice F+1. Por lo tanto, en los cálculos realizados en la etapa de procesamiento FP7, la suma del índice f en la ecuación (4) se realiza entre 1 y F+1, y la tabla tetradimensional W[7] tiene las siguientes dimensiones (1..L, 1..M, 1..F+1, 1..F). En las etapas de procesamiento FP8-FP11, las tablas tetradimensionales W[s] tienen las siguientes dimensiones (1..L, 1..M, 1..F, 1..F).
La figura 10 ilustra la última etapa de procesamiento o salida FP12 del procesador, según una realización. La etapa de salida FP12 comprende un componente de filtrado FT[12] que procesa las imágenes filtradas FI[11 ,f] para generar la imagen de salida PI.
La figura 11 ilustra la operación realizada por el componente de filtrado FT[12] de la etapa de salida FP12 según una realización. Cada píxel Pl[i,j] de la imagen de salida PI es calculado por el componente de filtrado FT[12] realizando una combinación lineal de píxeles vecinos seleccionados en las proximidades del píxel calculado en la tabla de píxeles tridimensional PX[i,j,f] que comprende todos los píxeles de las F imágenes filtradas FI[11 ,f]. De nuevo, puede añadirse un coeficiente de sesgo al resultado de la combinación lineal y el resultado de la suma puede transformarse mediante la misma función unitaria lineal rectificada (ReLU).
Según una realización, cada píxel PX[i,j] de la imagen de salida PI es calculado por el componente de filtrado FT[11], según la siguiente ecuación:
donde W[12] es una tabla tridimensional de coeficientes de ponderación W[12,l,m,f] que tiene las siguientes dimensiones (1..L, 1..M, 1..F], teniendo los coeficientes de ponderación W[12,l,m,f] valores reales, y a, b, c son coeficientes que pueden ajustarse con los parámetros L, M y U para seleccionar en la tabla PX los píxeles vecinos que están implicados en las combinaciones lineales. En el ejemplo de la figura 11, L y M se establecen a 3, y a y b se establecen a 3. Por lo tanto, cada píxel PX[i,j] de la imagen de salida PI se calcula a partir de 9 píxeles vecinos en cada una de las imágenes FI[11,f] (f = 1 a F).
La elección del número F de componentes de filtrado FT[s,f] en las etapas FPs (s = 1-3, 5, 7-11) resulta de una compensación recíproca entre la eficiencia de la red neuronal para eliminar el ruido de las imágenes de entrada y el tiempo de cálculo. En general, el número F se establece entre 10 y 100. En los ejemplos anteriores, F se establece a 64. Por lo tanto:
W[0] incluye 2 x 3 x 3 x 64 (= 1162) coeficientes de ponderación y B[0,f] incluye 64 coeficientes de sesgo,
W[s] y B[s,f], con s = 1, 2, 3 y 5, incluyen respectivamente 2 x 3 x 3 x 64 x 64 (= 73728) coeficientes de ponderación y 64 coeficientes de sesgo,
W[7] incluye 3 x 3 x 64 x 65 (= 37440) coeficientes de ponderación y B[7,f] incluye 64 coeficientes de sesgo,
W[s] y B[s,f], con s = 8, 9, 10 y 11, incluyen respectivamente 3 x 3 x 64 x 64 (= 36864) coeficientes de ponderación y 64 coeficientes de sesgo, y
W[12] incluye 3 x 3 x 64 (= 576) coeficientes de ponderación, y B[12] representa un coeficiente de sesgo.
Según una realización, todos los coeficientes W y B se determinan por medio de una fase de aprendizaje usando datos de entrenamiento. Según una realización, los datos de entrenamiento comprenden millones o miles de millones de series temporales de imágenes que se generan a partir de simulaciones de movimiento del suelo usando modelos elásticos simples de fallamiento. La serie temporal de imágenes se genera, por ejemplo, generando aleatoriamente una falla que se desliza lentamente con latitud, longitud, profundidad, ángulo de rumbo, ángulo de buzamiento y anchura aleatorios. Después, las series temporales de deformación del suelo se corrompen con una diversidad de señales de ruido. Por ejemplo, un ruido gaussiano correlacionado espacialmente simula retardos introducidos en las imágenes de InSAR por las turbulencias atmosféricas de diversas escalas de longitud. Además, puede usarse una diversidad de señales espurias adicionales para corromper adicionalmente la serie temporal de deformación del suelo. Por ejemplo, pueden añadirse señales de fallamiento transitorias para simular el escenario más desfavorable de señales relacionadas con el clima muy marcadas que parecen señales de fallamiento, pero sin correlación a lo largo del tiempo. También pueden introducirse píxeles o parches de píxeles erróneos que simulan errores de desenvolvimiento y descorrelación de píxeles en la serie temporal de imágenes.
Cada serie temporal de imágenes de entrenamiento está asociada con la imagen final que va a ser proporcionada por la red neuronal CNN. El entrenamiento de la red neuronal comprende hallar una solución a un sistema de ecuaciones en el que las incógnitas son los coeficientes W y B, conociéndose los píxeles PX[i,j] de la serie temporal de imágenes
RI y de la imagen resultante PI a partir de un caso de entrenamiento. Según una realización, una solución se calcula usando un método de minimización de descenso basado en gradiente iterativo usando una función de coste de modelo, seleccionándose los valores iniciales para los coeficientes W y B aleatoriamente entre -1 y 1, tal como para tener una distribución uniforme.
Según una realización, la función de coste de modelo MCF usada es la norma de error de la reconstrucción de deformación, por ejemplo, la suma de los valores absolutos de los errores de reconstrucción para cada píxel:
donde CPX[i,j] representa un píxel calculado en la imagen PI proporcionada por la red neuronal CNN a partir de una serie temporal de imágenes de entrenamiento, y MPX[i,j] representa el píxel correspondiente en la imagen final que va a proporcionarse, asociada con la serie temporal de imágenes de entrenamiento PI. Los coeficientes W y B se ajustan iterativamente para minimizar el resultado de la función de coste de modelo MCF sobre grupos de datos de entrenamiento. Cuando se usa una función de coste de este tipo que no es convexa, es casi seguro que nunca alcanzará un mínimo global, sino solo mínimos locales dependiendo de los valores iniciales. De este modo, se asegura alcanzar un conjunto singular de coeficientes W, B dependiendo de los datos de entrenamiento.
Según una realización, los coeficientes W y B se actualizan en cada iteración, basándose en el gradiente de la función de coste calculado a partir de un lote de datos, siguiendo la regla de ADAM, según la siguiente ecuación:
m
w = w - wstp —¡Tó----v e(7)
en donde:
w es uno de los coeficientes W y B, wstp y £ son valores constantes, establecidos respectivamente a 0,001 y 10-6,
m = Pi m g(1 - Pi), v = P<2>v g2(1 - £<2>), m y v se establecen inicialmente a 0, siendo P<1>establecidos respectivamente a 0,9 y 0,999, y siendo t el número de iteración,
g = gradiente (x, y), siendo x el conjunto de series temporales usado en una iteración, y es la deformación de salida verdadera correspondiente y g puede establecerse al valor medio del gradiente de la función de coste como una función de los parámetros del modelo,
Es evidente para los expertos en la técnica que pueden emplearse igualmente otros algoritmos de descenso de gradiente para ajustar los coeficientes W y B del modelo.
Una ventaja de las redes neuronales puramente convolucionales es que el tamaño de las imágenes que van a procesarse no depende del tamaño de las imágenes de entrenamiento: todas las ecuaciones anteriores (1) a (5) calculan un único píxel y se usan para calcular cada píxel de las imágenes y, por lo tanto, no dependen del número de píxeles de estas imágenes.
Según una realización, la red neuronal se entrena usando series temporales de entrenamiento sintéticas que comprenden nueve imágenes de 40 x 40 píxeles generadas a partir de deformaciones del suelo modeladas y ruido modelado.
Para probar la eficiencia de la red convolucional CNN entrenada, esta se usa para procesar series temporales de imágenes de InSAR de dos regiones que ya han sido analizadas “ manualmente” por expertos: la falla de Anatolia Septentrional en Turquía (datos de COSMO-SkyMed) y la falla de Chaman en la frontera entre Afganistán y Pakistán
(datos de Sentinel 1). Las imágenes de InSAR procedentes de Chaman comprenden series temporales de 9 imágenes de 7024 x 2488 píxeles que cubren un área de aproximadamente 180.000 km2. Tanto en Turquía como en Chaman, sigue habiendo en los datos señales más potentes que las identificadas por los expertos como fallamiento incluso después de correcciones atmosféricas convencionales. En ambos casos, ha sido necesario un conocimiento a priori de la ubicación de la falla y una inspección manual de los datos para producir series temporales de deformación de fallas. La red neuronal CNN entrenada usando datos de entrenamiento sintéticos y sin ningún ajuste fino adicional sobre datos reales, aísla y recupera automáticamente señales de deformación limpias en Turquía y Chaman, donde el análisis por expertos también halló señal. En la serie temporal de la falla de Anatolia Septentrional, la red neuronal CNN halló un deslizamiento en la línea de visión de 1,5 cm, prácticamente sin quedar señal de ruido en ningún otro lugar que en la falla, sin intervención humana y sin tener el conocimiento de la ubicación de la falla.
Debería observarse que todos los cálculos realizados cuando se ejecuta la red neuronal CNN para extraer la deformación del suelo a partir de la serie temporal de imágenes procedentes de Chaman pueden hacerse en un procesador gráfico de un ordenador personal convencional, con un tiempo de cálculo mínimo, por ejemplo, aproximadamente 4 minutos en un único procesador gráfico RTX 6000 de Nvidia® para extraer la deformación del suelo a partir de una serie temporal de 9 adquisiciones de 7024 x 2488 píxeles.
La figura 12 ilustra componentes de un ordenador personal convencional. El ordenador PC comprende al menos un procesador PRC y, acoplados operativamente al procesador, interfaces de usuario DSP, CM, memorias MEM, un procesador gráfico GP y circuitos de comunicación NIT. Las memorias MEM almacenan un sistema operativo y aplicaciones. Las interfaces de usuario DSP, CM pueden incluir, por ejemplo, pero sin limitación, un teclado numérico CM y una pantalla de visualización de ordenador DSP conectada a la tarjeta gráfica GP. Los circuitos de comunicación NIT posibilitan que el procesador PRC se acople operativamente a una red de comunicación electrónica NT, tal como Internet. La memoria MEM también incluye componentes de código de programa legibles por ordenador que implementan la red neuronal CNN. Por ejemplo, cuando estos componentes de código de programa legibles por ordenador son procesados por el procesador PRC, los componentes de código de programa se configuran para provocar la ejecución del método para procesar series temporales de imágenes, como se ha descrito anteriormente.
La descripción anterior de diversas realizaciones de la presente invención se proporciona con fines de descripción a un experto en la técnica relacionada. No pretende ser exhaustiva o limitar la invención a una única realización descrita. Numerosas alternativas y variaciones de la presente invención serán evidentes para los expertos en la técnica de la enseñanza anterior. En consecuencia, aunque se han analizado específicamente algunas realizaciones alternativas, otras realizaciones serán evidentes o podrán ser desarrolladas con relativa facilidad por los expertos en la técnica.
A este respecto, es evidente para un experto en la técnica realizar las operaciones realizadas por las etapas de procesamiento descritas en la figura 1 en otros órdenes, y el número de etapas de procesamiento puede hacerse variar. En particular, las etapas de procesamiento SL4 y SL6 podrían colocarse entre otras etapas de filtrado FPs. Las operaciones de combinación CB[s,f] realizadas por las etapas de procesamiento SL4, SL6 podrían realizarse de muchas otras formas aparte de seleccionar el valor de píxel máximo, tal como usando una función de suma o cálculo de promedio. El número de etapas de procesamiento S<l>y el número de imágenes procesadas por cada uno de los componentes de combinación CB[s,f] se establecen como una función del número de imágenes en la serie temporal de entrada RI, con la condición de que cada uno de los componentes de combinación CB[s,f] de la última etapa de procesamiento SLs en la red neuronal CNN proporcione una única imagen FI[s,f] para cada serie temporal. Además, la imagen de modelo EM podría introducirse en el procesamiento en otra etapa en la que ya no esté presente la componente temporal (el índice t), y también puede hacerse variar el número de etapas de filtrado FPs antes y después de la introducción de la imagen de modelo EM dependiendo del tipo de imágenes que procesar y, más particularmente, del tipo de señales de ruido que afectan a las imágenes que procesar.
Además, el número y las posiciones respectivas de los 35 píxeles vecinos y el número de imágenes adyacentes en la serie temporal tomadas en cuenta para calcular las combinaciones lineales también pueden hacerse variar dependiendo del tipo de imágenes que procesar y el tipo de señales de ruido que afectan a las imágenes que procesar. El ajuste de los parámetros L, M y U y los coeficientes a, b, c que selecciona los píxeles adyacentes implicados en las combinaciones lineales FT[s,f] puede realizarse como una función de la amplitud de las perturbaciones que afectan a las imágenes que van a procesarse. En general, estos parámetros y coeficientes se seleccionan entre 1 y 4.
La adición de los coeficientes de sesgo B[s,f] en las operaciones de filtrado FPs también es opcional, aunque esta introduce otros grados de libertad en el diseño de la red neuronal CNN. La eliminación de este coeficiente de sesgo podría compensarse fácilmente en términos de grado de libertad añadiendo una o más operaciones de filtrado FP, lo que introduce un número enorme de coeficientes de ponderación W y, por lo tanto, un número grande de grados de libertad.
La aplicación de la función unitaria lineal rectificada RL también es opcional y depende del tipo de imágenes que van a procesarse. En los ejemplos anteriores, y más particularmente cuando las imágenes que van a procesarse son imágenes de InSAR dañadas por turbulencia atmosférica, es deseable dar más importancia a los valores positivos que a los valores negativos de los píxeles. En otras aplicaciones, podrían adaptarse más otras funciones unitarias lineales rectificadas. Además, en lugar de usar una función de este tipo, podrían definirse diferentes coeficientes de ponderación W como una función del valor de píxel por el que se multiplica el coeficiente de ponderación.
Además, la fase de entrenamiento para calcular los coeficientes (coeficientes de ponderación W y coeficientes de sesgo B) de la red neuronal CNN no necesita realizarse antes de cada procesamiento de imagen verdadero, debido a que estos coeficientes solo dependen del tipo de imágenes y señales que procesar. En consecuencia, un ordenador que implementa una red neuronal diseñada para procesar un tipo de imágenes no necesita implementar el método de entrenamiento y el método para generar los datos de entrenamiento.
Además, el método descrito anteriormente no es aplicable exclusivamente a imágenes de InSAR del suelo adquiridas por satélites de observación terrestre. Este método puede adaptarse fácilmente para procesar series temporales de imágenes con ruido en las que se conocen modelos de ruido, para extraer señales de movimiento entre las imágenes de las series temporales.
Claims (15)
- REIVINDICACIONES i. Un método para procesar series temporales (RI[t]) de imágenes de una misma área, sometidas a ruido, comprendiendo el método: ejecutar por una red neuronal convolucional (CNN) una pluralidad de etapas de filtrado (FPs, s = 0 3, 5, 7-12) sucesivas que comprenden una primera etapa de filtrado (FP0) que recibe una serie temporal de imágenes de entrada (RI[t]), una última etapa de filtrado (FP12) que proporciona una imagen de salida (PI) y etapas de filtrado intermedias (FPs, s = 1-3, 5, 7-11), transformando, cada una de la primera etapa de filtrado y las etapas de filtrado intermedias, un conjunto (FI) de series temporales de imágenes filtradas (FI[s,t,f]) que incluyen inicialmente la serie temporal de imágenes de entrada, transformándose el conjunto generando un número (F, 1) de series temporales de imágenes filtradas combinando mediante combinaciones lineales (FT[s,t,f]) cada píxel (PX[i,j,t]) de cada imagen de cada serie temporal del conjunto con píxeles vecinos seleccionados en la imagen y en una imagen adyacente de las series temporales de imágenes del conjunto, usando cada combinación lineal un conjunto respectivo de coeficientes de ponderación y dando como resultado un píxel (PX[i,j,t,f]) de una de las imágenes filtradas generadas; y ejecutar por la red neuronal convolucional una o más operaciones de combinación (SL4, SL6), realizándose cada operación de combinación entre dos sucesivas de las etapas de filtrado intermedias, para reducir un número de imágenes en cada una de las series temporales de imágenes filtradas del conjunto, combinando imágenes de subconjuntos de imágenes adyacentes en cada una de las series temporales de imágenes filtradas del conjunto, reduciendo una última (SL6) de las operaciones de combinación cada serie temporal de imágenes filtradas del conjunto a una única imagen filtrada.
- 2. El método de la reivindicación 1, que comprende además, después de la última operación de combinación (SL6), introducir una imagen de modelo (EM) del área en el conjunto (FI) como una imagen filtrada adicional, yendo seguida la introducción de la imagen de modelo por una o más de las etapas de filtrado (FPs, s = 7 12).
- 3. El método de la reivindicación 1 o 2, en donde cada una de las combinaciones lineales (FT[s,t,f]) incluye un coeficiente de sesgo que se añade a un resultado de la combinación lineal.
- 4. El método de una de las reivindicaciones 1 a 3, en donde: cada uno de los coeficientes de ponderación tiene un valor que depende de un signo, positivo o negativo, de un valor de píxel por el que se multiplica este, o el resultado de cada una de las combinaciones lineales (FT[s,t,f]) se transforma mediante una función unitaria lineal rectificada.
- 5. El método de una de las reivindicaciones 1 a 4, que comprende además: generar la primera serie temporal de imágenes a partir de un modelo de deformación del suelo usando parámetros seleccionados aleatoriamente; generar series temporales de imágenes de entrenamiento, a partir de la primera serie temporal de imágenes, usando modelos de diferentes señales de ruido; y usar la serie temporal de entrenamiento generada para ajustar los valores de los coeficientes de ponderación.
- 6. El método de la reivindicación 5, en donde los valores de los coeficientes de ponderación se ajustan iterativamente aplicando un método de minimización de descenso basado en gradiente iterativo usando las series temporales de imágenes de entrenamiento y una función de coste de modelo (MEB), para minimizar el resultado de la función de coste de modelo.
- 7. El método de una de las reivindicaciones 1 a 6, en donde la generación del conjunto de series temporales de imágenes filtradas (FI[0,t,f]) a partir de la serie temporal de imágenes de entrada (RI[t]) se realiza usando la siguiente ecuación:en donde PX[i,j,t,f] es un píxel de una de las imágenes filtradas de la serie temporal f de imágenes filtradas, PX[i,j,t] es un píxel de una imagen de la serie temporal de imágenes de entrada, W[l,m,u,f] es uno de los primeros coeficientes para la serie temporal f de imágenes filtradas, B[f] es un coeficiente de sesgo para la serie temporal f, y LR() es una función unitaria lineal rectificada.
- 8. El método de una de las reivindicaciones 1 a 7, en donde cada combinación lineal de primeras operaciones de filtrado (FPs, s = 1-3, 5) de las operaciones de filtrado aplica la siguiente ecuación:en donde PX[i,j,t,f] es un píxel de una de las imágenes filtradas de la serie temporal de imágenes filtradas f, W[s,l,m,u,f,f] es uno de los segundos coeficientes de ponderación para la serie temporal f de imágenes filtradas para la operación de filtrado s, B[s,f] es un coeficiente de sesgo para la serie temporal f y la operación de filtrado s, y LR() es una función unitaria lineal rectificada con fugas.
- 9. El método de una de las reivindicaciones 1 a 8, en donde cada combinación lineal de segundas operaciones de filtrado (FPs, s = 7-11) de las operaciones de filtrado aplica la siguiente ecuación:en donde PX[i,j,f'] es un píxel de una de las imágenes filtradas f, W[s,l,m,f,f'] es uno de los segundos coeficientes de ponderación para las imágenes filtradas f para la operación de filtrado s, B[s,f’] es un coeficiente de sesgo para las imágenes filtradas f y la operación de filtrado s, y LR() es una función unitaria lineal rectificada con fugas.
- 10. El método de la reivindicación 9, en donde cada píxel (PX[i,j]) de la imagen de salida (PI) se calcula mediante la siguiente ecuación:en donde PX[i,j] es un píxel de la imagen de salida, PX[i,j,f] es un píxel de la imagen filtrada f, W[l,m,f] es uno de los terceros coeficientes, B es un coeficiente de sesgo y LR() es una función unitaria lineal rectificada con fugas.
- 11. El método de una de las reivindicaciones 7 a 10, en donde la función unitaria lineal rectificada con fugas LR es de tal modo que LR(x) = x si x > 0 y LR(x) = 0,5x si x < 0.
- 12. El método de una de las reivindicaciones 1 a 11, en donde cada una de las operaciones de combinación de imágenes (SL, s = 4, 6) aplica la siguiente ecuación:en donde PX[i,j,t,f] es un píxel de una de las imágenes filtradas de la serie temporal de imágenes filtradas f, y MÁX (PX[i,j,t+u,f]) es una función que proporciona el valor máximo de entre los valores de píxel PX[i,j,t+u,f] con u = -1, 0 y 1.
- 13. Un ordenador que comprende un procesador (PRC), una memoria (MEM), estando configurado el procesador para llevar a cabo el método de una de las reivindicaciones 1 a 12.
- 14. El ordenador de la reivindicación 13, que comprende además un procesador gráfico (GP) controlado por el procesador (PRC), estando configurado el procesador para configurar el procesador gráfico para llevar a cabo algunas de las operaciones del método.
- 15. Un producto de programa informático cargable en una memoria informática y que comprende porciones de código que, cuando son llevadas a cabo por un ordenador, configuran el ordenador para llevar a cabo el método de una de las reivindicaciones 1 a 12.
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| EP20157709.5A EP3866105B1 (en) | 2020-02-17 | 2020-02-17 | Method for processing insar images to extract ground deformation signals |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| ES3013922T3 true ES3013922T3 (en) | 2025-04-15 |
Family
ID=69631475
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| ES20157709T Active ES3013922T3 (en) | 2020-02-17 | 2020-02-17 | Method for processing insar images to extract ground deformation signals |
Country Status (5)
| Country | Link |
|---|---|
| US (1) | US20230072382A1 (es) |
| EP (1) | EP3866105B1 (es) |
| JP (1) | JP2023526713A (es) |
| ES (1) | ES3013922T3 (es) |
| WO (1) | WO2021165197A1 (es) |
Families Citing this family (13)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| SE544261C2 (en) | 2020-06-16 | 2022-03-15 | IntuiCell AB | A computer-implemented or hardware-implemented method of entity identification, a computer program product and an apparatus for entity identification |
| WO2023214915A1 (en) * | 2022-05-06 | 2023-11-09 | IntuiCell AB | A data processing system for processing pixel data to be indicative of contrast. |
| CN114594479B (zh) * | 2022-05-07 | 2022-07-26 | 中国测绘科学研究院 | 一种全散射体FS-InSAR方法及系统 |
| CN114966692B (zh) * | 2022-07-19 | 2022-11-08 | 之江实验室 | 基于Transformer的InSAR技术冻土区多变量时序形变预测方法及装置 |
| CN115453534B (zh) * | 2022-09-19 | 2023-05-16 | 中山大学 | 一种顾及解缠误差的序贯InSAR时序形变解算方法 |
| CN116299245B (zh) * | 2023-05-11 | 2023-07-28 | 中山大学 | 一种时序InSAR形变速率结果自适应镶嵌校正方法 |
| CN116659429A (zh) * | 2023-08-01 | 2023-08-29 | 齐鲁空天信息研究院 | 一种多源数据高精度时序地表三维形变解算方法和系统 |
| CN117148340B (zh) * | 2023-08-08 | 2024-10-18 | 长安大学 | 一种基于多源对地观测的活动性地质灾害探测方法与系统 |
| CN117213443B (zh) * | 2023-11-07 | 2024-03-19 | 江苏省地质调查研究院 | 一种天地深一体化地面沉降监测网建设与更新方法 |
| CN117671531B (zh) * | 2023-12-05 | 2024-08-02 | 吉林省鑫科测绘有限公司 | 一种无人机航测数据处理方法及系统 |
| CN118037563B (zh) * | 2024-01-19 | 2024-11-22 | 中国民用航空飞行学院 | 一种机场跑道沉降预测方法及系统 |
| CN118482685B (zh) * | 2024-05-11 | 2024-12-27 | 广东工业大学 | 融合时间序列与时空形变的地下开采沉降监测方法及系统 |
| CN118584488B (zh) * | 2024-08-02 | 2024-11-12 | 中国地震局第一监测中心 | 一种大气延迟误差修正的方法、装置、存储介质及电子设备 |
Family Cites Families (11)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US4975704A (en) * | 1990-01-26 | 1990-12-04 | The United States Of America As Represented By The Administrator Of The National Aeronautics And Space Administration | Method for detecting surface motions and mapping small terrestrial or planetary surface deformations with synthetic aperture radar |
| US7450054B2 (en) * | 2007-03-22 | 2008-11-11 | Harris Corporation | Method and apparatus for processing complex interferometric SAR data |
| EP2304466B1 (en) * | 2008-07-04 | 2015-10-28 | Telespazio S.p.A. | Identification and analysis of persistent scatterers in series of sar images |
| IT1394733B1 (it) * | 2009-07-08 | 2012-07-13 | Milano Politecnico | Procedimento per il filtraggio di interferogrammi generati da immagini sar acquisite sulla stessa area. |
| US8942454B2 (en) * | 2009-12-03 | 2015-01-27 | The United States Of America, As Represented By The Secretary, Department Of Health And Human Services | Signal to-noise enhancement in imaging applications using a time-series of images |
| CN104459633B (zh) * | 2014-12-01 | 2016-08-17 | 中国科学院电子学研究所 | 结合局部频率估计的小波域InSAR干涉相位滤波方法 |
| US10495750B1 (en) * | 2015-11-19 | 2019-12-03 | National Technology & Engineering Solutions Of Sandia, Llc | Spectral replacement to mitigate interference for multi-pass synthetic aperture radar |
| GB201709525D0 (en) * | 2017-06-15 | 2017-08-02 | Univ Nottingham | Land deformation measurement |
| US10430913B2 (en) * | 2017-06-30 | 2019-10-01 | Intel Corporation | Approximating image processing functions using convolutional neural networks |
| CN109541596B (zh) * | 2018-11-28 | 2022-05-20 | 中国电子科技集团公司电子科学研究院 | 基于深度学习算法的InSAR图像处理方法及装置 |
| WO2020240720A1 (ja) * | 2019-05-29 | 2020-12-03 | 日本電気株式会社 | 合成開口レーダの信号処理装置および信号処理方法 |
-
2020
- 2020-02-17 EP EP20157709.5A patent/EP3866105B1/en active Active
- 2020-02-17 ES ES20157709T patent/ES3013922T3/es active Active
-
2021
- 2021-02-15 JP JP2022549434A patent/JP2023526713A/ja active Pending
- 2021-02-15 WO PCT/EP2021/053628 patent/WO2021165197A1/en not_active Ceased
- 2021-02-15 US US17/904,454 patent/US20230072382A1/en active Pending
Also Published As
| Publication number | Publication date |
|---|---|
| EP3866105B1 (en) | 2024-12-25 |
| US20230072382A1 (en) | 2023-03-09 |
| EP3866105C0 (en) | 2024-12-25 |
| EP3866105A1 (en) | 2021-08-18 |
| WO2021165197A1 (en) | 2021-08-26 |
| JP2023526713A (ja) | 2023-06-23 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| ES3013922T3 (en) | Method for processing insar images to extract ground deformation signals | |
| Samsonov et al. | Application of DInSAR-GPS optimization for derivation of fine-scale surface motion maps of Southern California | |
| Sousa et al. | Persistent Scatterer InSAR: A comparison of methodologies based on a model of temporal deformation vs. spatial correlation selection criteria | |
| Nesterov et al. | GNSS radio tomography of the ionosphere: The problem with essentially incomplete data | |
| Karnesis et al. | Assessing the detectability of a Stochastic Gravitational Wave Background with LISA, using an excess of power approach | |
| Shirzaei et al. | Estimating the effect of satellite orbital error using wavelet-based robust regression applied to InSAR deformation data | |
| Angling et al. | Assimilation of radio occultation measurements into background ionospheric models | |
| Okura et al. | A method for weak-lensing flexion analysis by the HOLICs moment approach | |
| Camphuis et al. | SPT-3G D1: CMB temperature and polarization power spectra and cosmology from 2019 and 2020 observations of the SPT-3G Main field | |
| Chen et al. | ARU-net: Reduction of atmospheric phase screen in SAR interferometry using attention-based deep residual U-net | |
| Benoist et al. | Accounting for spatiotemporal correlations of GNSS coordinate time series to estimate station velocities | |
| Kitching et al. | Gravitational lensing accuracy testing 2010 (GREAT10) challenge handbook | |
| Falabella et al. | On the Phase Nonclosure of Multilook SAR Interferogram Triplets. | |
| Girard et al. | Sparse representations and convex optimization as tools for LOFAR radio interferometric imaging | |
| Zawadzki et al. | Regularized maximum likelihood image synthesis and validation for ALMA continuum observations of protoplanetary disks | |
| Hu et al. | Isolating orbital error from multitemporal InSAR derived tectonic deformation based on wavelet and independent component analysis | |
| Ayasso et al. | A variational Bayesian approach for unsupervised super-resolution using mixture models of point and smooth sources applied to astrophysical map-making | |
| Hobiger et al. | On the importance of accurately ray-traced troposphere corrections for Interferometric SAR data | |
| Sadeghi et al. | Monitoring land subsidence in a rural area using a combination of ADInSAR and polarimetric coherence optimization | |
| CN110286374B (zh) | 基于分形布朗运动的干涉sar影像仿真方法 | |
| Gallagher et al. | Two-dimensional drift velocities from the IMAGE EUV plasmaspheric imager | |
| Arminjon | Equations of motion for the mass centers in a scalar theory of gravitation | |
| Cornely | Flexible prior models: three-dimensional ionospheric tomography | |
| Kravchenko et al. | Orthorectification of Sich-2 satellite images using elastic models | |
| Kouba | Comparison of polar motion with oceanic and atmospheric angular momentum time series for 2-day to Chandler periods |