СПОСОБ ОПРЕДЕЛЕНИЯ ПАРАМЕТРОВ СОСТОЯНИЯ ПОЧВЕННО-РАСТИТЕЛЬНОГО ПОКРОВА ПО ДАННЫМ МНОГОСПЕКТРАЛЬНОГО АЭРОКОСМИЧЕСКОГО ЗОНДИРОВАНИЯ Российский патент 2011 года по МПК G01V9/00 

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

Известен способ определения объема фитомассы пробных площадок (в поле, в лесу и т.п.). В полевых условиях скашивается травянистая растительность в пределах выбранной площадки (2×2 метра) и взвешивается сырой вес образцов растительности, а при необходимости определяется и сухой вес. Так находится объем сырой и/или сухой фитомассы травянистой растительности. Для отдельных деревьев объем фитомассы определяется по таким эмпирическим параметрам, как диаметр ствола на уровне груди. Для выбранных тестовых участков при наземных обследованиях также существуют эмпирические связи между значениями зеленой фитомассы и общей биомассы древесины (Выгодская Н.Н., Козодеров В.В., Косолапов B.C. и др. Опыт использования аэрокосмических измерений для количественной оценки массы листвы лесных насаждений. Журнал «Лесоведение», 1995, №5, с.37-47).

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

Известен способ определения состояния почвы, включающий аэрокосмическую съемку поверхности открытых участков почвы в ранний весенний или поздний осенний период, преобразование полученного изображения, оценку фототона этих площадей с последующим выделением контуров и вычислением площадей выделенных контуров с привязкой обработанного изображения к почвенной карте и его трансформацию (Патент РФ №2265839).

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

Наиболее близким аналогом к предлагаемому изобретению является способ обработки многоспектральных изображений, который заключается в том, что после приема и регистрации на магнитном носителе цифровых данных аэрокосмического зондирования производится их географическая привязка, учет геометрических, радиометрических, радиационных, спектральных искажений данных, их отображение в определенной проекции карты. Сазонов Н.В., Дедаев Ю.Н., Вандышева Н.М. Назад, в будущее - к 30-летию создания автоматизированной информационно-управляющей системы «АИУС-Агроресурсы» // Информационный бюллетень, 2009, №2(69) и/или на сайте по адресу http://www.gisa.ru/53086.html?searchstring=сазонов.

На базе этого способа был создан автоматизированный комплекс, известный, как АИУС-Агроресурсы (Автоматизированная информационно-управляющая система - Агроресурсы), позволявший осуществлять обработку аэрокосмических изображений, полученных с помощью измерительных приборов, установленных на самолеты и советские спутники «Ресурс-Ф» и «Ресурс-О».

Основной недостаток известного способа состоит в том, что вместо реальных количественных оценок параметров состояния природных объектов рассчитываются лишь упрощенные эмпирические характеристики, известные как "вегетационные индексы". На практике при построении тематических карт вегетационных индексов или близких к ним "информационных продуктов" используется лишь два спектральных канала, охватывающих полосу поглощения хлорофилла (основного пигмента фотосинтезирующей растительности) и область максимума спектральной отражательной способности вегетирующей растительности в ближней инфракрасной части спектра. Одним из указанных информационных продуктов является "вегетационный индекс нормализованной разности" (NDVI/Normalized Difference Vegetation Index), который в соответствии с названием вычисляется как разность значений яркости в указанных выше двух каналах, деленной на их сумму. Такой подход не рассматривает все множество возможных решений обратной задачи оценки состояния растительности и не позволяет судить о неопределенности восстановления интересующих специалистов-биологов таких параметров состояния лесной растительности, как, например, объем фитомассы листвы/хвои. Именно этот параметр служит основным показателем нормального или стрессового состояния экосистем, которое возникает под влиянием дефицита увлажнения, загрязнений окружающей среды и других факторов.

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

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

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

Сущность предлагаемого способа заключается в том, что при определении параметров состояния почвенно-растительного покрова по данным многоспектрального аэрокосмического зондирования сначала производят прием и регистрацию на магнитном носителе цифровых данных аэрокосмического зондирования. Затем эти данные обрабатывают путем проведения их географической привязки, учета геометрических, радиометрических и спектральных искажений данных с последующим отображением обработанных данных в определенной проекции географической карты исследуемой местности. Кроме того, при обработке данных производится их тематическая обработка, заключающаяся в распознавании образов наблюдаемых объектов. При тематической обработке изображений для каждого элемента разрешения производится распознавание объектов аэрокосмического мониторинга. Для этого выделяют следующие классы: «облачность», «водоемы», «почвогрунты», «растительность». Класс растительность включает в себя лесные, луговые, болотные, сельскохозяйственные и другие экосистемы. Для элементов разрешения, относящихся к классу «растительность», производят количественную оценку параметров состояния этих объектов. Для лесной растительности определяется породный состав: лиственные, хвойные, смешанные экосистемы. Также определяется тип подстилающей поверхности (почва, трава, кустарники, болота и т.д.) и продукционные характеристики (скорость изменения содержания углерода в соответствующих экосистемах).

Перечень чертежей и краткое их описание.

Фиг.1. Результаты обработки многоспектрального изображения аппаратуры MODIS для района г.Москвы. Верхний - восстановление биомассы зеленых фитоэлементов (т/га); нижний - фотоизображений из системы Google Earth. На горизонтальной и вертикальной осях представлены приблизительные значения долготы и широты в градусах.

Фиг.2. Обработка многоспектрального изображения аппаратуры ЕТМ + для района г.Тверь. Верхний - восстановление биомассы зеленых фитоэлементов (т/га); нижний - фотоизображений из системы Google Earth. На горизонтальной и вертикальной осях представлены приблизительные значения долготы и широты в градусах.

Фиг.3. Пример определения продукционных характеристик леса для района г.Тверь по данным аппаратуры ЕТМ+. Красно-зеленой цветовой гаммой обозначена чистая первичная продукция (г/м²/год) лесной растительности, черным цветом обозначена растительность, не классифицированная как лес, серый - открытые почвогрунты, темно-серый - антропогенные объекты, синий - вода, фиолетовый - пропуски данных измерений.

Фиг.4. Принцип поиска частного решения обратной задачи.

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

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

где j - номер канала приемной аппаратуры k-го типа с функцией чувствительности аппаратуры R; Ω - телесный угол визирования; i - тип леса; r - тип межкроновой растительности; В - биомасса листвы/хвои; θ, φ (θ', φ') - зенитный и азимутальный углы визирования (падающего излучения), соответственно; D_canopy - сомкнутость полога; D_crown - плотность крон; δ₁ и δ₂ - доля теней в межкроновой области и на кронах деревьев; ρ₁, ρ₂, ρ₃ - спектральные коэффициенты отражения межкроновой области, кроновой области и многократного рассеяния внутри кроновых просветов; U_S и H_S - спектральные интенсивности падающего прямого и рассеянного солнечного излучения.

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

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

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

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

Частные решения обратной задачи находятся в процессе поиска соответствия многоканальных данных дистанционного зондирования в координатах "плотность лесного полога (D_canopy) - ажурность крон деревьев (A_crown)" наиболее близким численным значениям модельных яркостей прямой задачи (в дальнейших рассуждениях будет использоваться параметр плотности кроны D_crown=1-А_crown). На Фиг.4 приведен иллюстративный пример поиска частных решений для пары спектральных каналов. В каждом спектральном спутниковом канале модельные яркости, которые представляют собой поверхность L=L(D_canopy,D_crown), секутся плоскостью на уровне, соответствующем измеряемым яркостям. Частным решением является точка пересечения полученных таким образом линий уровня в плоскости аргументов D_canopy и D_crown. Каждой такой точке соответствует значение биомассы листьев B=B(D_canopy,D_crown), которая определяется в соответствии с известной параметризацией материалов аэрокосмической съемки.

Поскольку рассматриваемая обратная задача относится к классу математических задач, некорректных по Адамару, ее решение может не существовать, быть не единственным или существенно неустойчивым. Таким образом, для получения оптимального решения необходимо применять методы регуляризации. Очевидно, что расчетная база данных не описывает всех вариаций спектров растительных объектов, которые могут быть получены в результате реальных измерений. Таким образом, поиск решений рассматриваемой обратной задачи осуществляется с некоторой заданной точностью |L_измер.-L_{модельн.}|<Δ, (параметр Δ будем называть областью поиска), значение которой может изменяться в процессе поиска. Важно также установить разумное начальное значение Δ. Для грубого начального приближения точность решения будет невелика, при малых затратах компьютерного времени. Если же уменьшить начальное значение области поиска решения, то компьютерное время расчета будет большим, но при этом возрастет и точность. В расчетной процедуре область поиска увеличивается до тех пор, пока не находятся частные решения во всех парах каналов. В случае неединственности из полученного ансамбля решений выбирается наиболее вероятное.

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

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

1. Тип объекта (растительность, водные объекты, облачность, открытые почвогрунты, антропогенные объекты).

2. Прозрачность атмосферы. Обратная задача решается только для значения прозрачности атмосферы Р>0.7 (при углах визирования в надир). Приняты следующие градации: 1 - мутная Р~0.75, 2 - средне-прозрачная Р~0.79, 3 - легкая дымка Р~0.83, 4 - прозрачная Р ~ 0.88.

3. Объем фитомассы листвы/хвои лесной растительности. Возможные значения варьируются от нуля до 35 т/га.

4. Дисперсии ошибки объема зеленой фитомассы растительности в (т/га)².

5. Тип леса. Выделяется 11 классов породного состава от полностью лиственных до полностью хвойных пород.

6. Тип межкроновой растительности. Выделяется до 12-ти классов в соответствии с исходными данными модельных представлений о спектральных образах: яркая трава; более темная трава; болотная растительность: более темная и более светлая; мох; отсутствие растительности; кустарники; преобладание открытой водной поверхности между кронами деревьев.

7. Изрезанность верхней границы лесного полога.

8. Сомкнутость полога. Параметр означает относительную площадь покрытия элемента разрешения кронами деревьев и принимает значения от нуля до единицы. Значение 0 соответствует отсутствию лесной растительности, 1 - полному покрытию лесом или кустарниковой растительностью.

9. Ажурность крон деревьев. Параметр характеризует плотность листвы для отдельного дерева и имеет значения от единицы до нуля. 1 означает отсутствие листвы (голые ветки), 0 соответствует максимально плотной кроне.

10. Площадь проективного покрытия. Параметр означает относительную меру заполнения поверхности почвы фитоэлементами растений при рассмотрении сверху. 0 - отсутствие растительности, 1 - полное покрытие элемента разрешения фитоэлементами лесной растительности.

11. Доля поглощенной фотосинтетически активной радиации. Параметр означает долю солнечной радиации в диапазоне 380-710 нм, поглощенной листьями деревьев при фотосинтезе. Формально параметр принимает значения от 0 до 1.

12. Чистая первичная продукция. Параметр означает продукцию годовую углерода лесной древесной растительностью в единичный период времени. Единицы измерения - г/м²/год.

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

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

Следующие примеры иллюстрируют данный способ.

На Фиг.1 приведен пример обработки изображения аппаратуры MODIS спутника Terra на дату съемки 30.05.2002, соответствующей г.Москва и близлежащих районов Московской области. Данные MODIS имеют достаточно грубое пространственное разрешение (500 м), что в определенной степени затрудняет классификацию почвогрунтов, небольших водоемов и антропогенных объектов. Таким образом, изображение классифицировалось на облачность (белый цвет), водоемы (синий цвет), плотная городская застройка (серый цвет) и растительность (красно-зеленая цветовая гамма).

Валидация данного изображения может быть проведена на качественном уровне с использованием фотографических изображений из известной системы Google Earth, в которой данный регион представлен с высоким пространственным разрешением. Видно, что удается распознать крупные водоемы - Клязьминское, Пироговское, Учинское, Пестовское и Икшинское водохранилища. Река Москва практически везде в черте города имеет ширину 80-250 м, и таким образом она классифицируется только в районе Серебряного Бора, станции метро Кожуховская и Борисовских прудов. Всем крупным лесопарковым зонам внутри города соответствует значительный объем биомассы. Так, например, в северной части можно видеть парк Сокольники, Лосиный остров, Главный Ботанический Сад РАН, парк Сельскохозяйственной Академии им. К.А. Тимирязева, лесопарки Покровкое-Глебово, Серебряный Бор, Химкинский, Красногорский и Хлебниковский. В южной части хорошо видны Битцевкий и Бирюлевский лесопарки, в Восточной части - Измайловский и Кузьминский лесопарки.

На Фиг.2 приведен пример обработки изображения аппаратуры ЕТМ + спутника Landsat-7 на дату съемки 02.09.1999, соответствующей восточной части г.Тверь и ближайших окрестностей. Данные ЕТМ + имеют пространственное разрешение 30 м, что позволяет существенно улучшить качество классификации. На обработанном снимке выделены следующие классы: растительность (красно-зеленая цветовая гамма), водоемы (синий цвет), облачность (белый цвет, на приведенном примере отсутствует), открытые почвогрунты (светло-серый цвет) и антропогенные объекты (темно-серый цвет). Фиолетовым цветом выделены пропуски в данных. В приведенном примере пропуски возникали из-за солнечных бликов на водной поверхности. Данное изображение можно также сравнить с фотографическими снимками системы Google Earth.

На обработанном снимке можно видеть, что наряду с крупными водными объектами, такими как р. Волга и песчаный карьер, воспроизводятся более мелкие объекты, соответствующие городским отстойникам и портовой зоне. Ложного распознавания водных объектов не выявлено. Также хорошо распознаются открытые грунты (например, песчаные берега карьера) и промышленные зоны. Можно видеть, что наибольшие значения биомассы соответствуют лесным массивам. Наряду с крупными загородными лесными угодьями можно увидеть лесопарковые зоны внутри г.Тверь, например Бобачевскую рощу.

Кроме параметров, характеризующих состояние растительного покрова, модель позволяет оценивать некоторые продукционные характеристики, такие как доля поглощенной фотосинтетически-активной радиации и чистая первичная продукция. Под фотосинтетически-активной радиацией (ФАР) понимают ту часть солнечной радиации, которая поглощается растением при фотосинтезе. В российских работах - это диапазон длин волн 380-710 нм, именно такое определение используется в нашей программе. В зарубежных работах ФАР соответствует диапазону 400-700 нм.

Для расчета доли поглощенной ФАР в нашей модели применяется следующее соотношение FAPAR=A*F*[1+R*(1-(1-T)*F)], где А - коэффициент поглощения фитоэлементами деревьев различного типа, F - значение площади проективного покрытия, рассчитанное в нашей модели, R - коэффициенты отражения межкроновой растительностью различного типа (функция типа межкроновой раститительности), Т - коэффициент пропускания фитоэлементами деревьев различного типа (функция состава лесной растительности).

Чистой первичной продукцией называют массу углерода, поглощенного растениями из атмосферы на единичной площади в единичный период времени за вычетом дыхания растений. Устоявшаяся единица измерений - г/м²/год. В нашей методике используются параметризации - зависимости чистой первичной продукции от биомассы листвы, построенные по данным Московскго государственного университета леса. Параметризации адаптированы к растительности восточно-европейской части России. Туда вошли данные по сосне, ели, березе, клену, осине и ясеню. Соотношение для вычисления чистой первичной продукции имеет вид

где Y - чистая первичная продукция (текущий прирост углерода в г/м²/год), Х - биомасса (в т/га), W - процент хвойного леса в лиственном, А и В - параметры, оцениваемые по статистическим данным отдельно для лиственных и хвойных лесов (А<1).

Изобретение относится к области аэрокосмических методов исследования Земли и может быть использовано при оценке состояния почвенно-растительного покрова по данным дистанционного зондирования. Сущность: производят прием и регистрацию на магнитном носителе цифровых данных аэрокосмического зондирования. Обрабатывают эти данные путем проведения их географической привязки. Учитывают геометрические, радиометрические и атмосферные искажения данных. Отображают обработанные данные в определенной проекции географической карты исследуемой местности. Производят преобразование данных. Рассчитывают спектральную интенсивность уходящего излучения. Производят тематическую обработку с выделением класса «растительность». Определяют объем фитомассы для каждого элемента разрешения, относящегося к классу «растительность». Определяют параметры состояния, характеризующие тепло-, влаго- и энергообмен между пологом растительного покрова и атмосферой. Определяют биопродуктивность по классу «растительность». Технический результат: повышение эффективности обработки многоспектральных изображений. 4 ил.

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