Математические методы аэрокосмического мониторинга воздействия на лесные массивы аварий на магистральных трубопроводах

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

В настоящее время на околоземных орбитах летает целый комплекс искусственных спутников Земли (ИСЗ), обеспечивающих получение информации в интересах различных отраслей науки и экономики. Получение информации о земной поверхности, в том числе о лесах и древесно-кустарниковой растительности, обеспечивают различные беспилотные космические аппараты (КА) многоцелевого назначения, созданные специально для исследования природных ресурсов Земли и ОПС, оснащенные преимущественно нефотографическими электронными съемочными средствами.

Для съемки земной поверхности в оптическом диапазоне в целях исследования природных ресурсов Земли используются различные ИСЗ, запускаемые разными странами на околоземные орбиты. Это российские ИСЗ серии “Ресурс-Ф”, “Ресурс-01”, “Океан”, “Метеор”. В 1970-1980-х гг. съемку земной поверхности производили с ИСЗ серии “Космос”, “Метеор-Природа” и др. В 2005 г. осуществлен запуск нового отечественного ИСЗ нового поколения - “Монитор-Э”, в ближайшем будущем намечается запустить ИСЗ “Ресурс-ДК- 1” и др. Зарубежные аэрокосмические системы мониторинга включают американские ИСЗ “Landsat”, NOAA, EOS (Terra), IKONOS, QuickBird, французский SPOT, индийский IRS и др.

При космических съемках часто используют нефотографические системы в виде сканеров (линейчатых или матричных) или лазеров. Предпочтение нефотографическим системам отдается в связи с тем, что они обеспечивают оперативную (в реальном масштабе времени) доставку информации на Землю по радиоканалам в цифровом виде, позволяющем производить ее автоматизированную обработку. Информацию, полученную при сканерной съемке, передают в виде цифрового изображения по радиоканалу на приемный пункт или записывают на борту на магнитный носитель. Высокая чувствительность сканеров позволяет выполнять съемку в узких (несколько десятков нанометров) спектральных интервалах, в пределах которых различия между некоторыми объектами ОПС более четко выражены. Кроме того в цифровых данных, полученных с помощью сканеров, отсутствуют шумы, которые неминуемо появляются при фотосъемке и фотолабораторной обработке съемочных материалов.

Применение лазеров - привело к разработке активных оптических съемочных систем. С помощью лазера облучают снимаемую поверхность. Отраженный от нее сигнал принимает оптическая система. В результате съемки получают трехмерное цифровое изображение. Особенностью лазерной съемки является возможность проникать через лиственный покров деревьев. Примером лазерной системы является аппаратура канадской фирмы Optech Inc. - ALMT-1020 (Laser Terrain Mapping System), которая разработана с целью трехмерных топографических съемок поверхности и создания крупномасштабных планов местности.

Следующие характеристики аэрокосмической съемки являются важными для целей анализа воздействий аварийных выбросов объектов ПНГХК и МТП на лесные массивы: число спектральных зон; спектральные диапазоны съемки; разрешение на местности; ширина полосы захвата.

Современные космические системы (IKONOS, QuickBird, IRS и др.) обеспечивают получение изображений с пространственным разрешением от 0,6 до Юм. Панхроматические или зональные космические изображения высокого разрешения получают при съемке в видимой и ближней инфракрасной частях электромагнитного спектра.

Формы и размеры крон деревьев входят в число важных признаков дешифрирования и классификации лесных массивов в зонах техногенного воздействия. Анализ космических изображений масштаба 1:5000 и 1:10000 показал, что по ним можно уверенно идентифицировать конусовидные, овальные (выпуклые) и плоские формы крон деревьев.

Космические снимки с пространственным разрешением от 10-15 до 30-50 м получают с помощью космических систем, наиболее типичные из которых “Ресурс О” МСУ-Э (Россия), Landsat ЕТМ+ и “TERRA” ASTER (США), SPOT (Франция), IRS (Индия) и др. Эти спутниковые изображения относят к снимкам среднего разрешения, а съемку обычно производят в видимой, ближней и средней инфракрасной частях спектра.

Следует отметить, что космические снимки среднего разрешения обладают повышенным спектральным разрешением и могут быть получены на одну и ту же территорию с известной периодичностью траектории ИСЗ. Это является преимуществом среднемасштабной съемки, так как создает дополнительные возможности для динамической классификации и мониторинга состояния лесного покрова с учетом фенологического состояния лесов.

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

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

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

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

Космические снимки низкого разрешения с пространственным разрешением больше 150 м получают с помощью многих космических систем, например “Ресурс О” МСУ-СК (Россия), “TERRA” MODIS, NOAA AVHRR (США), SPOT Vegetation (Франция) и др. В ряде случаев пространственное разрешение космических систем низкого разрешения характеризуют размером пикселя как в метрической, так и в градусной системе. Так, детальные данные радиометра AVHRR, например, представлены размером пикселя 1x1 км или 30"х30", а обобщенные - размерами пикселя 8x8 км (4'х4') или 15x15 км (5'х 5' по широте и долготе). При низком разрешении космических снимков каждый пиксель изображения может отличаться комбинацией как отдельных компонент, в первую очередь, растительности и почв, так и объединением разнородных экосистем, а также категорий земель лесного фонда. Эта информация является удобной для анализа крупномасштабных аварий и пожаров и их последствий. Чаще всего для анализа техногенных изменений в растительном покрове применяются комбинации спектральных каналов в виде так называемых вегетационных индексов.

Ниже приведены краткие описания и формулы вегетационных индексов, широко используемых в оценках воздействия аварийных ситуаций на газо- и нефтепроводах на лесные массивы. В приведенных ниже формулах B_},j = 1...7

обозначают спектральные каналы Landsat 7.

Нормализованный разностный индекс растительности (NDVI), который характеризует чистую продукцию растительного покрова и транспирацию влаги

Относительный индекс растительности (RVI), который разделяет различные виды растительности, а также участки растительности одного вида, но с различной биомассой

Индекс растительности, который аналогично RVI способен разделять различные виды растительности, а также участки растительности одного вида, но с различной биомассой

Трансформированный индекс растительности (TVI) чувствительный к варьированию растительной биомассы

Исправленный трансформированный индекс растительности (CTVI). Характеризует изменение растительной биомассы

Трансформированный индекс растительности (TTVI). Характеризует изменение растительной биомассы

Нормализованный относительный индекс растительности (NRVI) - это показатель различных типов растительности

Перпендикулярные индексы растительности (PVI) - это показатели различия спектральных свойств почвы и растительности. Первый исправленный перпендикулярный индекс растительности (PVI1) учитывает различия спектральных свойств почвы, растительности и воды

где а и b - параметры почвенной линии.

Второй исправленный перпендикулярный индекс растительности (PVI2) - учитывает различия спектральных свойств почвы, растительности и воды

где а и b - параметры почвенной линии.

Третий исправленный перпендикулярный индекс растительности (PVI3) - учитывает различия спектральных свойств почвы, растительности и воды

Дифференциальный индекс растительности (DVI) - учитывает наклон почвенной линии

Дифференциальный индекс молодняков (AVI) - характеризует различия спектральных свойств молодняков и спелых древостоев

Скорректированный на почву индекс растительности (SAVI) - минимизирует почвенный фон, для чего в знаменатель индекса вводится корректирующий параметр

где L - корректирующий множитель, учитывающий влияние почвы.

Первый трансформированный скорректированный на почву индекс растительности (TSAVI1) - учитывает параметры почвенной линии и минимизирует почвенный фон

Второй трансформированный скорректированный на почву индекс растительности (TSAVI2) - учитывает параметры почвенной линии, индекс растительности скорректированный на почву и минимизирует почвенный фон

Первый модифицированный скорректированный на почву индекс растительности (MSAVI1) минимизирует почвенный фон

где: LO = 1 - 2b*NDVI*WDVI.

Второй модифицированный скорректированный на почву индекс растительности (MSAVI2) минимизирует почвенный фон

Взвешенный дифференциальный индекс растительности WDVI - учитывает наклон почвенной линии

Показатель зеленого цвета растительности (GVI), который чувствителен к концентрациям хлорофилла.

Модифицированный показатель зеленого цвета растительности MGVI

Относительный показатель зеленого цвета растительности (GreenNDVI) - очень чувствителен к концентрациям хлорофилла

Нормализованный разностный снежный индекс (NDSI), который очень чувствителен к мощности снега и льда равен

Нормализованный разностный водный индекс (NDWI) характеризует содержание воды в зеленой биомассе

Индекс, характеризующий активность хлорофилла

Индекс, характеризующий влажность поверхности

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

Индекс, который характеризует общую яркость

Индекс, характеризующий влажность растительности

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

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

(КМР), которая для повышения эффективности анализа должна быть использована совместно с вегетационными индексами.

Масштабы нефтедобычи поражают своим размером и хорошо видны из космоса. Одной из наиболее эксплуатируемых в России природных территорий является район Самотлора. На снимке из космоса, приведенном на рис.2.1.1, хорошо просматривается густая сетка нефтепроводов.

Рис.2.1.1. Густая сеть нефтепроводов в районе Самотлора (космический снимок со спутника Landsat-7 19 сентября 1999г., аппаратура ЕТМ+, пространственное разрешение 30м)

На рис.2.1.1 приведен скорректированный по контрасту третий спектральный канал космического снимка со спутника Landsat-7 от 19 сентября 1999г., полученный с помощью сенсора ЕТМ+ (пространственное разрешение 30м). Космический снимок имеет семь спектральных каналов. Спектральные диапазоны сенсора приведены в табл.2.1.1.

Таблица 2.1.1. Спектральные диапазоны сенсора

ЕТМ+ спутника Landsat7

Третий спектральный канал представляет красную область видимого спектра, в которой хорошо контрастирует растительный покров и сеть нефтепроводов.

Самотлорское нефтяное месторождение (Самотлор)— крупнейшее в России и одно из крупнейших в мире месторождений нефти. Расположено в Ханты-Мансийском автономном округе, вблизи города Нижневартовск, в районе озера Самотлор (61° с.ш., 76° в.д.). Географическая карта месторождения представлена на рис.2.1.2.

Рис.2.1.2. Географическая карта района нефтяного месторождения Самотлор

Нижневартовск — административный центр Ханты-Мансийского автономного округа (60°57’ с.ш., 76°33’ в.д.). Население 240 тыс. чел. Город был основан в 1909 на правом берегу реки Оби как пристань Нижневартовская. Крупнейшие предприятия нефтегазового комплекса: «Самотлорнефтегаз», «Нижневартовское нефтегазодобывающее предприятие», «ТНК- Нижневартовск».

Месторождение открыто в 1965 году. Максимум добычи нефти (около 150 млн т в год) пришёлся на 1980-е годы. В 21-м веке в связи с применением современных способов интенсификации нефтедобычи выработка нефти увеличилась. Всего за годы эксплуатации месторождения на нём было пробурено 16700 скважин и добыто более 2,3 млрд т нефти. В настоящее время разработку основной части месторождения ведёт ОАО «Самотлорнефтегаз», принадлежащее компании «ТНК-ВР».

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

Рис.2.1.3. Увеличенный фрагмент второго спектрального канала самотлорского месторождения

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

Наиболее часто используется методы градиентного выделения контуров в сочетании с пороговым ограничением. В этих методах используются, как линейные, так и нелинейные операторы. Это могут быть простые аппроксимации градиента или же операторы-шаблоны, использующие различные шаблоны в разных направлениях. Широкое применение получили методы Превита, Робертса, Кирша, Лапласса и др. методы выделения контуров. Для большинства методов найдены оптимальные значения пороговых критериев. Наилуч- шие результаты по выделению контуров на компьютерном изображении были получены для методов Кирша, Собеля и Превита. Однако, универсальных методов выделения контуров не существует. Каждый метод хорошо работает лишь в своей индивидуальной области применения. Так, например, использование оператора Лапласса даже для гладкой границы в результате приводит к изрезанному контуру. Это свойство оператора Лапласса можно объяснить на основе анализа его Фурье-спектра:

где В{х,у), В(о)_х,о)_х) - яркостное поле и его Фурье-образ, со*, (Оу- Фурье частоты.

Из формулы (2.1.24) следует, что оператор Лапласса подавляет низкие и усиливает высокие частоты. Для случая выделения широких плавных контуров следует пользоваться критерием:

Глаз человека хорошо выделяет участки изображения, для которых |AS| = max. Этот критерий может быть использован для получения оператора выделения границ, который сочетает в себе свойства градиента и оператора

Лапласса:

Для выделения сложных контуров на компьютерных изображениях природных объектов мозаичной структуры наилучшим образом подходит совместное применение процентильного фильтра и оператора слабого дифференцирования. Алгоритм процентильного фильтра имеет следующий вид:

где К - размер Sw - окрестности, K_mi„ и К_тах - минимальное и максимальное значения размеров S_kN - окрестности, D = —^(B-B)² - дисперсия яркост-

N i

ного поля внутри скользящего окна, N - размер аппертуры, D_min и D_max - минимальное и максимальное значения дисперсии, которые являются настроичны- ми параметрами алгоритма, В₀ - центральный элемент окна, р - номер центрального элемента S_kN - окрестности в вариационном ряду яркостей, {BJ - Sw - окрестность наиболее близкая по яркости к центральному элементу В₀, В’ - преобразованная яркость центрального элемента. Минимальный размер Sw - окрестности К_тШ должен быть меньше размеров тех яркостных элементов, которые необходимо сохранить при обработке. Максимальный размер К_тах должен быть равен размеру окна анализа.

После обработки компьютерного изображения с использованием оператора процентильного фильтра применялся оператор слабого дифференцирования, алгоритм которого имеет следующий вид:

где Dj - дисперсия вдоль одного из восьми радиусов, выходящих из центральной точки.

В настоящее время для цифрового анализа аэрокосмических изображений разработано большое количество комплектов программ. Это и профессиональные специализированные пакеты, такие как Idrisi, Envi, Erdas Imagin, a также многочисленные пакеты компьютерного дизайна, такие как Photoshop, Corel Draw, Jasc Paint Shop, Macromedia Flash и др. Кроме этих специализированных инструментов для обработки изображений используются также универсальные инструментальные пакеты, такие как MathCad и MatLab. Библиотека обработки изображений пакета MatLab позволяет решать практические все профессиональные задачи, как в пространственной, так и в частотной области. При этом в MatLab доступны такие важные библиотеки, как библиотека вейвлет-преобразований, многомерной статистики, нейронных сетей и пр.

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

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

Рассмотрим некоторые примеры, иллюстрирующие возможности классических методов цифровой обработки аэрокосмических изображений.

Важнейшую роль при обработке любых изображений играет гистограмма. Гистограмма представляет собой график распределения пикселей изображения по полному диапазону яркостей. На рис.2.1.4 представлен фрагмент третьего спектрального канала спутникового изображения района месторождения Самотлор. Соответствующая гистограмма представлена на рис.2.1.5. Как следует из гистограммы, градации яркости изображения занимают интервал (20 - 100) из 256 градаций полного диапазона яркостей, что делает изображение однотонным и с трудно различимыми деталями.

Рис.2.1.4. Фрагмент третьего спектрального канала спутникового изображения района месторождения Самотлор

Рис.2.1.5. Гистограмма космического снимка, представленного на рис.2.1.4

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

где bj,a. - новое и старое значение яркости пикселя, ^_min^_max - границы нового и а ,а - старого диапазонов.

min max ^г

В пакете MatLab для контрастирования используется библиотечная процедура imadjust. Далее приведен пример программной реализации методов контрастирования изображений в пакете MatLab.

%IMG_STRETCH

% нормализуем яркость изображения

function img_stretch()

п=3;

fname=[f num2str(n) '.tif]; fnamel=['fs' num2str(n) '.tif];

F = imread(fname, ’tif); figure, imshow(F);

%--------Построение гистограмм -----------

histob(F);

%------Растягиваем яркостный диапазон--------

% Fs=imadjust(F,,); % 1-й канал % Fs=imadjust(F,,); % 2-й канал Fs=imadjust(F,,); % 3-й канал % Fs=imadjust(F„); % 4-й канал %Fs-imadjust(F„); % 5-й канал %oFs=imadjust(F„); % 7-й канал imwrite(Fs,fnamel, 'tif); figure, imshow(Fs); return

% ======= Построение гистограммы ==============

function histojb(img) h=imhist(img); v-1:10:256; hv=h(v); figure, bar(v,hv); axis(); return

Результат работы программы для компьютерной обработки изображения (рис.2.1.4) представлен на рис.2.1.6.

Рис.2.1.6. Результат коррекции яркостного диапазона третьего спектрального канала спутникового изображения нефтяного месторождения Самотлор, представленного на рис.2.1.4

На следующем этапе анализа сцен и ситуаций используют разнообразные фильтры. В пакете MatLab фильтры реализованы в библиотечной процедуре edge. В этой процедуре применяется инструмент скользящего окна совместно с лапласианом (2.1.24). При этом используются маски с различными коэффициентами. При выделении границ линейных объектов широко применяются маски, которые получили названия фильтров Собеля, Превитта, Робертса и пр.

Прежде всего следует изучить на практике возможности выделения контуров с помощью процедуры edge при различных комбинациях входных параметров. Одним из основных параметров процедуры является величина порога фильтрации. Порог фильтрации задает величину лапласиана, при котором пиксель относиться к классу границы

где Т -порог фильтрации.

Далее приведена программа, которую можно использовать при выполнении лабораторных работ.

%oIMG FILTERS

% фильтры подчеркивания контуров

function imgJiltersQ fname-fs3.tif;

F = imread(fname, 'tif); figure, imshow(F);

%--------Пространственные фильтры -----------

% filter_b(F, 'sobel', 0.1);

% filter_b(F, 'prewitt', 0.1);

% filter_b(F,'roberts',0.1); filter b(F,'log',0.01,0.1); return

% ======= Фильтры =============

function filter_b(Ffilt,T,sig)

%o F - изображение, filt - фильтр, T - порого % sig - стандартное отклонение гауссиана error(nargchk(1,4, nargin)); if nargin < 4, sig = 2; end if nargin < 3, T = 0.1; end if nargin < 2, filt - 'sobel'; end = edge(Ffilt,T);

I - imcomplement(B); figure, imshow(I); return

В качестве примера на рис.2.1.7 приведен результат цифровой обработки фрагмента снимка, представленного на рис.2.1.6, с помощью фильтра Робертса. Преимущество полученного результата (рис.2.1.7) по сравнению с исходным изображением рис.2.1.6 состоит в том, что результирующее изображение является бинарным, а его много проще, чем исходное аналоговое изображение, разделить на отдельные объекты.

Рис.2.1.7. Результат цифровой обработки фрагмента снимка, представленного на рис.2.1.6, с помощью фильтра Робертса

Для выделения объектов в библиотеке пакета Matlab имеется специальная процедура bwmorph, которая позволяет выполнять различные преобразования бинарных объектов. Результат применения одной из этих операций к бинарному изображению, представленному на рис.2.1.7, показан на рис.2.1.8. Преобразование получено с помощью применения следующей простой цепочки операторов:

Различные варианты применения функции bwmorph позволяют выделять объекты на аэрокосмических изображениях.

Рис.2.1.8. Результат выделения нефтепроводов в районе месторождения Самотлор с помощью применения цепочки операторов (2.1.31) к изображению, представленному на рис.2.1.7

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