Evaluation of Sea Ice Drift in the Arctic Marginal Ice Zone based on Sentinel-1A/B Satellite Radar Measurements

全文:

Введение

Мониторинг дрейфа морского льда является важной составляющей исследований геофизических процессов в Арктике. Получаемые при этом данные можно использовать при решении широкого круга теоретических и практических задач. На сегодняшний день одним из наиболее перспективных подходов для решения этой задачи является анализ спутниковых снимков в видимом, инфракрасном и микроволновом диапазонах. Последовательности спутниковых изображений, относящиеся к одному и тому же фрагменту морской поверхности и узкому временному интервалу, могут быть использованы для восстановления динамики морской поверхности, в частности, скорости и направления дрейфа морского льда. При этом для решения описанной задачи часто применяются методы расчета оптического потока. По данному подходу существует обширная литература, в которой описаны как сами алгоритмы, так и особенности их использования на практике в зависимости от специфики исходных данных ^{1, 2} [1–7]. В силу объема имеющихся на сегодняшний день архивов спутниковых снимков актуальной является проблема автоматизации расчетов. Успешное решение этой задачи открывает возможность получения регулярных и продолжительных наблюдений за динамикой дрейфа морского льда.

Данное направление получило широкое развитие в последнее десятилетие, результаты исследований изложены, напр., в работах [8–12]. При этом используются различные алгоритмы трекинга применительно к данным MODIS, Sentinel-1, Sentinel-2, Landsat-8 и ряду других преимущественно в оптическом диапазоне. Важными факторами являются разрешение исходных снимков, возможность эффективной оптимизации вычислений с целью обработки больших массивов данных, автоматизация расчетов. В настоящей работе описан автоматизированный подход, позволяющий оценивать динамику морского льда по радиолокакционным (РЛ) изображениям, покрывающим площадь порядка нескольких сотен тысяч квадратных километров, с разрешением 40 м. По информации авторов, практические результаты по использованию автоматической системы обработки данных с такими параметрами публикуются впервые.

После запуска спутниковых радиолокаторов с синтезированной апертурой (РСА) Sentinel-1A/B появилась возможность получать пары РЛ-изображений (РЛИ) арктического региона с небольшой временной задержкой, что открыло возможность оценивать кинематические характеристики различных мезо- и субмезомасштабных процессов верхнего слоя арктических морей [13–17]. В приложении к оценке полей скорости дрейфа льда авторами работы [14] была показана возможность эффективно рассчитывать такие поля в областях интенсивной вихревой динамики в прикромочной ледовой зоне (ПЛЗ) с пространственным разрешением Ο(100 м) на основе метода максимальной кросс-корреляции. Дальнейшее развитие это направление получило в работах [18, 19], где авторы применили методы машинного обучения для идентификации вихрей ПЛЗ в данных спутниковых РСА.

Целью настоящей работы является описание общей методики и демонстрация эффективности расчетной схемы автоматической обработки спутниковых данных Sentinel–1A/B с целью восстановления полей дрейфа морского льда в ПЛЗ. Предложенная методика была применена к трехмесячному массиву данных, относящихся к району пролива Фрама. Использованный подход позволил в автоматическом режиме обработать статистически значимый объем данных и получить качественные результаты, которые подтверждают эффективность используемого подхода для оперативного мониторинга динамики дрейфа ледяных полей в прикромочной зоне Арктики.

Данные и методы

В качестве исходных данных использованы РЛИ Sentinel–1A/B, полученные для района пролива Фрама в 2017 и 2018 гг. В частности, использовались продукты L1 Ground Range Detected с размером пикселя 40 × 40 м и пространственным разрешением 93 × 87 м в режиме сьемки Extra Wide Swath (сверх-широкозахватный). Спутниковые РЛИ были получены с сайта https://scihub.copernicus.eu из архивов системы Европейского центра морских прогнозов Copernicus Open Access Hub.

Всего было отобрано 63 пары изображений, полученных в сентябре 2017 г., а также в июле и сентябре 2018 г. Для каждой пары РЛИ временной диапазон между последовательными измерениями составил от 48 до 52 мин. Радиолокационные изображения выбирались исходя из соображений информативности, определяемой выраженностью поверхностных структур в поле дрейфующего льда. Данные интерполировались на регулярную сетку с разрешением 40 × 40 м, после чего выделялись пересечения сканируемых областей на каждом РЛИ. Дополнительно изображения сглаживались при помощи медианного фильтра со скользящим окном размером 20 × 20 пикселей. На рис. 1 показана пара РЛИ Sentinel-1A/B за 26 июля 2018 г., полученных таким образом.

 

Р и с. 1. Пара РЛ-изображений Sentinel-1A (a) и Sentinel-1B (b) за 15 сентября 2018 г. с наличием дрейфующего льда в прикромочной ледовой зоне пролива Фрама

F i g. 1. A pair of Sentinel-1A (a) and Sentinel-1B (b) radar images for September 15, 2018 with the drifting ice present in the marginal ice zone of Fram Strait

 

Основным алгоритмом, использовавшимся при обработке данных, является алгоритм нормализированной максимальной кросс-корреляции (МКК) ¹ [1, 3, 5–7]. Суть метода заключается в нахождении максимально коррелирующих фрагментов на анализируемых изображениях, смещение которых друг от друга не превышает теоретически возможного. Для каждого пикселя на обоих снимках выделяются окрестности, для которых строится нормализированная кросс-корреляционная функция. Отклонение положения максимума этой функции от центральной точки принимается за искомый сдвиг, т. е. смещение фрагмента текстуры от снимка к снимку. Формально подход можно описать следующим образом: обозначим матрицы изображений как $I_1$ и $I_2$. Рассмотрим некоторый пиксель с индексами i и j. Пусть окрестность имеет квадратную форму и размер N × N, где N – нечетно;

$k=\left\{i-\left(N-1\right)/\frac{2}{i}+\left(N-1\right)/\mathrm{2,}l=\left\{j-\left(N-1\right)/\frac{2}{j}+\left(N-1\right)/2\right\}\right\},$

$M_1=\frac{1}{N^2}\sum _{k\mathrm{,1}}{I}_1\left(k,l\right)иM_2\left(u,v\right)=\frac{1}{N^2}\sum _{k\mathrm{,1}}{I}_2\left(k+u,l+v\right).$

Тогда нормализированная кросс-корреляционная функция будет иметь вид

$F\left(u,v,i,j\right)=\frac{{\displaystyle \sum _{i,k}\left(I_1\left(k,l\right)-M_1\right)\left(I_2\left(k+u,l+v\right)-M_2\left(u,v\right)\right)}}{\sqrt{\left({\displaystyle \sum _{l,k}{\left(I_1\left(k,l\right)-M_1\right)}^2}\right)\left({{\displaystyle \sum _{l,k}\left(I_2\left(k+u,l+v\right)-M_2\right)}}^2\right)}}.$

Здесь $-u_{\max }\le u\le u_{\max },-v_{\max }\le v\le v_{\max },$ где $u_{\max }\text{\hspace{0.33em}и\hspace{0.33em}}{v}_{\max }$ – заданные максимальные возможные сдвиги по обеим координатам. Искомые значения $u_0\text{\hspace{0.33em}и\hspace{0.33em}}{v}_0$ для пикселя c индексами i и j таковы, что $F\left(\left(u_0,v_0,i,j\right)=\max \left(F\left(u,v,i,j\right)\right)\right),$ где u и v принимают значения из указанного выше промежутка. Найдя описанные величины для всех $\left(N-1\right)/2<i<S_1-\left(N-1\right)/2+1$ и $\left(N-1\right)/2<j<-\left(N-1\right)/2+1$, где $\left(S_1,S_2\right)$ – размер изображений, получим матрицы U и V размером $\left(S_1-N-u_{\max }+l,S_2-N-v_{\max }+l\right)$.

На рис. 2 продемонстрированы обе компоненты $\left(u_0\text{\hspace{0.33em}и\hspace{0.33em}}{v}_0\right)$ векторного поля, построенного таким образом с параметрами N = 80 $u_{\max }=v_{\max }=50$ по данным, показанным на рис. 1.

 

Р и с . 2 . Пиксельные сдвиги, соответствующие компонентам полей скорости 0 u (a), 0 v (b) и полученные 26 июля 2018 г. на основе алгоритма МКК с использованием РЛ-изображений Sentinel-1A/B

F i g. 2. Pixel shifts corresponding to the velocity field components 0 u (a) and 0 v (b), and obtained using the MCC algorithm for the Sentinel-1A/B radar images for July 26, 2018

 

Здесь обращает на себя внимание соответствие фрагментов, содержащих полезные данные, фрагментам, где явно просматривается наличие ледяных скоплений на исходных изображениях. Остальная часть сцены либо не информативна, либо содержит шум, являющийся проявлением эффекта, который называют «ложными корреляциями». Основным общепринятым подходом для фильтрации соответствующих фрагментов является ограничение снизу на минимальное значение кросс-корреляционной функции и отсечение пикселей, не удовлетворяющих данному условию. Кроме этого, очевидно, эффективным подходом здесь мог бы стать морфологический анализ полученного поля. Однако авторы решили пойти другим путем, приняв во внимание то, что оператор $\left(I_1,I_2\right)\to \left(U,V\right)$ не является симметричным.

Рассмотрим некоторый пиксель с индексами (i, j). Обозначим индексы его окрестности как $k=\left\{i-\left(N-1\right)/\frac{2}{i}+\left(N-1\right)/2\right\},$ $l=\left\{j-\left(N-1\right)/\frac{2}{j}+\left(N-1\right)/2\right\}$ . Пусть $\left(u_0,v_0\right)$ – рассчитанный сдвиг и для любых смещений фрагмент $I_1\left(k,l\right)$ максимально коррелирует с фрагментом $I_2\left(k+u_0,l+v_0\right)$. Тогда корреляция $I_1\left(k-u_0,l+v_0\right)$ и $I_1\left(k-u_0,l+v_0\right)$ также должна быть высокой. Следовательно, если провести описанный расчет, взяв $I_2$ в качестве первого изображения, а $I_1$ – в качестве второго, результаты должны быть близкими при изменении знаков у обеих компонент. Назовем такой расчет обратным. Для оценки близости можно использовать косинус угла между векторами двух таких полей. На рис. 3, а показаны значения этой величины. Здесь явно просматривается соответствие фрагментов с высокими значениями косинуса фрагментам с информативными данными, что подтверждает эффективность описанного подхода. Установив оптимальное пороговое значение, можно выделить маску, для которой рассчитанные поля скорости являются релевантными. Для расчетов использовался порог, равный 0,92. Кроме этого, для удаления остаточных шумов можно применить градиентный фильтр, выделяющий области высокой пространственной неоднородности. Для каждого пикселя выберем максимальное значение вариации компонент u и v по окрестности размером 3 × 3 пикселя (рис. 3, b). Все фрагменты с превышением порогового значения 2,51 отсекаются.

 

Р и с. 3. Пространственное распределение косинуса угла между векторами рассчитанных сдвигов, полученных при прямом и обратном расчетах (a); оценка однородности компонент сдвигов по окну размером 3 х 3 пикселя (b)

F i g. 3. Spatial distribution of the cosine of angle between the vectors of calculated shifts resulted from direct and inverse calculations (a); assessment of the homogeneity of shift components over the 3 х 3 pixel window (b)

 

После осуществления описанной процедуры фильтрации в получаемой маске остается ряд небольших по площади отдельных фрагментов, информативность которых представляется сомнительной. Для их удаления маска разбивается на связные компоненты, для каждой из которых рассчитывается число содержимых пикселей. Если это число меньше заданного порога, компонента вырезается. В качестве оценочного значения выбрано отношение площади компоненты к площади всей сцены, в качестве порога – значение 0,25%.

Суммируя все вышеизложенное, можно кратко описать алгоритм обработки пары РЛ-изображений в виде следующих шагов:

– распаковка и репроекция исходных данных на регулярную сетку;

– прямой и обратный расчеты по методике нормализированной максимальной кросс-корреляции;

– фильтрация нерелевантных фрагментов;

– сжатие сцены с учетом границ полученной на предыдущем шаге маски;

– удаление мелких отдельных компонент маски;

– преобразование сдвигов в поле скорости (м/с).

Результат всего процесса обработки данных показан на рис. 4. Видно, что диапазон наблюдаемых значений скорости дрейфа льда составляет 0,1–0,8 м/с, среднее значение скорости ~ 0,2–0,3 м/с, доминирующее направление дрейфа льда – южное, юго-восточное. Максимальные значения скорости дрейфа наблюдаются в северо-западном секторе, а также на юго-восточной периферии показанной области.

 

Р и с. 4. Конечный результат расчета поля скорости (м/с) дрейфа льда в проливе Фрама за 26 июля 2018 г.

F i g. 4. Final result of calculation of the ice drift velocity (m/s) field in Fram Strait for July 26, 2018

 

Результаты

По итогам расчетов получено 63 векторных поля скорости дрейфа льда, описывающих горизонтальное перемещение скоплений морского льда в ПЛЗ пролива Фрама. Для пары изображений размером 10000 × 10000 пикселей время работы использованной реализации алгоритма МКК с размером скользящего окна 80 × 80 и максимальными сдвигами $u_{\max }=\mathrm{50,}{v}_{\max }=50$ при расчете в два потока на машине с процессором Intel i5 12500H под управлением операционной системы Debian GNU/Linux 12 составляет около двух часов.

Объективная оценка эффективности описанного метода является затруднительной в силу отсутствия эталонных данных о поверхностной морской динамике. Наиболее перспективным подходом здесь представляется сопоставление восстановленных по спутниковым данным полей скорости с натурными данными о направлении и скорости дрейфа льда, полученными на основе установленных на льдинах стационарных геолокационных датчиков (GPS-треке-ров). На рис. 5 показана траектория дрейфа такого датчика на фоне рассчитанного по спутниковым данным поля скорости.

 

Р и с . 5 . Сопоставление векторов скорости дрейфа льда по спутниковым изображениям с натурными измерениями, полученными по данным установленных на лед GPS-трекеров. Синей кривой показана траектория движения льдин по натурным данным. Красной стрелкой выделен фрагмент натурных измерений, совпадающий по времени со спутниковой съемкой. Моменты съемки Sentinel-1: 26.04.2019 в 05:02 и в 05:51. Траектория датчика – с 26.04.2019 (00:10) по 27.04.2019 (11:30)

F i g. 5. Comparison of the ice drift velocity vectors derived from satellite data with the field measurements obtained from the data of GPS trackers installed on ice. Blue curve shows the trajectory of ice floes based on the field data. Red arrow highlights a fragment of field measurements coinciding in time with satellite imagery. Time of Sentinel-1 acquisitions: at 05:02 UTC and at 5:51 UTC on 26.04.2019. GPS tracker is traced from 26.04.2019 (00:10) to 27.04.2019 (11:30)

 

К сожалению, авторы не обладают статистически значимым массивом натурных данных для проведения более детального сопоставления. Тем не менее из рис. 5 видно, что измерения по двум независимым источникам хорошо совпадают по направлению и по абсолютному значению. Кроме этого, определенное представление о релевантности подхода можно получить на основе экспертной оценки структуры восстанавливаемых полей и соответствия видимых смещений контрастных структур в текстурах изображений рассчитанным сдвигам. С этой точки зрения результаты расчетов демонстрируют высокое качество. Фрагменты с ярко выделяющимися рассогласованностями эффективно отсекаются с помощью процедуры автоматической фильтрации, при этом количество выброшенных пикселей пренебрежимо мало по сравнению с их общим числом в зонах с хорошо выраженными контрастами. Существенной избыточности или недостаточности этого алгоритма замечено не было. Мезомасштабные и субмезомасштабные динамические структуры (такие как отдельные вихри, вихревые диполи или ледовые филаменты), различимые на снимках, хорошо просматриваются также и в получаемых полях скорости.

Визуализация спутниковых продуктов с полями дрейфа льда

На следующем этапе продукты с восстановленными по квазисихронным спутниковым данным полями скорости дрейфа морского льда размещаются на веб-странице лаборатории морских полярных исследований (ЛМПИ) Морского гидрофизического института РАН по адресу http://polar-space.ru/arctic_currents. Для этого рассчитанные поля сначала интерполируются на сетку с шагом 40 м и записываются в виде netCDF-файлов. Для размещения в сети Интернет была разработана система визуализации данных на базе клиент-серверной архитектуры. На рис. 6 показана общая структура системы визуализации данного спутникового продукта.

 

Р и с. 6. Общая структура системы визуализации продуктов с полями дрейфа льда, представленная на сайте http://polar-space.ru/arctic_currents

F i g. 6. General structure of the system for visualizing the products with ice drift fields represented on the site http://polar-space.ru/arctic_currents

 

В рамках этой системы данные с полями дрейфа льда в формате netCDF хранятся в виде файлового архива на сервере. Для отображения скоростей дрейфа льда на географической карте был разработан специальный программный модуль на языке Python, преобразующий исходные данные в изображения с единой шкалой скорости для всех файлов. Название каждого файла рисунка содержит дату, что позволяет установить однозначное соответствие между запросом пользователя и отображаемым полем дрейфа льда.

Пользовательский интерфейс (ПИ) для выбора и отображения интересующего продукта был разработан с использованием javascript библиотек jQuery. Функции картографического сервиса реализованы с помощью библиотеки mapBox GL; ПИ позволяет выбирать поле дрейфа льда за интересующую дату, пролистывая массив дат или выбирая нужную дату из списка. На рис. 7 показан пример ПИ на сайте http://polar-space.ru/arctic_currents за выбранную дату 26 июля 2018 г.

 

Р и с. 7. Пользовательский интерфейс системы визуализации продуктов с полями дрейфа льда в прикромочной ледовой зоне пролива Фрама за 26 июля 2018 г. Слева – пример статического отображения полей дрейфа льда, справа – динамическое поле с использованием библиотеки Wind-JS

F i g. 7. User interface of the system for visualizing the products with ice drift fields in the marginal ice zone of Fram Strait for July 26, 2018. On the left is an example of a static display of ice drift fields, on the right – a dynamic field using the Wind-JS library

 

В рамках разработанного ПИ существует выбор отображения данных с полями дрейфа льда в виде не только статической картинки (рис. 7, слева), но и динамического поля (рис. 7, справа). Для реализации этой задачи используется библиотека Wind-JS с сайта https://github.com/Esri/wind-js, которая была адаптирована для отображения векторов дрейфа льда. Входными данными являются файлы в формате json, содержащие значения зональной и меридиональной компонент дрейфа льда. Для генерации этих файлов в автоматическом режиме был разработан специальный программный модуль на языке Python. К настоящему времени на веб-сайте ЛМПИ МГИ размещены все обработанные за 2017 и 2018 гг. данные. Кроме того, существующий интерфейс дополнен возможностью скачивания интересующих полей дрейфа льда в формате netCDF.

Заключение

В работе дано описание общей методологии обработки квазисихронных измерений спутниковых РСА Sentinel-1 A/B, работавших в тандеме в течение 5 лет с 2016 по 2021 гг. с целью получения регулярных по времени полей скорости дрейфа льда в арктической прикромочной зоне. Временное разрешение получаемых продуктов с полями скорости дрейфа льда составляет, как правило, одни сутки в европейском секторе Арктики, и от 1 до 5 сут в других арктических районах.

Ограниченная по набору исходных данных валидация получаемых полей скорости с натурными измерениями показала хорошее соответствие между ними. Тем не менее работа в этом направлении будет продолжена для определения объективных статистик качества рассчитываемых полей скорости дрейфа льда.

На фоне существующих продуктов по дрейфу льда в Арктике представленный подход и восстанавливаемые поля скорости отличаются значительно более высоким пространственным разрешением (40 м), относительной простотой используемых методов и нетребовательностью к вычислительным ресурсам.

Полученные продукты с полями скорости дрейфа льда в прикромочной зоне Арктики доступны на сайте http://polar-space.ru/arctic_currents и могут быть использованы специалистами для решения ряда практических задач, в том числе для оценки динамики и эволюции ПЛЗ и протекающих в ней процессов, а также для сравнения с иными спутниковыми продуктами и валидации моделей общей циркуляции океана более низкого разрешения.

Спутниковые данные Sentinel-1A/B были получены из архивов системы Европейского центра морских прогнозов Copernicus Open Access Hub с сайта URL: https://scihub.copernicus.eu.

 

¹ Barron J. L., Thacker N. A. Tutorial: Computing 2D and 3D Optical Flow. Manchester, 2005. 12 p. (Tina Memo No. 2004-012).

² Stark M. Optical Flow PIV: Improving the Accuracy and Applicability of Particle Image Velocimetry : Master Thesis. Zürich : ETH, Department of Mechanical and Process Engineering (D-MAVT), 2013. 57 p. https://doi.org/10.3929/ethz-a-009767070

作者简介

E. Plotnikov

Marine Hydrophysical Institute, Russian Academy of Sciences

编辑信件的主要联系方式.
Email: ev.plotnikov@ya.ru
ORCID iD: 0000-0003-4365-9369
SPIN 代码: 6601-1958

младший научный сотрудник, лаборатория морских полярных исследований, отдел дистанционных методов исследования

俄罗斯联邦, Sevastopol

I. Kozlov

Marine Hydrophysical Institute, Russian Academy of Sciences

Email: ik@mhi-ras.ru
ORCID iD: 0000-0001-6378-8956
SPIN 代码: 3267-6693

ведущий научный сотрудник, заведующий лабораторией морских полярных исследований, кандидат физико-математических наук

俄罗斯联邦, Sevastopol

E. Zhuk

Marine Hydrophysical Institute, Russian Academy of Sciences

Email: alenixx@gmail.com
ORCID iD: 0000-0002-4263-7734
SPIN 代码: 3814-6300

младший научный сотрудник, отдел океанографии

俄罗斯联邦, Sevastopol

A. Marchenko

University Centre in Svalbard

Email: aleksey.marchenko@unis.no
ORCID iD: 0000-0003-4169-0063
SPIN 代码: 9213-8763

доктор физико-математических наук, профессор, департамент арктических технологий

挪威, Longyearbyen

参考

补充文件