V. ПРИКЛАДНЫЕ СЮЖЕТЫ ГИС-МОДЕЛИРОВАНИЯ В ГЕОЭКОЛОГИИ

В реальных городах, разумеется, вместо представленных в Таблица 22.6 "чистых типов", обнаружим множество смешанных, для упорядочений которых необходимо проанализировать как (на самом деле) распределены характеризующие их базовые параметры: этажность, плотность и интенсивность застройки, а также количество и размерность (средняя площадь и периметр) домов в квартале. При этом следует учитывать то обстоятельство, что показатель средней этажности 5 весьма точно отражает действительность для квартала, застроенного двадцатью "пятиэтажками", но другой квартал с таким же средним значением может быть застроен несколькими трех-четырех-, девяти- и шестнадцати этажными домами. Поэтому "разброс" (стандартное отклонение) любого параметра имеет не меньшее значение, чем собственно индекс.

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

Для начала посмотрим, как распределяются кварталы по параметру Средней высоты домов (Hb) и Стандартному отклонению от средней высоты (Hdev). Инструмент ArcMAP10.x Анализ Группирования единственный, на сегодняшний день алгоритм кластерного анализа, "подсказывающей" пользователю оптимальное число классов (при условии что их менее 15), но как мы уже имели возможность убедиться далек от совершенства и пользоваться им лучше с "оглядкой" на другие возможности ГИС-анализа. Алгоритм использования для всех факторов будет аналогичным: первый запуск с максимально возможным числом групп производится для определения "оптимального" числа, второй с указанием рекомендованного числа групп (классов).

Мы "заказали" максимальное число групп - 15 и pdf-отчет, содержащий в последней части график Pseudo F-Statistic с указанием оптимального числа групп. В данном случае это число равно двум, что выглядит не очень правдоподобным для совокупности из более чем пяти тысяч кварталов.

Откроем таблицу с результатами группирования и используя инструменты геостатистического анализа построим гистограмму по полю SS_GROUP (номера групп) Geostatistical Analysis >> Explore Data >> Histogram:

Гистограмма позволяет обнаружить, что из 15-ти "заказанных" кластеров по меньшей мере 7-8 статистически значимы. Следовательно, мы можем повторить группирование для 7 классов.

Для интерпретации результатов группирования используем инструмент Summary Statistics ArcMAP10.x:

Посмотрим каким образом выделились классы: проще всего проводить такой анализ начиная с четко выделяющихся по экстремумам (максимальным или минимальным значениям) таксонов. Мы видим немногочисленный (всего 20 кварталов) класс (6) с максимально высокими домами-небоскребами выше 88 м, при разбросе чуть более 3 м (как мы уже знаем - это высота одного этажа). На противоположном конце шкалы высотности класс (5) - кварталы низкоэтажных (около 12 м) равных по высоте (разброс 1,8 м) домов: таких кварталов в Старой Москве по-прежнему очень много 2389; вероятно, к ним относятся морфотипы исторической застройки, двух-трехэтажные кварталы "немецкой" застройки и четырехэтажные "панельные" и кирпичные "хрущевки" и "брежневки".

Следующие два класса близки и по значениям высоты, и по разбросу: Класс (1) 1053 кварталов со средней высотой домов 16,7 м при разбросе 7,7 сформирован 3-4-хэтажными домами (с отдельными включениями более высоких домов; класс (4) 848 кварталов со средней высотой домов 21,7 м и разбросом 14,2 м - представляет морфотип преобладающих пятиэтажек, опять-таки с отдельными более высокими и более низкими строениями. Класс (2) включает кварталы с домами почти тридцатиметровой высоты и значительной дисперсией: следовательно, увеличение средней высоты сопровождается и увеличением разброса, что свидетельствует о том, что прирост высоты может обеспечиваться немногими строениями. Наконец класс (3) - сформирован очень разными строениями, сюда относятся и небоскребы 50-метровой высоты и обычные пятиэтажки, что определяет максимальный разброс 33 м.

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

Поскольку средняя площадь домов Москвы не связана линейной зависимостью с их высотой (коэффициент детерминации R² = 1,62%) можно ожидать иную кластеризацию для этого параметра. И действительно, Анализ Группирования ArcMAP10.x рекомендует 15 (!) кластеров, однако, гистограмма демонстрирует статистическую малочисленность по крайней мере половины таксонов (Рис. 22.28 )

Рассчитав суммарную статистику мы всегда можем выявить таксоны, представленные очень небольшим числом кварталов. В случае Москвы малочисленные классы при проверке оказались кварталами крупных торгово-развлекательных или культурно-деловых многофункциональных или коммунально-складских центров: они "выскочили" из общего ряда непомерно большими значениями площади строений, занимающих почти все пространство квартала. Таким образом, предварительное группирование может рассматриваться и как дополнительный метод обнаружения искажающих (outlet) значений. Например в модели, рассматривающей, в основном, жилые кварталы подобные "уникумы" лучше удалить. Запустив Анализ Группирования после коррекции данных получаем новое значение оптимального количества кластеров - 9 групп, для интерпретации которых проведем еще раз процедуру Суммарной Статистики.

Картограмма кварталов классифицированных по средней площади домов и разбросу значений демонстрирует преобладание двух классов, относящихся к морфотипу жилых кварталов. В классе 1 со средней площадью 560 м² преобладают небольшие дома - при ширине типового дома "хрущевки" около 12-15 м имеем дома длиной 38-46 м с 4 или 6 подъездами. В классе 7 со средней площадью 978 м² дома крупнее. Стоит ли сохранять редкие классы 4 и 5 (число кварталов соответственно 19 и 18) с очень крупными по площади строениями? Класс 4 сформировался за счет кварталов промзон, автотранспортных предприятий, заводоуправлений с крупными центральными строениями. 5 класс включает кварталы ритейлинговых строений, а также кварталы с немногочисленными, но крупными жилыми комплексами. Следовательно, целесообразно оставить их в нашей классификации.

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

Попытаемся получить комплексную классификации кварталов Москвы по совокупности разных морфологических признаков: выберем и признаки, характеризующие, норму (т.е., средние показатели по кварталу) и разнообразие (дисперсию значений, отклонения от среднего). Заметим попутно, что когда у исследователя в руках широкий набор факторов возникает неизбежный соблазн использовать их "полным списком", но опыт показывает, что такой подход редко бывает продуктивным. Для того чтобы убедиться в этом полезно поупражняться запуская Анализ Группирования несколько раз подряд с последовательно уменьшающимся числом параметров, (скажем, 10,8, 5...) и сравнивать результат.

Как можно убедиться (Рис. 22.32) результаты кластеризации по десяти и пяти факторам различаются (при условии одинакового числа групп - 15) не радикально: для выявления реальной картины дифференциации кварталов Москвы по морфологическим признакам зданий и плотности застройки вполне достаточно пяти факторов:

1. Средняя высота домов в квартале ,
2. Средняя площадь зданий в квартале ,
3. Стандартное отклонение площади всех поверхностей (стен и крыш) зданий в квартале ,
4. Плановая застроенность ("запечатанность): отношение общей площади домов (оснований) к площади квартала ,
5. Закрытость квартала .

Выбор оптимального числа групп - задача, которая ложится на плечи эксперта. В гистограмме частотности мы видим 10 групп, имеющих достаточное распространение - вспомним, что равно такое число классов предлагается в "идеальной" классификации Стюарта и Оке [Stewart, Oke, 2014]. Попытаемся получить аналогичную запустив Grouping Analysis с соответствующими опциями (10 групп, без пространственных ограничений).

Проанализируем итоговый результат кластерного анализа для кварталов Москвы

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

Самая малочисленная (всего 18) Группа (4) собрала кварталы со средней высотой 94 м - что ровно в два раза больше, чем у ближайшей по этому параметру группы (5), довольно значительной средней площадью домов и невысоким разбросом площади поверхностей зданий (что говорит об относительно гомогенности квартала). У кварталов этой группы относительно высокий показатель запечатанности (33%) и средние показатели закрытости - 24%. Произведя выборку по номеру группы SS Group = 4 в выходном слое процедуры Анализ Группирования мы можем просмотреть значения параметров для всех 18 кварталов, что дает возможность судить об их разбросе; таким образом, можно, например, убедиться, что размах высот строений в этом морфотипе составляет от 68 м до 148 м, а стандартное отклонение площади поверхностей для более чем половины кварталов равно нулю, т.е., перед нами варианты кварталов с высотной типовой застройкой. Подложив под выбранные кварталы снимок можно верифицировать это заключение (Рис. 22.36): в составе группы - новейшие жилые кварталы, занятые одним или несколькими высотными жилыми зданиями и бизнесцентрами, (так называемые "стилобаты" или "башни"), например - башня бизнесцентра "Домников" на ул. Маши Порываевой.

Группа (5) также немногочисленная (147) включает кварталы строений повышенной этажности (47 м), относительно небольшой средней площади, небольшой застроенности (24%), открытые (19%). Эта группа - новые кварталы, и поскольку архитекторы склонны экспериментировать с композицией "вертикальных форм", то вся группа более "пестрая" с большим разбросом значений по всем параметрам: по высоте (от 4 м3 до 73 м), по плотности застройки (6-24%), по закрытости (1-57%). Соответственно в кварталах группы (5) размещаются довольно разные объекты - жилые комплексы высоток различной конфигурации, деловые центры, гостиничные комплексы (Рис. 22.37).

Группа (9) собрала кварталы низкоэтажных 13 м строений и является второй по численности (1081 квартал). Два других ярких признака группы - минимальная средняя площадь основания и минимальная площадь поверхностей строений; таким образом, перед нами хорошо узнаваемый морфотип застройки конца 1950-х - начала 1960-х гг. Кроме жилых в группу попали и единичные кварталы сходной морфологии, но иного назначения: общественно-деловые, коммерческие, например - квартал студии Мосфильма (Рис. 22.38).

Близкая по высотности (14 м) группа (1) отличается максимальной численностью (1237 кварталов), чуть большей площадью строений, но значительно большим разбросом общей площади поверхностей (т.е., внутри квартала могут быть строения разной формы) и вдвое большей (33%) плотностью застройки. Группа (1) - это историческая и, зачастую,-ценная застройка с 3-6-этажными домами, многие из которых - памятники архитектуры. Значительная часть полузакрытых этого морфотипа расположена в пределах Садового Кольца (Рис. 22.39); за последние десятилетия внутренне пространство пополнилось современными жилыми и бизнес-комплексами более высокой этажности, выполненными по индивидуальным проектам.

Группа (6) отличается максимальной площадью строений в квартале (8272 м₂) при небольшой высоте (18 м) и низкой распространенности (всего 54 квартала) м высокой (2-й в рейтинге) запечатанности - 51%. Анализ картограммы показывает, что представители этой группы кварталов распространены спорадически, как в пределах Садового Кольца, так и в старых "примагистральных клиньях", единично - в новой застройке по соседству или в окружении группы 3 и группы 8. В центре Москвы морфотип (6) - это "кварталы одного строения": крупные торговые центры (Центральный детский магазин на Лубянке) а также жилые кварталы в стиле "сталинского ампира", в других местах - коммерческие и общественно-деловые кварталы на месте бывших промзон.

Группа (2) выделяется по параметру максимального стандартного отклонения площади поверхностей зданий и сооружений в квартале (18425 м²), при этом она малочисленна (82 квартала), характеризуется средней высотой (Hb_AVG = 24), выше, чем средняя плотностью застройки (37%) и средней закрытостью (22%). Кварталы группы разбросаны дисперсно по всем районам Москвы, встречаются как в центре, так и на периферии и сильно различаются по размерам. Проверка "представителей" группы (Рис. 22.41) подтверждает признак чрезвычайного разнообразия форм и размеров (следствием чего и является разброс площади поверхностей), иными словами это кварталы с самой разной архитектурной "начинкой", в том числе с единственным крупными домами индивидуального проекта ("Дома на Набережной"), советскими "высотками" с центральными башнями и "крыльями", а также - исторические промзоны с цехами в стиле конструктивизма 1930-х, наконец, реже прочих - жилые кварталы уникальной стилистики 1950-х гг.

Группа (10) отличается максимальной запечатанностью (55%) и закрытостью (55%), высотой близкой к минимальной (14 м) и также близкой к минимальной площадью строений. Абсолютное большинство кварталов данного морфотипа (Рис. 22.42) - исторические кварталы Москвы в пределах Садового кольца, с кирпичными строениями в среднем от трех до шести этажей по внешнему периметру (отсюда закрытость), которые изначально возводились по индивидуальным проектам, и, зачастую подвергались неоднократной перестройке. Сюда же попали исторические промышленные зоны с типичными образцами русской промышленной архитектуры конца XIX - начала XX вв.

Группа(3) выделилась наименьшей запечатанностью (17%) открытостью (13%) по внешнему периметру кварталов, при высоком разнообразии площадей поверхностей зданий. Высота домов - средняя, площадь домов - близка к минимальной. Большая часть объектов - крупноареальные кварталы за пределами Садового кольца.

Группа (8) достаточно многочисленна (969), средневысотная, с небольшой средней площадью оснований зданий; характеризуется невысокой запечатанностью (16%) и закрытостью (16%). Данный тип широко распространен за пределами исторического центра Москвы и хорошо представлен на периферии в виде крупных по площади кварталов с самым разным расположением домов внутри: рядовым, "скошенным", краевым вдоль красных линий (Рис. 22.44).

Группа (7) - относительно немногочисленная (384 квартала) низкоэтажных (14 м), но со значительной площадью оснований; дополнительно характеризуются довольно высокой запечатанностью (39%). В группу попали в основном, промышленные (в том числе - унаследованные, с измененными функциями), коммунально-складские и транспортные кварталы. Отметим, что многие из них включают исторически ценную архитектуру русского "заводского" стиля конца XIX - начала XX вв., а также постройки 1930-х гг. советской эпохи в стиле конструктивизма (Рис. 22.45)

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

Подведем предварительные итоги (полученные результаты обобщены в таблице 22.7).

1. Пяти признаков оказывается достаточным для дифференциации кварталов по морфологическим признакам даже для такого сложного города, каковым является Москва.
2. Мы выбрали количество классов (10) ориентируясь на "универсальную" классификацию, предложенную американскими специалистами, но очевидно, что как и все универсальные таксономические системы она плохо проходит проверку на локализацию. В условиях Москвы - с ее резко отличавшимися друг от друга историческими эпохами застройки и перестройки кварталов - показатели разброса значений по каждому параметру играют едва ли не большую роль, чем собственно параметры, и эта особенность в дальнейшем (по мере реновации старых кварталов, точечной застройки, появление многих индивидуальных проектов) будет только усиливаться. Поэтому 10 классов (групп) - навряд ли достаточно для отображения реальной картины.
3. Для получения более детального результата следует ввести ко всем (а не отдельным) параметрам дополнительные показатели разброса значений (в особенности для таких параметров как средняя высота и площадь строений в кварталах), и увеличить число классов (в терминологии ESRI - "групп").
4. Для получения более корректного результата целесообразно классифицировать морфологию кварталов внутри отдельных типов функциональных зон (жилых, общественно-деловых, торгово-коммерческих, производственных).

Общее описание морфотипа	№	Распространение	Высотность	Разнообразие форм зданий	Плотность застройки	Закрытость
Исторические жилые и старопромышленные средние по площади кварталы Центра и примагистральных лучей	1		Среднеэтажные	Значительное разнообразие форм	Средняя	Закрытые
Исторические жилые и старопромышленные маленькие кварталы Центра (в пределах Садового Кольца)	10		Низкоэтажные	Близкие по форме исторические типы*	Высокоплотная	Закрытые
Спорадически распространенные маленькие и небольшие по площади общественно-деловые и жилые кварталы	6		Среднеэтажные	Разнообразные по форме	Высокоплотная	Полузакрытые
Кварталы жилой и старопромышленной застройки Центра и ближней к Садовому Кольцу периферии	2		Среднеэтажные	Разнообразные по форме	Среднеплотная	Полуоткрытые
Различные по размеру (небольшие и крупные) кварталы старопромышленной застройки Центра, ближней и дальней периферии	7		Низкоэтажные	Разнообразные по форме	Высокоплотная	Полуоткрытые
Крупноареальные новые кварталы жилой застройки	8		Среднеэтажные	Однотипные (серийные) по форме	Низкоплотная	Открытые
Крупноареальные новые и новейшие кварталы жилой застройки	3		Среднеэтажные	Разнообразные по форме	Низкоплотная	Открытые
Спорадические распространенные на средней и дальней периферии кварталы жилой и общественно-деловой застройки	5		Высокоэтажные	Невысокое разнообразие форм	Низкоплотная	Полуоткрытые
Широко распространенные вдоль МКАДа и магистральных лучей различные по размерности кварталы жилой и общественно-деловой застройки	9		Низкоэтажные	Типовые и серийные форм	Низкоплотная	Открытые
Редко и спорадически распространенные маленькие кварталы из нескольких высотных зданий новейшей застройки	4		Высотная (Стилобаты)	Типовые и серийные форм	Среднеплотная	Полузакрытые

22.5. Моделирование индексов формы

Для расчета Микроклиматических Индексов и Индексов Формы (см. Таблица 22.3) необходимо сначала получить общую растровую поверхность кварталов, которая будет складываться из высоты земной поверхности (DSM - Digital Terrestrial Surface) и высоты зданий. Предварительно трансформируем полигональный файл кварталов Tracts в растр по полю Hb (высота строений) в качестве Value (переменной). Величину ячейки (cell size) следует выбрать такой, чтобы точность растра соответствовала размерности полигонов строений (для слоя Building OSM это приблизительно 4 ячейки). Далее используем инструмент Cell Statistics набора Local группы Spatial Analyst Tools ArcMAP10.x.

Input rasters or constant values: DSM_oldMSCW, Hb_DSM
Output raster: BuildDsmH
Overlay statistic (optional): SUM
Environment Settings >> Raster Analysis >> Cell Size = Same as a layer Hb (4 cell)

Рис. 22.46 Совмещенная растровая модель высоты зданий и сооружений с высотой земной поверхности

К Индексам формы относятся ряд показателей, так или иначе отображающих открытость (и наоборот - закрытость/затененность) различных поверхностей внутри квартала: крыш и стен строений, земной поверхности - обстоятельство, контролирующее способность получать и/или отражать солнечную радиацию. Инструменты для расчета этих и других микроклиматических индексов разработаны в SAGA GIS; как известно, все последние релизы пакета содержат специальные упаковки плагинов для имплементации в Toolbox ArcMAP10.x; тем не менее работа соответствующих скриптов в "родной" программе происходит и быстрее, и надежнее (без "вылетов").

Первый индекс Topographic Openness|Открытость Поверхности рассчитывается с помощью инструмента Topographic Openness набора Lighting, Visibility, относящегося к группе Terrain Analysis SAGA GIS и основан на алгоритме Йокоямы [Yokoyama et al., 2002]

Topographic Openness
Input:
Grid System: ...
Elevation = BuildDsmH
Output:
Positive Openness: create
Negative Openness: create
Option:
Radial limits (*радиус расчетного окна в единицах карты) = 10000 (default)
Direction: all (*открытость со всех сторон)
Number of Secto: 8 (*число румбов расчета)
Method: Lane Tracing
Units: Radian
Dufference from Nadir: checked (* т.е., открытость определяется сверху вниз от надира, в противном случае это будет видимость по сторонам горизонта в плане)

Рис. 22.47 Топографическая Открытость Поверхности (фрагмент)

Открытость Поверхности показывает насколько свободно (т.е., без визуальных барьеров и помех) может быть наблюдаем тот или иной объект предметно-пространственной среды (участок поверхности, стены, крыши) безотносительно к позиции наблюдателя и освещения (т.е., положения источника света, например - Солнца). Инструмент можно использовать для определения открытости с какой-либо из сторон горизонта; в этом случае в позиции Direction выставляется угол в градусах (например, если моделируется открытость солнечным лучам с юга, указываем Direction = 180⁰).

Рис. 22.48 Топографическая Открытость Поверхности c южной стороны (фрагмент)

Остальные Индексы Формы рассчитываются инструментом Sky View Factor (O.Conrad, 2008) набора Lighting, Visibility группы Terrain Analysis SAGA GIS. Подробное описание алгоритма приведено в публикациях [Oke, 2000; Anders et al., 2009;, Boehner, Antonic, 2009].

Sky View Factor
Input:
Grid System: ...
Elevation: BuildDsmH
Output:
Visible Sky: create
Sky View Factor: create
Terrain View Factor: create
Average View Distance: create
Option:
Maximum Search Radiu: 10000 (*default)
Number of Sectors: 8 (*число румбов расчета)
Method: Sell Size

Таким образом алгоритм Sky View Factor позволяет получить несколько параметров. Параметр Visible Sky|Видимость небесной полусферы внешне схож с показателем Топографическая Открытость Поверхности, но у него другой физический смысл, поскольку на выходе - видимая часть небесной полусферы выраженная в процентах, т.е., открытость в направлении снизу-вверх.

Рис. 22.49 Visible Sky|Видимость небесной полусферы

Sky View Factor (SVF)|Фактор Видимости Неба - основной результат, выражаемый в долях единицы

Рис. 22.50 Фактор Видимости Неба: (фрагмент)

Terrain View Factor - Фактор Видимости Поверхности: Доля открытой (для солнечной радиации) земной поверхности

Рис. 22.51 Фактор Видимости Поверхности (фрагмент)

Shadow coefficient|Коэффициент Затененности - более сложный показатель, демонстрирующий Aбсолютную затененность поверхности безотносительно к сторонам горизонта (Рис. 22.52).

Рис. 22.52 Коэффициент Затененности (фрагмент)

Average View Distance|Средняя Дистанция Видимости - Средняя дистанция до экранирующих взгляд препятствий или линии горизонта.

Рис. 22.53 Средняя Дистанция Видимости (фрагмент)

22.6.Моделирование локальных климатических параметров морфотипов

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

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

Potential Incoming Solar Radiation инструмент набора Lighting Visibility группы Terrain Analysis SAGA GIS рассчитывает потенциальную прямую солнечную радиацию на основе растра Высоты искусственной поверхности и растра Фактора видимости неба. Алгоритм позволяет сравнить излучение получаемое объектами предметно пространственной среды с идеальной плоской горизонтальностью поверхностью.

ет

Potential Incoming Solar Radiation
Input:
Grid System: ...
Elevation: BuildDsmH
Sky View Factor: Sky_View_F_MSCW
Output:
Direct Insolation: create
Diffuse Insolation :create
Compare with Flat Terrain :create (*излучение попадающее на плоскую поверхность)
Option:
Local Sky View Factor (подается на вход)
Units: kWh/m²
Shadows: Fat (*плотная - поскольку стены бросают сплошную тень)
Location : Constant Latitude (*выбор способа учета географического положения)
Latitude degree: 55; (*широта Москвы, градуса)
Min: 45 (*минуты)
Sec: 20 (*секунды)
Time Span: 4, 21 (*время захода и восхода Солнца в Июле на широте Москвы)
Resolution: 1,0 (*шаг моделирования равен одному часу)
Atmospheric effects: Lumped Atmospheric Transmittance (*Сосредоточенное Атмосферное Пропускание)

Рис. 22.54 Прямая солнечная радиация (фрагмент)

Кроме основного растра прямой солнечной радиации алгоритм рассчитывает Diffuse Insolation|Рассеянную радиацию а также для сравнения радиацию для идеальной плоскости Compare to Flat Terrain.

Рис. 22.55 a) Рассеянная радиация и b) радиация смоделированная для идеальной плоскости: (фрагменты)

Инструмент учета "шероховатости" поверхности, оказывающей влияние на ветровой поток, называется Wind Effect Windward Leeward Index|Ветровой Эффект со стороны господствующих ветров ; находится в наборе Morphometry группы Terrain Analysis SAGA GIS. Ветровой Эффект - это безразмерный показатель, значения ниже 1 указывают на области, находящиеся в ветровой тени, тогда как значения выше 1 указывают на области, подверженные воздействию ветра, все с учетом выставленного в опциях направления ветра ([Boehner, Antonic, 2009; Gerlitz, Conrad, Boehner, 2015]) Второй показатель - Высота эффективного ветрового потока в метрах от нулевой поверхности. /p>

Wind Effect (Windward / Leeward Index)
Input:
Grid System: ...
Elevation: BuildDsmH
Wind Direction: not set
Output:
Wind Effect: create
Effective air Flow Heigh: create (*высота эффективного ветрового потока с учетом шероховатости квартала)
Option:
Search Distance :3 km (*с этим показателем следует поэкспериментировать, он зависит от величины города и размерности кварталов
Constant Wind Direction: 330 (*направление в градусах, в котором дует ветер: если летом 2021 г. Москве преобладали юго-западные ветры, то ветровой поток был направлен на ССЗ)
Acceleration: 1.5 (*величина потенциального ускорения ветрового потока, изменяется от 1.0 до 1.5)
Elevation Averaging: not checked (*не учитываем конвективный подъем с увеличением расстояния)

Рис. 22.56 Ветровой Эффект со стороны господствующих ветров (фрагмент)

Картограмма (Рис. 22.56) Ветрового эффекта демонстрирует усиления ветра с наветренной стороны на уровне + 1,2 м/сек в том время как с подветренной наблюдается торможение до - 0,8.

Рис. 22.57 Высота Эффективного Ветрового Потока (фрагмент)

Высота эффективного ветрового потока изменятся от 250 м в местах защищенных от ветра (ветровой тени) и 360 м на наветренных сторонах крышах невысоких зданий до 460 м над высокими (Рис. 22.57).

Wind Exposition Index|Ветровая Экспозиция инструмент набора Morphometry группы Terrain Analysis) SAGA GIS, который вычисляет открытые ( > 1) и затененные < 1) области для всех направлений с использованием углового шага.

Wind Exposition Index
Input:
Grid System: ...
Elevation: BuildDsmH
Wind Exposition Index: create
Output:
Wind Effect: create
Option:
Search Distance :3 km
Angular Step Size (Degree): 15 (*угловой шаг расчета в градусах)
Acceleration:1.5 (величина потенциального ускорения ветрового потока, изменяется от 1.0 до 1.5)
Elevation Averaging = not checked (не учитываем конвективный подъем с увеличением расстояния)

Рис. 22.58 Ветровая экспозиция (фрагмент)

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

Рассчитанные Индексы Формы и Микроклиматические Индексы определены для объемной поверхности предметно-пространственной среды кварталов. Чтобы выявить значения для кварталов (как операционно-территориальных единиц необходимо использовать либо Zonal Statistic ArcMAP10.x, либо (не выходя из SAGA GIS, в которой получены индексы) инструмент Grid Statistics for Polygons набора Shapes-Grid Tools группы Shapes.

Grid Statistics for Polygons
Grids:
Grid System: ...
Grid: BuildBsmH, Visible Sky, Positive Openness, Positive Openness_180, Terrain View Factor, Sky View Factor, Average View Distance, Shadow Coefficient, Direct Insolation, Diffuse Insolation, Wind Effect, Effective air Flow Heights, Wind Exposition (*список гридов, из которых извлекаются значения)
Shapes:
Polygons: TractsOldMSCW
Statistic :create
Mean: Checked (* все остальные статистики в положении unchecked, поскольку мы вычисляем среднее значение по кварталам)
Option:
Field Naming: grid name
Method: simple and fast

22.7. Моделирование ландшафтных и температурных индексов подстилающей поверхности

Для полной характеристики морфотипов необходимо добавить к уже существующему показателю Vegetation Coverage (долю зеленых насаждений от площади квартала) два ландшафтных индекса, получаемых на основе снимков Landsat: NDVI Нормализованный Разностный Вегетационный Индекс и NDBI Нормализованный Разностный Индекс Застройки .

Для определения доли зеленых насаждений можно использовать слой Vegetation OSM, который по Москве оцифрован с неплохой полнотой и/или растровые данные по древесному покрову Global Forest Change, хранящиеся на сайте университета Мэриленда [Hansen_et_al_2013]. Растровый слой Tree Cover отображает только древесные насаждения, в полигональном слое OSM Vegetation оцифрованы (хотя и с разной степенью детальности) и другие ярусы растительного покрова: в том числе кустарники и обычные, а также партерные газоны. Для исследования микроклиматических свойств кварталов может оказаться полезным учитывать оба параметра.

Рис. 22.59 Сравнение покрытий, отражающих растительный покров a) растровый слой Tree Cover b) полигональный слой Vegetation OSM (фрагменты на район Воробьевых гор)

Очевидно, что слой Vegetation OSM точнее в передаче элементов озеленения, находящихся на учете муниципальных предприятий осуществляющих уход ("горзеленхозов") - газонов, партеров, древесно-кустарниковых групп, аллей, фрагментов парковых и примагистральных насаждений, в то время как растр Tree Cover лучше отражает очертания мелкомассивных насаждений, зарослей "самосева" и древесных групп придомовых пространств - важные элементы для рассматриваемой модели, поскольку внутриквартальные элементы зеленой инфраструктуры непосредственно влияют на формирование микроклимата.

Для расчета озелененности кварталов при работе с векторными данными OSM используем инструмент Tabulate Intersection, при работе с растрами Tree Cover - инструмент Grid Statistics for Polygons SAGA GIS или Zonal Statistic ArcMAP10.x.

Рис. 22.60 Картограммы доли зеленых насаждений от площади кварталов, 7 классов "natural breaks" a) по данным Tree Cover b) по слою Vegetation OSM <

Итоговые картограммы доли зеленых насаждений от общей площади кварталов различаются весьма существенно: слой Vegetation чуть лучше отражает аллеи и примагистральные полосные насаждения, газоны и цветочные партеры исторического Центра Москвы, но слой Tree Cover ближе к действительности в отображении древесных придомовых насаждений и также многочисленных самосевных насаждений и просто зеленых массивов крупных кварталов в новых районах

Для расчета ландшафтных индексов (а в дальнейшем - температуры поверхности и моделирования тепловых зон) нам понадобится снимки Landsat 8 на летний период, которые можно скачать на одном из широко используемых сайтов NASA. Подбор снимка удобнее всего провести на GloVis, выставляя диапазон даты Data Range = 01.07.2021 to 31.07.2021; Cloud Cover = 0 or empty to 5%. В довольно дождливое лето 2021 года снимков, удовлетворяющих таким условиям, оказалось всего два, дальнейший выбор осуществляем экспертно - оценивая качество материала.

Снимок перекрывает территорию Старой Москвы, поэтому после загрузки и визуализации (Композитные сочетания Landsat) каналов необходимо вырезать интересующую часть полигоном границ, что может быть сделано в ArcMAP10.x Clip by Mask или в SAGA GIS с помощью алгоритма Clip Grid with Polygon набора Shapes Grid Tools группы Shapes. Преимущество последнего варианта в возможности одновременной (за один прогон инструмента) и почти моментальной обработки всех каналов.

Grid Statistics for Polygons
Input:
Grid System: ...
Grid: LC08_L1TP_179021_20210713_20210721_01_T1_B1, LC08_L1TP_179021_20210713_20210721_01_T1_B1, ***__B3, ***__B4, ... ***__B11 (***список растров - каналов Landsat 8)
Shapes:
Polygons: OldMSCW_Boundary
Target Extent: polygons

Рис. 22.61 Landsat 8, композит 4-3-2 каналов (в "натуральных цветах") a) исходный снимок, b) фрагмент полученный Clip Grid with Polygon для территории Старой Москвы SAGA GIS

Рассчитаем два необходимых индекса NDVI Нормализованный Разностный Вегетационный Индекс и NDBI Нормализованный Разностный Индекс Застройки с помощью Калькулятора Растров в ArcMAP10.x или SAGA GIS. (Напомним, что в SAGA GIS существует специальный инструмент Vegetation Index набора Tools группы Imagery).

Рис. 22.62 Нормализованный Разностный Вегетационный Индекс (NDVI) для территории Старой Москвы по снимку Landsat 8 от 13.07.2021

Нормализованный индекс различий застройки|Normalized index of building рассчитывается по каналам SWIR1 (6 - Коротковолновый инфракрасный 1) и NIR (5 - Ближний инфракрасный) Landsat 8 по формуле:

NDBI = (SWIR1 - NIR) / (SWIR1 + NIR)

Рис. 22.63 Нормализованный индекс различий застройки (NDBI) для территории Старой Москвы по снимку Landsat 8 от 13.07.2021

Land Surface Temperature (LST)|Температура Подстилающей Поверхности - еще один параметр, рассчитываемый на основе снимков Landsat с использованием специальных тепловых каналов (band 10 и band 11). Существует несколько вариантов расчета LST-индекса: по формуле, с помощью калькулятора растра, и с помощью специальных плагинов для QGIS - QGIS SPT или LSTE. Работа с плагинами предоставляет оперативное удобство пользователю и занимает меньше времени, но использование формул позволяет понять суть и особенности алгоритма (который не идеален), поэтому целесообразно использовать оба подхода, что даст дополнительную возможность верификации результатов.

A) Расчет по формуле. Расчет Температуры Подстилающей Поверхности LST с помощью калькулятора растра на основе публикации [Alipour, Saradjian, Esmaeily, 2003] подробно описан на сайте GIS LAB в статье Оценка температуры поверхности из снимка Landsat-8 при помощи Land Surface Temperature QGIS Plugin.

Для расчета используем все тот же снимок Landsat 8, зафиксировавший состояние поверхности на полдень июльского антициклона в Москве. Потребуется канал 10 и два относящихся к нему параметра RADIANCE_MULT_BAND_10 и RADIANCE_ADD_BAND_10, значения которых выбираются из метафайла (txt) стандартного набора Landsat 8. Метафайл можно открыть в Блокноте или (при использовании) файлового менеджера Total Commander простым нажатием F3. Сведения в метафайле представлены в виде строк длинного столбца объединенных в группы; параметры RADIANCE_MULT_BAND и RADIANCE_ADD_BAND находятся в группе RADIOMETRIC_RESCALING (радиометрическая калибровка или коррекция).

Рис. 22.64 10-й канал Landsat 8 от 13.07.2021 для территории Старой Москвы

Теперь можно определить величину Lλ излучение земной поверхности, воспринимаемое на орбите Земли дистанционным сенсором, (при условии нерассеивающей атмосферы) по формуле:

Lλ = RADIANCE_MULT_BAND_10 ∗ BAND10 + RADIANCE_ADD_BAND_10
из метафайла MTL.txt имеем
RADIANCE_MULT_BAND_10 = 3.3420E-04;
RADIANCE_ADD_BAND_10 = 0.10000
в калькуляторе растра набираем
0.00033420 * "LC08_L1TP_179021_20210713_20210721_01_T1_B10.tif" + 0.10000

Рис. 22.65 Lλ - излучение земной поверхности территории Старой Москвы 13.07.2021

Далее для расчета параметра TOA brightness|Яркостная Температура, характеризующая интенсивность излучения, необходимо предварительно извлечь две термические константы из метафайла MTL.txt Landsat 8:

K1 b K2 - термические константы, выбираются из GROUP = TIRS_THERMAL_CONSTANTS метафайла MTL.txt
K1_CONSTANT_BAND_10 = 774.8853
K2_CONSTANT_BAND_10 = 1322.0789

Теперь мы можем рассчитать фотометрический параметр Яркостная Температура по формуле:

TB = K2 / ln [(K1 / Lλ) + 1]
в калькуляторе растра набираем
1322.0789 / Ln((74.8853 / "Lλ") + 1)
где:
TB-TOA яркостная температура (K)
Lλ-TOA спектральное излучение (Watts/(m2*srad*µm))
K1-Band-удельная константа теплового преобразования из метаданных (K1_CONSTANT_AND_x, где x - номер теплового канала)
K2-Band - удельная константа теплового преобразования из метаданных (K2_CONSTANT_AND_x, где x - номер теплового канала)

Рис. 22.66 Яркостная Температура для территории Старой Москвы 13.07.2021

На последнем шаге определяем Land Surface Temperature по формуле:

LST=TB/[1+(λ∗TB/c2)∗ln(emissivity)]
где
λ- длина волны излучения, для 10 канала: 10.8
c2=h∗c/s=1.4388∗10−2 m K = 14388 µm K
h, постоянная Планка = 6.626∗10−34 J s
s, постоянная Больцмана = 1.38∗10−23 J/K
c, скорость света = 2.998∗108 m/s

Для расчета LST необходимо еще знать величину так называемого Emissivity|Коэффициента излучения - безразмерной величины изменяющейся для большинства естественных поверхностей Земли в пределах от 0,6 до 1,0. Авторы алгоритма полагают [Alipour, Sarajianb, Esmaeily, 2003], что для застроенных поверхностей эта величина близка к 0,912. Вводя значение Emissivity и добавляя перевод в градусы Цельсия, получаем:

LST=(TB/(1+(10.8∗TB/14388)∗ln(0.98))-273.15
в калькуляторе растра
LST=("TB" / (1 +10.8 * "TB" / 14388) * Ln(0.98)) - 273.15

Те же авторы [Alipour, Sarajianb, Esmaeily, 2003] предлагают и альтернативный вариант расчета с заменой Коэффициента излучения на показатель NDVI, т.е., использовать следующее выражение:

LST=(TB/(1+(10.8∗TB/14388)∗ln(NDVI))-273.15
в калькуляторе растра
LST=("TB" / (1 +10.8 * "TB" / 14388) * Ln(NDVI)) - 273.15

Рис. 22.67 a) LST территории Старой Москвы - расчет с использованием значения Emissivity b) LST территории Старой Москвы расчет с использованием значения NDVI

Как можно убедиться (Рис. 22.67) расчет с использованием значения NDVI дает более тонкозернистый результат, но в целом закономерности изменения температуры поверхности по территории Старой Москвы в жаркий июльский антициклон аналогичны.

B. Расчет Температуры Подстилающей Поверхности LST с помощью плагина Land Surface Temperature Estimation для QGIS.

После знакомства с последовательностью вычисления алгоритм, используемый в плагине (подобное описание Land Surface Temperature Estimation Plugin User Manual) выглядит для пользователя понятным: собственно модуль состоит из четырех этапов, отраженных на соответствующих вкладках диалогового окна [Ndossi, Avdan, 2016]. До начала работы загружаем в проект QGIS-Next-GIS необходимые каналы Landsat 8: band 10, band 11, band 4 "red", band 5 "NIR". Далее приступаем к последовательной работе с диалоговыми окнами - вкладками плагина.

Рис. 22.68 Общий вид модуля Land Surface Temperature Estimation, запущенного в QGIS

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

Radiance
Sensor Type: Landsat
Band Number: 10
Thermal Band: "LC08_L1TP_..._T1_B10.tif"
Output Raster: Radiance
Metadata File: LC08_L1TP_... _T1_MTL.txt
Offset Calibration Factor: 0.10000 (* из метафайле GROUP = RADIOMETRIC_RESCALING RADIANCE_ADD_BAND_10)
Output Format: GTIFF
Add Result to Project: checked

После заполнения всех позиций необходимо запустить плагин; выходной файл используется в следующих операциях. Плагин LST QGIS-Next-GIS обрабатывает площадь в виде прямоугольника, построенного по касательной к краям сцены, поэтому для получения истинных значений, следует вырезать границами сцены необходимую область из растра: в QGIS-Next-GIS инструмент набора GDAL Обрезать Растр по Маске; только после этого появляются реальные значения параметра.

Рис. 22.69 Первый расчетный параметр Radiance для территории Старой Москвы: a) выходной растр, b) растр, вырезанный по границам объекта

Вторая вкладка модуля LST позволяет рассчитать Brightness Temperature|Яркостную Температуру по шкале Кельвина. На этой вкладке пользователь указывает используемый датчик, номер полосы, растр яркости и выходной растр.

Brightness Temperature
Sensor Type: Landsat TIRs
Band Number: 10
Radiance Raster: RadianceR
OutPut Raster: BrightTemp
OutPut Format: GTFF
Add Resalt to Project: checked

На финальном этапе плагин LST позволяет пользователю выбрать различные алгоритмы расчета поверхностной температуры, один из них предполагает использование параметра Emissivity|Коэффициента излучения, другой - основан на использовании индекса NDVI. Поэтому следующая вкладка позволяет расчитать данный индекс по загруженным каналом Red и и NIR Landsat 8.

NDVI
Red: LC08_L1TP_ ..._01_T1_B4
NIR: LC08_L1TP_ ... _01_T1_B5
OutPut Raster: NDVI
OutPut Format: GTFF
Add Resalt to Project: checked

Рис. 22.70 Значения a) Brightness Temperature и b) NDVI для территории Старой Москвы 13.07.2021

Альтернативная ветвь алгоритма (вместо NDVI) использование Коэффициента излучения|Emissivity, для расчета которого предлагается несколько "рецептов" - апробирование и выбор остаются "на совести" эксперта.

Land Surface Emissivity Algorithm: выбираем Zhang, Wang et al Algorithm
Normalised Differenced Vegetation Index: NDVI
OutPut Raster: Emissivity
OutPut Format: GTFF
Add Resalt to Project: checked

На последнем этапе рассчитывается собственно Температура Подстилающей Поверхности (LST), пользователь может заполнять эту вкладку дважды, используя поочередно параметр Emissivity или NDVI

Land Surface Temperature Algorithm
Sensor Type: Landsat TIRs
Band Number: 10
Brightness Temperature Raster: BrightTemp
Land Surface Emissivity Raster: Emissivity (* или NDVI)
OutPut Raster: LST
Unit: Celsius
OutPut Format: GTFF
Add Resalt to Project: checked

Как можно видеть (Рис. 22.71) температура поверхности объектов на территории города в жаркий летний "антициклональный" день изменяется в довольно широких пределах от 24 C^o над кронами насаждений крупных зеленых массивов до 47 C^o на искусственных поверхностях плотно застроенных городских кварталов.

Рис. 22.71 Итоговый результат моделирования температуры поверхности с использованием плагина QGSI LST: a) полная сцена - территория Старой Москвы и b) центральная часть города

C. Расчет с помощью плагина QGIS SPT Альтернативный плагин для QGIS-Next-GIS - QGIS SPT - предлагает "одноступенчатый" алгоритм c загрузкой четырех каналов Landsat 8 - 4, 5, 10, 11 и вводом дополнительных данных. Пользователь может использовать два варианта плагина - упрощенный и усложненный, различающиеся набором данных. В первом случае вводятся диапазон ожидаемых температур и коэффициент испаряемости, остальные параметры на вкладке Advanced Setting принимаются "as is"; во втором - вкладка Advanced Setting отмечается как редактируемая и далее в каждую позицию последовательно вводятся параметры: значения NDVI для почв и растительного покрова, значения Emissivity для почв и растительного покрова в обоих тепловых каналах TIR-1 и TIR-2, таким образом для этого варианта исполнения модуля необходимо заранее рассчитать NDVI сцены.

SPT Data satellite: Landsat 8
Output Temperature: Celsius Range Temperature: 10 - 40⁰
Tot. Water vapored: 0.5 - 3.0 (*default может быть изменен в диапазоне 0.5 - 3.0)
Band Red: LC08_L1TP_ ... _01_T1_B4 (* 4 красный)
Band NIR: LC08_L1TP_ ... _01_T1_B5 (*5 ближний инфракрасный) Band TIR-1: LC08_L1TP_..._01_T1_B10
Band TIR-2: C08_L1TP_... _01_T1_B11
Advanced Setting: checked (* усложненный инвариант)
At- transmittance: (* по умолчанию указано значение для US 1976, можно поменять на Mid Lat Summer, т.е., летнее для средних широт)
NDVI Soil: 0,17 (* выбрано по участку реконструкции с открытым грунтом на Воробьевых Горах)
NDVI Vegetation: 0, 53
TIR-1 Emis.Soil: 0,89
TIR-1 Emis.Vegetation: 0,98
TIR-2 Emis.Soil: 0,91
TIR-2 Emis.Vegetation: 0,99
Geometrical: 1,0

Для территории Москвы дефолтные значения дают результат (Рис. 22.72), превосходящий и результат, полученный расчетом по формуле и результат с использованием плагина LST. Поскольку плагин работает быстро, пользователь может экспериментировать с вводимыми дополнительными значениями: результат расчета наиболее чувствителен к изменению значений двух переменных Tot. Water vapored целесообразно также изменить позицию At-transmittance.; при введении собственных коэффициентов результат оказался несколько ближе к полученным двумя первыми способами (т.е., к расчету по формуле и с помощью QGIS LST модуля).

Рис. 22.72 Итоговый результат моделирования температуры поверхности с помощью плагина QGSI SPT: a) полная сцена - территория Старой Москвы, b) центральная часть города

Таким образом, мы получили достаточно близкие значения поверхностной температуры для всех использованных методов (расчет по формулам в калькуляторе растра, расчет по двум плагинам): в середине летнего дня в условиях июльского антициклона значения температуры поверхностей изменяются от 24,3 ^oС в пределах зеленых массивов до 47,8 ^oС или 49,8 ^oС в пределах наиболее "разогретых" застроенных кварталов Старой Москвы. Оставляя в стороне несущественные различия в результатах мы можем утверждать, что дифферент по параметру поверхностной температуры (т.е., температуры незакрытого и задернованного грунта, заасфальтированных участков, стен зданий и сооружений, наконец, крыш домов, выполненных из разных материалов) достигает 25 ^oС.

Разумеется, нами был использован единственный снимок с "идеальными" параметрами (подходящий день по погодным условиям, подходящее время суток и минимальная облачность), характерный для антициклона июля 2021 г. Можно предполагать, что эффект "острова тепла" выглядит иначе при несколько иных условиях (длительности установления антициклона, облачности, и т.д.). В этом смысле представляет интерес определение типичных "тепловых зон" по совокупности разновременных летних снимков с помощью инструмента ISODATA|Изокластерного анализа [Грищенко, 2014].

Инструмент Iso Cluster Unsupervised Classification ArcMAP10.x (Spatial Analysis Tools >> Multivariate) позволяет анализировать разновременные растры как набор каналов и выделять на этой основе типы кластеров. На эксперта ложится задача выбора числа кластеров и параметра Minimum class size, определяющего размерность минимального ареала, который образует самостоятельный кластер.

Iso Cluster Unsupervised Classification
Input raster bands: LC08_L1TP_2016..._01_T1_B10, LC08_L1TP_2017..._01_T1_B10, LC08_L1TP_2020..._01_T1_B10, LC08_L1TP_2022..._01_T1_B10 (* всего 10 растров летнего сезона с 2016 по 2021 гг.)
Output classified raster: Tempr_Claster_20
Number of classes: 16 (* см. ниже)
Minimum class size: 20 (* default, можно менять в обе стороны уменьшая или увеличивая дробность дифференциации)
Sample interval: 10 (default)

Из всех переменных диалогового окна Iso Cluster Unsupervised Classification решающее значение имеет число кластеров: в данном случае мы исходим из числа определенных выше морфотипов застроенной части Старой Москвы - 10; можно предположить также, что остальные типы "городской ткани", не относящиеся к выбранным нами кварталам, также генерируют не менее 6 качественно различных ареалов "теплового поля", тогда общее число ожидаемых кластеров - 16.

Рис. 22.73 Итоговый результат моделирования температуры поверхности с помощью плагина QGSI SPT: a) полная сцена - территория Старой Москвы, b) территория между лесопарками Измайлово и Кусково

Интерпретация результатов приведена в таблице; инструмент Зональная Статистика позволяет извлечь статистические показатели температуры поверхности и представить значения (среднее, разброс) для каждого типа ареала тепловой структуры города в условия летнего антициклона.

Табл. 22.8 Температурные характеристики различных "тепловых зон" Старой Москвы при летнем антициклоне июля 2021 г.

№	G-CODE	Вид "городской ткани"	AREA	MEAN	MAX	RANGE
1	2	Высотная застройка, в том числе в окружении открытых пространств	86	32,3	37,2	11,1
2	3	Долины малых в пределах зеленых массивов	34	27,2	30,8	4,5
3	6	Зеленые массивы с преобладанием широколиственных пород	3185	28,9	36	11,3
4	7	Центральные части акваторий озер и крупных рек	1051	28,3	33,8	9
5	8	"Ядра"и внутренние части городских лесопарков с насаждениями высокой плотности	6612	27,1	35,7	11,3
6	9	Внешние части городских лесопарков и насаждения низкой плотности с прогалинами и полянами	3185	28,9	36	11,3
7	12	Внутренняя часть периферии лесопарков и зеленых насаждений	1411	29,1	36,8	10,6
8	13	Периферийные части акваторий рек и озер	3218	30,1	38,5	12,9
9	14	Внешняя опушечная полоса лесопарков и зеленых насаждений	4463	31,5	37,9	11,2
10	15	Индустриальные и коммунально-складские площадки в том числе с открытыми пространствами и газонами	11119	34,3	44,4	17,5
11	20	Зоны отчуждения железных дорог, строительные площадки, реконструкция, открытые грунты	11757	37,6	47,8	19,9
12	16	Жилые кварталы с рыхлой и среднеплотной застройкой со сложившимся озеленением придомовых пространств	66	32,1	35,6	8,5
13	17	Кварталы со средней степенью озеленения придомовых пространств (рыхлая застройка)	103	32,9	39,3	14,2
14	18	Набережные с искусственными покрытиями, застроенные склоны долин	10123	33,5	40,3	13,2
15	19	Плотно застроенные и закрытые кварталы с невысокой и низкой степенью озеленения придомовых пространств	32730	35,5	43,9	16,3
16	20	Рыхлая застройка различного назначения (коммунально-складская, торговые центры) с преобладанием искусственных запечатанных поверхностей	11757	37,6	47,8	19,9

Как можно судить по картограмме [Рис. 22.73] и таблице территория Старой Москвы ожидаемо разделилась на три группы "тепловых зон", отчетливо различаемых по основным характеристикам: средней и максимальной температуре поверхности а также - по разнице максимальных и минимальных значений в пределах типа (Range). К первой группе, которую условно можно назвать перегретой тепловой зоной относятся:

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

В пределах ареалов перегретой тепловой зоны средняя температура поверхности выше 35 ^oС, а максимальная изменяется от 43 ^oС до 47 ^oС и более.

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

"Ядра"и внутренние части городских лесопарков с насаждениями высокой плотности;
Внутренняя часть периферии лесопарков и зеленых насаждений;
Внешняя опушечная полоса лесопарков и зеленых насаждений;
Внешние части городских лесопарков и насаждения низкой плотности с с прогалинами и полянами;
Зеленые массивы с преобладанием широколиственных пород;
Долины малых в пределах зеленых массивов;
Периферийные части акваторий рек и озер;
Центральные части акваторий озер и крупных рек.

В прохладной тепловой зоне средние температуры изменяются в диапазоне 27-31 ^oС, а максимальные в диапазоне 30-38 ^oС. Примечательно, что прослеживается четкое повышение температур при перемещении из внутренних ядерных частей зеленых массивов к их периферии и опушкам.

Наконец третья группа представлена кварталами жилой застройки различными по высоте и плотности, возрасту и озелененности: ее можно считать умеренно нагретой.

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

Средняя температура в этой группе 32-34 ^oС, но как можно видеть все ее кластеры характеризуются высокой степенью разброса (показатель RANGE) 14-16 ^oС, что означает присутствие большого разнообразия внутри каждого отдельного класса. Как отмечает М.Ю. Грищенко [Грищенко, 2014], для выявления дифференциации тепловых характеристик кварталов необходимо использовать большее число кластеров, либо запускать процедуру кластеризации предварительно вырезав сцену без крупных зеленых насаждений.

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

22.8. Регрессионный анализ температуры подстилающей поверхности города как зависимой переменной

Для анализа причинно-следственных связей необходимо включить зависимую переменную в набор параметров, характеризующих кварталы. Для этого можно использовать инструмент Grid Statistics for Polygons (Shapes >> Shapes-Grid_Tools) SAGA GIS: на входе растр температуры поверхности и шейп-файл кварталов, в качестве переменных статистик выберем MEAN, MAX, Variance.

Рис. 22.74 Температура поверхности по кварталам Старой Москвы: a) среднее значение b) максимальное значение

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

В нашем распоряжении - пять групп переменных:

Индексы Подстилающей Поверхности,
Индексы Компоновки,
Индексы формы,
Микроклиматические индексы,
Ландшафтные Индексы.

Если привлечь всю совокупность факторов мы почти наверняка получим "перегруженную модель", поэтому целесообразно провести предварительный отбор наиболее "многообещающих" переменных из каждой группы. Для этого построим точечные графики (scatterplot) для зависимой переменной и различных параметров в ArcMAP10.x (View >> Grafs >> Create scatterplot Matrix), либо в SAGA GIS (R(C на слое Tracts >> Atributes >> Scatterplot). В открывшейся вкладке в позицию Y назначаем зависимую переменную - T_MEAN, а в позиции X поочередно указываем различные независимые переменные.

Рис. 22.75 Диалоговое окно Scatterplot SAGA GIS

Вне зависимости от выбора программы (ArcMAP10.x или SAGA GIS) имеет смысл строить именно парные графики - это займет больше времени, чем построение множественного графика, но зато позволит лучше оценить переменные в качестве потенциальных кандидатов на объяснение феномена. Построим такие графики линейной регрессии и сравним коэффициенты ковариации для средних|MEAN и максимальных|MAX значений температур поверхности.

Рис. 22.76 Графики линейной регрессии для температуры поверхности (зависимая переменная) и предикторов: a) доли зеленых насаждений Percent_Tree_Cover b) стандартного отклонения средней площади домов Ab_DEV в кварталах; верхний ряд для средних значений и нижний ряд - максимальных значений

Несмотря на достаточно тесную зависимость между средними и максимальными значениями температуры поверхности (65%) рейтинг факторов обнаруживает значительные различия.

Рейтинг факторов для средних значений поверхностных температур

1. Percent_Tree_Cover - 21,3%
2. Ab_DEV - 9,6%
3. Percent_Veg - 7,1%
4. Visible_Sky - 5,2%
5. Positive_Opennes - 5,2%
6. Hb_AVG - 4,8%
7. Terrain_View_Factor - 4,7%
8. Hb_DEV - 4,6%
9. Abs_AVG - 4,4%
10. Diffuse_Insolation - 3,8%
11. Sky_View_Factor - 3,7%
12. Effective_Wind - 3,5%

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

1. Ab_DEV - 16,6%
2. Ab_SUM - 13,3%
3. Percnt_Tree_Cover - 12,6%
4. Pb_DEV - 8,1%
5. Pb_SUM - 7,5%
6. Af_SUM - 6,7%
7. Percent_Veg - 5,5%
8. Ac_DEV - 5,5%
9. Number_MAX - 5,5%
10. N_Dist_DEV - 5,3%
11. Visible_Sky_f - 4,4%
12. Diffuse_Insolation - 4,4%

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

Прежде чем обратиться к собственно регрессионному анализу мы должны оценить наборы факторов на предмет "перегруженности" модели, которая возникает при привлечении взаимосвязанных объясняющий переменных, используя процедуру Анализ Главных Компонент (PCA). В ArcMAP10.x инструмент PCA предназначен для обработки растров, и, чтобы воспользоваться им нам пришлось бы предварительно построить тематические растры для каждого поля - переменной полигонального файла квартала, поэтому обратимся к возможностям инструмента Principal Component Analysis SAGA GIS (Table >> Calculus).

Principal Component Analysis
Tables:
Input Table: Tracts
Attribute: Percent_Tree_Cover, Ab_DEV, ... (* отмечаем необходимые поля с переменными по списку для средних температур)
Output:
Principal Component Analysis = create
Options:
Method: variance-covariance matrix
Number of component: 3 (*default)

Располагая результатами Анализа Главных Компонент, можно освободиться от тесно связанных переменных, для которых коэффициент корреляции >0,65, при этом из каждой такой пары мы оставляем переменную с большим процентом детерминации относительно зависимой переменной (в нашем случае - температуры поверхности). Среди связанных переменных:

Группа микроклиматических факторов:
Sky View F - Terrain_View_Factor 0,93
Visible_Sky - Sky View F 0,82
Visible_Sky - Positive Openness - 0,95
Visible_Sky - Terrain_View_Factor 0,94
Оставляем Visible_Sky.

Группа индексов формы:
Ab_SUM - Af_SUM - 0,75
Ab_SUM - Pb_SUM - 0,83
Оставляем Ab_SUM
Pb_DEV - Ac_DEV - 0,70
Оставляем Pb_DEV.

Список факторов для регрессионного анализа средних значений поверхностных температур с указанием "объясняемой" доли множества кварталов для отдельных факторов:

1. Percent_Tree_Cover - 21,3%
2. Ab_DEV - 9,6%
3. Percent_Veg - 7,1%
4. Visible_Sky - 5,2%
5. Hb_AVG - 4,8%
6. Hb_DEV - 4,6%
7. Abs_AVG - 4,4%
8. Diffuse_Insolation - 3,8%
9. Effective_Wind - 3,5%

Список факторов для регрессионного анализа максимальных значений поверхностных температур с указанием "объясняемой" доли множества кварталов для отдельных факторов:

1. Ab_DEV - 16,6%
2. Ab_SUM - 13,3%
3. Percnt_Tree_Cover - 12,6%
4. Pb_DEV - 8,1%
5. Percent_Veg - 5,5%
6. Number - 5,5%
7. Visible_Sky_f - 4,4%
8. Diffuse_Insolation - 4,4%

Для выявления предполагаемых причинно-следственных связей используем Exploratory Regression ArcMAP10.x последовательно для средних и максимальных температур в качестве зависимых переменных.

Exploratory Regression
Input Features: Tracts
Dependent Variable: T_MEAN (* в варианте для максимальных значений T_MAX) Candidate Explanatory Variables: Weights Matrix File (optional): blanc
Output Report File (optional): Exploratory_MEAN_8f.txt

Располагая предварительными оценками попарных регрессионных уравнений полученных в SAGA GIS мы не можем ожидать достаточное количество проходящих моделей (Passing Model) с высокими значениями коэффициентов корреляции и их производных, поэтому есть смысл несколько "смягчить" пороговые значения различных показателей в "подвале" диалогового окна Exploratory Regression:

Search Criteria (RCl)
Maximum Number of Explanatory Variables (optional): 8 (по числу факторов)
Minimum Number of Explanatory Variables (optional): 3 (начальное число факторов в первой по очереди модели);
Minimum Acceptable Adj R Squared: 0,3 (наименьшее скорректированное значение R-квадрат для проходящих моделей; значение по умолчанию равно 0,5);
Maximum Coefficient p value Cutoff: 0,1 (значение отсечки для p-оценки, по умолчанию равно 0,05, указывающее, что проходящие модели будут содержать только объясняющие переменные, коэффициенты которых статистически находятся на уровне 95% достоверности, чтобы ослабить этот пороговый барьер для проходящих моделей, вводим большее значение p-значения);
Maximum VIF Value Cutoff: 7,5 (коэффициент инфляции дисперсии отражает допустимую степень избыточности (мультиколлинеарности) среди объясняющих переменных модели, поскольку мы попытались убрать избыточность на этапе подбора факторов - оставляем дефолтное значение);
Minimum Acceptable Jarque Bera p value: 0,05 (Минимально допустимое значение p-оценки возвращаемое диагностическим тестом Jarque Bera указывающее нормально ли распределены остатки модели, по умолчанию p равно 0,1, для смягчения этого критерия введем меньшее минимальное значение p);
Minimum Acceptable Spatial Autocorrelation p value: 0,05 (Минимально допустимое значение пространственной автокорреляции p-оценки, отражает результат проверки на остатки модели для пространственной кластеризации с использованием глобального Морана I Проходящие модели должны иметь большие значения p для этого диагностического теста -по умолчанию равно 0,1, для ослабления вводим 0,05).

Рассмотрим полученный TXT.File Отчета (см. Исследовательская регрессия). Первый раздел демонстрирует варианты таблиц значений контрольных коэффициентов для трех лучших моделей (Adjusted R-Squared Results) начиная от минимального числа переменных (3 из 9) до максимального (9 из 9).

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

Обращаем внимание на следующие обстоятельства:

1. Число "проходящих моделей"|Passing Models - ни одна модель не соответствует всем (даже смягченным) критериям;
2. Величина максимального достигнутого в моделях Скорректированного коэффициента детерминации AdjR2 - 0,36 (т.е., наша модель способна объяснить только 36% выборки):
3. Минимальное число вовлеченных переменных, после которого коэффициент детерминации AdjR2 больше не увеличивается - 6 (т.е., 8-й и 9-й не прибавляют "объясняющей силы" к модели) и могут быть удалены;
4. Перечни шести или семи "работающих на модель" переменных;
5. Знак связи для каждой из значимых переменных ("+", "-").

В следующем разделе Общее резюме Исследовательской Регрессии выясняем каким именно контрольным коэффициентам не соответствует построенная модель.

Рис. 22.78 Второй раздел Отчета линейной регрессии: Exploratory Regression Global Summary

Модель прошла смягченные пороговые значения трех контрольных коэффициентов Min Adjusted R-Squared, Max Coefficient p-value, Max VIF Value, т.е., удовлетворяет условиям минимально значимой детерминации и мультиколлинеарности, но ни один из 465 построенных инвариантов не прошел тест на нормальность ошибок (Jarque–Bera - минимальное значение 0,10), и тест на отсутствие пространственной автокорреляции (так называемая "нулевая гипотеза" Глобального Индекса Морана). Эти два теста - косвенные свидетельства ненормального распределения переменных и наличие собственной зависимости значения переменных от положения в пространстве.

Третий раздел Отчета Exploratory Regression Global Summary демонстрирует значимость предикторов и характер их связи с зависимой переменной - здесь все ожидаемо: шесть лидирующих факторов, способных в совокупности объяснить поверхностную температуры для 36,56% кварталов Старой Москвы.

Рис. 22.79 Третий раздел Отчета линейной регрессии: Exploratory Regression Global Summary

К факторам с наибольшими коэффициентами детерминации относятся в круглых скобках указан знак связи (+) - прямая, (-) - обратная:

1. Доля древесных насаждений в квартале 100% (-);
2. Доля "зеленки" (по OSM с учетом газонных покрытий и всех форм озеленения) 100% (-);
3. Разнообразие (стандартное отклонение) площади домов в квартале 100% (+);
4. Разнообразие (стандартное отклонение) высоты домов в квартале 89,02% (-);
5. Видимость неба 88,62%,(73,98 (+) и 26,02 (-);
6. Средняя высота домов в квартале 85,37% (+);
7. Рассеянная инсоляция 83,33% (+);
8. Средняя застроенность (запечатанность) кварталов 80,49% (+).

Тем не менее, и этот рейтинг не может пока считаться идеальным сочетанием "объясняющих" факторов для модели. Четвертый раздел Отчета Summary of Multicollinearity|Выводы по мультиколлинеарности показывает сколько раз близкие переменные, обеспечивающие одинаковый процент вариации включались в модель.

Рис. 22.80 Четвертый раздел Отчета линейной регрессии: Summary of Multicollinearity

В таблице Summary of Multicollinearity обнаруживаются две переменные с высоким VIF variance inflation factor|Коэффициентом инфляции дисперсии: VISIBLE_SKY (VIF=10,22) и DIFFUSE_INSOLATION (VIF=9,28), которые включались в модель 87 и 51 раз соответственно и которые имеют коэффициент ковариации 40,48%. Эти же две переменные присутствуют в последнем пятом разделе отчета в таблице Summary of Residual Spatial Autocorrelation (SA) во всех трех моделях с наибольшим показателем по данному признаку. Следовательно несмотря на значимость обеих переменных от одной из них лучше освободиться и запустить алгоритм еще раз - с шестью переменными, избавившись от слабого фактора WIND EFFECT, фактора с высокой мультиколлинеарностью DIFFUSE_INSOLATION и фактора с неоднозначным знаком детерминации Abs_AVG.

Рассмотрим вкратце результат оптимизации модели с шестью факторами. Инвариант "6 из 6" первого раздела текстового отчета показывает нам максимальный Скорректированный коэффициент детерминации AdjR² равный 0,35, следовательно, избавившись от трех переменных мы получили очень незначительное снижение объема объясненной выборки изучаемого феномена - 35% против 36% при девяти переменных. Возросла значимость оставшихся факторов (до 100%) и появилась однозначная связь для всех участвующих переменных - кроме одной (стандартное отклонение высоты зданий в квартале)

Рис. 22.81 Четвертый раздел Отчета линейной регрессии: Summary of Multicollinearity

Несмотря на то, что ни одна из 42 построенных инвариантов модели не прошла тест на нормальность ошибок (Jarque–Bera) и тест на отсутствие пространственной автокорреляции, нам удалось избавиться от мультиколлинеарности и получить "условно идеальную" модель, которая отображается и в таблице dbf при сортировке значения Акайке|AICc (от минимума):

1. Percent_Tree_Cover - 100%
2. Ab_DEV - 100%
3. Percent_Veg - 100%
4. Visible_Sky - 100%
5. Hb_AVG - 100%
6. Hb_DEV - 84,6%

Осуществим процедуру Exploratory Regression для зависимой переменной Максимальные температуры поверхности в квартале с несколько иным (по результатам Анализа Главных Компонент) набором факторов. Не повторяя всех этапов оптимизации (мы можем заранее исключить Diffuse_Insolation) обратимся сразу к итоговым результатам. Эта модель также не проходит тест Жака-Бера на нормальность ошибок и тест на отсутствие пространственной автокорреляции, но позволяет выделить 5 факторов со 100% значимостью (знак связи указан в скобках):

1. Ab_SUM Сумма площадей оснований зданий (+),
2. Ab_DEV Среднее квадратическое отклонение площади зданий (+),
3. Percnt_Tree_Cover Доля древесных насаждений - tree cover (-),
4. Percent_Veg Общая доля озеленения по OSM (-),
5. Visible_Sky Фактор видимости неба (+).

Используя данный список переменных для зависимой переменной Максимальная температура поверхности обратимся к алгоритму Geographically Weighted Regression|Географически Взвешенной Регрессии ArcMAP10.x. Предварительно необходимо получить Матрицу Пространственных весов с помощью инструмента Generate Spatial Weights Matrix.

Generate Spatial Weights Matrix Input Feature Class: Tracts.shp
Unique ID Field=ID
Output_Spatial_Weights_Matrix_File: SW_OLS5factor.swm
Conceptualization_of_Spatial_Relationships="K_NEAREST_NEIGHBORS"
, Distance_Method: "EUCLIDEAN",
Row_Standardization: "ROW_STANDARDIZATION",

GWR ArcMAP10.x дает три основных результат: таблицу диагностики, выходной класс объектов (шейп-файл), содержащий все локальные оценки, растровые поверхности, отображающие каждый коэффициент отдельно.

ArcToolbox >> Spatial Statistic Tools >> Modeling Spatial Relationships >> Geographically Weighted Regression
Input Feature Class: Tracts.shp
Dependent Variable: T_MAX
Explanatory Variable: Ab_SUM, Ab_DEV, Percnt_Tree_Cover, Percent_Veg, Visible_Sky (*5 переменных)
Output Feature Class: .....GWR5values.shp
Kernel Type: ADAPTIVE
Bandwidth method :BANDWIDTH_PARAMETER
Weights: T_MAX
Coefficient raster workspace: ...MRPHO
Leave all other fields blank or as filled by default
OK

Рис. 22.82 Таблица диагностики DBF-файла географически взвешенной регрессии

Модель GWR объясняет около 87,5% вариации зависимой переменной, что значительно лучше чем 37% обычной модели линейной регрессии OLS. Ликвидация двух переменных позволила значительно оптимизировать модель: в частности, критерий Акайке (AICc) уменьшился с 19205,0 до 57,1.

Сравним отображение стандартизированных остатков для Географически взвешенной регрессии GWR и линейной регрессии Метода наименьших квадратов OLS.

Рис. 22.83 Распределение остатков a) для OLS, b) GWR

Как можно видеть (Рис. 22.83) распределение остатков в картограмме GWR содержит меньше кварталов в крайних значения стандартного отклонения ( -2.5 Std.Dev. : +2.5 Std. Dev.) чем в картограмме OLS и, что важнее они распределены случайно и не образуют скоплений (красных или синих кварталов).

Итак, модель Взвешенной географической регрессии с пятью независимыми переменными объясняют 87% изменчивости признака - максимальных температур поверхности города в условиях летнего антициклона. Насколько такой результат можно считать удовлетворительным? С одной стороны, построенная модель не может претендовать на исчерпывающее объяснение феномена: явно упущены некие важные факторы, которые (допустим) могут оказаться связанными с более тонко локализованными причинами - материалом стен (бетон, кирпич), типом (плоская, двух- или четырехскатная) и материалом крыш, экспозицией стен и крыш (ориентированные в южном направлении нагреваются сильнее). Мы вправе также предположить, что кварталы в целом для данного анализа не вполне подходящие операционно-территориальные единицы, поскольку реальная ситуация может быть в значительной степени дифференцирована внутри квартала. С другой стороны, можно достаточно уверенно утверждать, что формирование "острова тепла" в Москве в условиях антициклона достаточно четко связано как с индексами формы, так и с индексами подстилающей поверхности. Риски "перегрева" возрастают с увеличением площади домов и сложности их формы (фактор разброса значений площади), а также - с увеличением индекса видимости неба. Охлаждают город зеленые насаждения, причем как древесные насаждения (кроны деревьев), так и кустарниковая растительность и газонные покрытия.

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