Рекурсивное миграционное преобразование радарограммы в системе подповерхностной радиолокации Текст научной статьи по специальности «Физика»

ИНФОРМАТИКА, ВЫЧИСЛИТЕЛЬНАЯ ТЕХНИКА И УПРАВЛЕНИЕ

УДК 621.376.3

РЕКУРСИВНОЕ МИГРАЦИОННОЕ ПРЕОБРАЗОВАНИЕ РАДАРОГРАММЫ В СИСТЕМЕ ПОДПОВЕРХНОСТНОЙ РАДИОЛОКАЦИИ

© 2016 А. Б. Борзов, Л. В. Лабунец, А. В. Зотов, Р. Е. Копейкин

Московский государственный технический университет имени Н.Э.Баумана

Статья поступила в редакцию 24.05.2016

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

ВВЕДЕНИЕ

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

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

Борзов Андрей Борисович, доктор технических наук, профессор, заведующий кафедрой СМ-5. E-mail: borzov@rambler.ru

Лабунец Леонид Витальевич, доктор технических наук, старший научный сотрудник, профессор кафедры СМ-5. E-mail: labunets@bmstu.ru

Зотов Алексей Владимирович, студент кафедры СМ-5. E-mail: aleksei.zotov.ml@gmail.com Копейкин Роман Евгеньевич, ассистент кафедры СМ-5, аспирант. E-mail: r.kopeykin@gmail.com

ДИСКРЕТНАЯ МОДЕЛЬ ВОЛНОВОГО ПОЛЯ

Натурные измерения выполняют с помощью однопозиционного георадара, излучатель и приемник которого практически совмещены в пространстве. Регистрируемый массив данных описывает рассеянное средой трехмерное волновое поле х, г), где t - время, х и г - пространственные координаты, направленные вдоль и в глубину исследуемой поверхности. Соответствующую дискретную модель поля, заданного на ортогональном растре, удобно представить в виде трехмерной матрицы

1(п1, п2, п3) = Цп1 М, п2 Дх, п3 Дг),

0 < пг < Ыг -1, г = 1,2,3. (1)

Здесь Дt и Дх , Дг - интервалы временной и пространственной дискретизации, согласованные по теореме отсчетов Котельникова - Шеннона [1, с. 58] с верхними граничными частотами Ш., г = 1,2,3 (частотами Найквиста).

Фурье-спектр поля вычисляют по формуле

ф(ю, кх, кг )= | dt | | / (, х, г)ехр (- у ю t -/кхх - ]к2г ,(2)

если ю < Ж1, |кх < ж2, |кг <Ж3 и ф(ю ,кх, к г) 0 в противном случае. Здесь, в отличие от геофизической транскрипции [2, с. 33; 3, с. 89], принята стандартная радиофизическая форма записи непрерывного преобразования Фурье. Также следует отметить, что ширина полосы частот спектра (2) определяется значениями 2Ш., г = 1,2,3.

Если выполняются ограничения Дt <%/ Ж1 , Дх <%! Ж2 , Дг <я/ Ж3, то периодическое продолжение спектра (2) в частотном пространстве (ю, к х, к г), как результат процедуры дискретизации данных (1), не приводит к эффекту наложения спектральных копий [1, с. 57]. В случае равенств в указанных выше ограничениях опорная область спектра (2) и ее периодические копии образуют плотную упаковку в частотном пространстве,

что обеспечивает наименьшую плотность узлов ортогонального растра в исходном пространстве (t, x, z). Это, в свою очередь, позволяет реконструировать значения поля

N j-1 N 2-1 N 3-1

f(l, X, z) = ¿ ¿ ¿ f(n 1, n 2, n 3 )(l - n jAí, x - n 2Ax, z - n 3Az)

nj = 0 n 2 = 0 n 3 = 0

(3)

с помощью оптимальной интерполирующей функции

( ) = sin ( 11)sin (2 X)sin ((3 z) (4)

ф \ 1 , x, z Г — (4)

V ' W11 W 2 x W 3 z

с нулевой или приемлемо малой ошибкой e(l,x,z) — f(l,x,z)- f (l,x,z).

МИГРАЦИОННОЕ ПРЕОБРАЗОВАНИЕ

В ряде работ [4, 5] отмечалось, что наиболее эффективные алгоритмы миграции сигналов, сформированных элементарными отражателями однородной среды, основаны на суммировании амплитуд сигналов каждого рассеивателя вдоль дифракционной гиперболы [6 - 8] или методе Стол-та восстановления геофизического изображения зондируемой среды [2, с. 44; 3, с. 95]. Указанные алгоритмы устойчивы к шумам измерений, не вносят значимых интерференционных искажений в процессе преобразования и требуют относительно небольших вычислительных затрат.

Радарограмма представляет собой двумерное сечение матрицы (1) r(t, x) = f (t, x, z = 0), т. е. содержит результаты экспериментальных измерений поля на поверхности среды z = 0. Миграционное преобразование радарограммы r(t, x) позволяет восстановить волновое поле в вертикальном разрезе среды p(t, z) = f (t = 0, x, z) для начального момента времени t = 0. Основой такого рода реконструкции является скалярное волновое уравнение [3, с. 105] в частотном пространстве (ю, kx)

dF (ю , k\z )

—г-.—V- — jkzdz.

F (ю , kx\z )

Здесь F (ю , kx|z) - спектр Фурье двумерного сечение волнового поля f (t, x, z) для фиксированного значения заглубления z > 0 ; волновое число k z удовлетворяет соотношению

k1 +k2 =

4 ю

2

(5)

V" (х, г )

V (х, г) - скорость распространения ЭМВ в среде. Решение указанного выше однородного линейного дифференциального уравнения с разделяющимися переменными имеет вид дисперсионного соотношения

р(® ,кх\г)=Я(ю ,кх)ехр{е(ю ,кх\г)}, (6)

е(ю ,kx|z)— j

4ю^ - к 2 J| 2( \ kx

0 1 v (x, u)

\1/2

du, (7)

где Я(ю , кх)= F(ю , кх|г = о) - спектр Фурье радарограммы.

В соответствии со свойствами преобразования Фурье полученное решение показывает, что волновое поле / (, х, г) на глубине г получают с помощью сдвига гармоник радаро-граммы г (г, х ) = / (, х, 2 = 0) на величину накопленной фазы е ( ю , кх | г ). Миграционное преобразование методом фазового сдвига [2, с. 42; 3, с. 107] предполагает вычисление обратного непрерывного преобразования Фурье с помощью дисперсионного соотношения (6). С учетом тождества ехр ((ю г) = 1 для t = 0 получим

ш 2

Р (х , г ) = / ( = 0, х , г )= I Я (кх = г )еХР {кхх }кх ,

-Ш 2

Ш1

Я(кх ,2)= I Я(ю ,кх )ехр{е(ю ,кх|г)}ю. (8) -ш1

Иными словами, волновое поле в вертикальном разрезе среды р(х,г) представляет собой обратное непрерывное преобразование Фурье спектрограммы Я (к х , г ) по волновому числу кх. Спектрограмма, в свою очередь, является результатом фазовых сдвигов спектра Фурье радарограммы для текущих глубин г и последующего интегрирования по частоте ю.

ОТОБРАЖЕНИЕ СТОЛТА

В однородной среде скорость распространения волнового фронта не зависит от координат вертикального разреза V (х, г ) = V. В этом случае накопленная фаза (7) увеличивается линейно в зависимости от глубины

е(ю > кх\г )= ^ ^ ^ю) - кх = кг (ю > кх )г .

Соответствующее интегральное преобразование Фурье, основанное на дисперсионном соотношении (6), приобретает вид

W2 Wj

p (x, z)— j dk x j R ( ю , kx )exp {jkxx + jkz ( ю , kx )z }dco .

-W2 -W,

Замена переменных в полученном равенстве

с помощью подстановки ю = vл^' к^ + к^/2 приводит к обратному преобразования Фурье

Ш 2 Шз

р (х, г) = 3 (кх, кг ) | dk х | Я{ю(кх , кг ), кх }ехр{{кхх + jkzz }кг .

- W2 -W3

для спектра R(ro, kx) радарограммы, рассчитанного на нелинейном растре |ra,(kx , kz), kx} . Здесь

i \ vk z J(,kz)= z

k 2 + k 2

якобиан преобразования. Последовательность вычислений, реализующая алгоритм миграционного отображения Столта, представлена в Приложении 1. Результаты применения алгоритма для тестового изображения, хранящегося в базе данных на сайте [9] обучающего курса профессора Guy G. Drijkoningen факультета геоинженерных технологий университета Дельфта, иллюстрирует рис. 1. Вычисления выполнялись для следующих параметров: v =1550 м/с, At = 0,004 с, Ax = 50м, Az = 6,2м, Nl = 1001, N2 =141, N3 =1001. В данном примере алгоритм Столта отображает дифракционные гиперболы исходных сейсмических данных в два точечных источника волнового поля геофизического разреза однородной среды.

С методической точки зрения радиофизическая интерпретация нелинейного отображения Столта состоит в преобразовании радарограммы двумерным фильтром. Передаточная функция такого фильтра в частотном пространстве (ю, kx) имеет вид [2, с. 34; 3, с. 107]

H (ю , kx|z)

= exPi jz

где фиксированное значение заглубления г > 0 трактуется как параметр.

ОТОБРАЖЕНИЕ ГАЗДАГА

Рассмотрим многослойную среду, в которой скорость распространения ЭМВ не зависит от горизонтальной координаты х и изменяется локально постоянно от слоя к слою, т. е.

'Ух, 0 < 2 < 2!

'(x, z ) = v (z )

V2 ' z 1 ^ z < z 2

, (9)

VL > zL-1 ^ z < zL

Рис. 1. Миграционное преобразование Столта: а - исходные данные; б - мигрированные данные

Аб(ю , kx|z )=б(ю , kx|z + A z )-

-Э(ю , kx\z ) = kz (ю , kx ) A z, (10)

kz

(11)

где V, и - скорость распространения ЭМВ и заглубление нижней границы 1-ого слоя, I = 1, 2, ..., Ь.

Изменение скорости распространения зондирующего излучения в среде естественным образом учитывает миграционное отображение Газдага [2, с. 42]. В рамках этого подхода дисперсионное соотношение (6) представляют в рекуррентной форме

^ (со , кх\г + Д 2 )= ^ (со , кх\г )ехр{у Дб(со , кх\г )},

где с учетом формулы (7) фазовый сдвиг де( га , кх\г) спектра Фурье нисходящего волнового поля при переходе от заглубления г к заглублению г + Дг имеет вид

Алгоритм Газдага прямого отображения по нисходящим волнам г (, х р (х, 2 ), представлен в Приложении 2.

Миграционное преобразование радаро-граммы с помощью алгоритма фазового сдвига предполагает наличие априорной информации о распределении скорости распространения ЭМВ по глубине. Такого рода информацию получают на этапе предварительной обработки данных с помощью современных методов структурного анализа радарограммы [10].

Результаты применения алгоритма для сейсмических измерений, хранящихся в базе данных на сайте [9], демонстрирует рис. 2. Вычисления выполнялись для следующих параметров: Дt = 0,004 с, Дх = 20 м, Дг = 4 м, Ы1 = 3001, Ы2 = 1024, Ы3 = 1046. Представленное на рис. 3 локально постоянное распределение скоростей (9) по Ь = 14 слоям среды загружалось из отдельного файла базы данных [9].

Рис. 2. Миграционное преобразование Газдага: а - исходные данные; б - мигрированные данные

х, м

Рис. 3. Распределение скорости распространения волн по слоям неоднородной среды

РЕКУРСИВНОЕ ОТОБРАЖЕНИЕ СТОЛТА

Вычисление интеграла (7) с учетом модели (9) дает кусочно-линейную зависимость накопленной фазы от глубины

е(ю

* (2 >( *(L )(

< Z < Z 1

-де(2)(

ю ,к.

Z J < z < z 2

L-1

-Z Ае1

i =1

С+1)

< Z < Z ,

Ю .

)z },

(12)

где г0 = 0, Де(г+1)(с , кх )= {к(+1)(с , кх)-к() (с , кх )} 2/ - коррекция набега фазы за счет изменения скорости распространения ЭМВ в среде. Значения волнового числа кг находят по формуле (11) для скорости V (2 )= у.

Кусочно-линейная зависимость (12) позволяет представить дисперсионное соотношение (6) в следующем виде

р(с , кх|2)= в(с , кх\21-1 )ехР{к2)( I = 1, 2, ...,Ь. (13)

Здесь спектр Фурье В (со , к х | 2/-1) восходящего волнового поля для верхней границы 2 /-1 текущего I -ого слоя вычисляют рекурсивно по слоям среды. В первом слое 2 0 < 2 < 21 спектр восходящего волнового поля на поверхности среды представляет собой радарограмму

в(с,кх\2о )= Я(с,кх )

В последующих слоях г /-1 < 2 < 21, I = 2, ..., Ь спектр восходящих волн на глубине 21-1 рассчитывают по предыдущим спектральным оценкам нисходящего поля

в (с , к х | 2 /-1 ) =

= F (со , кх | 2 /-1) ехр {-)к(/)(со , кх ) 2 /-1}. (14)

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

Анализ рекурсии Столта выполнялся для тестовой радарограммы в виде суммы двух дифракционных гипербол и гауссовского шума (рис. 4.а). Скорости распространения волн (равные тангенсам углов наклона асимптот) в двух слоях среды принимали значения 1 и 2 соответственно. Миграцию гипербол в точечные источники волнового поля геофизического разреза двухслойной среды иллюстрирует рис. 4.б.

Результаты тестирования рекурсивного алгоритма Столта демонстрируют экономию вычислительных затрат по сравнению с алгоритмом Газдага пропорциональную количеству слоев среды.

■

. ... •. ' . v. .-

•V w^/ - . v.. _

Рис. 4. Рекурсивное миграционное преобразование Столта: а - исходные данные; б - мигрированные данные

Приложение 1. Миграционное преобразование Столта.

Шаг 1. Вычислить двумерное дискретное преобразование Фурье (ДПФ) радарограммы r(nl5 n2)= r(n1At, n2Ах) по времени t и координате x.

R ( ) V V ( ) I ■ 2 я . 2 я I

R{m 1, m 2) = Z Z r\n 1> n ^exP |-J^n 1 m 1 - J^n 2 m 2 I,

n1= 0 n 2 = 0

. 2 я

2 я

N1

N 2

(П.1)

где Я(т1, т2)= Я(т1 Дю, т2Дкх),

0 < т1 < -1, I = 1,2,

- дискретная модель с пектра Фурье Я(0,, кх), рассчитанная на ортогональном растре с интервалами дискретизации

2Wl 2 л

А®:

А к.

N1 -1 (N1 -1 )А t

2W 2 = 2 я ' " N2 -1 (N2 -1)

в частотном домене 0<Q<2W1 =

0 < кх < 2W2 =

2 я

А х

2 я

д7

Шаг 2. Выполнить дискретизацию пространственной частоты кг в интервале 0 < к2 < 2W3 = 2 л/ Д2 с шагом

и

2W 3 2%

A kz —-— 7-1-.

z N3 -1 (N3 -1) Az

(П.2)

Отобразить узлы ортогонального растра ( т 2 Д кх , т3 Д кг) на пространство частот {ю(т2,т3), т2Дкх} с помощью нелинейного соотношения (5).

(т 2, т 3 )) 1(т 2 Д кх )2 +(т 3 Д кг )2-(П.3)

ю

Преобразовать частоты ю ( m 2, m 3 )

косновно-

2 %

му домену 0 < О < — ДПФ радарограммы

Д

О,(т2,т3) = {ю(т2,т3)}},

Дг

0 < тг < N -1, I = 2,3, (П.4)

где функция |а}Ь означает остаток от деления величины а на целое число порций Ь.

Шаг 3. Выполнить интерполяцию отсчетов ДПФ радарограммы на нелинейном растре

2, т 3), т 2 Д кх }

Я{(т 2, т 3 ), т 2 }=

N >-1

р(m 2,m 3 ^

vm3Ak z

2Л] (m2Akx)2 +(m3 Akz)

=R7{ü(m2 ,m3), m2}.

к и k.

x z

p(n2,n3) = p(n2Ax,n3Az)

N

1 N 2-1 N 3-1 |

— E E Pim2, m^exP| j

N 3 m2—0 m3 — 0 ^

2 %

2 %

Приложение 2. Миграционное преобразование Газдага.

Шаг 1. Вычислить по формуле (П.1) двумерное ДПФ ^(ш1, т2) радарограммы.

Шаг 2. Инициализировать цикл по глубине г. п3 = 0, Дг = гь/((3 -1).

Вычислить значения дискретной модели волнового поля на поверхности объекта локации

г = 0, р(п2,п3 ) = р(п2Дх,п3Дг)= г(0,п2Дх)= г(0,п2). Шаг 3. Цикл по глубине г. п 3 = п 3 + Ь г = г + Дг.

Шаг 4. Инициализировать цикл по волновому числу к .

m 2 — -1, A kx —

2W 2

2 %

N2 -1 ((2 -1 )Ax :

k x — A k x.

1

E R (m!, m 2 )p{o(m 2, m 3 )- m1Aю},

да1—0

В соответствии со свойством дуальности ДПФ [1, с. 93] оптимальная интерполирующая функция p(D) в частотном пространстве имеет вид, аналогичный интерполирующей функции (4) в исходном пространстве

/„ч sin(т R Q/ 2 ) sin (%Q/ Aю2

P(Q)— V n/o n/л , (П.5)

т R Q/ 2 %Q/Aю

где тR — (Nj -1) Al - длительность отраженных сигналов радарограммы.

Шаг 4. Вычислить оценку ДПФ Р ( m 2 , m 3) волнового поля p(x, z) в вертикальном разрезе среды с помощью операции масштабирования

Шаг 5. Цикл по волновому числу кх.

т 2 = т 2 +1, к х = к х + Дк х.

Обнулить начальное значение спектрограммы

Я (т 2, п 3 )= Я (кх , г )= 0.

Шаг 6. Инициализировать цикл по частоте ю .

п л 2Ш1 2 п 0

т 1 = 0, Дю= —-7 = ^-ГП^ ю = 0.

1 N1 -1 (Ж 1 -1 )Д г

В соответствии с выбранным алгоритмом численного интегрирования выражения (8) сформировать массив квадратурных коэффициентов С(т 1), 0<т 1 <N1 -1. Вычислить начальный

спектр Фурье F( 0 , кх| г)волнового поля для текущих значений глубины г = п 3 Дг и волнового

числа k x — m 2 Ak x

F (m 1 , m ^ n 3 ^F (ю , k^z )=R (m 1 , m 2 ).

Обновить значение спектрограммы

R (m 2 , n 3 ) = R (m 2 , n 3 ) + Aro C (m 1)) (m 1 , m 2 n 3 ).

Шаг 7. Цикл по частоте ю.

m 1 — m 1 +1, ю — ю + Aю.

Вычислить по формулам (9) и (10) текущий фазовый сдвиг

Ae(m 1

m 2 n 3 —

)=A^ j2"^} -(m 2 Akx )2 .

Шаг 5. Реконструировать дискретную модель волнового поле в вертикальном разрезе среды с помощью обратного ДПФ по волновым числам

.. n 2 m 2 + j — n 3 m 3 I,

N 2 N 3

(П.6)

0 < nt < Ni -1, i — 2,3.

Вычислить текущий спектр Фурье F (ю , k x | z) волнового поля

F(ni! , m2 | n3 )=F( m¡ -1, m21 n3 )exp {Ae( 1П1 , m2 n3 )}.

Обновить значение спектрограммы R (m 2,« 3 ) = R (m 2, и 3 ) + Aro C (m 1)) (m 1 , m 2 « 3 ).

Шаг 8. Критерий продолжения цикла по частоте ю . Если m 1 < N1 -1, то идти к Шагу 7.

Шаг 9. Критерий продолжения цикла по волновому числу k x. Если m 2 < N 2 -1, то идти к Шагу 5.

Шаг 10. Критерий продолжения цикла по глубине г. Если п 3 < N 3 - 1, то идти к Шагу 3.

Шаг 11. С помощью обратного ДПФ спектрограммы по волновому числу к х реконструировать дискретную модель поля в вертикальном разрезе среды.

Р (п 2, п 3 )= р (п 2 Д х, п 3Дг )=

N 2 -1

^ X R (m 2 , n 3 )exP

N 9 „, _ _п

(

2 m 2 = 0

2 %

/-n 2 m 2

N 2 2 2

Л

0 < n 2 < N 2 -1.

нисходящего волнового поля для дискретных значений глубины г ^- < т 3 Дг < г ^

F(т2,т3) = В(т2,тъ\ гм)ехр({3Дкг г};

• масштабировать отсчеты ДПФ нисходящего волнового поля для дискретных значений глубины

г1 _1 < т 3 Д г < г 1

р (m 2,m3) =

V,m 3 А к^

2Л] (m 2 Akx)2 + (m 3Akz)

= F (m 2, m 3 ).

Приложение 3. Рекурсивное миграционное преобразование Столта.

Шаг 1. На ортогональном растре с интервалами дискретизации Дю и Д кх вычислить

по формуле (П.1) двумерное радарограммы.

Шаг 2. Инициализировать цикл по слоям среды. Выбрать шаг дискретизации Дг глубины. Сформировать массивы V. и г, I = 1, 2, ..., Ь для скоростей распространения ЭМВ и заглублений нижних границ слоев. По формуле (11) для скорости V(г ) = ¥1 сформировать массив волновых чисел

к(г1 \т 1, т 2)= к(г1)(т 1 Дю, т 2 Дкх), 0 < т1 < N. -1, 1 = 1,2.

Вычислить ДПФ восходящего волнового поля на поверхности среды I = 0, г0= 0.

В(т 1,т2] г0 )= Я(т 1,т2),

0 < т1 < Ni -1, I = 1,2.

Шаг 3. Цикл по слоям среды I = 1+1. Шаг 4. В текущем слое г ^- < г < г ^ реализовать нелинейное преобразование Столта. Для этого:

• выполнить дискретизацию пространственной частоты кг с шагом Дкг в интервале 0 < к2 < 2Ш3 = 2п/ Дг по формуле (П.2);

отобразить узлы ортогонального растра

(т 2 Д кх , т 3 Д кг) на пространство частот

{ю(т2,т3), т2 Дкх} по формуле (П.3) для скорости распространения ЭМВ V = V1;

• преобразовать частоты ю (т 2,т 3) к основ-

2 п

ному домену 0 < о < — по формуле (П.4);

• для дискретных значений глубины г1 -1 < т 3 Дг < г1 с помощью функции (П.5) выполнить интерполяцию отсчетов ДПФ восходящего волнового поля на нелинейном растре

{(т2,т3), т2Дкх}

N

В(т 2 , т э| г/-1 )= Е В (m1, т 2 г/-1 )Р{о(т 2 , т 3 )- т1Дю};

т^=0

по формуле (13) вычислить отсчеты ДПФ

Шаг 5. Рекурсия спектральных оценок в пространстве частот (р, kx) на ортогональном растре {mjAw, m2Akx}:

• по формуле (13) вычислить ДПФ нисходящего волнового поля для нижней границы z, текущего слоя среды

F[mi ,m^ )= B(mi ,m2I zi-i )exp{./k!')(ml,m2)z/

• по формуле (11) для скорости v(z ) = Vl+1 распространения ЭМВ в последующем слое среды сформировать массив волновых чисел

к(z(+1 )(m 1, m2) 0 < mt < Nt - L i = 1,2;

• по рекурсивной формуле (14) вычислить ДПФ восходящего волнового поля для верхней границы z, последующего слоя среды

B(m1 ,mi\zi )= F(m1 ,mi\zi )exp{-/k(z(+1 ^ m1 ,m2 )z(}.

Шаг 6. Критерий продолжения цикла по слоям. Если (< L, то идти к Шагу 3.

Шаг 7. С помощью обратного ДПФ спектра Р(m2,m3) по формуле (П.6) реконструировать дискретную модель p(n2 , n3)волнового поля в вертикальном разрезе среды.

СПИСОК ЛИТЕРАТУРЫ

1. Даджион Д., Мерсеро Р. Цифровая обработка многомерных сигналов: Пер. с англ. М.: Мир. 1988. - 488 с.

2. Клаербоут Д. Ф. Сейсмическое изображение земных недр: Пер. с англ. М.: Недра. 1989. 407 с.

3. Claerbout J.F. Basic Earth Imaging. Cambridge: Free Software Foundation. 2010. 218 p.

4. Зеленков А.В., Зеленков С.А. Выбор метода миграции для обработки сигналов подповерхностной радиолокации // Электроника и электротехника. Каунас: Технология, 2004.№ 3(52). С. 41 - 46.

5. Klokov A., Sato M. Application of 3D- migration algorithm to GPR on an irregular ground surface // Proc. of the IGARSS, 2011. P. 870 - 873.

6. Yilmaz Oz., Doherty S.M. Seismic Data Analysis: processing, inversion and interpretation of seismic data. V. 1. Society of Exploration Geophysics, 2001.

7. Линников О. Н., Сосулин Ю. Г., Толмазов Б. Б. Обработка сигналов и изображений в импульсном радиолокаторе подповерхностного зондирования // Радиотехника и электроника. 2004, Т. 49. № 7. С. 824 - 839.

8. Юфряков Б.А., Суриков Б.С., Сосулин Ю.Г., Линников О.Н. Метод интерпретации данных радиолокаторов

подповерхностного зондирования // Радиотехника и электроника. 2004. Т. 49. № 12. С. 1436 - 1451. 9. Practical Course Introduction to Reflection Seismology [Электронный ресурс] - URL: http://geodus1. ta.tudelft.nl/PrivatePages/G.G.Drijkoningen/intro_

seismology (дата обращения 15.04.2016).

10. Лабунец Л.В. Симаков В.В. Структурный анализ сигналов в системах подповерхностной радиолокации // Электромагнитные волны и электронные системы. 2013. Т. 18. № 8, С. 49 - 68.

RECURSIVE MIGRATION RADARGRAM TRANSFORM IN SUBSURFACE RADAR SYSTEMS

© 2016 A.B. Borzov, L.V. Labunets, A.V. Zotov, R.E. Kopeykin

Bauman Moscow State Technical University

A computationally effective recursive algorithm of migration transform of radargram in subsurface radar system is proposed. The algorithm is based on the analysis of descending and ascending wave fields in subsequent layers of inhomogeneous surface via a combination of Stolt and Gazdag migrations. Keywords: subsurface radar, digital signal processing, migration transform.

Audrey Borzov, Doctor of Engineering Sciences, Professor,

Head at the SM-5Department. E-mail: borzov@rambler.ru

Leonid Labunets, Doctor of Engineering Sciences, Senior

Research Fellow, Professor at the SM-5 Department.

E-mail: labunets@bmstu.ru

Alexei Zotov, Student of SM-5 Department.

E-mail: aleksei.zotov.ml@gmail.com

Roman Kopeykin, Assistant Lecturer at the SM-5 Department.

E-mail: r.kopeykin@gmail.com