Трехмерное геологическое моделирование верхней части подземного пространства и его использование на примере исследования в районе Тунчжоу, г. Пекин, Китай

ВВЕДЕНИЕ

Быстрое освоение городского подземного пространства (ГПП) требует понимания его объемного строения, особенно в верхней части. Его геологическая среда включает дисперсные и скальные грунты, воду и газы [1]. Сложность и изменчивость геологических условий являются важнейшими источниками неопределенностей и потенциальных рисков при будущем строительстве в ГПП. Исследования геологических условий верхней части городского подземного пространства постоянно развивались в течение последних 20 лет [2–6]. Однако детальное определение его литолого-стратиграфического строения по-прежнему остается сложной задачей из-за неоднородности свойств неглубоко залегающих геологических слоев, которые имеют отношение к строительству городской инфраструктуры, в том числе инженерных сетей [7, 8]. Хотя методы геологических и геотехнических исследований позволяют получать подробные данные в отдельных точках [8], свойства грунтов могут сильно изменяться в вертикальном направлении. Эффективным методом описания изменчивости свойств грунтов в горизонтальном и вертикальном направлениях является трехмерное (3D) геологическое моделирование. Построить 3D геологическую модель для городской территории довольно сложно.

В последние десятилетия 3D геологическое моделирование быстро развивается в качестве эффективного инструмента визуализации ГПП [8–20]. Этот метод может улучшить представление о строении подземной среды [21–26] и описать свойства, которые неоднородно распределены внутри геологических тел [17, 27–29].

Кроме того, для оценки ресурсов [6, 9], изучения геологических условий [30–33] и геотехнических свойств грунтов [7, 10, 34, 35] было построено множество 3D геологических моделей, способствующих пониманию устройства верхней части подземного пространства. Так, исследователи из Британской геологической службы [2, 5, 12, 36, 37] построили множество 3D геологических моделей всех масштабов (от площадок до городов и до сухопутной и шельфовой частей Великобритании), используя различные программные инструменты и подходы к моделированию [38]. Хоу и др. [9] создали 3D геологическую модель для описания ГПП города Фошань в китайской провинции Гуандун. Йе и др. [33] построили 3D модель осадок грунтов для моделирования смещений системы подземных водоносных горизонтов в центре г. Шанхая (Китай). Хёйер и др. [39] разработали 3D геологическую модель для территории г. Самсё (Дания) с высоким разрешением для обновления оценки рисков для свалки Пиллемарк. Андерсон и др. [3] построили комплексную 3D геологическую модель для территории г. Вайле (Дания) в целях поддержки градостроительства путем предоставления карт планирования. Чэнь и др. [40] создали комплексную 3D модель, разработанную на основе методов многоточечной статистики, для территории эстуария реки Миньцзян (Китай) в целях более точной визуализации геологической среды. В целом, для улучшения понимания геологии в различных масштабах 3D геологические модели используются все чаще.

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

Основная цель данного исследования – разработка эффективного и практичного процесса 3D геологического моделирования для условий речных отложений на урбанизированной территории и обсудить возможности ее применения. В качестве примера был взят район Тунчжоу г. Пекина, для которого были проанализированы данные по 40 новым инженерно-геологическим скважинам в сочетании с 1506 файлами данных по ранее выполненным изысканиям в целях обобщения информации по литологическим условиям – необходимого этапа реализации метода вертикальных разрезов. Была подробно рассмотрена процедура использования этого метода. И наконец, авторы разработали процесс моделирования и комплексную 3D модель подземного пространства района Тунчжоу до глубины 50 м. Эта новая модель была использована в качестве инструмента для визуализации изменчивости геологических единиц по глубине и геометрии и для характеристики множества свойств (например, модуля всестороннего сжатия, анализируемого в этой статье) каждой единицы, а затем для оценки геологических условий. Полученные результаты в сочетании с картой изопахит (равных мощностей) каждой геологической единицы, полученной на основе созданной модели, могут улучшить понимание геологических условий ГПП и предоставить ценную информацию для снижения рисков при строительстве в будущем. Еще одним ключевым направлением использования достигнутых результатов является создание динамической модели мониторинга, учитывающей время, однако для ее разработки необходимы дополнительные исследования.

ЧЕТВЕРТИЧНАЯ ГЕОЛОГИЯ

Верхняя часть подземного пространства под Пекинской равниной включает отложения нескольких огромных аллювиальных конусов выноса и осадочных впадин. Район Тунчжоу г. Пекина расположен на юго-востоке Пекинской равнины, причем его рельеф является немного выше на северо-западе и ниже на юго-востоке. На северо-запад текут три реки (Вэньюй, Чаобай и Ляншуй), которые сливаются в юго-восточной части исследуемого района (рис. 1). По всей территории широко распространены четвертичные отложения, которые, как правило, залегают на коренных породах. Под урбанизированной частью голоценовые и позднеплейстоценовые отложения залегают на неогеновых толщах и локально – на палеогеновых (в Юго-Восточной впадине). Район Тунчжоу является зоной пересечения двух крупных аллювиальных конусов выноса рек Юндин и Чаобай [41]. Наиболее ранним и крупным является аллювиальный конус выноса реки Юндин. Его вершина расположена в районе Шицзиншань, простирается он в направлении СЗ-ЮВ, а его литолого-стратиграфическое строение меняется от однослойных структур, сложенных гравием, на западе до многослойных структур, сложенных пылеватыми глинами, глинами, песками и гравием, на востоке [41] (здесь и далее названия грунтов переведены дословно, поскольку они были даны авторами в соответствии с унифицированной системой классификации грунтов USCS, которая является более градуированной и несколько смещенной по сравнению с классификацией по российскому ГОСТ 25100-2020. – Ред.).

Рис. 1. Карта распределения аллювиальных конусов выноса на территории Пекинской равнины (адаптировано по [41]) (а); распределение ранее существовавших и новых скважин на исследуемой территории (б); скважины на изучаемой территории, выбранные для моделирования (в)

Чаобай, вторая по величине река в г. Пекине, берет начало в горах Яньшань и течет на юг, образуя обширный аллювиально-пролювиальный конус выноса (см. рис. 1). Непрерывное поднятие гор и интенсивное опускание равнин привели к наличию мощных осадков. Максимальная мощность четвертичных отложений наблюдается в южной части, где она превышает 450 м [42]. Это объясняет большую изменчивость последовательностей напластования по всей этой территории.

ИСХОДНЫЕ ДАННЫЕ И МЕТОДЫ ИССЛЕДОВАНИЯ

Геологические данные

В общей сложности для исследуемой территории было собрано 1506 отсканированных файлов с данными исследований по скважинам и с соответствующими геотехническими отчетами. Из них 614 файлов поступили из архивов городского строительства района Тунчжоу, а остальные 892 файла были получены в Институте архитектурных обследований и проектирования района Тунчжоу. Эти скважины по большей части были неглубокими (30 м или менее) и располагались вокруг местных строительных площадок в северной части рассматриваемой территории. Некоторые из них были относительно плохо литологически описаны. Для получения дополнительной информации по литологии и определения типичных геологических слоев для моделируемой территории было пробурено 40 новых более глубоких скважин в южной части района Тунчжоу (см. рис. 1 а, б). Максимальная глубина некоторых из них достигла 100 м (например для скважины 1-11 на рисунке 1, a). Всего для моделирования были отобраны данные по 476 скважинам (см. рис. 1, в), которые обычно имели глубину более 30 м и были подробно описаны с предоставлением информации о местоположении, идентификации, глубине, высотных отметках, координатах, а также о кровле, подошве и литологии слоев.

В процесс моделирования в качестве важного фактора также были включены разломы. Всего было рассмотрено пять разломов, верхние точки которых располагались не глубже 50 м, – разломы Нанькоу-Суньхэ (ns на рисунке 1, а), Наньюань-Тунсянь (nt), Чжанцзявань (zjw), Сядянь (xd) и Яосиньчжуан (yxz) [43–45]. Данные по разломам были взяты из технического отчета проекта под названием «Исследование основных геологических проблем в субцентре района Тунчжоу», который был выполнен в 2016–2017 годах Пекинским институтом геологических исследований (изысканий). Однако мало где были точно задокументированы точные местоположения верхних точек разломов в неглубоко залегающих четвертичных отложениях (менее 50 м), поэтому плоскость сместителя моделировалась на основе средней амплитуды вертикального смещения по разлому, которая составляла около 1–2 м.

Методы

Процесс трехмерного моделирования

В данной работе для построения 3D геологической модели были использованы как явные, так и неявные методы. Моделирование проводилось в программном пакете SKUA-GoCAD, применявшемся как к тем, так и к другим методам, с использованием последовательности операций Structure and Stratigraphy и Reservoir Properties («Структура и стратиграфия» и «Характеристика мощностей и свойств слоев»), а также в программе Creatar Xmodelling для явного моделирования четвертичных отложений до небольшой глубины (рис. 2). Явное моделирование, например метод вертикальных разрезов (см. следующий параграф), подразумевает использование геологами геологических концепций и их взаимодействий, в то время как неявный подход относится к применению программ со встроенными математическими функциями, основанными на геологических концепциях [38, 46].

Рис. 2. Процесс моделирования

Метод вертикальных разрезов

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

Рис. 3. Исходные скважины и обобщенное интерполированное литолого-стратиграфическое строение по линии разреза WE08 (названия грунтов переведены дословно, поскольку они были даны авторами в соответствии с унифицированной системой классификации грунтов USCS, которая является более градуированной и несколько смещенной по сравнению с классификацией по российскому ГОСТ 25100-2020. – Ред.)

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

1. Преобразование данных в цифровую форму. Информация по уже существующим скважинам была оцифрована на основе исходных бумажных копий документов. Сбор данных был выполнен в соответствии с «Нормами оцифровки текстово-графических геологических данных в сканированном виде» (SZ1999001-2000), выпущенными Информационным центром Министерства земель и природных ресурсов Китайской Народной Республики (ICMLR), и с «Техническими требованиями к созданию национально значимой базы данных по геологическим скважинам», составленными Геологической службой Китая (CGS). Для точного выполнения оцифровки данных необходимо было выполнить следующие шаги:

• проверить исходные бумажные копии документов и новую информацию по скважинам; описать конкретное расположение различных слоев, включая глубину залегания их кровли и подошвы, координаты (x, y, z) и т.д.;
• сравнить каждый случай зафиксированной литологической информации с данными по другим скважинам, расположенным вокруг, особенно в непосредственной близости; если искомый тип грунта, зафиксированный в одной скважине, совпадает с таковым в соседних скважинах на аналогичной глубине, то для искомого слоя следует использовать название этого грунта или рассмотреть возможность его залегания в виде линзы;
• выполнить проверку атрибутивных данных – дополнительную проверку геотехнических параметров, которые по большей части относятся к типам грунтов; свойства каждого типа грунта обычно имеют определенный диапазон значений атрибутов, например природное содержание воды в пылеватой глине обычно составляет от 15 до 40%, а предел пластичности этого грунта составляет от 12 до 23%; поэтому величины показателей свойств можно использовать для дополнительной проверки соответствия установленной литологии; зарегистрированные данные, которые показались в той или иной степени ненадежными, были помечены и затем отброшены.

2. Обобщение данных. Зарегистрированная информация по типам грунтов иногда могла быть противоречивой, например в одном отчете был указан суглинок (sub-clay), а в другом – пылеватая глина (silty clay) или песчанистая глина (sandy clay); кроме того, были довольно распространены переходные типы грунтов, такие как пылеватый песок (silty sand) частично с пылеватой глиной (silty clay) (см., например, скважины на рисунке 3); поэтому литологическую информацию обобщали до того, как соединять между собой разрезы скважин.

В этом исследовании обобщение привело к выделению в общей сложности шести категорий грунтов, таких как насыпной грунт (artificial ground), глинисто-пылеватый грунт (clayey silt), пылеватая глина (silty clay), песок от пылеватого до мелкозернистого, песок от средне- до крупнозернистого/гравия. Такие типы грунта, как пылеватая глина (silty clay), тяжелая пылеватая глина (heavy silty clay) и суглинок (sub-clay), на основе их описаний были объединены в категорию «пылеватая глина» (silty clay); песчанисто-пылеватый грунт (песчанистый алеврит, sandy silt), пылеватый грунт (silt) и мелкозернистый песок – в категорию «песок от пылеватого до мелкозернистого»; гравелистый песок, крупнозернистый песок, песок от средне- до крупнозернистого и гравий – в категорию «песок от средне- до крупнозернистого или до гравия».

Затем стратиграфическим слоям были присвоены коды в виде уникальных номеров. Это было сделано путем сопоставления литологических данных по скважинам с известной информацией по региональной стратиграфии и «назначения» слоев в зависимости от мощности и последовательности напластования описанных отложений в каждой скважине. Например, в результате такого обобщения для скважины 043542-1 были «назначены» следующие слои грунта (сверху вниз): насыпной грунт (слой 10), глинисто-пылеватый грунт (слой 22), песок от пылеватого до мелкозернистого (слой 24), глинисто-пылеватый грунт (слой 32), пылеватая глина (слой 33), песок от пылеватого до мелкозернистого (слой 34), пылеватая глина (слой 43), песок от пылеватого до мелкозернистого (слой 44), пылеватая глина (слой 53) и две линзы (cl и fn) (рис. 4).

Рис. 4. Литологическое обобщение на примере скважины 043542-1

3. Соединение выделенных категорий слоев. Ключевым шагом в этой процедуре является выбор маркирующего слоя. В качестве типичного маркирующего для исследуемой территории был выбран слой крупнозернистого песка или гравия. Были выбраны следующие принципы соединения слоев в практической работе (они уже применялись при моделировании подземного пространства центральной части г. Пекина [47]).

• Соединение основных/главных слоев. Под основным/главным слоем понимается мощная единица, которая широко и непрерывно распространена в пределах рассматриваемой территории. Для изучаемого района между скважинами были соединены слои насыпного грунта, глинисто-пылеватого грунта, пылеватой глины, песка от пылеватого до мелкозернистого и песка от средне- до крупнозернистого или до гравия.
• Присвоение кодов слоям и линзам. Код слоя или линзы отражал порядок расположения и тип слагающего его грунта в модели и был уникальным.
• Одинаковая в вертикальном направлении последовательность одинаковых слоев, зарегистрированная для разных скважин, расположенных на одной и той же линии вертикального разреза (например, скважин по линии WE08 на рисунке 5), по существу, определяла одну и ту же группу слоев, при условии что уклон между скважинами для нее был небольшим.
• Выделение линз. Линзу выделяли, если ее мощность превышала 0,5 м и она выклинивалась с двух концов диаметра скважины или на 1/2 общего расстояния между скважинами. Тонкие слои мощностью менее 0,2 м, но с типичными стратиграфическими особенностями (например в последовательности слоев глин), имевшие значение для моделирования, также сохранялись в модели.

Рис. 5. Линии разрезов, разломы и пример геологических единиц для исследуемой территории: a – линии разрезов и скрытые активные разломы с глубиной верхней точки разлома менее 50 м (разломы Наньюань-Тунсянь (nt), Нанькоу-Суньхэ (ns), Яосиньчжуан (yxz), Чжанцзявань (zjw), Сядянь (xd)); б – геологические единицы, развитые вдоль разреза NS12

Затем были построены вертикальные разрезы по линиям, которые проходили примерно в ортогональных направлениях в соответствии с распределением скважин и осадочных фаций (см. рис. 5). Для исследуемой территории был создан «каркас» из 44 разрезов расположенных с шагом 1–3 км (иногда и больше) друг от друга. Из них 26 линий разрезов (от WE1 до WE26) проходили с востока на запад и 18 (от NS01 до NS18) – с севера на юг (см. рис. 5, а). На рисунке 5, б для примера представлен разрез NS12, проходящий с севера на юг через 35 скважин.

Кригинг

При пространственном моделировании применялся обычный (ordinary) кригинг – алгоритм оценки с взвешиванием по расстоянию, при котором используется вариограмма для характеристики изменчивости в пределах рассматриваемой территории и для оптимизации весов, присваиваемых точкам данных [48–50]. Обычный кригинг дает оценку переменной Z(x) в месте, где не было наблюдений, на основе ее средневзвешенного значения по результатам наблюдений в соседних местах. Соответствующая теория может быть описана путем рассмотрения внутренней случайной функции, обозначаемой Z(x_i), где x_i представляет собой локации выборки [51]. Оцененное средневзвешенное значение Z'(x₀), получаемое с помощью предиктора обычного кригинга для неисследованного (без измерений) места, определяется следующим образом:

где Z(x_i) – значение переменной Z в точке выборки x_i; λ_i – вес, назначенный для каждого значения Z(x_i); n – количество используемых значений.

Чтобы обеспечить несмещенность оценки, сумма весов всех значений при обычном кригинге равна единице:

До проведения обычного кригинга должна быть выполнена оценка вариограммы, описывающей пространственную структуру распределения данных с использованием геостатистической характеристики [49, 53, 54]. Такая вариограмма определяется как 1/2 ожидаемого квадрата разницы между парами значений в точках Z(x) и Z(x+h), находящихся на расстоянии h друг от друга [55, 56]:

Получить надежную вариограмму по всей рассматриваемой территории было невозможно из-за большой изменчивости свойств в различных геологических единицах. Поэтому в данном исследовании использовали множество вариограмм, чтобы лучше охарактеризовать неоднородно распределенные свойства путем анализа параметров вариограммы для каждой единицы по отдельности с помощью инструмента Variogram Analyzer («Анализатор вариограмм») в программном пакете SKUA-GoCAD.

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

Стратиграфия

Для построения разрезов было выделено шесть категорий грунтов (насыпной грунт, глинисто-пывлеватый грунт, глина, пылеватая глина, песок от пылеватого до мелкозернистого, песок от средне- до крупнозернистого или до гравия). Для построения 3D структурной модели (модели слоев) были интегрированы данные по 476 скважинам и 44 разрезам. В эту модель вошли 10 стратиграфических границ (маркеров) до глубины 50 м от поверхности земли (рис. 6–9). Сплошность единиц отложений локально прерывалась либо активными разломами, либо эрозионными разломами, либо эрозионными границами.

Рис. 6. Трехмерная структурная геологическая модель (модель слоев) для территории района Тунчжоу г. Пекина

Рис. 7. Распределение верхнего слоя песка от пылеватого до мелкозернистого и включенных в него линз

Рис. 8. Схема разрезов для первой циклотемы в пределах исследуемой территории

Рис. 9. Распределение и мощность геологических единиц в пределах исследуемой территории в виде схем разрезов для: а – второй циклотемы; б – третьей циклотемы

Отложения в районе Тунчжоу до глубины 50 м были в целом разделены на 10 основных/главных геологических слоев. Затем эти слои были сгруппированы в три циклотемы аллювиальных отложений (см., например, рис. 4). Трехмерная модель показала пространственные структуры, включающие циклические последовательности из переслаивающихся глин, пылеватых грунтов, песков, гравия и линз (см. рис. 7–9). (Далее будут описаны смоделированные слои с названиями, обобщенными так, как было описано ранее. – Ред.)

Отложения примерно до глубины 17 м от поверхности (то есть первая циклотема) включают слой насыпного грунта и переслаивающиеся единицы пылеватых грунтов и глин, песков и гравия. Мощность гравия, локально встреченного на севере изучаемой территории, составляет около 3 м, но он отсутствует в центральной и южной частях этого района. Выше гравия по всей территории наблюдается слой пылеватого или мелкозернистого песка мощностью около 9 м, включающего полевой шпат, кварц, слюду, оксиды железа и совсем немного органических материалов. В юго-восточной части имеется прослой пылеватой глины (см. рис. 7). Песок в некоторых местах перекрывается последовательностью слоев пылеватой и тяжелой пылеватой глины. Выше залегает верхний слой глинисто-пылеватого грунта мощностью до 7 м (в среднем 3 м), содержащий слюду, оксиды железа и бирюзу, а внутри него локально встречаются линзы пылеватой глины. Еще выше по всей площади имеется слой насыпного грунта мощностью 1–2 м (рис. 8).

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

Третья циклотема характеризуется широким распространением по всей исследуемой территории единиц пылеватых глин и мелкозернистых песков (см. рис. 9, б). Первые в основном включают коричневато-желтые или серовато-желтые пылеватые глины от среднеплотных (medium) до плотных/полутвердых (stiff) мощностью от 1 до 10 м (в среднем 4,7 м), часто содержащие слюду, бирюзу, а также оксиды железа и органические материалы. А вторые включают коричнево-желтые или серые мелкозернистые пески, плотные и водонасыщенные, содержащие полевой шпат, кварц, слюду, а также органические материалы, со средней толщиной около 9,6 м. В этой мощной единице также изредка встречаются линзы гравия, крупнозернистого песка или глинисто-пылеватого грунта. На рисунке 10 показано распределение слоев вдоль разреза NS12.

Рис. 10. Распределение слоев на основе построенной 3D геологической модели: а – 3D схема разрезов; б – распределение слоев вдоль разреза NS12

Трехмерная модель свойств

Для построения моделей свойств (таких физических и геотехнических параметров, как природная влажность w(%), начальный коэффициент пористости/пустотности e и модуль всестороннего сжатия E_s) на основе рассмотренной выше 3D структурной модели использовался инструмент Reservoir Properties («Характеристика мощностей и свойств слоев») в программном пакете SKUA-GoCAD. Для примера на рисунке 11 показано распределение значений предела текучести на основе построенной 3D модели свойств.

Рис. 11. Распределение значений предела текучести грунтов Wl на основе трехмерной атрибутивной геологической модели для района Тунчжоу г. Пекина

Влажность w(%) и начальный коэффициент пристости/пустотности e

Влажность грунтов до глубины 50 м в пределах исследуемой территории (по результатам лабораторных исследований образцов) варьирует от 6,5 до 53,1% при среднем значении около 24,8% без каких-либо трендов, связанных с глубиной или литологией (рис. 12). На рисунке 12, в показаны все результаты измерений влажности грунтов. Видно, что значения w в основном сосредоточены в диапазоне 20–35%, что соответствует также и тому, что представлено на рисунках 12, а, б. А начальный коэффициент пористости/пустотности варьирует по вертикали от 0,4 до 1,2 (рис. 13) при среднем значении около 0,8. На рисунке14 показано распределение значений w и e для трех единиц пылеватой глины на основе 3D модели свойств. При этом имеется следующий тренд: пылеватая глина на большей глубине характеризуется более низким начальным коэффициентом пористости/пустотности. Тогда как ее влажность во всех трех единицах колеблется в сходных интервалах.

Рис. 12. Влажность грунтов w вдоль разреза WE18 в направлении З-В (а) и вдоль разреза NS12 в направлении С-Ю (б) на основе 3D модели; результаты всех измерений величин w в зависимости от глубины (в)

Рис. 13. Начальный коэффициент пористости/пустотности грунтов вдоль разрезра WE18, проходящего в направлении З-В (а), и разреза NS12, проходящего в направлении С-Ю (б); результаты всех измерений величин e в зависимости от глубины (в)

Рис. 14. Распределение влажности грунтов w, % (а) и начального коэффициента пористости/пустотности e (б) для трех слоев пылеватой глины в пределах района Тунчжоу

Модуль всестороннего сжатия E_s

На рисунке 15 показано распределение величин модуля всестороннего сжатия E_s (модуля объемной упругости). Видно, что значения E_s постепенно увеличиваются с глубиной. Для определения особенностей распределения величин E_s в зависимости от глубины далее будут отдельно рассмотрены различные ее интервалы (0–10, 10–17, 17–35 м и глубже 35 м) на основе локальных особенностей циклических последовательностей отложений, описанных выше.

Рис. 15. Распределение величин модуля всестороннего сжатия Es (модуля объемной упругости) в отложениях района Тунчжоу до глубины 50 м: а – трехмерная модель; б – схема разрезов

Для трех единиц пылеватой глины значения E_s варьируют в пределах:

• 2–27 МПа при среднем значении 6,4 МПа на глубине 0–10 м;
• 4–33 МПа при среднем значении 9,6 МПа на глубине 10–17 м;
• 5–47 МПа при среднем значении 11,8 МПа на глубине 17–35 м;
• 6–54 МПа при среднем значении 16,8 МПа на глубине 35–50 м.

Как показано на рисунке 16, слой пылеватой глины на небольшой глубине (0–10 м) с величиной E_s в основном меньше 7,5 МПа был определен как единица со сжимаемостью от средней до высокой (E_s=7,5 в соответствии с документом DBJ 11-501-2009, разработанным Пекинским геотехническим институтом и Пекинским институтом архитектурного проектирования в 2017 г.).

Слои пылеватой глины в интервалах глубины 10–17 и 17–35 м со значениями E_s, обычно сосредоточенными в диапазонах 4–7,5 и 7,5–11 МПа соответственно (см. рис. 16), были определены как единицы со сжимаемостью от средней до высокой и от низкой до средненизкой соответственно, а слой на глубине 35–50 м – как низкосжимаемая единица. Но все же, как видно из рисунка 17, на котором отражены величины E_s для пылеватой глины на разной глубине, значения E_s пылеватой глины в пределах единиц распределены неоднородно.

Рис. 16. Диаграммы значений Es для пылеватой глины в пределах разных интервалов глубины (1–10, 10–17, 17–35 и более 35 м) (Авторы не обозначили вертикальные оси на столбчатых диаграммах – должно быть, это количество измерений. – Ред.)

Рис. 17. Распределение величин Es для пылеватой глины на разной глубине (около, 10, 23, 35, 45 м) на основе 3D модели свойств

ПРИМЕНЕНИЕ

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

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

Изменчивость по глубине

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

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

Как указывалось выше, в районе Тунчжоу было выявлено десять основных/главных единиц. Из построенной модели можно экспортировать карты их мощности – например, карту мощности первого слоя песка (единицы 24), который непрерывен, имеет среднюю мощность около 9 м и широко распространен в районе Тунчжоу (рис. 18).

Рис. 18. Карта мощности первой единицы песка (единицы 24), полученная на основе 3D геологической модели для исследуемой территории

Кроме того, в качестве одного из геологических показателей для оценки возможности использования городского подземного пространства (его пригодности для строительства) выступает его литолого-стратиграфическое строение. В этом отношении авторы настоящей статьи путем обобщения выделили три типа структур (однослойную, двуслойную и многослойную, см. таблицу 2) при работе над проектом «Стандарт для геологической оценки ресурсов городского подземного пространства в центральной части г. Пекина» (в соответствии с документом DB11/T 1895-2021, разработанным Пекинским институтом геологических исследований в 2022 г.). Затем был рассчитан вес этого показателя с использованием нечеткого математического метода (основанного на нечеткой логике), как это делалось и для других базовых геологических показателей, таких как геотехнические свойства. Однако для других ограничений, например таких как активные разломы, оседание, проседание, разжижение песка, карстовые провалы и т.д., веса определяются с помощью метода анализа иерархий [58]. Такой подход согласуется с моделью для оценки, описанной в работе [59].

Таблица 2. Типы литолого-стратиграфического строения как один из показателей для оценки городского подземного пространства (адаптировано по [59, 60])

По 3D модели подземного пространства исследуемой территории до глубины 50 м было интерпретировано многослойное строение при локальной двуслойности или локальной неоднородности в отдельных интервалах глубины с относительно более низкой точностью количественной оценки, чем для однослойной зоны (см. таблицу 2).

Распределение значений E_s

Еще для применения построенной 3D модели важен ее потенциал для характеристики большого количества свойств в пределах выделенных единиц. В оценке городского подземного пространства большую роль играет анализ инженерно-геологических свойств, неоднородно распределенных внутри геологических единиц [6, 9, 61]. В качестве одного из показателей для этой оценки, как уже говорилось, приняты геотехнические свойства. Модуль всестороннего сжатия E_s (модуль объемной упругости) с учетом его влияния на инженерно-геологические условия включен в «Стандарт геологической оценки ресурсов городского подземного пространства в центральной части г. Пекина» (документ DB11/T 1895-2021, разработанный Пекинским институтом геологических исследований в 2022 году). Результаты количественной оценки этого модуля, выполненной авторами настоящей статьи, приведены в таблице 3. Они разделены на три интервала, содержащие от 6 до 10, от 3 до 6 и от 0 до 3 баллов. Эти интервалы соответствуют простым, средним и сложным геотехническим условиям соответственно. Затем на основе 3D модели свойств по значениям E_s рассчитываются три пространственные зоны. Две единицы пылеватой глины с величинами E_s более 15 интерпретируются как попадающие в итервал 6–10 баллов (рис. 19).

Таблица 3. Значения E_s и литолого-тратиграфическое строение, рассматриваемые в качестве критериев оценки пригодности городского подземного пространства для строительства

Рис. 19. Распределение значений Es, превышающих 15, в верхней и нижней единицах пылеватой глины (единицах 23 и 43 соответственно) на основе 3D модели свойств

Динамическая модель мониторинга

Динамическую модель мониторинга получают путем включения в 3D геологическую модель результатов измерений (например, с помощью распределенных волоконно-оптических систем мониторинга), выполняемых в течение некоего временного интервала. На рисунке 20 показаны деформации сжатия геологических единиц (песка и пылеватой глины) в интервале глубины 30–40 м в течение периода с 18.07.2019 по 22.09.2020. Эти слои демонстрируют низкую сжимаемость при средней вертикальной деформации 0,97 мм с июля 2019 года по сентябрь 2020 года (по данным ежегодного отчета по мониторингу геологической безопасности ресурсов подземного пространства на посту мониторинга в районе Тунчжоу г. Пекина за 2020 год), что значительно ниже, чем у верхнего слоя пылеватой глины.

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

На основе приведенных выше данных по E_s и по вариациям мощности можно сделать вывод, что для подземного строительства на исследуемой территории подходящими с точки зрения вертикальной устойчивости представляются единицы песка и пылеватой глины на глубине от 30 до 40 м. Этот вопрос будет дополнительно изучен при продолжающихся исследованиях с использованием 3D геологических моделей в качестве основы для построения модели динамического мониторинга для района Тунчжоу.

ЗАКЛЮЧЕНИЕ

Основная цель данного исследования состояла в том, чтобы представить эффективную и практичную процедуру 3D геологического моделирования четвертичных осадочных отложений на урбанизированной территории и обсудить ее потенциальное применение. В качестве типичного примера был выбран район Тунчжоу г. Пекина. Чтобы по-новому взглянуть на геологические условия до небольшой глубины, были использованы данные по 476 инженерно-геологическим скважинам (436 старым и 40 новым более глубоким). Затем эти сведения были использованы для 3D геологического моделирования грунтовых условий и оценки пригодности городского подземного пространства для строительства.

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

1. Была использована высококачественная информация по старым и новым инженерно-геологическим скважинам для определения литолого-стратиграфического строения грунтовой среды и построения вертикальных разрезов для территории района Тунчжоу г. Пекина;

2. С помощью полученной 3D геологической модели было представлено множество деталей грунтовых условий до глубины 50 м. В модель вошли 10 главных единиц, сгруппированных в три циклотемы (три циклические последовательности переслаивающихся глин, пылеватых грунтов, песков и гравия с переменным количеством линз).

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

4. Был выполнен анализ динамической модели мониторинга на основе полученной 3D модели. Он показал, что геологические единицы (песка и пылеватой глины) на глубине от 30 до 40 м со средней вертикальной деформацией 0,97 мм с июля 2019 года по сентябрь 2020 года подходят для подземного строительства на рассматриваемой территории с точки зрения вертикальной устойчивости. Этот вопрос будет дополнительно изучен при продолжающихся исследованиях с использованием построенной 3D геологической модели в качестве основы для динамической модели мониторинга подземного пространства территории района Тунчжоу.

-

Данная работа была выполнена при финансовой поддержке проекта «Трехмерное моделирование городских подземных ресурсов и окружающей среды в районе Тунчжоу г. Пекина» (0747-1761SITCN070).

Авторы хотели бы поблагодарить своих коллег из Пекинского института геологических исследований за помощь и поддержку, редактора и рецензентов за конструктивные предложения и замечания, а также Чао Гао (Chao Gao) за ценное обсуждение.

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

ИСТОЧНИК ДЛЯ ПЕРЕВОДА

He H., Xiao J., He J., Wei B., Ma X., Huang F., Cai X., Zhou Y., Bi J., Zhao Y., Wang Ch., Wei J. Three-dimensional geological modeling of the shallow subsurface and its application: a case study in Tongzhou District, Beijing, China // Applied Sciences. MDPI, 2023. Vol. 13. Article 1932. DOI: https://doi.org/10.3390/app13031932, URL: https://www.mdpi.com/2076-3417/13/3/1932

СПИСОК ЛИТЕРАТУРЫ, ИСПОЛЬЗОВАННОЙ АВТОРАМИ ПЕРЕВЕДЕННОЙ СТАТЬИ

1. Wang, C.S.; Zhou, C.H.; Peng, J.B.; Fan, J.; Zhu, H.H.; Li, X.Z.; Cheng, G.H.; Dai, C.S.; Xu, N.X. A discussion on high-quality development and sustainable utilization of China’s urban underground space in the new era. Earth Sci. Front. 2019, 26, 1–8, (In Chinese with English abstract).
2. Aldiss, D.T.; Black, M.G.; Entwisle, D.C.; Page, D.P.; Terrington, R.L. Benefits of a 3D geological model for major tunnelling works: An example from Farringdon, east-central London, UK. Q. J. Eng. Geol. Hydrogeol. 2012, 45, 405–414. [CrossRef].
3. Andersen, T.R.; Poulsen, S.E.; Pagola, M.A.; Medhus, A.B. Geophysical mapping and 3D geological modelling to support urban planning: A case study from Vejle, Denmark. J. Appl. Geophys. 2020, 180, 104130. [CrossRef].
4. Costa, A.L.; Sousa, R.L.; Einstein, H.H. Probabilistic 3D alignment optimization of underground transport infrastructure integrating GIS-based subsurface characterization. Tunn. Undergr. Space Technol. 2017, 72, 233–241. [CrossRef]. 
5. Culshaw, M.G. From concept towards reality: Developing the attributed 3D geological model of the shallow subsurface. Q. J. Eng. Geol. Hydrogeol. 2005, 38, 231–284. [CrossRef]. 
6. Zhou, D.K.; Li, X.Z.; Wang, Q.; Wang, R.; Wang, T.D.; Gu, Q.; Xin, Y.X. GIS-based urban underground space resources evaluation toward three-dimensional land planning: A case study in Nantong, China. Tunn. Undergr. Space Technol. 2019, 84, 1–10. [CrossRef].
7. Dong, M.; Neukum, C.; Hu, H.; Azzam, R. Real 3D geotechnical modeling in engineering geology: A case study from the inner city of Aachen, Germany. Bull. Eng. Geol. Environ. 2015, 74, 281–300. [CrossRef].
8. Maio, R.D.; Paola, C.D.; Forte, G.; Piegari, E.; Urciuoli, G. An integrated geological, geotechnical and geophysical approach to identify predisposing factors for flowslide occurrence. Eng. Geol. 2020, 267, 105473. [CrossRef].
9. Hou,W.; Yang, L.; Deng, D.; Ye, J.; Clarke, K.; Yang, Z.; Zhuang,W.; Liu, J.; Huang, J. Assessing quality of urban underground spaces by coupling 3D geological models: The case study of Foshan city, South China. Comput. Geosci. 2016, 89, 1–11. [CrossRef].
10. Rienzo, F.D.; Oreste, P.; Pelizza, S. Subsurface geological-geotechnical modelling to sustain underground civil planning. Eng. Geol. 2008, 96, 187–204. [CrossRef].
11. Rienzo, F.D.; Oreste, P.; Pelizza, S. 3D GIS Supporting Underground Urbanisation in the City of Turin (Italy). Geotech. Geol. Eng. 2009, 27, 539–547. [CrossRef].
12. Royse, K.R.; Rutter, H.K.; Entwisle, D.C. Property attribution of 3D geological models in the Thames Gateway, London: New ways of visualising geoscientific information. Bull. Eng. Geol. Environ. 2008, 68, 1–16. [CrossRef].
13. Breunig, M. An approach to the integration of spatial data and systems for a 3D geo-information system. Comput. Geosci. 1999, 25, 39–48. [CrossRef].
14. Lemon, A.M.; Jones, N.L. Building solid models from boreholes and user-defined cross-sections. Comput. Geosci. 2003, 29, 547–555. [CrossRef].
15. Wu, Q.; Xu, H.; Zou, X.K. An effective method for 3D geological modeling with multi-source data integration. Comput. Geosci. 2005, 31, 35–43. [CrossRef].
16. Maccormack, K.; Arnaud, E.; Parker, B.L. Using a Multiple Variogram Approach to Improve the Accuracy of Subsurface Geological Models. Can. J. Earth Sci. 2017, 55, 786–801. [CrossRef].
17. Touch, S.; Likitlersuang, S.; Pipatpongsa, T. 3D geological modelling and geotechnical characteristics of Phnom Penh subsoils in Cambodia. Eng. Geol. 2014, 178, 58–69. [CrossRef].
18. Turner, A.K. Challenges and trends for geological modelling and visualisation. Bull. Eng. Geol. Environ. 2006, 65, 109–127. [CrossRef].
19. Zhang, F.; Zhu, H.H.; Ning, M.X. Modeling method of 3D strata suitable for massive data. Chin. J. Rock Mech. Eng. 2006, 25, 3305–3310, (In Chinese with English abstract).
20. Zhu, L.F.; Zhang, C.J.; Li, M.J.; Pan, X.; Sun, J.Z. Building 3D solid models of sedimentary stratigraphic systems from borehole data: An automatic method and case studies. Eng. Geol. 2012, 127, 1–13. [CrossRef].
21. Caumon, G.; Gray; Antoine, C.; Titeux, M.-O. Three-Dimensional Implicit Stratigraphic Model Building From Remote Sensing Data on Tetrahedral Meshes: Theory and Application to a Regional Model of La Popa Basin, NE Mexico. IEEE Trans. Geosci. Remote Sens. 2013, 51, 1613–1621. [CrossRef].
22. Caumon, G.; Collon-Drouaillet, P.; Veslud, C.L.C.D.; Viseur, S.; Sausse, J. Surface-Based 3D Modeling of Geological Structures. Math. Geosci. 2009, 41, 927–945. [CrossRef].
23. Frank, T.; Tertois, A.L.; Mallet, J.L. 3D-reconstruction of complex geological interfaces from irregularly distributed and noisy point data. Comput. Geosci. 2007, 33, 932–943. [CrossRef].
24. Jшrgensen, F.; Hшyer, A.-S.; Sandersen, P.B.E.; He, X.L.; Foged, N. Combining 3D geological modelling techniques to address variations in geology, data type and density – an example from Southern Denmark. Comput. Geosci. 2015, 81, 53–63. [CrossRef].
25. Pellerin, J.; Lйvy, B.; Caumon, G.; Botella, A. Automatic surface remeshing of 3D structural models at specified resolution: A method based on Voronoi diagrams. Comput. Geosci. 2014, 62, 103–116. [CrossRef].
26. Wang, Z.G.; Qu, H.G.;Wu, Z.X.; Yang, H.J.; Du, Q.L. Formal representation of 3D structural geological models. Comput. Geosci. 2016, 90, 10–23. [CrossRef].
27. Gallerini, G.; Donatis, M. 3D modeling using geognostic data: The case of the low valley of Foglia river (Italy). Comput. Geosci. 2009, 35, 146–164. [CrossRef].
28. Travelletti, J.; Malet, J.P. Characterization of the 3D geometry of flow-like landslides: A methodology based on the integration of heterogeneous multi-source data. Eng. Geol. 2012, 128, 30–48. [CrossRef].
29. Zhu, L.F.; Li, M.J.; Li, C.L.; Shang, J.G.; Chen, G.L.; Zhang, B.;Wang, X.F. Coupled modeling between geological structure fields and property parameter fields in 3D engineering geological space. Eng. Geol. 2013, 167, 105–116. [CrossRef].
30. Collon, P.; Steckiewicz-Laurent, W.; Pellerin, J.; Laurent, G.; Caumon, G.; Reichart, G.; Vaute, L. 3D geomodelling combining implicit surfaces and Voronoi-based remeshing: A case study in the Lorraine Coal Basin (France). Comput. Geosci. 2015, 77, 29–43. [CrossRef].
31. Price, S.J.; Terrington, R.L.; Busby, J.; Bricker, S.; Berry, T. 3D ground-use optimisation for sustainable urban development planning: A case-study from Earls Court, London, UK. Tunn. Undergr. Space Technol. 2018, 81, 144–164. [CrossRef].
32. Rodrнguez, A.A.; Kuhlmann, U.; Marschall, P. 3D modelling of the Excavation Damaged Zone using a Marked Point Process technique. Geomech. Energy Environ. 2018, 17, 29–46. [CrossRef].
33. Ye, S.; Luo, Y.;Wu, J.; Yan, X.;Wang, H.; Jiao, X.; Teatini, P. Three-dimensional numerical modeling of land subsidence in Shanghai, China. Hydrogeol. J. 2016, 24, 695–709. [CrossRef].
34. May, M.E.; Dlala, M.; Chenini, I. Urban geological mapping: Geotechnical data analysis for rational development planning. Eng. Geol. 2010, 116, 129–138. [CrossRef].
35. Rohmer, O.; Bertrand, E.; Mercerat, E.D.; Rйgnier, J.; Alvarez, M. Combining borehole log-stratigraphies and ambient vibration data to build a 3D Model of the Lower Var Valley, Nice (France). Eng. Geol. 2020, 270, 105588. [CrossRef].
36. Burke, H.; Mathers, S.J.;Williamson, J.P.; Thorpe, S.; Ford, J.; Terrington, R.L. The London Basin Superficial and Bedrock LithoFrame 50 Model; British Geological Survey: Nottingham, UK, 2014.
37. Ford, J.R.; Burke, H.; Royse, K.R.; Mathers, S.J. The 3D Geology of London and the Thames Gateway: A Modern Approach to Geological Surveying and Its Relevance in the Urban Environment; British Geological Survey: Nottingham, UK, 2008.
38. MacCormack, K.E.; Berg, R.C.; Kessler, H.; Russell, H.A.J.; Thorleifson, L.H. 2019 Synopsis of Current Three-Dimensional Geological Mapping and Modeling in Geological Survey Organizations; AER/AGS Special Report 112; Alberta Energy Regulator/Alberta Geological Survey: Edmonton, AB, Canada, 2019.
39. Hшyer, A.S.; Klint, K.E.S.; Fiandaca, G.; Maurya, P.K.; Christiansen, A.V.; Balbarini, N.; Bjerg, P.L.; Hansen, T.B.; Mller, I. Development of a high-resolution 3D geological model for landfill leachate risk assessment. Eng. Geol. 2019, 249, 45–59. [CrossRef].
40. Chen, Q.; Liu, G.; Ma, X.; Li, X.; He, Z. 3D stochastic modeling framework for Quaternary sediments using multiple-point statistics: A case study in Minjiang Estuary area, Southeast China. Comput. Geosci. 2019, 136, 104404. [CrossRef].
41. Cai, X.M.; Luan, Y.B.; Guo, G.X.; Liang, Y.N. 3D Quaternary geological structure of Beijing plain. Geol. China 2009, 36, 1021–1029, (In Chinese with English abstract).
42. Zhi, C.; Chen, H.; Li, P.; Ma, C.; Zhang, J.; Zhang, C.;Wang, C.; Yue, X. Spatial distribution of arsenic along groundwater flow path in Chaobai River alluvial–proluvial fan, North China Plain. Environ. Earth Sci. 2019, 78, 259. [CrossRef].
43. Che, Z. A study of the fault activity in the capital circle. North China Earthq. Sci. 1993, 11, 23–34.
44. Jiang, W.L.; Hou, Z.H.; Xie, X. Research on paleoearthquakes in Jiuxian trenches across Nankou-Sunhe fault zone in Changping County of Beijing plain. Sci. China 2002, 45, 160–173. [CrossRef].
45. Zhang, S.M.;Wang, D.D.; Liu, X.D.; Zhang, G.H.; Zhao, J.X.; Luo, M.H.; Ren, J.J.;Wang, R.; Zhang, Y. Using borehole core analysis to reveal Late Quaternary paleoearthquakes along the Nankou-Sunhe Fault, Beijing. Sci. China Ser D Earth Sci. 2008, 51, 1154–1168. [CrossRef].
46. Huang, M.; Ninic, J.; Zhang, Q. BIM, machine learning and computer vision techniques in underground construction: Current status and future perspectives. Tunn. Undergr. Space Technol. 2021, 108, 103677. [CrossRef].
47. He, H.H.; Jing, H.; Xiao, J.Z.; Zhou, Y.X.; Liu, Y.; Li, C. 3D geological modeling and engineering properties of shallow superficial deposits: A case study in Beijing, China. Tunn. Undergr. Space Technol. 2020, 100, 103390. [CrossRef].
48. Michel, D.; Roger, A. Geostatistical Ore Reserve Estimation; Elsevier Scientific Pub. Co.: Amsterdam, The Netherland, 1977.
49. Isaaks, E.H.; Srivastava, R.H. An Introduction to Applied Geostatistics; Oxford University Press: New York, NY, USA, 1989.
50. Krige, D.G. A statistical approach to some basic mine valuation problems on theWitwatersrand. J. South Afr. Inst. Min. Metall. 1952, 52, 119–139.
51. Triantafilis, J.; Odeh, I.; Mcbratney, A.B. Five Geostatistical Models to Predict Soil Salinity from Electromagnetic Induction Data Across Irrigated Cotton. Soil Sci. Soc. Am. J. 2001, 65, 869–878. [CrossRef].
52. Lloyd, C.D.; Atkinson, P.M. Assessing uncertainty in estimates with ordinary and indicator kriging. Comput. Geosci. 2001, 27, 929–937. [CrossRef].
53. Kravchenko, A.N. Influence of Spatial Structure on Accuracy of Interpolation Methods. Soil Sci. Soc. Am. J. 2003, 67, 1564–1571. [CrossRef].
54. Zimmerman, D.L.; Zimmerman, M.B. A Comparison of Spatial Semivariogram Estimators and Corresponding Ordinary Kriging Predictors. Technometrics 1991, 33, 77–91. [CrossRef].
55. Liu, D.; Wang, Z.; Zhang, B.; Song, K.; Duan, H. Spatial distribution of soil organic carbon and analysis of related factors in croplands of the black soil region, Northeast China. Agric. Ecosyst. Environ. 2006, 113, 73–81. [CrossRef].
56. Webster, R.; Oliver, M.A. Geostatistics for Environmental Scientists, Second Edition; Wiley: Chichester, UK, 2008.
57. Schokker, J.; Sandersen, P.; Beer, H.D.; Eriksson, I.; Kallio, H.; Kearsey, T.; Pfleiderer, S.; Seither, A. 3D Urban Subsurface Modelling and Visualisation—A Review of Good Practices and Techniques to Ensure Optimal Use of Geological Information in Urban Planning. COST European Cooperation in Science and Technology Sub-Urban WG2 Report. 2017. Available online: https://www.google.com.hk/url?sa=t&rct=j&q=&esrc=s&source=web&cd=&cad=rja&uact=8&ved=2ahUKEwi9 suvtxPH8AhVfr1YBHSDnDOIQFnoECAkQAQ&url=https%3A%2F%2Fwww.researchgate.net%2Fpublication%2F317304360 _3D_urban_subsurface_modelling_and_visualisation_-_a_review_of_good_practices_and_techniques_to_ensure_optimal_ use_of_geological_information_in_urban_planning&usg=AOvVaw2M8yqT1ejj3boaKZ6frxeT (accessed on 26 January 2023).
58. Saaty, T.L. Decision making with the analytic hierarchy process. Int. J. Serv. Sci. 2008, 1, 83–98. [CrossRef].
59. He, J.; Zhou, Y.X.; Zheng, G.S.; Wang, J.M.; Liu, Y. Research on the Geological Suitability Evaluation System of Underground Space Resource Utilization in Beijing. Chin. J. Undergr. Space Eng. 2020, 16, 955–966, (In Chinese with English abstract).
60. He, J.; He, H.H.; Zheng, G.S.; Liu, Y.; Zhou, Y.X.; Xiao, J.Z.; Wang, C.J. 3D geological modelling of superficial deposits in Beijing City. Geol. China 2019, 46, 244–254, (In Chinese with English abstract).
61. Thierry, P.; Prunier-Leparmentier, A.-M.; Lembezat, C.; Vanoudheusden, E.; Vernoux, J.-F. 3D geological modelling at urban scale and mapping of ground movement susceptibility from gypsum dissolution: The Paris example (France). Eng. Geol. 2009, 105, 51–64. [CrossRef].