Введение

В основе георадиолокации лежит регистрация отражений электромагнитных волн от границ, где электрические параметры среды (удельная проводимость или диэлектрическая проницаемость) претерпевают изменение [1, с. 24–28, 47–67]. Сверхширокополосные импульсы формируют электромагнитные волны, которые, проходя через толщу грунта, встречают на пути различные включения и слоистые неоднородности; отраженные от них сигналы улавливаются приемной антенной, после чего усиливаются, оцифровываются и поступают на компьютерную обработку. Каждое отражение образует на мониторе трассу, а их совокупность – георадиолокационный разрез (радарограмму). Таким образом, на экране формируется изображение, по которому судят о расположении и характере выявленных объектов.

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

В настоящее время для электродинамического моделирования часто применяют метод конечных разностей во временной области (FDTD) [2; 3]. Есть и другие методы моделирования [4–6]. Отдельно можно выделить метод матриц линий передач (TLM) [7] и метод импедансного аналога пространства [8]. Однако упомянутые методы не учитывают или делают затруднительным учет дисперсии и поляризационных потерь в неоднородных средах, к тому же дополнительно требуют значительных вычислительных ресурсов, а использование рекурсивных методов расчета значительно увеличивает время вычислений. Указанных выше недостатков можно избежать, если применить дифференциальные ABCD-матрицы. Данный метод базируется на разбиении среды на элементарные фрагменты, каждый из которых описывается передаточной матрицей, сформированной на основе решений уравнений Максвелла. Сама вычислительная процедура реализуется через метод прямого счета [9, с. 162].

В 1986 г. сотрудниками лаборатории георадиолокации ИГДС СО РАН были проведены исследования стоянки древнего человека «Диринг-Юрях» в Центральной Якутии (возраст 0,2–3 млн лет), где на данных георадиолокации наблюдался эффект отклонения отраженного сигнала от нулевой оси. Тогда же В. Н. Федоровым было выдвинуто предположение, что это эффект накопления заряда на границе раздела талых и мерзлых пород.

Как было показано Н. О. Кожевниковым в публикации [10], аналогичные явления регистрировались при МПП-исследованиях (методом переходных процессов) в зонах многолетней мерзлоты (Аляска, Канада, Якутия). При этом наблюдались сигналы с немонотонным спадом и сменой знака ЭДС. Такое поведение переходных процессов было впоследствии обозначено [11] как эффект быстрой индукционно-вызванной поляризации (ВПИ) вмещающей среды.

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

Материалы и методы исследования

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

Математическая модель дисперсионной среды. Основным параметром, характеризующим физическую среду для электромагнитных волн (ЭМВ) при подповерхностной георадиолокации, является комплексная диэлектрическая проницаемость, когда относительная комплексная магнитная проницаемость среды μr = 1 [1]:

, (1)

где – относительная комплексная диэлектрическая проницаемость, ε′– действительная (статическая) часть, связанная с накоплением энергии в диэлектрике приложенного поля, а ε″– мнимая часть, связанная с потерями энергии в среде. Эти потери на частотах ниже 1 ТГц вызваны потерями на релаксацию εd и проводимость σcm [12, с. 108–118]:

, (2)

где ω = 2πf – круговая частота, ε0 – абсолютная диэлектрическая проницаемость.

Аномальная дисперсия почв, грунтов и горных пород связана с существованием нескольких релаксационных процессов [10] с различными временами релаксации. Чаще всего это связано c содержанием воды, которая и определяет частотную зависимость диэлектрических проницаемостей, диссипативных и поляризационных потерь [13].

Диэлектрическую релаксацию в мерзлых породах можно описать моделью Дебая [12, с. 114] с временем релаксации τрел ≅ 100 мкс [10]:

, (3)

где εcm – диэлектрическая проницаемость среды на низких частотах; ε∞ – значения диэлектрической проницаемости среды на высоких частотах.

Можно предположить, что быстро устанавливающаяся поляризация в мерзлых породах представляет собой диэлектрическую релаксацию льда, содержащегося в порах мерзлой породы [12, с. 142–157].

Эти явления можно объяснить эффектом Максвелла – Вагнера, когда в гетерогенных средах под действием электрического поля, направленного перпендикулярно границам между областями, отличающимися отношением диэлектрической проницаемости к удельной электропроводности, накапливаются заряды [10; 14]. Модель гетерогенной геологической среды, образованная непрерывной вмещающей фазой и овальными включениями приведена на рис. 1.

Рис. 1. Модель гетерогенной геологической среды, образованная непрерывной вмещающей фазой и овальными включениями Примечание: составлен авторами на основе источников [10; 14]

Чтобы учесть в модели тепловые потери, потери на отражение и излучение, в выражении (2) введена проводимость σcm. При этом следует отметить, что при релаксационных потерях особое значение будет иметь проводимость среды на высоких и ультравысоких частотах – σ∞ [15]. В таком случае представить время релаксации τрел в соотношении Дебая можно как

, (4)

где εст и ε∞ – значения диэлектрической проницаемости среды на низких и высоких частотах соответственно, σ∞ – проводимость среды на высоких частотах. Тогда получается, что если учесть релаксационные потери, то полную проводимость горной породы Ẏ можно записать как

. (5)

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

В этой модели принято допущение о том, что элементарный объем однородной среды ∆x∆y∆z стремится к бесконечно малым линейным размерам. Электродинамическая задача решается с помощью метода ABCD-матриц, где для каждой ячейки строится своя матрица передачи по стандартной схеме [16]. Чтобы реализовать этот подход, используют комплексные величины проводимости Ẏ и сопротивления Ż среды, определяемые формулами

. (6)

Тогда уравнения Максвелла для гармонических ЭМВ запишутся как [9]:

(7)

Постоянная распространения и волновое сопротивление , а затухание α и фазовая скорость V запишутся как

дБ/м

м/с (8)

Рассмотрим плоскую ЭМВ или Т-волну в однородной среде. Пусть Т-волна распространяется в направлении , а вектор электрического поля Ė параллелен орте . Учитывая, что для Т-волн

получим из (7) систему из двух однородных дифференциальных уравнений (ОДУ):

(9)

Система ОДУ (9) описывает распространение Т-волны с линейной поляризацией в одномерной однородной среде.

Переходя к конечным разностям dz→Δz, систему ОДУ (9), описывающую распространение ЭМВ в направлении увеличения координаты z, можно переписать следующим образом:

(10)

Соотношения (10) пересчитывают поля на элементарном участке слева направо, в направлении увеличения z (от z к z+Δz). Это дифференциальная матрица передачи ЭМВ участка среды мощностью Δz. Однако при построении итоговой матрицы передачи удобнее вести обратный пересчет – от правой границы к левой. Поэтому далее матрица участка вводится для пересчета полей справа налево (от z+Δz к z). Тогда, в соответствии с теорией многополюсников [17], обратная дифференциальная ABCD-матрица участка среды мощностью Δz запишется как

. (11)

При таком определении матрица (11) является обратной и позволяет последовательно пересчитывать поля от правой границы слоя (z = L) к левой (z = 0). Величина Δz должна быть много меньше длины волны исследуемого диапазона частот. Устойчивость и условие сходимости дифференциальной ABCD-матрицы (11) рассмотрены в работе [9].

Неоднородную одномерную среду мощностью h можно представить в виде каскадного соединения N ABCD-матриц слоев c различными Ẏj и Żj, j = 1,2,…,N. Для определения матрицы передачи [a]Σ такой среды достаточно перемножить ABCD-матрицы:

, (12)

где .

Тогда результирующую матрицу передачи многослойной неоднородной среды, связывающую поля на правой границе среды (z = L) с полями на левой границе (z = 0), можно записать как

, (13)

где Ėx(0)Ėx(L) – амплитуда вектора электрического поля слева и справа от среды соответственно; Ḣy(0)Ḣy(L) – амплитуда вектора магнитного поля слева и справа от среды соответственно.

а) б)

Рис. 2. Зондирующий импульс Eзонд георадара серии «ОКО»: а) временная характеристика; б) амплитудно-частотная характеристика Примечание: составлен авторами по результатам данного исследования

Была разработана программа, моделирующая процесс распространения видеоимпульса георадара серии «ОКО» с центральной частотой 250 МГц в слоистой дисперсионной среде (рис. 2).

Условия распространения Т-волн определяются граничными условиями слева и справа от среды:

, (14)

где Eзонд – амплитуда зондирующего импульса; Z1,2 – волновые сопротивления сред слева и справа от исследуемой среды соответственно.

Решая совместно (13) и (14), находим амплитуды магнитного Hy(L) и электрического Ex(L) полей справа:

(15)

где – элементы результирующей матрицы [a]Σ. Остальные амплитуды электрического и магнитного полей находим из (13) и (14).

Результаты исследования и их обсуждение

Адекватность предложенной модели проверялась путем расчетного моделирования с последующим сравнением с эталонными лабораторными измерениями, выполненными на различных образцах четвертичных отложений Центральной Якутии методом георадиолокации [18]. В этих опытах изучалась зависимость скорости распространения ЭМВ от состояния мерзлости и влажности. Результаты показали, что все изученные породы демонстрируют значение коэффициента Nv, которое с относительным отклонением 4,3 % хорошо описывается линейным уравнением

, (16)

где VM – скорость ЭМВ в мерзлом образце (0,150 м/нс); VT – скорость ЭМВ в различных состояниях мерзлости, вплоть до полной оттайки образца; W – весовая влажность песка (3–22 %); c – поправочный коэффициент, зависящий от типа образца: для горного песка 0,3, для кварцевого – 11. Таким образом, скорость ЭМВ в талом горном песке снижается от 5 до 34 % в зависимости от влажности и от 15 до 45 % в речном песке.

Был рассчитан случай – речной песок, влажность 22 %. На рис. 3 приведены частотные зависимости затухания α и скорости распространения V ЭМВ в мерзлом и талом песках.

Расчеты показали, что с увеличением частоты потери α и фазовая скорость V ЭМВ в талых песках растет, при этом на низких частотах эти потери уменьшатся, как это показано на рис. 3, а. На рис. 3, б, приведена частотная зависимость скорости ЭМВ в талом и мерзлом песках. Наблюдается хорошее соответствие расчета с экспериментальными данными, что подтверждает эквивалентность математической модели и экспериментальных результатов.

а) б)

Рис. 3. Частотная зависимость затухания α (а) и скорости распространения V (б) ЭМВ в мерзлом и талом речном песке с влажностью 22 % Примечание: составлен авторами по результатам данного исследования

а) б)

Рис. 4. Отраженный сигнал U(x) георадара: а) в начале при x = 0 и в конце талого слоя при x = L; б) амплитудно-частотная характеристика сигнала в начале и конце талого слоя Примечание: составлен авторами по результатам данного исследования

Рис. 5. Георадарные данные сезонных наблюдений: в летнее время (а) и зимой (б) Примечание: составлен авторами по результатам данного исследования

Рассчитывалось распространение видеоимпульса в образце мерзлого песка, состоящего из сезонно-талого слоя мощностью 2 м и мерзлого слоя. На рис. 4, а, приведены отраженные сигналы U(x) георадара на границе воздух – талый грунт при x = 0 и на границе талый грунт – мерзлый грунт при x = L. На георадарной трассе наблюдается отклонение сигнала от нулевой оси в положительную сторону на границе талого и мерзлого слоев (в диапазоне 15–35 нс). На рис. 4, б, приведены амплитудно-частотные характеристики этих сигналов. На них можно заметить сильное поглощение высокочастотных составляющих зондирующего сигнала (синяя пунктирная линия свыше 100 МГц). Отклонение в положительную или отрицательную сторону по шкале амплитуд на рис. 4, а, определяется поляризацией вектора Е зондирующего импульса. Это отклонение вызвано более сильной поляризацией талого слоя под действием внешнего электрического поля, чем мерзлого слоя, что зависит от величин σ∞ и ε∞. Можно предположить, что в талом слое происходит накопление электрического заряда на время τрел. Отметим, что на скорость распространения отраженного сигнала влияет проводимость слоев, при этом σ∞ влияет больше, чем σст.

Для сравнения на рис. 5 представлены экспериментальные радарограммы, полученные георадаром «ОКО-3» с антенной 250 МГц. Георадиолокация выполнена вдоль борта карьера на месторождении песка «Хатынг-Юряхское» в Центральной Якутии в летнее и зимнее время.

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

Заключение

Моделирование данных георадиолокации, полученных на образцах четвертичных пород Центральной Якутии с разной степенью мерзлости, подтвердило наличие поляризационного эффекта между слоями горных пород, отличающихся по электропроводности. Данный эффект, обусловленный механизмом Максвелла – Вагнера, проявляется на радарограммах как смещение сигнала относительно нулевой линии при переходе ЭМВ из проводящей в диэлектрическую среду. Одновременно с этим в частотном спектре отраженного импульса преобладает низкочастотная компонента, что связано с диэлектрическими потерями релаксационного типа, описываемыми моделью Дебая. Сопоставление с натурными георадиолокационными измерениями на участке с сезонно-талым слоем показывает, что учет этого поляризационного явления при обработке и интерпретации радарограмм дает возможность достоверно определять криогенное состояние массива и локализовать участки с талыми породами. Также появляется перспектива мониторинга динамики криогенных процессов, активизирующихся при оттаивании под влиянием природно-климатических факторов или антропогенных нагрузок.