Моделирование работы георадара

  • Tutorial
Георадар (радиотехнический прибор подповерхностного зондирования, GPR, Ground Penetrating Radar), применяющийся в настоящее время весьма широко — от картирования нор кроликов и изучения ящериц до поиска мин, остается достаточно дорогим удовольствием.

image

Дисплей георадара (кадр из британской телепередачи «Команда времени»)

Но оценить его возможности и изучить различные аспекты взаимодействия электромагнитного поля георадара со средой можно и без приобретения/аренды «железного» устройства. В этом нам поможет пакет gprMax (gpr — от абреввиатуры GPR, Max — начальные буквы фамилии Джеймса Клерка Максвелла, заложившего основы электродинамики), распространяющийся под лицензией GNU GPL v3.
Авторы этого проекта, начатого в 1996 году — Крэйг Уоррен (Craig Warren) из Нортумбрийского университета и Антонис Джианопулос (Antonis Giannopoulos) из Эдинбургского университета. Пакет был изначально разработан на C, а затем полностью переписан на комбинации Python 3/Cython.

image

Установка пакета требует наличия установленных компилятора, поддерживающего OpenMP (Microsoft Visual C++ 2015 Build Tools (рекомендуется именно эта версия!) для Windows / gcc для Linux), библиотеки NumPy и компилятора Cython. Загрузив из репозитария на GitHub и распаковав исходный код проекта, переходим в корневую папку и выполняем команды:

python setup.py build
python setup.py install

В качестве «быстрого старта» рассмотрим работу с пакетом на несложном двумерном примере — передающая антенна T импульсного радара (impulse GPR) излучает электромагнитный импульс, часть энергии которого непосредственно достигает приемной антенны R в виде прямой волны (DW — direct wave), а часть проникает сквозь песок, отражается от поверхности проводящего цилиндра и достигает приемной антенны в виде отраженной волны (RW — reflected wave):

image

Формат входного файла
Создадим папку models в корневой папке проекта, в которую поместим текстовый файл hello.in, содержащий команды для выполнения моделирования (приводимые далее команды соответствуют актуальной (третьей) версии проекта).

Каждая команда имеет вид:

#команда: параметр_1 параметр_2 параметр_3 ...

На одной строке можно записать только одну команду, причем первым символом строки, содержащей команду, должен быть #.

Команды можно сопровождать комментариями:

##комментарий

Порядок команд важен для команд конструирования объектов — такие команды выполняются в порядке их следования во входном файле.

Форма импульса

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

image

  1. один период синусоиды (sine)
  2. гауссов импульс (gaussian)
  3. «мексиканская шляпа» -«mexican hat» (ricker) — пропорционален второй производной гауссовой функции, форма кривой такого импульса напоминает сомбреро (эта форма импульса была использована американским геофизиком Норманом Рикером (Norman Ricker) в 1953 году при изучении сейсмических сигналов)

Выбираем для нашего примера гауссов импульс (тип импульса — gaussian) с центральной частотой $f_с = 1 \ \ ГГц = 1 \cdot 10^9 Гц$, задаваемый командой:

#waveform: gaussian 1 1e9 pulse

(1 — условная амплитуда импульса, pulse — метка импульса)

В этом случае импульс, используемый при моделировании, описывается выражением:

$W(t) = e^{-2 \cdot {\pi}^2 \cdot {f_c}^2 {(t-{1 \over {f_c}})}^2} $


Модель среды и система координат

При 2D-моделировании исследуемая область разбивается на ячейки заданного размера, а система координат модели выглядит так — оси X и Y образуют расчетную плоскость (шириной $w$ и высотой $h$), протяженность модели по оси Z имеет величину, равную шагу дискретизации $\Delta$.

image

При выборе шага дискретизации можно следовать эмпирическому правилу — размер шага не должен превышать одной десятой наименьшей длины волны, исследуемой в модели («10 ячеек на длину волны»):

$\Delta \le 0,1 \cdot {\lambda}_{min} $


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

Скорость распространения электромагнитной волны в среде с относительной диэлектрической проницаемостью ${\epsilon}_r$, выраженная в сантиметрах в наносекунду — принятая в радиолокации единица скорости, определяется выражением:

$v = {30 \over {\sqrt {{\epsilon}_r}}}$


Длина волны в сантиметрах при этом определяется выражением:

$\lambda = {v \over f}$


($f$ — частота в ГГц).
Для вывода формы гауссова импульса и его спектра можно использовать команду:

python -m tools.plot_source_wave gaussian 1 1e9 5e-9 1e-12 -fft

(gaussian — тип импульса, 1 — амплитуда импульса, 1e9 — центральная частота (1 ГГц), 5e-9 — длительность отображения импульса (5 нс), 1e-12 — шаг по времени (1 пс))

image

По спектру гауссова импульса с $f_с = 1 \ \ ГГц$ определяем, что на уровне -40 дБ частота $f \approx 3 \ \ ГГц$.
В рассматриваемом примере средой, в которой находится зондируемый объект, выберем сухой песок с относительной диэлектрической проницаемостью ${\epsilon}_r = 3$.
Найдем скорость распространения электромагнитной волны в песке:

$v = {30 \over {\sqrt {3}}} = 17,3 \ \ см/нс$


Определим длину волны в песке:

$\lambda = {v \over f} = {17,3 \over {3}} = 5,8 см = 58 мм $


Исходя из этого, выбираем шаг, одинаковый для всех осей ($\Delta = {\Delta}_x = {\Delta}_y = {\Delta}_z$) и равный 2 мм = 0,002 м из соображений удобства (в 1 см укладывается целое число шагов):

#dx_dy_dz: 0.002 0.002 0.002

Ограничим моделируемую область прямоугольником шириной $w$, равной 80 см = 0,8 м, и высотой $h$, равной 60 см = 0,6 м:

#domain: 0.60 0.60 0.002

(для двумерной модели требуется указать толщину, равную одному шагу (0.002))

Размер области моделирования и величина пространственного шага определяют число ячеек модели и, соответственно, требования к памяти компьютера.

Опишем песок с удельной проводимостью $\sigma = 0,01 \ \ См/м$ и относительной диэлектрической проницаемостью ${\epsilon}_r = 3$ командой:

#material: 3 0.01 1 0 sand

(1 соответствует относительной магнитной проницаемости ${\mu}_r$, равной единице (магнитные свойства отсутствуют), 0 — отсутствию магнитных потерь, а sand — метка этого материала).

Заполним песком большую часть моделируемой области (от y = 0 до y = 38 см = 0,38 м):

#box: 0 0 0 0.80 0.38 0.002 sand

(0 0 0 — координаты левого нижнего угла, 0.80 0.38 0.002 — координаты правого верхнего угла (0.002 — шаг дискретизации))
Оставшаяся часть является свободным пространством (метка free_space), по свойствам практически эквивалентным воздуху (${\epsilon}_r = 1$, ${\mu}_r = 1$, ${\sigma} = 0$).
Границы области моделирования представляются как идеально поглощающий электромагнитное излучение материал (Absorbing Boundary Condition, ABC).

В качестве мишени выберем цилиндр из идеального проводника (полностью отражающий электромагнитное излучение) радиусом 6 см = 0,06 м с центром, расположенным в точке с координатами x = 25 см = 0,25 м и y = 10 см = 0,1 м:

#cylinder: 0.250 0.1 0 0.250 0.1 0.002 0.06 pec

(pec — идеально проводящий материал)

Антенны

Моделируемый георадар оснащен двумя антеннами — передающей и приемной.
В нашем учебном примере представим передающую антенну диполем Герца, длина которого равна шагу дискретизации (в трехмерном случае возможен выбор антенны из обширной библиотеки), лежащим на песке (работа в контакте с зондируемой средой) на расстоянии 5 см слева от середины области (x = 35 см = 0,35 м, y = 38 см = 0,38 м):

#hertzian_dipole: z 0.35 0.380 0.0 pulse

(z — ось поляризации диполя (для двумерного случая (2D TMz-режим) допустима только z), pulse — метка формы импульса, излучаемого антенной)

Приемная антенна располагается обычно на небольшом неизменном удалении от приемной, которое называется базой антенного блока (такой вариант взаимного расположения антенн называется «common-offset»). В качестве базы выбираем расстояние 10 см, таким образом, горизонтальная координата равна 35 + 10 = 45 см = 0,45 м (5 см справа от середины области):

#rx: 0.45 0.38 0.0

Интервал моделирования

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

Определим примерное время, требуемое сигналу в рассматриваемом случае, приняв расстояние от антенн до мишени $h = 18 \ \ см$:
$t \approx {2 \cdot h \over {v}} = { 2 \cdot 18 \over 17,3} = 2,1 \ \ нс$
Учитывая, что вершина гауссова импульса радара с центральной частотой 1 ГГц сдвинута относительно начала оси времени на 1 нс, выбираем временное окно 5 наносекунд:

#time_window: 5e-9

Моделирование

Таким образом, содержание входного файл таково:

hello.in
#domain: 0.80 0.60 0.002
#dx_dy_dz: 0.002 0.002 0.002
#time_window: 5e-9
#material: 3 0.01 1 0 sand
#waveform: gaussian 1 1e9 pulse
#hertzian_dipole: z 0.35 0.380 0.0 pulse
#rx: 0.45 0.38 0.0
#box: 0 0 0 0.80 0.38 0.002 sand
#cylinder: 0.40 0.1 0 0.40 0.1 0.002 0.06 pec


Запускаем процесс моделирования:

python -m gprMax models\hello.in

Для выполнения моделирования используется метод конечных разностей во временной области (FDTD, finite-difference time-domain) (базовый алгоритм метода был предложен Кэйном Йи (Kane Yee)), в честь которого ячейки, на которые разбивается модель, названы ячейки Йи (Yee). Численное решение получается во временной области путем решения уравнений Максвелла для каждой ячейки.

В двумерном случае (2D TMz-режим) рассчитываются только компонента $E_z$ электрического поля и компоненты $H_x$ и $H_y$ магнитного поля.

При превышении объема доступной памяти выдается сообщение о нехватке памяти для выполнения моделирования:

image

Если какой-то из объектов вышел за пределы области моделирования, то выдается сообщение об ошибке:

image

Для выполнения моделирования с описанной моделью потребовалось всего лишь около 56 МБайт оперативной памяти (при уменьшении шага в два раза — до 1 мм — требования к памяти повышаются до 99 МБайт).

После завершения моделирования в папке models появляется файл hello.out, содержащий результаты моделирования в формате HDF5, специально разработанном для хранения численных данных.

Построение трасс

Для визуального представления результатов построим трассы:

python -m tools.plot_Ascan models\hello.out

Каждая представленная в открывшемся окне трасса (A-scan) отображает временной график одной из составляющих электромагнитного поля в точке расположения приемной антенны:

image

На трассах видна прямая волна, пришедшая напрямую от передающей антенны (DW), и отраженная от мишени волна (RW).

Горизонтальная ось времени связана с глубиной мишени, отразившей сигнал, через скорость электромагнитной волны в веществе.

А что будет, если поставить во входном файле команду cylinder перед командой box?

image

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

Построение профиля

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

Изменим описание антенн, сместив их к началу горизонтальной оси:

#hertzian_dipole: z 0.10 0.38 0.0 pulse
#rx: 0.20 0.38 0.0

Зададим шаг перемещения антенн, равный 1 см = 0,01 м:

#src_steps: 0.01 0 0
#rx_steps: 0.01 0 0

Таким образом, содержание входного файл таково:

hello.in
#domain: 0.80 0.60 0.002
#dx_dy_dz: 0.002 0.002 0.002
#time_window: 5e-9
#material: 3 0.01 1 0 sand
#waveform: gaussian 1 1e9 pulse
#hertzian_dipole: z 0.10 0.38 0.0 pulse
#rx: 0.20 0.38 0.0
#src_steps: 0.01 0 0
#rx_steps: 0.01 0 0
#box: 0 0 0 0.80 0.38 0.002 sand
#cylinder: 0.40 0.1 0 0.40 0.1 0.002 0.06 pec


Запускаем моделирование в пакетном режиме:

python -m gprMax models\hello.in -n 50

(50 — число шагов перемещения радара).

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

image

После окончания моделирования в папке модели присутствуют 50 файлов hello1.out… hello50.out.
Объединяем эти файлы в файл hello_merged.out командой:

python -m tools.outputfiles_merge models/hello

Строим профиль:

python -m tools.plot_Bscan models\hello_merged.out Ez

(Ez — составляющая электромагнитного поля, по которой строим профиль — компонента, непосредственно преобразуемая в напряжение)

image

Как видим, на профиле сверху отображается горизонтальная полоса, обусловленная прямой волной, а ниже — характерная гипербола, обусловленная отраженной волной и отображающая изменение расстояния от мишени до георадара при его перемещении.

Легенда справа показывает соответствие цветов и напряженности поля $E_z$ (показанной на трассе):
красный — $E_z > 0$
белый — $E_z = 0$
синий — $E_z < 0$
Анализируя такие профили, можно делать выводы о глубине, размерах и даже форме мишеней, причем для этого используются и нейронные сети.

image

Трассы и профиль на дисплее георадара (кадр из британской телепередачи «Команда времени»)

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

Создадим в нижней части модели второй слой песка с диэлектрической проницаемостью ${\epsilon}_r = 9$:

hello.in
#domain: 0.8 0.6 0.002
#dx_dy_dz: 0.002 0.002 0.002
#time_window: 5e-9
#material: 3 0.01 1 0 sand_1
#material: 9 0.01 1 0 sand_2
#waveform: gaussian 1 1e9 pulse
#hertzian_dipole: z 0.10 0.38 0.0 pulse
#rx: 0.20 0.38 0.0
#src_steps: 0.01 0 0
#rx_steps: 0.01 0 0
#box: 0 0 0 0.80 0.38 0.002 sand_1
#box: 0 0 0 0.80 0.10 0.002 sand_2


image

Как видно, ниже «следа» прямой волны (DW) появился линейный сегмент отраженной от границы раздела двух слоев песка волны (RW).

За рамками этого примера остались такие возможности gprMax как трехмерное моделирование, сложные формы поверхности, детальные модели антенн, учет дисперсии электромагнитных волн… При этом пакет можно применять не только для моделирования работы георадара, но и для исследования распространения электромагнитных волн в различных средах, а для ускорения расчетов использовать технологию NVIDIA CUDA, обеспечивающую прирост скорости в десятки раз по сравнению с параллельными вычислениями на CPU посредством OpenMP. Для повышения гибкости моделирования во входном файле можно размещать блоки кода на Python.

Некоторые примеры использования пакета gprMax:


Полезные ссылки:

официальный сайт gprMax
Руководство пользователя gprMax
gprMax на YouTube
AdBlock похитил этот баннер, но баннеры не зубы — отрастут

Подробнее
Реклама

Комментарии 15

    0
    Интересно видит ли этот прибор песок с повышенным содержанием частиц золота?
    +1
    Спасибо за интересную статью! Давно интересовался GPRmax.
    Можно ли все-таки задавать материалы с дисперсией? Ведь в реальности на такой широкой полосе даже у простой воды проводимость хорошо вырастает, а значит импульс начнет терять свои ВЧ-составляющие в отличие от модели. Кстати все оценки скоростей также будут приблизительными для материалов с дисперсией (особенно если на пути сигнала несколько слоев).

    Антенна явно моделируется неверно, и выглядит скорее как идеальный источник. Обычный диполь это резонансная антенна, которая работает хорошо на очень узкой полосе частот, а значит при попытке использовать ее для передачи импульсов, у которых полоса огромная, форма импульса будет неизбежно искажаться. Только некоторые типы сверхширокополосных антенн могут похвастать способностью передавать импульс без значительных фазовых искажений. Вот на картинке пример подачи monocycle pulse на разные типы антенн и что в итоге излучается в пространство:

    То есть в вашем примере с диполем, вместо идеального импульса, мы должны увидеть примерно то, что на картинке вверху слева. Надеюсь другие антенны задаются более реалистично.
      +1

      Спасибо!
      Я же и написал, что описываю простой учебный пример, в котором и среда без дисперсии, и антенна идеальная. В описании антенны я указал, что в пакете имеется библиотека реальных антенн, которой и можно воспользоваться уже при трехмерном моделировании. А в конце статьи я указал, что пример — простейший, и возможностей намного больше. Как задавать среды с дисперсией описано здесь — http://docs.gprmax.com/en/latest/input.html (смотрите add_dispersion_debye etc)

        +1
        Понял, надеюсь напишите статью и про 3d моделирование, было бы интересно!
        +1
        Судя по тому, что зарегистрированные импульсы остались ограниченными по времени (т.е. источник все-таки остался широкополосным), и по тому, что задача двумерная (т.е. бесконечно протяженная по направлению, перпендикулярному экрану), можно предположить, что в данном случае источником была бесконечная нить тока, которая в двумерном случае моделируется как источник в одной ячейке метода FDTD.
          +1

          Вот как представляется диполь Герца в gprMax:
          image

        +1
        У Вас принципиальная ошибка при моделировании распространения сигнала по трассе, а именно в предположении, что диэлектрическая проницаемость среды является вещественной константой.
        На самом деле диэлектрическая проницаемость является комлекснозначной функцией зависящей от частоты колебаний (минимальная модель такой среды может быть построена с использованием однокомпонентной модели Дебая).
        Значительная частотная дисперсия наблюдается уже на частотах выше 1 ГГц, поэтому приведенные Вами графики не будут соответствовать действительности.
        В частности, используемое Вами неявно, приближение групповой скорости, будет давать значительную ошибку в оценке расстояния до диэлектрической неоднородности.
        Если Вас это интересует, то можем связаться
          +1

          Спасибо за внимание к статье!
          Так я же и не претендовал на то, чтобы в учебном примере описать все возможности пакета — о чем и написано в статье. Дисперсия вполне себе моделируется посредством add_dispersion_*, но я умышленно не стал усложнять этот пример.

          0

          возможно, подскажете — а какие-то DIY аппаратные системы существуют в разумной ценовой категории(100-200евро, которые не жалко потратить на то, чтобы это попробовать вживую) для работы с gprMax?

            0

            gprMax — это симулятор георадара, непосредственно для работы с ним сам радар не требуется. gprMax используется, когда требуется, например, исследовать зависимость радарных сигнатур от тех или иных факторов. Это намного быстрее сделать в симуляторе, чем возиться с реальным радаром.
            Но DIY-проекты георадаров существуют — например, вот такой:
            https://hackaday.io/project/4440-open-ground-penetrating-radar

              0
              Спасибо за ответ.
              Я проверил (не вчера) не один из из DYI проектов георадаров (ибо давно интересуюсь этой тематикой) — вопрос был о разумном железе для хобби за 200 евро.

              И, возможно, знаком с симуляторами(порою, более сложными, чем георадар) намного лушче вас — поэтому о симуляторах и не спрашивал;)
                0

                Если бы Вы ещё и свои мысли ПОЛУШЧЕ выражать научились, было бы вообще прекрасно :-)))

            0
            Самой большой проблемой георадара является обычная вода — стоит намочить грунт и глубина зондирования становится смешной (первые десятки см).
            Кроме того точность измерения расстояния сильно плывет.
            С практической точки зрения — все применение георадара это поиск арматуры в ЖБ плитах.
            Если нужно без буровых работ узнать распределение слоев пород (построить 3D модель распределения сопротивления породы), то гарантированно рабочий метод это электроразведка становлением поля.
            Стоимость аппаратуры меньше, а диапазон глубин от метра до 10 км. Результат принимается геологией (в отличии от георадара)
            тут море картинок sibgeotech.ru/rezultaty
              0

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

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

            Самое читаемое