Программа для построения экологических карт

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

Отвлекаясь от сути, отмечу, что несмотря на всю несуразность сложения двух величин в килограммах и метрах, в зависимости от решаемой задачи такая операция может быть принята как вполне корретная. Например, сравниваются два набора измерений (вес1, длина1) и (вес2, длина2), то в качестве интегральной характеристики может быть предложено max(вес1, вес2) - min(вес1, вес2) + max(длина1, длина2) - min(длина1, длина2). Данная величина всегда неотрицательна и имеет одну важную интерпретацию — она обращается в нуль тогда и только тогда, когда измерения в обоих случаях (1 и 2) совпадают. С другой строны, она не имеет физического смысла — поскольку в ее формировании участвуют два разноразмерных показателя. Тем не менее, если нашей целью является отслеживание равенства наборов (вес1, длина1) и (вес2, длина2), ее использование вполне правомерно — она обращается в нуль только при их равенстве. Однако, гораздо приятней иметь величину, имеющую еще и вполне ясную физическую интерпретацию.

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

Что подразумевает термин —  "экологические карты"?  Каждый словарь трактует это понятие по своему. Достаточно общим, является следующее определение (из "Географического словаря"): экологическая карта отражает взаимодействие живых организмов (в т. ч. людей) со средой; в более широком смысле – взаимодействие социально-экономических и природных геосистем.

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

Алгоритм построения экологических карт (Б.И. Сёмкин)

На основе данных обучающей выборки, для всевозможных пар факторов, рассчитываются матрицы взаимной встречаемости. Если размерность факторного пространства \(m\), то количество вычисляемых матриц встречаемости равно \(m(m-1)/2\). Размерность каждой из матриц определяется количеством градаций соответствующих ей пары факторов. Элементами матриц являются частоты (эмпирические оценки условных вероятностей) взаимной встречаемости факторов, найденные по обучающей выборке.

Далее, каждой картируемой точке \(x\) с учетом вычисленного по обучающей выборке массива матриц взаимной встречаемости сопоставляется матрица \(S_{m\times m}(x)\). Элементами этой матрицы являются значения функции сопряженности \(f(a,b,c)\), вычисленные на значениях элементов матриц взаимной встречаемости. Поясним это следующим примером.

Пусть рассматривается трехфакторное пространство с факторами \(\Phi_1, \Phi_2, \Phi_3\), каждый из которых имеет соответственно \(s_1,s_2\) и \(s_3\) градаций. Таким образом, общее число матриц встречаемости, вычисленное по некоторой обучающей выборке, будет равно $3\cdot2/2=3$. Обозначим эти матрицы для определенности \(M(\Phi_1, \Phi_2)\), \(M(\Phi_2, \Phi_3)\) и \(M(\Phi_1,\Phi_3)\), т.е. \(M(\Phi_i, \Phi_j)\) матрица взаимной встречаемости градаций факторов \(\Phi_i\) и \(\Phi_j\) размерностью \(\dim M(\Phi_i, \Phi_j)=s_i\times s_j\).

Если картируемой точке \(x\) соответствует набор градаций факторов \((\Phi_1(x),\Phi_2(x),\Phi_3(x))\), то соответствующая матрица \(S_{3\times 3}(x)\) будет определяться следующим образом: $$
S_{3 \times 3}(x)  = \begin{array}{*{20}c}
   {} & {\Phi _1 } & {\Phi _2 } & {\Phi _3 }  \\
   {\Phi _1 } &  0  & {a_1 (x)} & {a_2 (x)}  \\
   {\Phi _2 } & {a_1 (x)} &  0  & {a_3 (x)}  \\
   {\Phi _3 } & {a_2 (x)} & {a_3 (x)} &  0   \\
\end{array}, \,\,\mbox{т.е.}\,\, S_{3 \times 3}(x)=\left(
\begin{array}{ccc}
     0  & {a_1 (x)} & {a_2 (x)}  \\
    {a_1 (x)} &  0  & {a_3 (x)}  \\
    {a_2 (x)} & {a_3 (x)} &  0
\end{array}
\right).
$$
Значение \(a_1(x)\) для пары факторов \(\Phi_1\) и \(\Phi_2\), вычисляется по элементам матрицы встречаемости \(M(\Phi_1,\Phi_2)\). Если \(i_{1,2}\) и \(j_{1,2}\) — номер строки и номер столбца матрицы \(M(\Phi_1,\Phi_2)\) соответствующие градациям \(\Phi_1(x)\) и \(\Phi_2(x)\), то $$
\begin{array}{l} a_1(x)=f(a,b,c),\,\,\mbox{где}\\
a=\sum\limits_{k = 1}^{s_2 } {M_{i_{1,2} k} (\Phi _1 ,\Phi _2 )},\\
b=\sum\limits_{k = 1}^{s_1 } {M_{kj_{1,2} } (\Phi _1 ,\Phi _2 )},\\
c=M_{i_{1,2} j_{1,2}} (\Phi _1 ,\Phi _2 ).
\end{array}
$$
Аналогичным образом определяются значения \(a_2(x)\) и \(a_3(x)\).

Функция сопряженности \(f(a,b,c)\) задается пользователем (по умолчанию в программе используется мера Жаккара, т.е. \(f(a,b,c)=c/(a+b-c)\)). В алгоритме  используются только симметричные меры сопряженности, поэтому матрицы \(S(x)\) являются симметричными по построению. Диагональные элементы матрицы \(S(x)\) не участвуют в дальнейших вычислениях и полагаются равными 0.

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

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

Реализация предложенного алгоритма выполнена в программе ECOMAP, разработанной на базе программного пакета научно-технических рассчетов MatLab.

Скачать программу ECOMAP (412,7 KB)

Для запуска программы ECOMAP необходим установленный программный пакет MATLAB версии 7 и выше.

Диалоговое окно программы ECOMAP


Особенности программы ECOMAP

Загрузка данных осуществляется из специально подготовленного mat-файла. Если исходные данные представлены в виде excel-таблицы, то средствами MatLab ее легко конвертировать в необходимый mat-файл.
Данные mat-файла (здесь, собственно, хранится обучающая выборка) имеют следующую структуру: 
  1. Переменная FactorNames — представляет массив названия факторов понятных пользователю и расположенных в той последовательности, в которой они даны в таблице Table.
  2. Переменная FactorTypes —  массив, описывающий возможные значения для категориальных переменных. Здесь же задается тип переменных. Порядок должен соответствовать порядку в переменной FactorNames. В каждой колонке указывается ее тип: 'cat' — категориальный, 'num' — числовой.
  3. Таблица данных Table. Содержит данные обучающей выборки. Ее колонки соответствуют порядку в переменных FactorNames, FactorTypes.

Процесс построения экологической карты в программе ECOMAP состоит из следующих этапов:

  • Верхнетреугольный блок (без главной диагонали) матриц факторных отношений \(S(x)\) вытягивается в вектор размерности \(m(m-1)/2\), где \(m\) — размерность факторного пространства. Данная операция реализуется для каждой матрицы факторных отношений; их общее число определяется количеством картируемых точек (\(N\)). В результате формируется матрица размерности \(m(m-1)/2\times N\), столбцы которой являются векторами, подлежащими кластеризации.
  • Для всевозможных пар полученных векторов вычисляется матрица взаимных расстояний в соответствии с выбранной (или непосредственно определенной) пользователем функцией расстояния между векторами \(x\) и \(y\): \(d(x,y)\). В программе реализовано 8 вариантов  функций \(d(x,y)\); по умолчанию используется евклидово расстояние, т.е. \(d(x,y)=((x-y)^T(x-y))^{1/2}\). При необходимости пользователь может самостоятельно задать вид функции расстояния \(d(x,y)\), предварительно выбрав во вкладке "Выберите функцию расстояния" опцию "Задано пользователем". Размерность матрицы взаимных расстояний (равна \(N(N-1)/2\)) определяется количеством \((N)\) картируемых точек.
  • В соответствии с определенным пользователем критерием кластеризации формируется иерархическое дерево, число листов которого равно \(N\).
  • По заданному пользователем количеству кластеров (\(P\)) (по умолчанию в программе \(P=3\)) проводится кластеризация векторов: каждому вектору назначается номер кластера, которому он принадлежит. Процедура кластеризации представляется собой последовательное перемещение от корня иерархического дерева, до момента выделения \(P\) или более кластеров.

Функция расстояния \(d(x,y)\). Пользователь может самостоятельно определить функцию расстояния \(d(x,y)\) или воспользоваться 8 вариантами ее реализации в программе. Следует иметь ввиду, что при самостоятельном определении функции \(d(x,y)\) время вычисления матрицы взаимных расстояний, а следовательно и время всего процесса кластеризации, существенно возрастет. Учитывая, что размерность этой матрицы велика и
определяется количеством картируемых точек, рекомендуется использовать предусмотренные в программе варианты функций расстояния.

При пользовательском определении функции расстояния следует записывать только тело функции, писать \(d(x,y)=\) не следует. При этом полагается, что переменные \(x\) и \(y\) являются векторами и операции с ними записываются как операции над векторами. Например, если мы хотим определить пользовательскую функцию расстояния \(d(x,y)=\sqrt{\sum\limits_i(x_i-y_i)^2}\), т.е. евклидово расстояние между векторами (на практике такой необходимости не возникнет, поскольку евклидово расстояние является штатной функцией расстояния программы ECOMAP), то в поле "Функция кластеризации \(d(x,y)=\)" необходимо записать: $$
\begin{array}{l}
\textsf{sqrt}((x-y)'*(x-y)),\,\,\mbox{или}\\
((x-y)'*(x-y))^{\wedge} 0.5,\,\,\mbox{или}\\
((x-y)'*(x-y))^{\wedge} (1/2),\\
\end{array}
$$
где \(\textsf{sqrt}\) — операция извлечения квадратного корня, "'" — символ транспонирования вектора (в соответствии с синтаксисом языка MatLab), "*"\, — операция умножения матриц (в данном случае вектор-строка умножается на вектор-столбец), "\(^{\wedge}\)" — операция возведения в степень.

В программе ECOMAP предусмотрены следующие штатные функции \(d(x,y)\) (они являются частью встроенной в MatLab возможности кластерного анализа):

  • Евклидово расстояние: \(d(x,y)=\sqrt{\sum\limits_i(x_i-y_i)^{2^{\phantom{Q}}}}\).
  • City-block расстоняние:  \(d(x,y)=\sum\limits_i|x_i-y_i|\).
  • Расстояние Минковского: \(d(x,y)=\left(\sum\limits_i(x_i-y_i)^p\right)^{1/p}\). Параметр \(p\) определяется пользователем по запросу в процессе работы программы. Следует иметь ввиду, что при \(p=2\) и \(p=1\) расстояние Минковского отождествляется соответственно с евклидовым и City-block расстояниями.
  • Расстояние Махаланобиса: \(d(x,y)=(x-y)^TV^{-1}(x-y)\), где \(V\) — выборочная оценка ковариационной матрицы.
  • Расстояние Хэмминга: \(d(x,y)\) — отношение количества неравных друг другу компонент векторов \(x\) и \(y\) к размерности вектора \(x\) (или \(y\)).
  • 1-Jaccard coef.:  \(d(x,y)\) — отношение количества неравных нулю и друг другу компонент векторов \(x\), \(y\) к общему числу ненулевых компонент в обоих векторах \(x\) и \(y\).
  • 1-Cosine coef.:  \(d(x,y)=1-\cos(\angle x,y)\), где \(\cos(\angle x,y)=\dfrac{x^Ty}{\sqrt{(x^Tx)(y^Ty)}}\).
  • Расстояние Чебышева: \(d(x,y)=\mathop {\max }\limits_i |x_i  - y_i |\).

При построении иерархического дерева используется функция расстояния между кластерами, определенная на элементах матрицы всевозможных
парных расстояний. Пользователь может выбрать четыре способа вычисления расстояния (не путать с расстоянием между векторами, определяемым функцией \(d(x,y)\)) между кластерами (критерия кластеризации), каждому способу будет соответствовать свой исход процесса кластеризации.

Для пояснения методов вычисления расстояния между кластерами введем обозначения:

  • \(X_{ri}\), \(Y_{sj}\) — i-ый и j-ый вектора, принадлежащие кластерам \(r\) и \(s\) соответственно.
  • \(n_r\), \(n_s\) — количество векторов в кластере \(r\) и \(s\).
  • \(D(r,s)\) — расстояние между кластерами \(r\) и \(s\).

В зависимости от выбора пользователя расстояние \(D(r,s)\) определяется следующим образом:

  • Минимальное расстояние: \(D(r,s)=\mathop {\min }\limits_{i,j} d(X_{ri} ,X_{sj} ),\,i=1\ldots n_r,\,j=1\ldots n_s\).
  • Максимальное расстояние:  \(D(r,s)=\mathop {\max }\limits_{i,j} d(X_{ri} ,X_{sj} ),\,i=1\ldots n_r,\,j=1\ldots n_s\).
  • Расстояние Евклида: \(D(r,s)=((\overline{X}_r-\overline{X}_s)^T(\overline{X}_r-\overline{X}_s))^{1/2}\), где \(\overline{X}_r=\dfrac{1}{n_r}\sum\limits_iX_{ri}\) и \(\overline{X}_s=\dfrac{1}{n_s}\sum\limits_jX_{sj}\).
  • Среднее растояние: \(D(r,s)=\dfrac{1}{n_r n_s}\sum\limits_{i,j}d(X_{ri},X_{sj})\).

После того, как  критерий кластеризации  выбран (определено расстояние \(D(r,s)\)), осуществляется процесс построения иерархического бинарного дерева, которое хранится в памяти ЭВМ в виде матрицы \(Z_{N-1\times3}\) (\(N\) — количество классифицируемых объектов). Формирование матрицы происходит по следующему принципу (изложен в документации к инструментарию Statistics Toolbox для среды MatLab).


Пример работы программы

Построим карту кластеров сходных факторных отношений для территории учебно-опытного лесхоза "Дальневосточный". Имеющийся массив данных — 1885 описаний (распределенных регулярно с шагом 500 м по всей территории лехоза) с географической привязкой к системе координат WGS-84. Рассмотрим ситуацию, когда картируемая территория, и территория, которой соответствует обучающая выборка, совпадают. При построении карты кластеров будем учитывать только факторы геоморфологической группы — высоту над уровнем моря (13 градаций), экспозицию склона (11 градаций) и крутизну склона (7 градаций).

Примем следующие сокращенные обозначения факторов: \(\Phi _1\)="Высота над уровнем моря", \(\Phi _2\)="Экспозиция склона", \(\Phi_3\)="Крутизна склона".

Факторы имеют следующие градации (целочисленный номер градации соответствует номеру столбца или строки матрицы взаимных встречаемостей): фактор \(\Phi_1\): 1 — 0–100 м,   2 — 101–150 м, 3 — 151–200 м, 4 — 201–250 м, 5 — 251–300 м, 6 — 301–350 м, 7 — 351–400 м, 8 — 401–450 м, 9 — 451–500 м, 10 — 501–550 м, 11 — 551–600 м, 12 — 601–650 м, 13 — 651–700 м; фактор \(\Phi_2\): 1 — равнина,   2 — пойма, 3 — север, 4 — северо-восток, 5 — восток, 6 — юго-восток, 7 — юг, 8 — юго-запад, 9 — запад, 10 — северо-запад, 11 — водороздел; фактор \(\Phi_3\): 1 — менее \(1^{\circ}\), 2 — от
\(1^{\circ}\) до \(6^{\circ}\), 3 — от \(6^{\circ}\) до \(10^{\circ}\), 4 — от \(10^{\circ}\) до \(15^{\circ}\), 5 — от \(15^{\circ}\) до \(20^{\circ}\), 6 — от  \(20^{\circ}\) до \(25^{\circ}\), 7 — более \(25^{\circ}\).

Найденные по обучающей выборке матрицы взаимной встречаемости градаций факторов имеют вид: $$
M(\Phi_1,\Phi_2)=\left[ \begin {array}{ccccccccccc}
10&16&0&3&0&0&6&4&0&0&0\\10&44&36&13&10&14&23&19&13&22&1\\20&38&60&28&14&24&40&32&21&28&0\\14&34&49&49&19&29&36&45&22&35&0\\7&19&71&41&19&25&51&34&25&41&0\\13&8&56&38&22&17&54&47&27&31&0\\2&12&21&23&15&14&31&23&22&30&8\\0&6&16&7&9&7&8&13&9&13&9\\0&2&6&5&0&2&4&1&6&8&10\\0&0&1&0&0&0&2&0&2&4&0\\0&0&0&0&0&0&1&0&0&1&0\\0&0&1&0&0&0&0&1&1&1&0\\0&0&0&1&0&0&0&0&0&0&0\end
{array} \right]
$$ $$
M(\Phi_2,\Phi_3)=\left[ \begin {array}{ccccccc}
21&55&0&0&0&0&0\\9&90&42&29&6&3&0\\0&108&112&69&22&3&3\\0&74&73&55&6&0&0\\0&42&45&20&1&0&0\\0&56&39&30&6&1&0\\0&72&87&63&26&7&1\\1&74&82&47&11&3&1\\0&76&42&27&2&1&0\\0&97&65&42&9&1&0\\8&20&0&0&0&0&0\end
{array} \right]
$$$$
M(\Phi_1,\Phi_3)=\left[ \begin {array}{ccccccc}
8&26&3&1&1&0&0\\5&131&45&15&7&2&0\\5&153&84&47&12&3&1\\6&117&126&64&12&5&2\\0&81&132&96&22&2&0\\5&88&99&96&18&5&2\\6&91&52&43&7&2&0\\2&43&35&14&3&0&0\\2&31&5&3&3&0&0\\0&2&3&2&2&0&0\\0&0&1&0&1&0&0\\0&0&2&1&1&0&0\\0&1&0&0&0&0&0\end
{array} \right]
$$

Приведем в виде матрицы (\(F\)) фрагмент (первые 10 записей) используемого массива описаний. При учете всех 1885 описаний данная матрица будет иметь 1885 строк. Каждая строка матрицы \(F\) соответствует отдельному описанию и содержит градации факторов \(\Phi_1\), \(\Phi_2\) и \(\Phi_3\) (нижний индекс соответствует номеру столбца \(F\)). $$
F=\left[ \begin {array}{ccc}
3&7&3\\3&10&2\\2&8&2\\3&2&2\\3&4&2\\2&3&2\\2&2&3\\3&8&3\\4&4&3\\3&4&2\end
{array} \right]
$$
Для построения карты кластеров используем следующие дополнительные условия эксперимента: мера сопряженности факторов — коэффициент
Жаккара, функция расстояния между векторами \(x\) и \(y\): \(d(x,y)\) — "CityBlock", межклассовое расстояние — "Максимальное расстояние", количество кластеров равно 4.
Соответствующая матрице \(F\) матрица парных сопряженностей факторов (\(Q\)) будет иметь вид: $$
Q=\left[ \begin {array}{cccccccccc} {\frac {40}{521}}&{\frac
{28}{491}}&{\frac {19}{405}}&{\frac {19}{223}}&{\frac
{28}{485}}&{\frac {2}{27}}&{\frac {11}{85}}&{\frac {8}{123}}&{\frac
{49}{491}}&{\frac {28}{485}}\\{\frac
{21}{202}}&{\frac {153}{916}}&{\frac {131}{838}}&{\frac
{153}{916}}&{\frac {153}{916}}&{\frac {131}{838}}&{\frac
{5}{83}}&{\frac {21}{202}}&{\frac {126}{793}}&{\frac
{153}{916}}\\{\frac {29}{252}}&{\frac
{97}{881}}&{\frac {74}{909}}&{\frac {90}{853}}&{\frac
{37}{449}}&{\frac {108}{973}}&{\frac {21}{362}}&{\frac
{41}{362}}&{\frac {73}{722}}&{\frac {37}{449}}\end {array} \right]
$$

Первая строка матрицы \(Q\) содержит коэффициенты  Жаккара для факторов \(\Phi_1\) и \(\Phi_2\), вторая строка — для факторов \(\Phi_1\) и \(\Phi_3\), третья строка — для факторов \(\Phi_2\) и \(\Phi_3\). Столбцы матрицы \(Q\) соответствуют описаниям, представленным строками матрицы \(F\). Например, в пятой строке матрицы \(F\) определены градации \((\Phi_1,\Phi_2,\Phi_3)=(3,4,2)\). Тогда, соответственно, в пятом столбце матрицы \(Q\) будет \(\frac{28}{485}\) — коэффициент Жаккара для пары градаций \((\Phi_1,\Phi_2)=(3,4)\), \(\frac{153}{916}\) — коэффициент Жаккара для пары градаций \((\Phi_1,\Phi_3)=(3,2)\), \(\frac{37}{449}\) — коэффициент Жаккара для пары градаций \((\Phi_2,\Phi_3)=(4,2)\).

В соответствии с выбранными насройками процедуры кластеризации столбцам матрицы \(Q\) будут присвоены номера кластеров $$\left[
\begin {array}{cccccccccc} 2&2&1&2&2&2&4&2&2&2\end {array}
\right].$$

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

Сглаживание (последний рисунок) выполнено методом k-ближайших соседей при  k=40. Интерфейс программы ECOMAP позволяет изменять параметр сглаживания, обеспечивая подбор ситуации с наилучшим визуальным разбиением области на кластеры.

Экологическая карта лесхоза Дальневосточный Экологическая карта лесхоза Дальневосточный Экологическая карта лесхоза Дальневосточный

Публикации

 

blog comments powered by Disqus