CC BY
55
Серия «Науки о Земле» 2013. Т. 6, № 2. С. 114-124 Онлайн-доступ к журналу: http://isu.ru/izvestia
И З В Е С Т И Я Иркутского государственного университета
УДК 004.93+14:902.2
Дешифрирование мультиспектральных космических
снимков для решения задач поиска и идентификации
*
археологических памятников
Е. П. Крупочкин ([email protected]) А. Р. Батуев ([email protected])
Аннотация. В статье представлены результаты разработки и апробации алгоритмов автоматизированного поиска и идентификации археологических памятников по данным мультиспектральной космической съемки. Рассмотренные примеры опираются на экспериментальный материал, полученный по репрезентативным площадкам Алтая (бассейн р. Чуя). Для решения поставленных цели и задач предлагается технология алгоритмизированного дешифрирования, основанная на двух подходах. Первый базируется на математическом описании пространственного распределения и корреляции яркостей элементов космических снимков, для реализации второго предложен алгоритм многоступенчатой классификации. В ходе экспериментов выявлены проблемы, которые снижают возможности использования разработанных алгоритмов и модулей: сезонность съемки и высокая дифференциация спектральной кривой по индексным изображениям, невысокое пространственное разрешение при сравнительно малых размерах объектов поиска, отсутствие спектральных библиотек заданных районов поиска и др. Вместе с тем полученные результаты свидетельствуют о целесообразности применения разработанных подходов и алгоритмов, что подтверждается стабильным результатом при значительной вариации параметров и условий поиска.
Ключевые слова: мультиспектральная космическая съемка, данные дистанционного зондирования, дешифрирование объектов археологии, алгоритмы компьютерного распознавания объектов поиска, спектральные свойства растительности, индикационные дешифровочные характеристики археологических памятников.
Введение
Один из наиболее острых вопросов сибирской археологии заключается в отсутствии детальной археологической карты как на территорию Сибири, так и на территорию Алтая. Для ее разработки приоритетное значение имеет решение не только прикладных, но и теоретических вопросов, связанных:
• с развитием концепции геоархеологического картографирования;
* Исследования выполнены при финансовой поддержке РФФИ, грант № 13-05-00995а.
• с созданием сквозной технологии, позволяющей разрабатывать цифровые археологические карты требуемого масштаба и территориального охвата;
• с внедрением в археологические изыскания единых стандартов и форматов данных и др.
На современном уровне эволюции методов геоинформатики и дистанционного зондирования стало возможным оперативное изготовление археологических карт с высокой точностью, детальностью и наглядностью. С учетом явных научных потребностей археологии с одной стороны и новых возможностей в области информационных технологий с другой актуальной задачей является разработка новых методов поиска, идентификации и картографирования памятников археологии с использованием ГИС и мультиспектральных данных дистанционного зондирования.
Цель исследований заключается в разработке алгоритмов автоматизированного дешифрирования мультиспектральных космических снимков для поиска и последующего ГИС-картографирования памятников археологии на примере репрезентативных площадок Алтая (бассейн р. Чуи). Район исследований в физико-географическом плане расположен в Юго-Восточном Алтае, в территориально-административном - Кош-Агачском и Улаганском районах Республики Алтай.
Для достижения цели мы предлагаем технологию алгоритмизированного дешифрирования, основанную на двух подходах: 1) базируется на математическом описании пространственного распределения и корреляции яркостей элементов космических снимков, 2) реализован по алгоритму многоступенчатой классификации Decision Tree.
В первом подходе использована взвешенная мера количественной оценки близости текстурных образцов, построенных с применением струк-турно-яркостных характеристик. Для определения оптимальной меры количественной оценки близости текстурных образцов проведены соответствующие эксперименты. При этом особое внимание уделено анализу взвешенного евклидового расстояния.
Основой предлагаемого метода является теория марковских случайных полей с гиббсовским распределением вероятности (МСП/ГРВ). Согласно данному подходу цифровое изображение представляет собой совокупность элементов (в случае цифрового изображения - пикселей), взаимодействующих между собой. В МСП/ГРВ-модели предполагается, что состояние каждого элемента растра (значение сигнала в этом элементе) зависит от состояний его соседей по окрестности. Из этого следует, что в зависимости от текстуры меняется характер взаимодействия элементов растра. Соответственно, для определения границ однородных текстурных участков изображения необходимо произвести анализ взаимодействия элементов растра между собой.
С позиции формализованного, т. е. алгоритмического описания предлагаемого метода, его основой является совокупность элементов рас-
тра, так называемая клика. Так, кликой пары (Я, ■) является подмножество к с Я такое, что:
1) к состоит из одного элемента растра (т, п);
2) к - это такое подмножество элементов растра, в котором любая пара элементов является соседней (или парой элементов для одной окрестности). Таким образом, если для любого (т,п) Ф (к,1) (т,п) е к и (к,1) е к, то (т, п) е ■(к,1) и (к,1) е ■(т,п).
Разработка алгоритма поиска значимых клик
Эксперименты с генерацией и сегментацией искусственных и натуральных визуальных образцов с помощью МСП-моделей позволили найти эвристическую стратегию поиска, дающую удовлетворительные по релевантности результаты в описании структур взаимосвязей однородных текстур:
• найти гистограммы различий оттенков серого (ГРОС) для всех семейств клик в окне поиска;
• вычислить расстояния между нормализованными ГРОС и распределением НСП;
• выбрать семейства клик, чьи расстояния превышают заданное пороговое значение.
Предлагаемый алгоритм может служить основой для получения разметок элементов растра как результата сегментации снимков высокого и сверхвысокого разрешения. Технология представляет собой итеративный процесс расчёта значимых клик участков, представленных на анализируемом изображении. Расположение и размер областей расчёта клик определяются координатами и размером движущихся расчётных прямоугольных окон. Физически алгоритм реализуется на основе прохода окна как по горизонтали (построчно), так и по вертикали (по столбцам).
Работа алгоритма по выделению однородных участков кусочно-текстурного изображения на основе применения модели эскиза текстурного образца включает несколько стадий.
Стадия 1. Нахождение значимых клик в двух соседних расчётных окнах. Численное выражение полученного эскиза текстуры записываем в специальную матрицу, далее называемую матрицей значений клик. Элемент матрицы Хц равен 0, если клика типа I, взятая в направлении ц, не является значимой, и 1 - если она является таковой. Таким образом, каждый элемент матрицы значений клик ||хц|| является численным выражением степени значимости клики типа I, взятой в направлении Ц.
Стадия 2. Формирование матрицы различий и анализ данных матрицы различий. Предполагается, что участки матрицы ||Оу||, соответствующие однородным участкам кусочно-текстурного изображения (под однородными участками в данном случае понимается участок изображения, содержащий не более один текстурный образец), должны содержать элементы, либо равные между собой, либо имеющие минимальное различие. Напротив, участки матрицы, относящиеся к окрестностям элементов границ
между различными участками анализируемого изображения, должны содержать элементы, существенно отличающиеся по значениям от остальных.
Стадия 3. По индексам элементов, содержащих существенно отличающиеся значения, восстанавливаются соответствующие координаты пикселей анализируемого изображения. Найденные пиксели визуализируются, представляя полученные границы между различными текстурными участками анализируемого кусочно-текстурного изображения.
Выбор граничных элементов в матрице различий осуществляется на основе предположения о том, что десятипроцентная погрешность является достаточной для того, чтобы утверждать факт сходства или различия двух выборок [3].
Для решения задачи идентификации археологических объектов на снимках было разработано приложение Object Detection (рис. 1).
Рис. 1. Алгоритмизированная блок-схема, описывающая работу модуля Object Detection в среде IDL
Первая фаза тестирования алгоритма, проведенная ранее [2], показала в целом положительный результат. Эксперименты показали целесообразность его применения при обработке снимков высокого и сверхвысокого разрешения. По нашим данным, такое разрешение должно составлять как минимум 'Л от среднего размера объектов.
Дальнейшая апробация разработанного продукта выполнялась по имеющимся мультиспектральным снимкам QuickBird и Ikonos. Особенность данного подхода - предварительный визуальный анализ снимков с созданием текстурных образцов ROI (регионов интереса). Теоретически ROI могут включать как точки, так и линии и полигоны. Тем не менее для объектов, имеющих незначительные (2-3 пикселя или - 4 ~ 6,5 м) размеры, следует устанавливать точечный тип ROI. Если объекты имеют большой диаметр или площадь, следует использовать полигональный (площадной) тип ROI-образца. После записи и сохранения ROI непосредственно запускается сам модуль.
Далее задается 2 ключевых параметра - размер окна поиска (Window Size) и шаг сдвига оператора (Offset).
Проведенные эксперименты показали, что если размеры объектов поиска (Ri) приблизительно известны, в этом случае следует установить как минимум двойное значение Window Size, т. е. 2*R¿ пикселей. С другой стороны, значительное увеличение окна оператора не приводит к улучшению достигнутого результата. Для параметра Offset оптимальными являются значения в диапазоне от 1/2 до 1/3 от размера окна.
Приведем пример обработки снимка QuickBird с разрешением 2,4 м. Визуально на снимке с большим трудом можно обнаружить группы памятников или их отдельные элементы. Однако применение машинного анализа приведенного снимка позволяет получить более убедительную картину. Так, к примеру, использование алгоритма классификации Isodata по отношению к обработанным модулем Object Detection данным позволило распознать на снимке группу объектов (рис. 2), а также отдельно расположенные объекты - тюркские оградки, межевые знаки, каменные кольца и т. п.
Оперирование настройками разработанной программы позволило улучшить результат поиска, особенно для объектов, выстроенных в цепочку. В таких случаях больший интерес представляют обнаруженные группировки объектов, например цепи курганных насыпей пазырыкской культуры.
Второй подход, как было отмечено выше, реализован на основе алгоритма Decision Tree. Особенность алгоритма заключается в возможности многоступенчатой классификации, при этом классификатор состоит из ряда бинарных решений, которые используются, чтобы правильно определить категорию для каждого пикселя. Решение введенного выражения на каждой ступени классификации делит пиксели изображения на два класса. Каждый новый класс может быть разделен еще на два класса по другому выражению и так до получения желаемого результата.
Рис. 2. Результат классификации изображения, обработанного модулем Object Detection
Отличительная особенность алгоритма - возможность использования в одном классификаторе данных из различных источников и файлов. Главным инструментом анализа является язык IDL, использующий синтаксис, аналогичный сервису ENVI Band Math. Выражения включают математические функции, операторы отношений, функции - IDL (сложение (+); вычитание (-); умножение (*); деление (/); sin(x); cos(x); tan(x); arcs - asin(x), acos(x), и atan(x); sinh(x), cosh(x), и tanh(x); операторы отношений (LT, LE, EQ, NE, GE, GT) и логические операторы (AND, OR, NOT, XOr) и др.). Таким образом, в ходе обработки создается новый бинарный файл с двоичными значениями 0 и 1. Результат 0 посылается по ветви «Нет», результат 1 - по ветви «Да».
Как видно из приведенного примера (рис. 3), процедура включает три основных шага: 1) вычисление Ndvi с построением индексного файла, 2) переклассификацию полученного массива с заданным условием, 3) формирование бинарного файла с требуемым результатом.
Начало алгоритма
//Шаг 1
Вычисление индекса NDVI
(float(b4)-float(b3)/float(b3)-float(b4))
//Шаг 2
Если значение Ndvi > 0. 16
Сохранить пиксели, удовлетворяющие условию в массиве «Class 1»
Иначе
Обнулить пиксели из массива «Node 1»
Если значение Ndvi < 0. 3
Сохранить пиксели, удовлетворяющие условию в массиве «Class 1»
Иначе
Обнулить пиксели из массива «Node 1»
//Шаг 3
Цикл по всем значениям массива «Class 1» образовать в новый класс и присвоить значение 1
По всем нулевым значениям массива «Node 1» образовать новый класс «Class 0» Конец цикла
Конец алгоритма
Рис. 3. Алгоритм классификации мультиспектрального снимка, реализованный на языке IDL
Приведем анализ одного из результатов обработки (рис. 4). Согласно заданной в алгоритме последовательности вначале выполняется расчет вегетационного индекса NDVI (Normalized Difference Vegetation Index). Его вычисление базируется на двух наиболее стабильных участках кривом спектральной отражательной способности сосудистых растений - красной (0,62-0,75 мкм) и ближнем инфракрасном (0,75-1,3 мкм) участках спектра. Его значения тем больше, чем больше зеленая фитомасса. Однако на его значения также оказывают влияние видовой состав растительности, ее сомкнутость, состояние, экспозиция и угол наклона поверхности, цвет почвы под разреженной растительностью и др. [4].
По данным полевых описаний проективных покрытий и проведенных ранее тестов на примере репрезентативных площадок археологического комплекса «Юстыт», мы обнаружили, что наибольшее число объектов сосредоточено в интервале NDVI от 0,16 до 0,3. Аналогичные значения были заданы и в алгоритме поиска для района «Курайка». Как видно из полученного результата, большая часть объектов сосредоточена именно в классе - 1, что объясняется близкими физико-географическими условиями. Таким образом, объекты распределены на границе между участками почвы (без сосудистых растений) и разреженной растительностью (см. табл.). Классы, не попавшие в данный интервал или превышающие его, обнули-лись до фоновых. В сравнении с методами кластеризации мы получили более стабильный результат, учитывающий физико-географические особенности района поиска.
Рис. 4. Фрагмент бинарного снимка Ikonos, полученного в результате обработки по алгоритму Decision Tree (участок долины р. Курайки, бассейн р. Чуи)
Таблица
Значения NDVI для зеленой растительности (по данным ИТЦ «СканЭкс» [1])
Тип объекта Коэффициент отражения* в красной области спектра Коэффициент отражения* в ближней инфракрасной области спектра Значение NDVI
Густая растительность 0,1 0,5 0,7
Разреженная растительность 0,1 0,3 0,5
Открытая почва 0,25 0,3 0,025
Облака 0,25 0,25 0
Снег и лед 0,375 0,35 -0,05
Детальный анализ изображения позволил также выявить ряд объектов, относящихся к классу - 0. Как выяснилось, для последних характерно расположение в лесных колках, что затрудняет процесс индикационного дешифрирования. При этом необходимо учитывать средние размеры археологических объектов, которые соответствуют 3-5 м. В совокупности ряд негативных факторов - незначительные размеры объектов поиска,
плохая выраженность в рельефе, высокая лесистость отдельных участков района исследований и др. - понижают достоверность результатов дешифрирования на отдельных участках до 10 %.
Обсуждение результатов исследований
Практическая реализация первого подхода, реализованного по данным космической съемки в оптическом диапазоне, показала удовлетворительные результаты с достоверностью поиска 60-87 %. Как выяснилось, точность алгоритмов дешифрирования определяется морфометрической сложностью объектов и их размерами. На основе анализа обработанных снимков установлена следующая закономерность - чем меньше размер объекта, тем менее достоверен результат.
В силу известных причин мы не всегда имеем возможность тестировать новые алгоритмы по снимкам сверхвысокого разрешения. Тем не менее проведенные эксперименты (разрешение 1,0 м и 2,4 м) показали целесообразность использования предлагаемых разработок. Расчеты показали -для оптимизации затрат и получения максимального результата пространственное разрешение снимка должно составлять как минимум 'Л от среднего размера объектов. Таким образом, при диаметре курганных насыпей от 4 до 8 м и разрешении 2 м результат будет более чем убедительным.
Второй подход основан на реализации метода многоступенчатой классификации. Его преимущество заключается в возможности включения в алгоритм поиска одновременно нескольких условий и показателей. В процедурах предобработки и постобработки согласно условиям алгоритма могут быть использованы дополнительные сцены (фрагменты снимков). Кроме того, в сравнении с другими методами анализа многозональных изображений, метод многомерной классификации выдает результат в виде бинарного файла. Таким образом, результатам поиска всегда соответствует только один класс, при этом нули образуют фон. Такая форма представления данных наиболее удобна, поскольку исключает необходимость неоднозначной интерпретации полученных результатов. Еще одним немаловажным обстоятельством является скорость и точность векторных преобразований бинарных файлов, что также подчеркивает целесообразность применения апробированного нами метода при решении поставленных задач.
Наряду с положительными результатами отметим ряд вопросов, проработка которых необходима в ближайшей перспективе. Так, к примеру, многие объекты слабо выражены в трёхканальном изображении, соответственно модуль поиска Object Detection не может их распознать. Отсюда необходимо дальнейшее совершенствование алгоритма поиска с ориентацией на мультиспектральную съемку. Кроме того, проблема, связанная с существенными временными затратами на процедуру обработки ДДЗ с большой площадью покрытия, решена частично. Поэтому вполне очевидно - формализация и совершенствование структурной схемы разработанного модуля является актуальной задачей.
Как видно из представленных результатов, общей проблемой для всех апробированных методов и алгоритмов дешифрирования является сезонность съемки и высокая дифференциация спектральных отражательных свойств подстилающей поверхности. Таким образом, для повышения эффективности алгоритмов поиска и идентификации археологических объектов необходимо продолжение цикла экспериментов с использованием ДДЗ разного пространственного и радиометрического разрешения. Для совершенствования данной процедуры и повышения дешифровочных свойств снимков будут разработаны спектральные библиотеки для заданных районов поиска.
Исследования выполнены при финансовой поддержке РФФИ, грант № 13-05-00995а.
Список литературы
1. Использование карт вегетационных индексов, полученных по космическим снимкам: NDVI в сельском хозяйстве // Офиц. сайт инженер.-технол. центра «СканЭкс». - URL: http://www.scanex.ru/ru/monitoring/default. asp?submenu= farm-ing&id=ndvi. (дата обращения: 07.10.2013).
2. Крупочкин Е. П. Археологическое картографирование в России: состояние и пути развития / Е. П. Крупочкин, К. А. Боенко // Изв. Алт. гос. ун-та. - 2010. -№ 3-2(67). - С. 110-119.
3. Кук Д. Компьютерная математика. / Д. Кук, Г. Бейз. - М. : Наука, Гл. ред. физ.-мат. лит., 1990. - 384 с.
4. Черепанов А. С. Вегетационные индексы // Геоматика. - 2011. - № 2. -С. 98-102.
Interpretation of multispectral satellite imagery to meet the challenges of search and identification of archaeological sites
A. R. Batuev, E. P. Krupochkin
Annotation: The paper presents the results of the development and testing of algorithms for automated search and identification of archaeological sites using multispectral imagery. These examples are based on the experimental data obtained on representative sites of Altai (watershed Chuya). To achieve the goals and objectives of the proposed technology algorithmized interpretation, based on two approaches. The first is based on the mathematical description of the spatial distribution and correlation of luminance elements of satellite images for the implementation of the second algorithm of multistage classification. The experiments reveal problems that reduce the possibility of using the developed algorithms and modules - seasonal shooting and high differentiation of the spectral curve of the images on the index, low spatial resolution at a relatively small size objects of the search , the lack of spectral library search target areas , etc. However, the results demonstrate the feasibility of using the developed approaches and algorithms, as evidenced by a stable result in significant variations in parameters, and search terms.
Key words: multi-spectral satellite imagery, remote sensing data interpretation archaeological sites, object recognition algorithms for computer search, spectral properties of vegetation, interpretive Indicative characteristics of archaeological sites.
Крупочкин Евгений Петрович, кандидат географических наук, доцент Алтайский государственный университет 656049, г. Барнаул, пр-т Ленина, 61 тел.: +7-983-381-3453
Krupochkin Evgeny Petrovich Ph. D. in Geography, docent Altai State University 61, Lenin prosp., Barnaul, 656049 tel.: +7-983-381-3453
Батуев Александр Раднажапович, доктор географических наук, профессор Институт географии им. В. Б. Сочавы СО РАН
664033, г. Иркутск, ул. Улан-Баторская, 1 тел.: (3952) 42-27-19 Иркутский государственный университет 664003, Иркутск, ул. К. Маркса, 1 тел.: (3952)52-10-71
Batuev Alexander Radnazapovich, Doctor of Geography, professor, V. B. Sochava Institute of Geography SB RAS
1, Ulan-Batorskaya Str., Irkutsk, 664033 tel.: (3952) 42-27-19; The Irkutsk State University 1, K. Marx str., Irkutsk, 664003 tel.: (3952) 52-10-71
