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

Камбалин Игорь Олегович аспирант; географический факультет; Московский Государственний Университет 119991, Россия, г. Москва, ул. Ленинские Горы, 1 Kambalin Igor Olegovich Postgraduate student; Faculty of Geography; Lomonosov Moscow State University 119991, Russia, Moscow, Leninskie Gory str., 1 igorkambalin@gmail.com Другие публикации этого автора     Кошурников Андрей Викторович ORCID: 0000-0001-6160-7795 доктор геолого-минералогических наук преподаватель; кафедра гляциологии и криолитологии; Московский государственный университет им. М.В. Ломоносова 119192, Россия, г. Москва, ул. Столетова, 9, кв. 12 Koshurnikov Andrey Viktorovich Doctor of Geology and Mineralogy Lecturer; Department of Glaciology and Cryolithology; Moscow State University 119192, Russia, Moscow, Stoletova str., 9, sq. 12 koshurnikov@msu-geophysics.ru Другие публикации этого автора

Балихин Ермолай Игоревич аспирант; геологический факультет; Московский государственный университет им. М.В. Ломоносова 119234, Россия, г. Москва, ул. Ленинские горы, 1 Balihin Ermolai Igorevich Postgraduate student; Faculty of Geology; Lomonosov Moscow State University 119234, Russia, Moscow, Leninskie gory str., 1 ermikus2@mail.ru Другие публикации этого автора

Введение

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

Параметризация среды в криолитозоне требует понимания направлений и характера изменчивости физических свойств пород. Повторно-жильные льды создают резкий контраст свойств по латерали, а пластовые льды обуславливают изменения по вертикали. Переходы между различными типами грунтов, например, от глины к скальным породам, также вносят вертикальные контрасты. Такие свойства, как влажность, засоленность и обводнённость – оказывают значительное влияние на физические характеристики среды ^[1]. Однако отсутствие единых подходов к параметризации физических свойств для вероятностного моделирования и интерполяции приводит к тому, что исследователи часто опираются на ограниченные или эмпирические данные ^[2-4]. Это снижает точность моделей и затрудняет экстраполяцию полученных результатов.

В таких условиях становится необходимым создание комплексных методов, объединяющих данные геофизических исследований и продвинутые статистические модели ^[5-7]. Современные технологии, включая методы трёхмерного статистического моделирования ^[8], позволяют более точно учитывать пространственную изменчивость данных. Однако успешность их применения зависит от правильного подбора параметров, что часто требует дополнительных исследований и тестирования на конкретных геологических условиях ^[9].

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

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

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

3. Изучить методы Байесовского кригинга и адаптировать их к специфике криолитозоны с формированием рекомендаций по их применению.

4. Проанализировать изменчивость физических свойств грунтов, выявив связь с геоморфологическими и литологическими особенностями.

5. Выполнить валидацию созданных моделей, используя данные буровых работ и полевых наблюдений ^[10].

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

В ходе работ проведена апробация метода трёхмерного Байесовского кригинга, применяемого для интерполяции геофизических данных. Исследования выполнены на территории шлаковых отвалов никелевого завода в г. Норильск, где ранее подобный метод не применялся. Использование метода трёхмерного Байесовского кригинга позволило выявить закономерности в пространственном распределении физических свойств многолетнемерзлых пород, а также изучить их изменчивость.

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

Обзор методов статистического моделирования

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

Методы статистического моделирования можно подразделять по характерным особенностям ^{[11, 12]}:

По типу входных данных:

- пространственные: кригинг, методы пространственной интерполяции (IDW - Inverse Distance Weighting, сплайн-функции);

- временные: анализ временных рядов;

- многомерные: кластерный анализ, регрессионные методы.

По целям анализа:

- интерполяция: кригинг, IDW, сплайны;

- классификация: Random Forest, кластерный анализ;

- регрессия: корреляционный и регрессионный анализ;

- оценка неопределённости: Bayesian Methods, Байесовский кригинг.

По алгоритмической сложности и требованиям к количеству данных для вычислений:

- простейшие: IDW, линейная регрессия;

- сложные: Gradient Boosting, Deep Learning, Байесовский кригинг.

Геофизические исследования в криолитозоне, осложненные неоднородностью среды (изменчивость структуры грунтов, вариации теплопроводности, сложные процессы замерзания-оттаивания, высокая электропроводность таликов вместе с повышенным электрическим сопротивлением льдистых пород и т.д.), наиболее интересны для статистического моделирования [Douglas, Zheng, Gruber, Jørgensen]. Это обусловленно в первую очередь большим количеством получаемых данных (в отличие от буровых), а также их неоднозначностью. То есть, полученные данные сами по себе могут описывать сразу несколько геокриологический ситуаций и для однозначной интерпретации необходимо производить дополнительные исследования.

Для геофизических исследований следует выделить следующие особенности и требования к статистическим моделям:

- геостатистические подходы должны учитывать анизотропию физико-механических свойств мерзлых грунтов ^{[3, 13]};

- методы прогнозирования должны адаптироваться для обработки данных с пропусками и высокой степенью шума;

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

Байесовский кригинг 3D, поставляемый в программном комплексе ArcGIS Pro, является наиболее подходящим методом для работы с геофизическими данными в криолитозоне, так как удволетворяет вышеуказанным требованиям.

В сравнении с другими геостатистическими методами, такими как IDW, сплайны или классический кригинг, Байесовский кригинг 3D обеспечивает более точное моделирование сложных геофизических процессов засчёт учета неопределенности ^[14].

В то же время, от методов сложного машинного обучения, таких как Random Forest или Gradient Boosting, Байесовский кригинг отличается встроенными механизмами для учёта пространственной корреляции и анизотропии данных. Методы машинного обучения также могут учитывать пространственные отношения входных данных, особенно алгоритмы Deep Learning, но это требует большого объёма обучающих данных и встраивания энергозатратных модулей, рассчет которых занимает больше времени.

Определение параметров построения модели методом Байесовского кригинга 3D

В ArcGIS Pro Байесовский кригинг в 3D предоставляет пользователю набор параметров, которые можно настроить для оптимизации итоговой модели. Основные настройки включают:

1. Выбор поля для учёта ошибки измерений – при измерении одного и того же параметра разным оборудованием, для которого известны ошибки измерения;

2. Модель семивариограммы, определяющей форму пространственной корреляции исходных данных:

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

- linear – линейная модель предполагает пропорциональное увеличение дисперсии с расстоянием. Применяется для данных с линейной корреляцией;

- thin Plate Spline – эта модель подходит для данных с высокой гладкостью. Часто используется для интерполяции сложных поверхностей;

- exponential – экспоненциальная модель хорошо работает с данными, где корреляция быстро спадает с увеличением расстояния;

- whittle – модель, которая используется для данных с более сложной пространственной структурой, включая анизотропию;

- K-Bessel – модель, применяемая для описания сложных корреляционных структур, которые включают цикличные или повторяющиеся процессы.

3. Тип трансформации – позволяет учитывать наличие выбросов (значений резко отличающихся от основного массива данных) в исходных данных:

- none – используется, если исходные данные не имеют выбросов или были предварительно обработаны;

- emprical – сглаживает массив данных, если в нём присутствуют выбросы;

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

- subset size – устанавливает размер подвыборки (локальной модели) данных, используемой для расчёта локальных параметров интерполяции. Увеличение размера подвыборки уменьшает количество подвыборок, что может повысить устойчивость модели, так как большее количество точек снижает влияние выбросов и шума. Однако это может сгладить локальные вариации. Уменьшение Subset Size делает модель более адаптивной к локальным особенностям, но увеличивает её чувствительность к шуму и нестабильность. Однако уменьшение размера подвыборки приводит к увеличению количества таких подвыборок, что может увеличить время вычислений.

4. Local model area overlap factor – определяет степень перекрытия между подвыборками. Каждая точка данных может попадать в несколько подвыборок. Более высокие значения overlap factor создают более сглаженную выходную поверхность, так как каждая точка обрабатывается несколько раз. Однако увеличение overlap factor приводит к росту времени вычислений. Значения параметра должны быть в пределах от 1 до 5. Реальное перекрытие обычно больше указанного значения, чтобы обеспечить одинаковое количество точек для каждой локальной модели.

5. Number of simulated semivariograms – количество семивариограмм строящихся для каждой локальной модели. Чем больше это число, тем более стабильной будет модель, однако это значительно увеличивает время расчета. Подбор этого параметра в первую очередь зависит от зашумленности данных и обычно хватает от 100 до 500 симуляций. Для того чтобы проверить, достаточно ли выбранное количество, можно в построенной модели оценить ошибку предсказания. Так, если дисперсия исходных данных больше среднеквадратичной ошибки модели на порядок, то часто можно считать модель верной.

6. Order of trend removal – позволяет удалить тренд изменения значений по вертикали. Такая операция часто требуется, если измерения по вертикали производятся значительно более часто, чем по горизонтали, а значит, тренды оказывают большее влияние на общее сглаживание модели по вертикали. Для анализа именно локальных изменений, такого рода тренды можно удалить, что позволит получить большую вариативность значений по вертикали:

- none – тренд не удаляется, этот параметр используется по умолчанию;

- first order – удаляет линейный тренд первого порядка в вертикальном направлении, что позволяет стабилизировать расчёты и уменьшить влияние сильных вертикальных изменений значений. Такой подход полезен для данных, которые изменяются по вертикали быстрее, чем по горизонтали.

7. Elevation inflation factor – определяет масштабирование значений в вертикальном направлении перед расчётом модели и интерполяцией. Этот параметр корректирует разницу в изменчивости значений по вертикали и горизонтали, делая один единичный шаг по высоте эквивалентным одному шагу по горизонтали. Более высокие значения параметра усиливают влияние вертикальной структуры данных, что стабилизирует расчёты и улучшает точность интерполяции. Значение по умолчанию рассчитывается автоматически с использованием метода максимального правдоподобия и обычно находится в диапазоне от 1 до 1000. Пользователь может задать собственное значение для уточнения модели на основе кросс-валидации.

8. Search neighborhood – определяет параметры поиска соседей для интерполяции. У этого параметра есть ряд настроек, которые также следует рассмотреть:

- max neighbors – максимальное количество соседних точек, используемых для расчёта значений. Увеличение этого параметра повышает устойчивость, но может сгладить локальные детали;

- min neighbors – минимальное количество соседей, необходимое для выполнения расчёта. Низкие значения повышают гибкость модели, но могут ухудшить стабильность в разреженных данных.

Sector type: Делит пространство вокруг точки на секторы (платоновы тела):

1) 1 Sector (Sphere): Все ближайшие соседи из любого направления. Простое решение для равномерных данных.

2) 4 Sectors (Tetrahedron): Деление на четыре сектора.

3) 6 Sectors (Cube): Деление на шесть регионов.

4) 8 Sectors (Octahedron): Деление пространства на восемь секторов.

5) 12 Sectors (Dodecahedron): Деление на 12 секторов для равномерного покрытия пространства.

6) 20 Sectors (Icosahedron): Деление на 20 секторов для более детализированного поиска соседей.

Рис. 1. Пример деления области поиска на 4 сектора [Esri. 3D Search Neighborhoods. URL: https://pro.arcgis.com/ru/pro-app/latest/help/analysis/geostatistical-analyst/3d-search-neighborhoods.htm (дата обращения: 10.12.2024)]

Во время настройки самой модели можно выбрать то, как будут обрабатываться точки с совпадающими координатами, что может быть полезно в случае проведения двойного измерения для проверки получаемых данных. В настройке среды Coincident Points можно выставить следующие параметры:

- mean of values at coincident locations (по умолчанию) – для всех совпадающих точек на одном местоположении вычисляется среднее значение;

- еxclude all coincident data – полностью исключает все совпадающие точки из анализа;

- minimum value at coincident locations – использует минимальное значение из всех совпадающих точек на данном местоположении;

- include all coincident data – все данные на совпадающих точках включаются в анализ без усреднения или исключения.

Подбор параметров для Байесовского кригинга 3D на примере данных частотного зондирования

Для отработки метода подбора данных используются данные об удельном электрическом сопротивлении, полученные методом частотного зондирования на полигоне близ террикона шлаковых отвалов города Норильск (рис. 2). Данные включают в себя 11 профилей, расстояние между которыми 100 метров, а частота съемки по профилю – каждые 2 метра. Суммарное количество точек – 69 076. В ходе обработки данных каждой точке были назначены координаты XY по встроенному в оборудование GPS трекеру, а координата Z назначалась как разница интегральной глубины замера и высоты точки на цифровой модели рельефа. Цифровая модель строилась по заранее классифицированным точкам, с исключением точек высокой и низкой растительности, а также водоёмов.

Рис. 2. Контур участка изысканий

В пределах исследуемой территории распространены повторно-жильные льды, преимущественный состав отложений представлен тундровыми глеево-торфянистыми и торфяными (глееземы торфянистые и торфяные тундровые) почвами, залегающими на глинистых и тяжелосуглинистых щебнистых породах ^{[15, 16]}.

На полигоне проведены исследования, направленные на картирование техногенного загрязнения металлическими частицами, сдуваемыми с террикона и переотлагающимися на участке ^[9].

Для того, чтобы построить пространственную модель распределения удельного электрического сопротивления, сформируем требования к искомой модели с учётом имеющихся данных о среде:

1. Измеряемые данные, исходя из знаний о среде и предварительному анализу точек на этапе камеральной обработки, подвержены значительной анизотропии свойств. Это объясняется в первую очередь тем, что электрическое сопротивление зависит как от литологического состава (и наличия металлического загрязнения), так и от фазового состояния воды в грунте. Если первое имеет наибольшую изменчивость по латерали, то второе обладает наибольшей изменчивостью с глубиной. Необходимо выбрать для моделирования тип вариограмм, подходящий для данных с выраженной анизотропией, а именно – Whittle.

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

3. Данные расположены достаточно близко друг к другу в пределах профиля, но не между профилями, и сеть не является равномерной. Примем значение масштабирования высот (Elevation Inflation Factor) равным 10, что отражает отношение расстояния между профилями к глубинности метода. В то же время, радиус поиска соседей необходимо взять больший, чем расстояние между профилями в 1-1.5 раза.

4. Данные не являются однородными во всех направлениях из-за не изометричной сети. Для учёта этого факта следует взять Sector type не меньший чем Cube (более 4 секторов).

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

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

Таблица 1

Таблица 2

На основе полученных 3D воксельных моделей распределения полей удельного электрического сопротивления создан комплект карт и разрезов (рис. 3-7), на которых полученные поля проинтерпретированны и апробированны при помощи заверочных буровых работ. На представленных ниже изображениях показаны результаты трёхмерного статистического моделирования распределения геоэлектрических свойств многолетнемерзлых пород с использованием различных параметров Байесовского кригинга.

Рисунок 3. Модель (№1) распределения геоэлектрических свойств многолетнемерзлых пород. Карта и разрез

Рисунок 4. Модель (№2) распределения геоэлектрических свойств многолетнемерзлых пород. Карта и разрез

Рисунок 5. Модель (№3) распределения геоэлектрических свойств многолетнемерзлых пород. Карта и разрез

Рисунок 6. Модель (№4) распределения геоэлектрических свойств многолетнемерзлых пород. Карта и разрез

Рисунок 7. Модель (№5) распределения геоэлектрических свойств многолетнемерзлых пород. Карта и разрез

Рассмотрим полученные результаты подробнее, сравнивая различные варианты моделей и их соответствие фактическим данным бурения по скважинам №5, 6 и 8.

Скважина №5 расположена на северном склоне естественного возвышения, что обеспечивает стабильное промерзание пород. В керне представлены лёгкие суглинки и тяжелая супесь с массивной криоструктурой, обусловливающие высокое электрическое сопротивление.

Модели №4 и №1 не обнаруживают талых пород до глубин 20 метров, хоть и позволяют судить о некотором понижении сопротивления в окрестностях.

Модели №5 и №2 хоть и отражают положение талой толщи достаточно близко к действительности – являются недостаточно достоверными. У модели №2 значения по мощности талого слоя оказались завышенными, и кровля талой толщи обнаружена на 0,4 м выше, чем в действительности. А в модели №5 границы получились сглаженными и отражающими более не естественное состояние среды, но переходы между отдельными областями.

Модель №3, в свою очередь, наиболее точно соответствует фактическим данным керна, показывая талый слой уже на глубине около 17,5 м, подтверждая высокую точность экспоненциальной семивариограммы и адекватность параметров поиска соседей (15 точек, 6 секторов). Более того, характер распределения поля в окрестностях скважины №5 более четко отображен именно в рамках модели №3, отражая структуру согласованную с экспозиционными условиями.

Среда в пределах скважины №8 также характеризуется устойчиво высоким сопротивлением почти на всей глубине, обусловленным наличием мерзлых легких суглинков и супесей. Однако, вследствие меньшей льдистости, контрастность среды около этой скважины ниже, чем около скважины №5. Ни одна из моделей не может демонстрировать абсолютную точность в полученных результатах. Однако, наиболее приближенно к действительности получилось оконтурить границу талого слоя в рамках модели №3 на глубине 18,4 м, что лишь немного не соответствует фактическому положению границы (18,7 м). Остальные модели демонстрируют либо заниженные значения – модель №1, либо завышенные – модели №4 и 2, отражая границу на глубине всего 16,0-16,5 м. Модель №5 галлюцинировала возникновение высокоомного слоя под мерзлым на глубинах 18,0-20,0 м.

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

При дальнейшем анализе выявлено, что бурый окрас отложений, и аномально низкое удельное электрическое сопротивление, говорит о захоронении в этом месте металлических отходов. Данная особенность чётко зафиксирована всеми моделями, однако модель №3 снова демонстрирует наиболее детальное и структурированное отображение низкоомного слоя, что важно для оценки масштабов загрязнения.

Преимущества использования модели №3:

1. Высокая точность отображения подошвы мерзлых пород и наличия талых слоев, подтверждаемая керновыми данными.

2. Хорошее соответствие фактическим наблюдениям по границам зон различной омности.

3. Устойчивость к влиянию шума и локальных аномалий.

Недостатки остальных моделей:

1. Модель №1 демонстрирует слабое соответствие фактическим данным в зонах резких границ.

2. Модели №2 и №4 менее детализированы и часто сглаживают важные детали.

3. Модель №5, несмотря на высокую детализацию, показывает повышенную чувствительность к шумам и нестабильность в отражении мелкомасштабных аномалий.

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

Более того, применение метода подбора параметров моделирования исходя из особенностей геокриологического строения действительно позволяет добиться наиболее точных результатов при моделировании, а неправильный выбор, обусловленный недостаточным изучением среды приводит к ошибкам, которые, потенциально, могут пагубно отразиться в ходе инженерно-геологических работ на территории ^{[17, 18]}.