Подход к моделированию изменений земной поверхности с использованием клеточных автоматов

стические характеристики этого режима течения, безусловно, должны иметь отличия от развитого турбулентного потока. Косвенным подтверждением изложенного является то, что при развитом турбулентном потоке в трубке (в нашем случае при числах Рейнольдса большее 2500, а для более крупных пузырей, исследованных в работе [1], при числах Ке>5000) одно- и двухфазные мелкопузырьковые и МГС потоки практически не отличимы друг от дру-

га и описываются в безразмерных координатах единой кривой (в данном случае законом Блазиуса).

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

СПИСОК ЛИТЕРАТУРЫ

1. Бурдуков А.П., Валукина Н.В., Накоряков В.Е. Особенности течения газожидкостной пузырьковой смеси при малых числах Рейнольдса // Журнал прикладной механики и технической физики. - 1975. - № 4. - С. 137-141.

2. Валукина Н.В., Козьменко Б.К., Кашинский О.Н. Характеристики мелкопузырьковой смеси при течении в вертикальной трубе // Инженерно-физический журнал. - 1979. - Т. 36. -№ 4. - С. 695-699.

3. Bosheniatov B.V. Measurement of parameters of two-phase flow with very small gas bubbles in horizontal tubes // Intern. Conf. on the Methods of Aerophys. Research: Proc. Pt. IV. - Novosibirsk, 2004. - P. 81-86.

4. Бошенятов Б.В. Гидродинамика микропузырьковых газожидкостных сред // Известия Известия Томского политехнического университета. - 2005. - Т. 308. - № 6. - С. 156-160.

5. Прантль Л., Титьенс О. Гидро- и аэромеханика. - М.-Л.: ОНТИ НКТП СССР, 1935. - Т. II. - 283 с.

УДК 504.064

ПОДХОД К МОДЕЛИРОВАНИЮ ИЗМЕНЕНИЙ ЗЕМНОЙ ПОВЕРХНОСТИ С ИСПОЛЬЗОВАНИЕМ КЛЕТОЧНЫХ АВТОМАТОВ

А.В. Замятин, Н.Г. Марков

Томский политехнический университет E-mail: zamyatin@tpu.ru

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

Введение

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

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

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

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

Теоретические основы моделирования изменений

земной поверхности

Один из наиболее перспективных и рациональных способов моделирования изменений земной поверхности основан на использовании статистического аппарата марковских цепей [1]. В этом случае вероятностная информация об изменениях земной поверхности содержится в матрице вероятностей переходов Р=[р], каждый элемент которой является вероятностью перехода типа земной поверхности а>1 в тип а, где у = 1,...,М, М- число типов земной поверхности на изучаемой территории, выявленное на этапе интерпретации АИ. Количественная информация об изменениях содержится в матрице фактических переходов территорий Мфакг=[т,:ф] и матрице ожидаемых переходов территорий Можид=[т™]. Матрица Мфакт содержит количество элементов изображения (пикселей), перешедших из типа земной поверхности в тип а,, а Можид - количество пикселей изображения, которые должны поменять свое состояние (значение) в результате процесса моделирования изменений земной поверхности.

Множество КА, с использованием которых осуществляется моделирование, формируется при последовательном сканировании исходного изображения.

При этом каждый КА можно представить в виде квадратной матрицы MCA=[cii] порядка d. Значение центрального элемента c'kh каждой такой матрицы зависит некоторым образом от значений всех элементов этой матрицы c'k=fcA(cuA2,---,Ckh,---,Cdd). Полученное при этом значение c'kh является элементом результирующего изображения. При традиционном подходе поиск функции fCA чаще всего основан на использовании какого-либо простого правила (например, используется принцип мажоритарности), найденная функция используется для всех КА моделируемого изображения. В то же время в предлагаемом подходе к моделированию изменений земной поверхности применяется такой способ формирования функции fCA, при котором для каждого КА используется три различных вероятностных компонента (рфжг, р/6", р.""01), что позволяет учесть особенности формирования правил взаимодействия каждого используемого КА. При этом каждый из компонентов представляет собой соответствующие вероятности р. перехода типа земной поверхности ю; в тип ю. Рассмотрим более подробно особенности расчета каждой такой вероятности.

Компонент, учитывающий факторы влияния. В настоящее время при моделировании изменений заселенных территорий (с помощью моделей, известных как Urban Growth Models) широко исполь-

Вероятностные карты факторов влияния

Слои исходной тематической карты

Фактор влияния F

Фжгор влияния 1 2 jr

Вероятностные карты факторов влияния

Вероятность для типа M

Вероятности переходов

Слои исходной тематической карты

Тип поверхности M

Тип поверхности 1

Рис. 1. Использование нейросети для формализации факторов влияния: а) обучение нейросети; б) получение отклика

2

зуется подход, основанный на применении так называемых карт соответствия [1]. Эти карты несут информацию о вероятности заселения каждой точки земной поверхности изучаемого фрагмента территории в зависимости от таких факторов влияния как расстояние до автомагистралей, расстояние до крупных торговых центров, расстояние до источников питьевой воды, величина уклона местности и т.п. Для применения подобного подхода при моделировании изменений земной поверхности необходимо построение таких карт соответствия для всех типов земной поверхности, представленных на тематической карте. Учитывая возможные сложные взаимозависимости между подобными факторами, задача построения подобных карт соответствия для каждого типа земной поверхности представляется сложной и трудно формализуемой.

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

В общем случае применяемый МСП имеет К=М¥ входных и М выходных нейронов, где М -число типов земной поверхности на разновременных тематических картах, ¥ - число факторов влияния. Процесс получения необходимой вероятностной информации с помощью нейросети, учитывающей выявленные факторы влияния, заключается в том, что каждый аксон /-ого выходного нейрона интерпретируется как вероятность перехода анализируемого пикселя типа щ исходного изображения в тип щ (р,факт).

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

окрестности, совпадающей с элементами текущего КА, вычисляется вероятность рувер=«/ру, /=1,2,...,т, где т - количество типов земной поверхности в анализируемой окрестности, а пу - количество элементов типа щ/ в анализируемой окрестности.

Компонент, учитывающий пространственную составляющую. Третий компонент формирования правил функционирования КА базируется на использовании пространственных характеристик каждого из типов земной поверхности, основанных на вычислении фактора насыщенности, предложенного в [3]. Использование факторов насыщенности происходит следующим образом. Для каждого щ формируется соответствующий вектор Р™^/!/;,.../}, г=1,...,М, где М - количество типов земной поверхности. Этот вектор Г,™ несет информацию о насыщенности того или иного типа земной поверхности на всем изображении в целом. Затем, для каждой точки изображения рассчитывается локальный вектор насыщенности р0", после чего вероятность рупрст определяется как р/ПРст~1/ ^(р нас, р л°к), где й - евклидово расстояние между двумя этими векторами насыщенности.

После определения каждого из вероятностных компонентов, результирующую вероятность перехода каждого из элементов исходного изображения в выбранной окрестности текущего КА можно представить как р/"н=/сА(р)~р/ер+рг,прс1+р/акг.

Немаловажной необходимой деталью при проведении моделирования изменений земной поверхности является не только определение вероятности перехода руф"н, но и порядок изменения элементов результирующего изображения. Этот порядок в рамках предлагаемого подхода определяется следующим образом. Самыми первыми изменяются элементы изображения, результирующая вероятность для которых руфин имеет наибольшее значение, а самыми последними - те элементы изображения, для которых эта вероятность наименьшая. Предварительные результаты исследования эффективности применения такого ранжирования показывают, что оно позволяет получать значительно более адекватные результаты моделирования изменений земной поверхности [4].

Исследование предлагаемого подхода с использованием модельных разновременных изображений

Исследование эффективности предлагаемого подхода с точки зрения точности было проведено с использованием модельных разновременных тематических карт. Модельные данные представляли собой 30 различных комбинаций из трех модельных изображений, имитирующих разновременные АИ, на которых изменение границ и площадей типов земной поверхности происходит аналогично изменениям, наблюдаемым на реальной земной поверхности. Моделирование осуществляется по двум изображениям, а точность оценивается по третьему изображению с помощью матрицы пере-путывания и интегрального критерия точности -каппа индекса согласия (КИС) [4]. Для получения

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

1

0,9 0,8 0,7 0,6 О 0,5 ^ 0,4 0,3 0,2 0,1 0

Рис. 2.

1

3

4

1,2->3 2,3->4 4,3->2 3,2->1

Пары разновременных изображений

Точность прогнозных карт, полученных с использованием различных составляющих в правилах функционирования КА: 1) р, 2) рвер, 3) ре и р¡прст, 4) точность карты, полученной с помощью модуля СА_МагО

В отличие от вероятностных характеристик pi¡ev и Pi|,vcт, рассчитываемых исключительно по разновременным тематическим картам, вероятность p¡fш определяется с использованием дополнительной априорной информации о процессах, происходящих на исследуемом фрагменте земной поверхности. Для простоты, исследования проводились при условии отсутствия влияния априорной информации, поэтому вероятность p¡fш не учитывалась. Отметим также, что проводились сравнительные исследования эффективности построения прогнозных карт, выполненные с помощью существующего модуля моделирования изменений земной поверхности CA_Markov, представленного в растровой геоинформационной системе (ГИС) Ип8132 (США) [7]. На рис. 2 в качестве примера представлена часть результатов проведенных исследований. На оси абсцисс показаны различные виды комбинаций четырех модельных изображений, имитирующих разновременные АИ. Так, запись 1,2—>3 означает, что расчет параметров и построение прогнозной карты проводилось на основе модельных изображений № 1 и 2, а точность полученной прогнозной карты оценивалась по изображению № 3.

Результаты проведенных исследований показали высокую значимость вероятностных характеристик piíщ> и p¡j,vcт при формировании правил функционирования КА, а также то, что созданное программное обеспечение (ПО), реализующее предложенный подход, позволяет строить значительно более точные прогнозные карты, чем модуль CA_Markov ГИС Шя32.

Применение предлагаемого подхода для анализа

динамики земной поверхности фрагмента территории Ханты-Мансийского автономного округа

В ГИС лаборатории ТПУ были программно реализованы предлагаемый подход к моделированию изменений земной поверхности и соответствующие алгоритмы, а также необходимые функции предварительной обработки и автоматизированной интерпретации АИ [8,9]. Апробация предлагаемого подхода проводилась при решении практически

важной задачи анализа динамики земной поверхности фрагмента территории Ханты-Мансийского автономного округа (ХМАО). Решение задачи позволит выявить тенденции в изменении земной поверхности в зоне магистрального нефтепровода, проходящего через исследуемую территорию.

Для решения указанной задачи были использованы снимки системы ДЗЗ Landsat (ETM+), полученные Центром приема и обработки данных ДЗЗ Югорского НИИ информационных технологий на 10 августа 1999 г., 6 сентября 2000 г и 3 июля 2002 г. Все снимки имели уровень обработки 1G и были переведены в проекцию UTM на эллипсоиде WGS84 по данным орбитальной привязки. Учитывая точность проведенной геометрической коррекции и разрешение исходных АИ, рабочий масштаб тематических карт будет соответствовать картам масштаба 1:200 000.

Для анализа эффективности предлагаемого подхода и созданного ПО в целом в качестве базовых АИ для проведения интерпретации и моделирования использовались только АИ 1999 и 2000 гг. В этом случае АИ 2002 г. можно принимать в качестве тестового (эталонного) для оценки качества результатов интерпретации и моделирования в виде прогнозной карты 2002 г.

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

Этап 1. Из имеющихся космических снимков, покрывающих территорию ХМАО, была составлена композиция, которая позволила выявить такие фрагменты для всех трех космических снимков, которые не закрыты облачным покровом и содержат общий фрагмент земной поверхности на исследуемой территории. Учитывая небольшие вариации космических снимков по фенологической фазе, имеется возможность построения обучающих выборок для одних и тех же типов земной поверхности, представленных на изучаемой территории в основном различными классами растительности.

Этап 2. Выбранные фрагменты космических снимков, которые являются общими для всех разновременных АИ, была подвергнуты неконтролируемой классификации по методу ISODATA. Для этого использовалось ПО системы ER Mapper 5.5, а количество кластеров при неконтролируемой классификации было экспериментально принято равным десяти. С учетом результатов этой классификации, а также имеющихся данных наземных наблюдений, для каждого из космических снимков были построены обучающие выборки для таких типов земной поверхности как пойменные участки, кустарники, мшистые болота и вода. Затем все космические снимки были подвергнуты тематической

2

интерпретации с помощью созданного ПО. Это позволило получить разновременные растровые тематические карты, на которых представлены тематические классы, соответствующие построенным обучающим выборкам (рис. 3).

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

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

Таблица 1. Матрица вероятностей переходов

Тип Пойменные участки Кустарники Вода Мшистые болота

Пойменные участки 0,375998 0,221505 0,001555 0,400942

Кустарники 0,101386 0,843558 0,000232 0,054824

Вода 0,094620 0,005566 0,812616 0,087199

Мшистые болота 0,043285 0,017094 0,002255 0,937366

Расчет основывается на тематических картах, полученных с космических снимков 1999 и 2000 гг. Рассчитывается также матрица ожидаемых переходов (табл. 3), необходимая для получения тематической карты на следующий момент времени -2001 г. Отметим, что матрица ожидаемых переходов и матрица фактических переходов формируются из пикселей соответствующих тематических карт (изображений), поэтому в табл. 2 и 3 приводятся элементы матрицы в пикселях.

Таблица 2. Матрица фактических переходов

Тип Пойменные участки Кустарники Вода Мшистые болота

Пойменные участки 130818 77247 554 139822

Кустарники 44624 375793 128 23614

Вода 1073 58 8560 867

Мшистые болота 38144 15356 2046 834497

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

Этап 5. После того, как вся необходимая вероятностная и другая количественная информация получена, проводится непосредственно моделирование изменений земной поверхности и получение прогнозной карты 2002 г. (рис. 4, б). В качестве промежуточного результата моделирования, позво-

£ У ^ ^ ¿Г .: :': ^

йЩЗЭ ЫМ-?

„«^ ь  -ч , V.

а б

Рис. 4. Фрагменты прогнозных карт, полученные с помощью ПО системы: а) 2001 г., б) 2002 г.

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

| [пойменные участки

ляющего оценить интересующие нас изменения земной поверхности 2001 г., получена прогнозная тематическая карта 2001 г. (рис. 4, а).

Таблица 3. Матрица ожидаемых переходов

Тип Пойменные участки Кустарники Вода Мшистые болота

Пойменные участки 62924 96774 1000 187743

Кустарники 55534 328251 337 60038

Вода 1315 348 6944 1949

Мшистые болота 51835 35857 3642 798709

Этап 6. Для того чтобы осуществить дополнительную (сравнительную) оценку полученной с помощью ПО системы прогнозной карты на 2002 г., была сформирована еще одна прогнозная карта на 2002 г. Эта карта была получена с использованием аналогичных входных данных (разновременных тематических карт), что и в первом случае, но с помощью модуля CA_Markov, представленного в ГИС Мп-8132. Также во втором случае в качестве дополнительного, промежуточного результата моделирования была получена прогнозная карта на 2001 г. Заметим, что при получении прогнозных карт использовались параметры модуля CA_Markov, задаваемые по умолчанию, а размер анализируемой окрестности КА был аналогичен размеру в первом случае и составлял 5x5.

Таблица 4. Матрица значений насыщенности исходной тематической карты на 2000 г.

Тип Пойменные участки Кустарники Вода Мшистые болота

Пойменные участки 0,293475 0,62907 3,618187 -0,181239

Кустарники -1,202167 1,113413 2,90909 -2,136832

Вода -5,148586 -4,03138 4,551272 -5,124147

Мшистые болота -0,70012 -0,823881 3,101082 0,596449

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

По результатам пространственного анализа тематической карты на 2000 г., можно сделать вывод, что пойменные участки исследуемого фрагмента земной поверхности занимали общую площадь около 257,6 км2, а к 2002 г. по данным сформированной прогнозной карты - около 205,9 км2, то есть общая площадь пойменных участков уменьшилась приблизительно на 51,6 км2. В то же время мшистые болота на 2000 г. занимали общую площадь того же фрагмента земной поверхности около 1198,5 км2, а к 2002 г. - около 1258,1 км2, то есть можно говорить об увеличении площади данного типа земной поверхности на 59,6 км2. Водные объекты исследуемой земной поверхности, а также границы кустарников практически не претерпели никаких изменений, что хорошо видно из табл. 2 и 3 и прогнозных карт (рис. 4).

Заключение

Предлагаемый подход к моделированию изменений земной поверхности, реализованный в ПО лаборатории ГИС ТПУ, использует в своей основе вероятностное формирование правил функционирования КА. Предварительные результаты исследований, полученные с использованием модельных разновременных изображений, показали перспективность предложенного подхода. Результаты апробации при решении задачи анализа динамики земной поверхности фрагмента территории ХМАО, подтвердили результаты исследований, полученные на модельных изображениях, и показали основные тенденции изменений ландшафтных типов изучаемой территории.

Исследования выполнены при поддержке РФФИ (грант № 03-07-90124) и Минобрнауки РФ (грант № 432).

СПИСОК ЛИТЕРАТУРЫ

1. Baker W.L. A review of models of landscape change // Landscape Ecology. - 1989. - № 2. - P. 111-133.

2. Yeh A.G.-O., Li X. Simulation of development alternatives using neural networks, cellular automata, and GIS for urban planning // Photogrammetric Engineering and Remote Sensing. - 2003. -V. 69. - № 9. - P. 1043-1052.

3. Verburg P.H. e.a. A method to analyse neighborhood characteristics of land use patterns // Computers, Environment and Urban Systems. - 2003. - V. 23. - P. 383-400.

4. Zamyatin A.V., Markov N.G. Designing Forecast Thematic Maps Using Time Series Remotely Sensed Images // ISPRS Proc. of The Intern. Archives of The Photogrammetry, Remote Sensing and Spatial Information Sciences, V. XXXV, part B. - Istanbul, Turkey, 2004. - P. 492-497.

5. Richards J.A., Xiuping Jia. Remote Sensing Digital Image Analysis: An Introduction. -Berlin: Springer, 1999. - P. 363.

6. Clarke K.C., Hoppen S., Gaydos L. A self-modifying cellular automaton model of historical urbanization in the San Francisco Bay area // Environment and Planning B: Planning & Design. - 1997. -№ 24. - P. 247-261.

7. Clark Labs. Geoinformation system Idrisi 32. - http:// www.clar-klabs.org/ IdrisiSoftware/ (01.05.2005).

8. Замятин А.В., Марков Н.Г., Напрюшкин А.А. Адаптивный алгоритм классификации с использованием текстурного анализа для автоматизированной интерпретации аэрокосмических изображений // Исследование Земли из космоса. - 2004. -№ 2. - С. 32-40.

9. Markov N.G., Napryushkin A.A., Zamyatin A.V., Vertinskaya E.V. Adaptive Procedure of RS Images Classification with Use of Extended Feature Space // Proc. of SPIE. - 2003. - V. 4885. - P. 489-500.