Определение площадей озер по данным со спутников серии Sentinel-2

В последние десятилетия под влиянием глобальных и региональных изменений в климатической системе и воздействием антропогенных факторов на территории Российской Федерации наблюдаются значительные изменения состояния и гидрологического режима водных объектов [1]. Известно [2, 3], что фундаментальное значение для понимания и оценки степени воздействия климатических изменений и человеческой деятельности на наземные водные ресурсы имеют морфометрические характеристики озер, такие как площадь водного зеркала, уровень и объем воды.

В последние годы для изучения площадей озер широко используются архивные и оперативные данные, получаемые со спутников серии Landsat с периодичностью 16 дней [3-5]. Однако пространственное разрешение мультиспектральных изображений, получаемых со спутников Landsat, составляет 30 м.

С 2016 г. на регулярной основе с периодичностью съемки в 10 дней стали доступны данные со спутника Sentinel-2А. Пространственное разрешение изображений Sentinel-2 в видимом и ближнем инфракрасном диапазонах спектра составляет 10 м, что значительно лучше разрешения изображений Landsat. Второй спутник с аналогичными характеристиками Sentinel-2B был запущен 7 марта 2017 г., что позволило улучшить периодичность съемки до 5 дней. Более высокое разрешение изображений дает возможность использовать их для мониторинга небольших по площади озер.

Точность определения площадей озер зависит не только от пространственного разрешения изображений, но и от используемых методов выделения водной поверхности. В настоящее время при выделении водной поверхности пиксели изображений, как правило, жестко разделяются на водные и неводные. При этом наличие смешанных граничных пикселей (mixels), которые захватывают как воду, так и сушу, может приводить к существенным ошибкам при определении площадей водоемов, а также к расхождениям между результатами различных методов выделения водной поверхности [6].

В данной работе проводится сравнительный анализ известных методов выделения водной поверхности и предлагается новый метод автоматического определения площади озер по спутниковым данным серии Sentinel-2 на основе субпиксельного анализа смешанных пикселей.

Экспериментальное сравнение методов выделения водных объектов

В настоящее время известен целый ряд методов выделения водных объектов по данным оптической спутниковой съемки [7, 8], которые условно можно разбить на несколько групп: 1) пороговые методы, основанные на применении определенного спектрального канала или комплексного спектрального признака (индекса); 2) методы, использующие пороговое дерево решений на основе заданных спектральных характеристик; 3) методы на основе классификации с обучением; 4) методы, использующие алгоритмы кластеризации, и др.

Самым простым способом выделения водных объектов является пороговое разделение пикселей по значениям яркости в заданном канале. Обычно используется спектральный канал из ближнего или среднего инфракрасного диапазона, потому что в этих диапазонах наблюда-– 2 – ется сильное поглощение излучения водой в отличие от растительности и сухой почвы. Однако этот метод подвержен влиянию шума и допускает ошибки, относя теневые области и часть антропогенных объектов к воде [9].

Помимо исходных спектральных каналов широко распространено применение водных индексов. Эти индексы усиливают контраст между водными поверхностями и другими объектами. Наиболее распространенными [3-5] являются нормализованный разностный водный индекс NDWI [10] и его модификация MNDWI [11]. Предполагается, что водным объектам соответствуют положительные значения этих индексов, а остальным – отрицательные. NDWI использует видимый зеленый (G) и ближний инфракрасный (NIR) каналы:

NDWI =

G - NIR

G + NIR

В свою очередь MNDWI отличается использованием канала SWIR вместо канала NIR. Однако необходимо отметить, что для снимков Sentinel-2 канал SWIR (№ 11) имеет разрешение 20 м, что снижает пространственную точность результатов по сравнению с использованием NDWI.

Многочисленные исследования показывают [10, 12], что зачастую оптимальный водный порог для индексов NDWI, MNDWI может быть отличен от нуля, что может ограничивать их применение в автоматическом режиме. Кроме того, методы на основе водных индексов допускают ошибки на теневых областях [9].

Третий водный индекс строится как отношение значений в синем канале (B) и канале SWIR: B/SWIR. Предполагается, что для водных объектов значения этого индекса больше 1. Данный индекс был выбран, так как он необходим для выделения воды при формировании продукта Sentinel-2 уровня обработки 2А [13].

Построение дерева решений, использующего заданные пороги по нескольким спектральным каналам и индексам (вегетационному, водному и др.), позволяет сократить количество ошибочно отнесенных к воде объектов. Например, с помощью нормализованного вегетационного индекса NDVI можно отделить «зеленую воду» (с повышенным значением NDVI) от сильно увлажненной почвы [7]. Как правило, дерево решений строится отдельно для каждой спутниковой платформы, а для высокого разрешения не применяется.

Для методов классификации с обучением необходимо задание представительных обучающих выборок, что требует привлечения экспертных знаний и значительных временных затрат. Методы кластеризации (например, широко распространенный алгоритм ISODATA) дает возможность автоматически разделить изображение на некоторое количество кластеров, из которых эксперт выбирает относящиеся к водным объектам [14]. Это обстоятельство, в частности, не позволяет использовать методы классификации и кластеризации для обработки спутниковых данных в автоматическом режиме [15].

В работе [8] предложен двухэтапный метод ECCA-Water для автоматического выделения водных объектов и сильно увлажненных почв. Он основан на применении непараметрического ансамблевого алгоритма кластеризации ECCA и дерева решений, построенного с применением средних значений спектральных характеристик кластеров. Алгоритм ECCA использует информацию в трех спектральных каналах: зеленом, красном и ближнем инфракрасном. В дереве – 3 – решений задействован также водный индекс NDWI и вегетационный индекс NDVI. Этот метод позволяет выделять водные объекты на изображениях не только низкого и среднего, но и высокого (2-10 м) пространственного разрешения.

Далее приводятся результаты, полученные в ходе применения автоматических методов определения площадей водного зеркала для Большого Ярового и Красиловского озер, расположенных на территории Алтайского края, по данным со спутника Sentinel-2.

Красиловское озеро (53°11’06”N, 84°21’35”E) представляет собой непроточный водоем со сложной конфигурацией акватории (рис. 1 a ). Озеро Большое Яровое (52°52’00”N, 78°36’50”E) – это бессточный водоем с соленой водой, имеющий овальную форму (рис. 2 a ). В табл. 1, 2 приведены значения площадей Красиловского и Большого Ярового озер за 2016 г., вычисленные по спутниковым данным Sentinel-2, с помощью метода ECCA-Water и пороговых методов на основе водных индексов NDWI, MNDWI и B/SWIR, а также их стандартные среднеквадратические отклонения (СКО).

a                               b                               c

Рис. 1. RGB-композит снимка Sentinel-2 оз. Красиловское за 26.08.2016 (каналы: Red, NIR, Blue) (a); результаты выделения акватории озера по снимку за 15.09.2016 методами ECCA-Water (b) и NDWI (c)

Fig. 2. RGB composite of Sentinel-2 image of lake Krasilovskoye taken on 26.08.2016 (bands: Red, NIR, Blue) (а); water extraction results from the image taken on 15.09.2016 with ECCA-Water (b) and NDWI (c) methods

b

c

Рис. 2. Результаты выделения акватории оз. Большое Яровое по спутниковому снимку Sentinel-2 за 14.05.2016 методами ECCA-Water (a, b) и MNDWI (c)

Fig. 3. Water extraction results from Sentinel-2 image of lake Big Jarovoye taken on 14.05.2016 with ECCA-Water (a, b) and MNDWI (c) methods

Таблица 1. Сравнение результатов определения площади оз. Красиловское по данным Sentinel-2 различными методами

Table 1. Comparison of the results of determining the area of lake Krasilovskoe according to Sentinel-2 data by different methods

	Площадь (км2)
Дата съемки	ECCA-Water	NDWI	MNDWI	B/SWIR
04.06.2016	0.653	0.68	0.743	0.65
26.08.2016	0.631	0.65	0.668	0.602
15.09.2016	0.635	0.877	0.992	0.853
05.10.2016	0.643	0.672	0.685	0.618
СКО (исключая результат за 15.09.2016)	0.011	0.016	0.039	0.024

Таблица 2. Сравнение результатов определения площади оз. Большое Яровое по данным Sentinel-2 различными методами

Table 2. Comparison of the results of determining the area of lake Big Jarovoye according to Sentinel-2 data by different methods

	Площадь (км2)
Дата съемки	ECCA-Water	NDWI	MNDWI	B/SWIR
01.05.2016	73.55	74.17	74.35	73.54
14.05.2016	73.58	74.17	74.37	73.57
13.06.2016	73.6	76.52	78.93	76.04
11.09.2016	73.47	74.1	74.25	73.48
18.09.2016	73.54	74.15	74.26	73.53
СКО (исключая результат за 13.07.2016)	0.047	0.033	0.061	0.037

Анализ полученных результатов показывает, что методы на основе водных индексов относят тени от облаков к водным объектам, что приводит к существенным ошибкам в случае, если тени пересекают границу анализируемого водоема (например, на снимках оз. Красиловское за 15 сентября и оз. Большое Яровое за 13 июня). В то же время метод ECCA-Water не допускает подобных ошибок (рис. 2). Кроме того, методы на основе индексов NDWI и MNDWI ошибочно относят золоотвал ТЭЦ г. Яровое, расположенный возле химического предприятия «Алтай-химпром», к акватории оз. Большое Яровое (рис. 2 b ).

Если исключить снимки, на которых выделению водоемов мешают тени от облаков, то представленные результаты вычисления площадей в течение года характеризуются низкой дисперсией для каждого из методов. При этом наименее устойчивые показатели демонстрирует метод на основе индекса MNDWI.

Таким образом, метод ECCA-Water является наиболее подходящим для автоматического вычисления площади озер по спутниковым данным Sentinel-2. При отсутствии теней на изображениях также эффективен метод, основанный на применении водного индекса B/SWIR, однако данный индекс использует канал SWIR с разрешением 20 м. Эти методы применимы для – 5 –

Рис. 3. График изменения площади оз. Красиловское в 2017 г., построенный по данным спутника Sentinel-2 (км2)

Fig. 3. The graph of the change of lake Krasilovskoe area in 2017 according to Sentinel-2 satellite data (in square kilometers)

мониторинга площадей озер. На рис. 3 приведена внутригодовая динамика площади оз. Краси-ловское в 2017 г., вычисленная по данным Sentinel-2 с помощью метода ECCA-Water.

Анализ табл. 1, 2 показывает устойчивые расхождения между результатами различных методов выделения водных объектов. Детальное рассмотрение выявило, что отличия наблюдаются на границах воды и суши и объясняются тем, что различные методы по-разному классифицируют смешанные граничные пиксели. Площадь, занимаемая смешанными пикселями, напрямую зависит от разрешения съемки и длины береговой линии. Например, при пространственном разрешении 10 м для оз. Большое Яровое (рис. 2), имеющего правильную овальную форму с периметром 40 км, площадь смешанных пикселей составляет ~0.4 км2 (0.5 % от общей площади). А для озера Красиловское (рис. 1), имеющего сложную конфигурацию береговой линии с периметром 10 км, эта площадь составляет ~0.1 км2 (14 % от общей площади). Таким образом, ошибка определения площади может оказаться выше для водоемов меньшего размера и характеризующихся большим коэффициентом извилистости береговой линии (отношение длины береговой линии к длине окружности круга, имеющего площадь, равную площади озера).

Далее предлагается новый метод определения площади водных объектов по данным Sentinel-2, основанный на методе ECCA-Water и субпиксельном анализе смешанных пикселей.

Определение площади озер с использованием субпиксельного анализа смешанных пикселей

Для обработки смешанных пикселей будем использовать линейную модель смеси [16]. Она подразумевает, что спектральные характеристики I_M смешанного пикселя определяются взвешенной суммой характеристик входящих в него классов в соответствии с занимаемыми ими площадями. Для случая двух классов ( W – вода и NW – суша) это можно выразить следующим образом:

I M I W ^X S W + I NW ^X S NW , где S W + S NW 1.

Здесь SW и SNW – доли пикселя, занимаемые соответствующими классами. Для определения соотношения площадей классов внутри пикселя необходимо знать эталонные значения – 6 – спектральных характеристик этих классов (IW и INW). В одномерном случае SW можно найти по формуле

S W = ( I M – I NW ) / ( I W – I NW ).

В качестве эталонов, как правило, выбираются характеристики отдельных «чистых» пикселей. Однако фиксированные значения эталонов не подходят при обработке границ водной поверхности [6], что связано с существенными локальными изменениями спектральных характеристик воды в прибрежной зоне. Исходя из предположения, что эталонные значения характеристик воды должны быть более близки к характеристикам соседних водных пикселей, будем выбирать эталонные пиксели в окрестности каждого рассматриваемого смешанного пикселя.

В предлагаемом методе определения площади водоемов можно выделить три основных этапа.

Этап 1. Выделение поверхности воды на изображении с помощью алгоритма ECCA-Water. При этом изображение разбивается на множество кластеров, часть из которых помечена как водные. Для каждого i -го кластера вычисляется среднее значение водного индекса NDWI, которое будем обозначать NDWI _i ^avg. Рассматриваемый водный объект, определяемый как пространственно-связная водная область, выделяется в отдельный кластер L .

Этап 2. Определение смешанных пикселей скользящим окном размером 5×5 пикселей. Если в рассматриваемом окне содержится хотя бы один водный пиксель из кластера L и при этом присутствует хотя бы один неводный пиксель, то центральный пиксель окна считается смешанным.

Этап 3. Для каждого смешанного пикселя определяется доля водной поверхности с помощью линейной модели по формуле (2). При этом в качестве спектральных характеристик используются значения водного индекса NDWI. Эталонная характеристика суши INW определяется как средний показатель NDWIavg для всех пикселей внутри окна, не принадлежащих к i водным кластерам. Для установления эталонной характеристики воды для всех пикселей внутри окна, принадлежащих к водному кластеру L, определяется минимальное значение водного индекса NDWImLin . Учет данной локальной характеристики осуществляется с весом k, определяемым как отношение числа водных пикселей в локальном окне к общему числу пикселей в окне. В итоге эталонная характеристика воды определяется по формуле

I_w = NDWI min X k + NDWI avg x ( 1 - k ) .

Этап 4. Вычисление площади водного объекта с учетом результатов обработки смешанных пикселей.

В табл. 3 представлены площади оз. Красиловское и оз. Улугколь (расположенного в Республике Хакасия), вычисленные по спутниковым снимкам Ресурс-П и IKONOS соответственно при варьировании пространственного разрешения изображений (методом осреднения). Изображение IKONOS с разрешением 0.8 м было получено с помощью процедуры паншарпенинг в программном пакете ENVI. Из данных табл. 3 видно, что при ухудшении разрешения определяемая площадь озер без учета смешанных пикселей уменьшается. При этом применение процедуры коррекции позволяет значительно повысить устойчивость получаемых результатов и снизить их зависимость от пространственного разрешения съемки.

Таблица 3. Площадь озер до и после коррекции смешанных пикселей при изменении пространственного разрешения изображений

Table 3. The area of the lakes before and after mixed pixels correction with varying spatial resolution

Площадь оз. Красиловское, км2 Снимок Ресурс-П			Площадь оз. Улугколь, км2 Снимок IKONOS
Разрешение, м	Алгоритм ECCA-Water	После коррекции	Разрешение, м	Алгоритм ECCA-Water	После коррекции
2.14	0.728	0.736	0.8	6.18	6.21
6.42	0.692	0.738	3.2	6.12	6.19
10.7	0.696	0.752	9.6	6.08	6.22
15	0.672	0.748	16	5.99	6.22
СКО	0.019	0.008	СКО	0.080	0.013

Таблица 4. Площадь озера Красиловское на 20.05.2017 до и после коррекции смешанных пикселей

Table 4. The area of lake Krasilovskoe at 20.05.2017 before and after mixed pixels correction

Метод вычисления площади	Площадь, км2 Снимок Канопус-В	Площадь, км2 Снимок Sentinel-2
ECCA-Water без коррекции	0.660	0.674
ECCA-Water с коррекцией	0.728	0.713

Для валидации предложенного метода был проведен эксперимент с использованием муль-тиспектрального и панхроматического снимков Канопус-В оз. Красиловское, а также снимка Sentinel-2, полученных 20 мая 2017 г. Эталонная площадь оз. Красиловское, вычисленная визуально-инструментальным методом на основе панхроматического снимка с разрешением 2.1 м, составила 0.723 км2. В эксперименте площадь озера вычислялась в автоматическом режиме по мультиспектральным изображениям Sentinel-2 и Канопус-В (разрешение 10.5 м) методом ECCA-Water без проведения коррекции смешанных пикселей и с коррекцией. Результаты представлены в табл. 4. Ее анализ показывает, что коррекция смешанных пикселей позволяет увеличить точность определения площади озера.

Заключение

Результаты проведенных экспериментальных исследований показали, что методы на основе водных индексов обеспечивают выделение водных объектов по спутниковым данным Sentinel-2 в автоматическом режиме. Однако они могут ошибочно относить к воде некоторые антропогенные объекты. Кроме того, эти методы относят к водным объектам тени от облаков, что может приводить к существенным ошибкам при определении площадей водоемов. В то же время рассмотренный метод ECCA-Water, основанный на применении непараметрического ансамблевого алгоритма кластеризации ECCA и дерева решений, не допускает подобных ошибок.

В работе показано, что для водоемов небольшого размера смешанные пиксели значительно влияют на точность определения их площади по спутниковым снимкам. При – 8 – ухудшении пространственного разрешения съемков наблюдается снижение определяемой площади.

Предложен новый метод автоматического определения площадей озер на основе субпиксельного анализа смешанных пикселей, который позволяет снизить влияние пространственного разрешения спутниковых изображений на точность определения площадей исследуемых водоемов. Это дает возможность расширить область использования данных Sentinel-2 для мониторинга даже небольших по площади озер.