Возможность применения карт элементарных геохимических ландшафтов для картографирования содержания тяжелых металлов и металлоидов в городских почвах

Город – сложная природно-антропогенная среда, ключевыми компонентами которой являются люди, промышленные объекты, селитебные кварталы, парки и разветвленная транспортная сеть. В результате деятельности человека помимо механической трансформации ландшафтов (выравнивание, искусственные холмы, бетонирование берегов рек) существенное влияние на здоровье людей оказывает химическое загрязнение тяжелыми металлами и металлоидами (ТММ), недоступное для глаза человека, но имеющее наиболее сильное негативное воздействие (Сает и др., 1990; Demetriades, Birke, 2015; Mapping…, 2011). Определение их содержания в почвах городских ландшатов имеет большое значение с точки зрения оценки рисков для здоровья населения. Накопление ТММ в городской среде снижает качество воздуха, почв, а также угрожает здоровью при попадании в организм человека через кожу, пищеварительную систему или при вдыхании. При этом из-за сложности, трудоемкости и дороговизны анализов шаг опробования почв обычно составляет порядка 500–700 м или больше. Однако для оценки городской территории с точки зрения загрязнения ТММ необходимо иметь возможность выполнять детальное картографирование и районирование городской среды. Для достоверного картографирования и последующей оценки рисков необходимо интерполировать значения, полученные в точках опробования, с учетом факторов, влияющих на распределение элементов в почвах. К таким факторам относятся: локальные особенности циркуляции атмосферы, определяющие направление и интенсивность разноса загрязнений от источников; рельеф местности, влияние которого в значительной степени определяет латеральный перенос вещества, механический и литологический состав грунтов, определяющий возможность вертикальной миграции и накопления загрязнений. При этом учесть влияние рельефа сравнительно проще, чем учесть влияние других факторов, с одной стороны, благодаря доступности данных (в виде цифровых моделей рельефа, ЦМР), и, с другой стороны, в результате развития методов геоморфометрического анализа и цифрового почвенного картографирования (Florinsky, 2016). Привлечение сведений о высотах местности в виде ЦМР потенциально может повысить детальность и достоверность моделирования загрязнения почв в условиях, когда отбор проб выполняется по относительно редкой сетке.

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

ОБЪЕКТЫ И МЕТОДЫ

Объектом исследования являются почвы г. Дархана, третьего по величине города Монголии, административного центра аймака Дархан-Уул. Население г. Дархан, по данным статистической службы Монголии, превышает 86 тыс. человек (Population…, 2020) . Город расположен в долине р. Хара (правый приток р. Ор-хон), в пределах крупной геоморфологической области - Центрально-Монгольской приподнятой равнины между Сибирской и Китайской платформами (Геоморфология,,,, 1982) . Климат резкоконтинентальный со значительными годовыми и суточными колебаниями температуры воздуха ( t_янв = -18...-25 °C, t_июн = +18...+20 °C), резкой сменой сезонов, значительной сухостью воздуха, малым количеством осадков (300-350 мм в год). В течение года преобладают ветры южного и северного направлений (Мурзаев, 1952) .

Основным источником загрязняющих веществ в районе г. Дархан является промышленная зона к югу от городской застройки. Здесь располагаются, в частности, Дарханская ТЭЦ, работающая на буром угле, крупнейший в Монголии металлургический завод и некоторые другие производства. Другой источник, специфический для городов Монголии – юрточные кварталы, где широко применяется отопление углем при помощи небольших бытовых печек (Kosheleva et al., 2019) . Загрязняющие вещества с выбросами поступают в атмосферу, а затем осаждаются и накапливаются в городских почвах.

Почвенно-геохимическая съемка территории проводилась летом 2011–2014 гг. Было отобрано 126 смешанных проб, включая три фоновые. Пробы отбирались из поверхностного горизонта (0– 10 см) по нерегулярной сетке с шагом примерно 500–700 м (Kosheleva et al., 2019; Timofeev et al., 2019) . Схема размещения точек опробования показана на рисунке 1.

Валовое содержание ТММ в пробах почв и растительного материала анализировалось масс-спектральным методом с индуктивно-связанной плазмой в ВИМС. Для подробного анализа выбраны типичные для данной территории загрязнители: As, Cd, Cr, Cu, Pb, Sb, W (Kosheleva et al., 2019) . При обработке полученных данных вычислялись коэффициенты концентрации EF l = C u /C b , где C u и C b - концентрация элемента в городских и фоновых образцах.

Рис. 1. Географическое положение г. Дархан ( а ) и размещение точек отбора проб на его территории ( b ).

Fig. 1. The location of Darkhan, Mongolia ( a ), and the distribution of soil samples over study area ( b ).

В настоящее время картографирование непрерывно изменяющихся свойств почв обычно выполняется с привлечением методов геостатистики (Геостатистика…, 2007). Для оценки возмож- ности картографирования был выполнен геостатистический анализ значений коэффициентов концентрации ТММ в точках опробования и построены вариограммы показателей. Вариограммы приведены на рисунке 2. Видно, что для всех коэффициентов концентрации, за исключением Cr и Cu, характерны примерно одинаковые величины разности вне зависимости от расстояния. Такой вид – так называемый “чистый эффект самородка”, pure nugget, – имеют вариограммы данных, в которых отсутствует пространственная корреляция (Демьянов, Савельева, 2010).

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

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

По набору и определенным соотношениям ведущих геохимических процессов Б.Б. Полынов (1956) различал элювиальные (автономные), супераквальные и субаквальные ландшафты. Эта классификация была дополнена М.А. Глазовской (2002) трансэлювиальными ландшафтами склонов, трансаккумулятивными нижних частей склонов, аккумулятивно-элювиальными замкнутых понижений на водораздельных поверхностях. Эти кате- гории ЭЛ в классификации А.И. Перельмана соответствуют уровню рода (Перельман, Касимов, 2000).

Рис. 2. Полувариограммы коэффициентов концентрации (DF) изучаемых ТММ.

Fig. 2. Semivariograms of DF concentration coefficients of considered HMMs.

Карты родов ЭЛ, как правило, составляются вручную на основе экспертных знаний и поэтому отличаются высокой степенью субъективности. Качество полученного результата напрямую зависит от личного опыта исследователя. Существуют единичные попытки использовать формальные критерии и геоинформацион-ные подходы (Богданова и др., 2012) . В настоящем исследовании мы используем автоматическое выделение ЭЛ уровня рода в соответствии с формальными критериями, приведенными в таблице 1 (Тимофеев, Энтин, 2018) .

Таблица 1. Критерии выделения родов элементарных геохимических ландшафтов. Источник: (Тимофеев, Энтин, 2018)

Table 1. Criteria for mapping elementary geochemical landscapes. Source: (Timofeev, Entin, 2018)

Расположение контуров	Крутизна, °	Кривизна, м-1	Элементарный ландшафт (ЭЛ)	Индекс
Водораздел	<2	≈ 0, >0	Автономный первого порядка (плоские водораздельные поверхности)	А1
Между поймой и водоразделом	<2	Любое значение	Автономный второго порядка (надпойменная терраса, четко выраженная в рельефе)	А2
Водораздел	<4.5	<0	Автономный второго порядка (надпойменная терраса, четко выраженная в рельефе)	А2
Внутри поймы	2–4.5	>0	Трансэлювиальный со слабым выносом	ТЭ1
Между поймой и водоразделом	2–4.5	Любое значение	Трансэлювиальный со слабым выносом	ТЭ1
Водораздел	2–4.5	≈ 0, >0	Трансэлювиальный со слабым выносом	ТЭ1
Между поймой и водоразделом, граничит с ТЭ1	4.5–10	Любое значение	Трансэлювиальный с умеренным выносом	ТЭ2
Между поймой и водоразделом, граничит с ТЭ1 или ТЭ2	10–20	Любое значение	Трансэлювиальный с активным выносом	ТЭ3
Между поймой и водоразделом, граничит с ТЭ1, ТЭ2 или ТЭ3	> 20	Любое значение	Трансэлювиальный с интенсивным выносом	ТЭ4
Внутри поймы	2–4.5	<0, ≈ 0	Трансаккумулятивный	ТА
Внутри поймы	<2	-, 0, +	Супераквальный	СА

Исходными данными для создания карты родов ЭЛ послужили цифровая модель рельефа (ЦМР) и векторные наборы данных об объектах гидрографической сети (р. Хара). В качестве ЦМР используется фрагмент глобального регулярно-сеточного набора данных о высотах земной поверхности SRTM (Farr et al., 2007). Шаг регулярной сетки (размер ячейки) этого набора составляет 1″, что примерно соответствует 30 м на местности. Фрагмент был трансформирован в систему координат UTM, а таже подвергнут фильтрации по методу Sun et al. (2007). Затем была выполнена передискретизация набора данных на сетку более высокого разрешения (10 м) – такая операция необходима для упрощения векторизации растра ЭЛ на заключительных этапах работы.

Набор пространственных данных об объектах гидрографической сети был создан путем визуального дешифрирования мозаик космических снимков сверхвысокого разрешения, полученных с космических аппаратов WorldView-1, WorldView-2, WorldView-3 и доступных в виде покрытия в сервисе Google Earth. Дешифрирование выполнено при помощи инструментария Google Earth, полученный набор содержал контуры объектов гидрографической сети полилинейной геометрии без дополнительной атрибутивной информации.

На основе ЦМР были рассчитаны растры крутизны и кривизны склонов. Расчет проводился в SAGA (Conrad et al., 2015) . Полученные значения были переклассифицированы в соответствии с таблицей 1. Следует заметить, что значение кривизны, рассчитанное по ЦМР, практически никогда не может быть строго равно 0 из-за особенностей машинных вычислений, а также вследствие наличия шума и артефактов интерполяции в исходных данных. Поэтому для отделения значений кривизны было выбрано пороговое значение, равное по модулю 0.0005 м^-1. Значения, меньшие по модулю, чем указанная пороговая величина, были приняты как приблизительно равные 0.

Границы поймы определялись путем выделения зоны подтопления по разности высот ячеек, соответствующих положению водотока, и прилегающих территорий. Разность высот была рассчитана в SAGA (инструмент Vertical Distance to Channel Network). Ячейки ЦМР, для которых эта разность оказалась ниже 2 м (известный по данным полевого обследования уровень подъема воды в р. Хара), маркированы как относящиеся к пойме. Для определения положений водоразделов была рассчитана общая водосборная площадь. Расчет выполнен в SAGA на основе алгорит- ма Multiple Flow Direction. Ячейки с низкими значениями водосборной площади помечались как потенциально относящиеся к водораздельным пространствам. Однако низкие значения водосборной площади, рассчитанной по ЦМР, характерны не только для водораздельных пространств, но и, вследствие особенностей процедуры расчета, к положительным формам рельефа в пойме реки. Поэтому в качестве дополнительного критерия использовалось максимальное значение крутизны склона в радиусе 5 ячеек ЦМР. Обычные значения крутизны склона для поймы не превышают 5°, поэтому ячейки, максимальная крутизна в окрестности которых превышала 5°, интерпретированы и промаркированы как водораздельные, для остальных ячеек эта маркировка снята. Полученные результаты положения пойм и водоразделов объединены в тематический слой с выделением трех классов: “пойма”, “промежуточное”, “водораздел”.

Для создания карты ЭЛ был выполнен оверлей четырех указанных выше классифицированных слоев. В результате был получен растр, отражающий уникальные комбинации возможных классов. Некоторые уникальные комбинации соответствуют одному и тому же роду ЭЛ: например, комбинации классов “внутри поймы” и крутизны менее 2° соответствуют супераквальным ландшафтам (СА) при любых возможных значениях кривизны. Поэтому результат оверлея был повторно переклассифицирован в соответствии с таблицей 1. После этого растр родов ЭЛ был конвертирован в векторный набор данных. На заключительном шаге полученный векторный набор был генерализован при помощи инструмента присоединения полигонов небольшой площади (порядка площади 1–2 ячеек исходного растра) к более крупным соседним контурам. Вся вышеописанная процедура проиллюстрирована на рисунке 3.

Для картографирования значений коэффициента концентрации использовалась площадная интерполяция (areal interpolation) по принципу “полигон-к-полигону” (Krivoruchko et al., 2011) .

Рис. 3. Методика построения карты родов элементарных геохимических ландшафтов (ЭЛ) и получения карт коэффициентов концентрации на ее основе.

Fig. 3. Map of elementary geochemical landscapes (genera) and map of concentration coefficients derivation methodic.

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

В настоящей работе в качестве исходных полигонов использовались полигоны Тиссена (диаграмма Вороного), построенные по точкам опробований; целевыми полигонами являлись контуры родов ЭЛ, рассчитанных на основе ЦММ.

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

Процедура повторяется для каждой точки в исходном наборе (Демьянов, Савельева, 2010) . При этом характеристики ошибок результатов интерполяции должны быть не хуже, чем аналогичные показатели для результатов интерполяции традиционными методами (например, ОВР или обычного кригинга) – которые, как было указано выше, заведомо непригодны для имеющегося набора данных. Поэтому одновременно с интерполяцией по предложенной методике была выполнена интерполяция традиционными методами и оценка этих результатов методом кросс-валидации.

РЕЗУЛЬТАТЫ И ОБСУЖДЕНИЕ

На рисунке 4 показана карта элементарных ландшафтов (ЭЛ), построенная на основе ЦМР. Видно, что контуры различных элементарных ландшафтов распределены по территории неравномерно. Супераквальные ландшафты (СА) образуют один непрерывный контур вокруг русла р. Хара. К нему примыкает несколько крупных контуров трансаккумулятивных ландшафтов (ТА). Автономные ландшафты (А1, А2) занимают небольшие площади вблизи водоразделов и не образуют связных, протяженных контуров. Также обращает на себя внимание относительно случайное чередование ТА и ТЭ1: небольшие по площади контуры одного из типов находятся внутри более крупного контура другого типа; граница между крупными контурами также довольно извилистая. Причина появления такой “мозаичной” картины состоит, вероят- нее всего, в сохранении на используемой ЦМР локальных особенностей поля высот, соответствующих формам рельефа разного размера и масштабного уровня.

Рис. 4. Карта элементарных геохимических ландшафтов (ЭЛ) уровня рода. Объяснения индексов приведены в таблице 1.

Fig. 4. Map of elementary geochemical landscapes (genera). See index explanations in table 1.

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

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

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

Видно, что, как и ожидалось, наиболее высокие концентрации загрязнителей наблюдаются в южной части изучаемой территории, ближе к основным источникам загрязнителей. Абсолютные значения полученных коэффициентов варьируют в разных пределах: наиболее высокие значения (более 5) наблюдаются для Cr, Cu и W, в то время как обычные значения для Cd и Pb варьируют от 1 до 2. В остальном же наблюдаемые различия показателей связаны, с различными источниками. As и Cr имеют достаточно однородное распределение со слабовыраженной аккумуляцией к северу от Дарханской ТЭЦ на склоне отрога хребта, отделяющего промышленную зону от селитебной. Cu, Sb, Cd имеют значительно более дифференцированное распределение. Кроме аномалии вблизи ТЭЦ, зоны аккумуляции сформировались в юрточных кварталах в северной части города. Сравнение карт As, Cr, Cu, Sb и Cd позволяет сделать вывод о том, что источником этих элементов являются выбросы золы при сжигании угля. Промышленное высокотемпературное сжигание высвобождает более широкий перечень загрязнителей, нежели низкотемпературное в юртах. В селитебной зоне с многоэтажной застройкой в центре города и в юрточном квартале на севере обнаружены высококонтрастные аномалии Pb вблизи дорог. Их основным источником являются продукты эксплуатации автотранспорта и в особенности этилированный бензин. Полученные данные подтверждают результаты ранее проведенных исследований (Kosheleva et al., 2019).

Оценка результатов интерполяции показателей, на основе которых были построены карты, представленные на рисунке 5, выполнялась методом кросс-валидации.

В таблице 2 приведены среднеквадратические ошибки (RMSE), рассчитанные для каждого EFi ТММ на основе методики интерполяции, описанной в настоящей работе. Для сравнения приводятся среднеквадратические ошибки, полученные для интерполяции того же набора данных детерминистическим методом (интерполяцией по метду обратных взвешенных расстояний, ОВР) и методом ординарного кригинга.

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

В одном случае (Sb) разница между рассчитанными значениями ошибок незначительна. Еще в одном случае (Cr) среднеквадратическая ошибка интерполяции, выполненной предлагаемым методом, оказалась хуже ошибки детерминистической интерполяции, но лучше ошибки интерполяции методом ординарного кригинга. Наконец, результат интерполяции As по формальным критериям уступает всем альтернативным подходам. Приведенные значения, на наш взгляд, не могут служить убедительным доказательством преимущества интерполяции с использованием данных о рельефе, но свидетельствуют об определенном потенциале такого подхода применительно к картографированию содержания ТММ.

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

Рис. 5. Карты коэффициентов концентрации ТММ.

Fig. 5. Maps of HMM concentration coefficients.

DFi

| более 6 ^J| 3 - 3,5 ■ 5.5-6  | 2,5-3

• ■    5-5,5  || 2-2,5

И 4,5-5     1,5-2

• ■    4-4,5     1-1,5

| 3,5-4 менее 1

О 1 2 3 км

Рис. 5 (продолжение). Карты коэффициентов концентрации ТММ.

Fig. 5 (continuation). Maps of HMM concentration coefficients.

Таблица 2. Среднеквадратические ошибки (RMSE) интерполированных значений коэффициентов концентрации. Полужирным начертанием выделены лучшие значения для каждого ТММ

Table 2. Root mean square errors (RMSE) of interpolated concentration coefficients values. Best value for each HMM is highlighted

ТММ	Предложенный метод	ОВР	Ординарный кригинг
As	1.00	0.91	0.88
Cd	1.12	1.33	1.15
Cr	1.13	1.07	1.19
Cu	1.09	1.19	1.15
Pb	2.53	2.61	2.76
Sb	1.32	1.33	1.31
W	1.77	1.78	1.71

В частности, процедура получения карты ЭЛ на основе ЦМР оставляет возможность для совершенствования за счет применения процедур генерализации и более дробного разделения крупных контуров. Для более точной оценки необходимо выполнение пробоотбора на содержание ТММ с шагом меньшим, чем используемый в настоящем исследовании.

ЗАКЛЮЧЕНИЕ

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

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