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

ВВЕДЕНИЕ

Л. Вердонк с соавторами в начале своей статьи [1] отметили, что при археологических исследованиях, как и в других областях, все чаще стали использовать многоканальные георадары с многоэлементными антеннами (антенными решетками) непрерывного излучения со ступенчатой перестройкой частоты (SFCW) или с импульсными многоэлементными антеннами. Такие системы позволяют быстро проводить георадиолокационную съемку на большой площади, причем с достаточно высокой частотой дискретизации (близкой к предписанной теоремой Найквиста). Эти так называемые антенные решетки в большинстве своем состоят из антенн, закрепленных на большой раме и расположенных таким образом, чтобы расстояние между линиями съемки было меньше физического размера антенн (например: каждая приемная антенна может быть предназначена для записи сигналов двух соседних передающих антенн; комбинации «передатчик – приемник» могут быть повернуты на 90 град. друг по отношению к другу и на 45 град. к направлению съемки).

Однако для приобретения больших многоканальных систем требуются значительные денежные затраты, которые не всегда возможны. Поэтому авторы работы [1] вспомнили об альтернативе, предложенной в 2005 году Лекебушем (Leckebusch), – о параллельном использовании нескольких пар одиночных георадиолокационных антенн. Такие модульные системы, конечно, несравнимы по скорости работы с многоканальными, позволяющими быстро осуществлять комбинированный сбор данных при высокой частоте дискретизации. Но с другой стороны, их можно постепенно расширять – и тогда затраты будут поэтапными. К тому же с модульными системами легче работать в ограниченных пространствах, на территориях с неровной поверхностью или при неподходящем буксирующем транспортном средстве, чем с крупными многоканальными многоэлементными антеннами.

В рассматриваемой статье [1] обсуждаются двумерные и трехмерные георадиолокационные исследования с полным (максимально высоким) разрешением (full-resolution survey) с использованием импульсной модульной георадарной системы, указываются требования к позиционированию, приводятся примеры фильтрации собранных данных разными методами для устранения полос шума и помех на горизонтальных (временных) срезах результатов трехмерной съемки и сравнивается их эффективность.

 

ВЫПОЛНЕНИЕ СЪЕМКИ

Вердонк и соавторы [1] выполнили пробные георадиолокационные измерения с помощью георадарной системы Spidar производства компании Sensors & Software с антеннами pulseEKKO PRO на 500 МГц при археологических исследованиях на территориях древнеримских городов Мариана (коммуна Луччана, регион Корсика, Франция) и Колония Ульпия Траяна (г. Ксантен, земля Северный Рейн-Вестфалия, Германия). Далее для краткости они называют вторую площадку Ксантеном.

Антенны для обеих съемок закреплялись на колесной раме, буксируемой квадроциклом. Корпус одиночный антенны на 500 МГц имел ширину примерно 0,23 м, расстояние между средними точками антенн составляло 0,25 м (рис. 1). При этом методики съемки на двух указанных территориях различались.

 

 

3D-съемка с полным разрешением в Мариане

В Мариане авторами работы [1] была проведена трехмерная съемка с полным разрешением с использованием трех пар антенн на 500 МГц на площади 85 м х 35 м в октябре 2010 года. Для определения местоположения использовался отслеживающий тахеометр (система наземного позиционирования).

Поскольку для выполнения 3D-съемки часто используется поперечный шаг между профилями не более четверти длины волны, для съемки в Мариане этот шаг должен был составить 0,05 м, учитывая среднюю подповерхностную скорость распространения радиоволны указанной частоты 0,1 м/нс. В результате поперечное (относительно линии движения при съемке) расстояние между антеннами в 5 раз превысило желаемый шаг между профилями. Поэтому, чтобы получить нужные 0,05 м, использовалась схема записи, показанная на рисунке 2, а: пять последовательных проходов выполнялись со сдвигом в поперечном направлении на 0,05 м, а после них следовали интервалы по 0,55 м.

 

 

Авторы статьи [1] напоминают, что для несмещенной трехмерной георадиолокационной съемки точность координат и точность записи должны быть лучше 1/8 длины волны, чтобы однозначно привязать каждое георадарное измерение к правильному положению на сетке с ячейками шириной в 1/4 длины волны. Для указанных выше частоты и скорости распространения радиоволны 1/8 ее длины равна 0,025 м. В соответствии с этим должна быть выполнена коррекция с учетом задержки по времени между измерением положения и его объединением с данными георадара и ошибок синхронизации между измерениями расстояния и угла. Кроме того, точность координат подразумевает, что во время съемки должно тщательно отслеживаться расстояние между соседними профилями в 1/4 длины волны. В случае системы с поперечным расстоянием между центрами антенн 0,20 м или больше, как это было использовано для съемок, рассматриваемых в статье [1], ее авторам пришлось пройти профили между ранее записанными, чтобы получить набор трехмерных данных с полным разрешением (см. рис. 2, а). В этом случае было очень важно как можно более точно следовать теоретически выбранным поперечным шагам. Для 3D-съемки в Мариане вдоль рулеток, протянутых поперек будущих профилей через каждые 40 м, с нужным шагом размещались колышки, которые временно удалялись для каждого прохода. Поскольку рельеф местности был относительно ровным, а съемка проводилась очень медленно (квадроцикл двигался со скоростью около 1 м/с, и два человека за 13 часов выполнили 234 прохода длиной по 85 м в одном направлении), можно было сохранять расположение колышков точно по теоретически выбранным линиям сетки и следовать по этим линиям с точностью до 0,03 м.

 

 

2D-съемка в Ксантене

В Ксантене авторы статьи [1] провели в июле 2011 года двумерную георадарную съемку с поперечным шагом между профилями 0,25 м на площади около 6600 кв. м. Они использовали шесть пар антенн на 500 МГц, закрепленных на раме, буксируемой квадроциклом. Поперечное расстояние между центрами антенн было равно желаемому поперечному шагу между профилями 0,25 м. Использованная схема сбора данных показана на рисунке 2, б.

Для позиционирования использовались направляющие тросы. Они были проложены вдоль линий движения на поперечном расстоянии друг от друга 6 м с использованием фиксированных топографически определенных точек. Между двумя соседними тросами можно было пройти четыре полосы сканирования шириной по 1,5 м. Линейные координаты измерялись с помощью колеса счетчика пройденного пути, запускающего георадарное зондирование через каждые 0,05 м, а координатные метки делались на записи через каждые 25 м. Для 2D-съемки, выполненной в зигзагообразном режиме, два человека сделали 47 проходов средней длиной по 95 м за 4 часа.

Как и в Мариане, георадиолокационные измерения в Ксантене были выполнены с использованием суммирования (накопления) в одну запись по 8 последовательных трасс (stack of 8) для улучшения соотношения сигнал/шум, временного интервала дискретизации (temporal sampling interval) 0,2 нс и временного окна (развертки сигнала, time window) 60 нс. При обоих исследованиях скорость движения составляла около 1 м/с, хотя с теми же настройками была возможна скорость до 2 м/с. Так что по сравнению со съемкой с одной парой антенн время записи в полевых условиях в целом было сокращено более или менее пропорционально количеству используемых пар антенн.

 

ОБРАБОТКА ДАННЫХ

Вердонк и соавторы [1] подвергли первичные результаты съемки этапам обработки, перечисленным на рисунке 3. После частотной фильтрации (с окном в 4 нс) было выполнено выравнивание начала отсчета времени (time-zero alignment) между разными антеннами. Для каждой радарограммы рассчитывалась средняя трасса. Для этих средних трасс были выровнены (aligned) отрицательные прогибы (отрицательные фазы зондирующего импульса, negative troughs) после начала воздушной (прямой) волны. Впоследствии была применена функция усиления, основанная на обратной средней огибающей амплитуд всех трасс результатов съемки, сглаженной с помощью фильтра скользящего среднего с длиной окна в 15 временных отсчетов (3 нс).

 

 

Этапы обработки результатов съемки в Мариане (см. рис. 3) включали в том числе привязку трасс к горизонтальным координатам системы наземного позиционирования, записанным центральной антенной, над которой была установлена призма электронного тахеометра (см. рис. 1). Поскольку этот тахеометр позволял получать координаты с частотой 5 Гц, данные о местоположении не были доступны для каждой трассы. Поэтому промежуточным трассам были присвоены координаты с помощью сплайн-интерполяции между двумя ближайшими результатами измерений тахеометра. Координаты двух боковых антенн x и y были рассчитаны с учетом направлений от центральной антенны на две ближайшие точки с уже известными координатами (с использованием угла между этими направлениями и теоретическим поперечным направлением, а также известного расстояния между серединами центральной и боковых антенн, равного 0,25 м). Затем эти данные были интерполированы на регулярную сетку с размером ячеек 0,05 м х 0,05 м с использованием триангуляции Делоне, включая линейную интерполяцию между амплитудами в углах окружающего треугольника. После ослабления полос шума и помех, которое обсуждается в следующих разделах, были введены статические поправки по высоте для исключения влияния рельефа. Чтобы получить плоскую поверхность для миграции методом фазового сдвига, к слегка неровной территории съемки была подобрана плоскость таким образом, чтобы минимизировать среднеквадратическую ошибку. Подповерхностная скорость радиоволны оценивалась путем применения 2D-миграции по алгоритму Столта к отдельным георадиолокационным профилям с использованием диапазона постоянных скоростей. Было обнаружено, что скорость уменьшалась с 0,11 м/нс для мелкозалегающих слоев до приблизительно 0,08 м/нс на глубине 1 м. Затем данные, полученные в Мариане, подверглись 3D-миграции методом фазового сдвига с использованием полученной функции скорости, которая также использовалась для преобразования временного разреза в глубинный.

Что касается съемки в Ксантене, то, как напоминают авторы статьи [1], она выполнялась вдоль направляющих тросов без системы наземного позиционирования. При этом для расширения или сужения промежутков между георадиолокационными трассами там, где это было необходимо, использовались координатные метки, чтобы скомпенсировать возможное несовершенство калибровки датчика пройденного пути (см. рис. 3). Для данных из Ксантена не потребовалось никакой топографической коррекции, так как территория съемки была практически плоской. Преобразование временного разреза в глубинный было основано на равномерной скорости радиоволны 0,075 м/нс, определенной с помощью тестов на миграцию с постоянной скоростью с использованием 2D-миграции по алгоритму Столта, что также использовалось для миграции отдельных профилей. Затем данные между 2D-профилями были проинтерполированы с помощью бикубического алгоритма для создания псевдотрехмерных горизонтальных срезов в точках сетки с размером ячеек 0,05 м х 0,05 м.

 

СИНФАЗНЫЕ ШУМ И ПОМЕХИ (LINEAR NOISE)

Вердонк с коллегами [1] столкнулись с тем, что на георадиолокационных профилях из Марианы и Ксантена были видны значительные неинформативные горизонтальные и субгоризонтальные полосы (рис. 4, а). И это не было связано с прямой волной, генерируемой передатчиком, принадлежащим другой паре антенн, поскольку в использованной при съемке системе Spidar только один передатчик может быть активен в заданном временном окне приемника, что устраняет интерференцию. Дело в том, что невозможно создать две антенны с абсолютно идентичными амплитудными и фазовыми спектрами, как бы ни старался производитель. Поэтому на разрезах, записанных с помощью разных антенн, картины шумов и помех всегда немного различаются, что вызывает появление полос вдоль линий съемки на горизонтальных срезах (рис. 5).

 

Рис. 4. Георадиолокационный профиль из Ксантена, записанный вдоль линии y = 39,5 м (ее местоположение показано стрелкой на рис. 15, a): а – необработанный профиль, на котором видна так называемая горизонтальная и субгоризонтальная полосчатость (синфазные шум и помехи); б – тот же профиль после удаления фона путем вычитания среднего по всем трассам; в – результат 1D дискретного вейвлет-преобразования строки 185 (на 37 нс) профиля «б», рассчитанный с использованием вейвлета Добеши 4 (db4), при этом «S» обозначает исходный сигнал, A8 – аппроксимацию, D1–D8 – уровни детализации (details); г – результат фильтрации профиля «б» путем 1D дискретного вейвлет-преобразования, обнуления коэффициентов, приуроченных к D7, D8 и A8 для каждой строки; д – результат τ-p-преобразования профиля «б» (представления профиля «б» в τ-p-области на основе преобразования Радона, где τ – время прохождения сигнала, p – величина, обратная кажущейся скорости) после удаления отображений воздушных и приземных (ground) волн-помех; е – отфильтрованный профиль, полученный следующим образом: из профиля «д» были сохранены только данные для p в области от минус 0,33 нс/м до плюс 0,33 нс/м (для τ>20 нс), а результат после обратного преобразования Радона был вычтен из соответствующей области профиля «б» для получения отфильтрованного профиля «е» [1]

 

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

 

Эти различные картины шумов и помех могут быть усилены за счет собственных шумов и помех от работающей аппаратуры (аппаратного «звона»), при проведении съемки на глинистых или влажных грунтах или при плохом либо изменяющемся в зависимости от места измерений контакте между антенной и грунтом. Поскольку антенна георадара касается земли, на ее входное  сопротивление (входной импеданс) сильно влияет грунт. Это может вызвать рассогласование между импедансом антенны и волновым сопротивлением питающего кабеля, в результате чего возникнут токи между фидером и концами антенны. Следовательно, передающая антенна будет генерировать несколько затухающих во времени импульсов вместо одного. Да и в приемной антенне может возникнуть такой остаточный отклик. Кроме того, во время рассматриваемых съемок Вердонк с коллегами [1] заметили, что рисунок горизонтальных полос изменился, когда немного переместили антенный кабель. Действительно, кабель вблизи антенны тоже может выступать в качестве источника излучения и генерировать или усиливать аппаратный «звон». При этом испытания авторов работы [1] показали, что в условиях без аппаратного «звона» (например, на песчаных грунтах) различия в амплитудах для разных антенн были менее значительными.

 

ИСПОЛЬЗОВАННЫЕ МЕТОДЫ ФИЛЬТРАЦИИ

Вердонк и соавторы [1] рассмотрели и опробовали несколько методов ослабления полос синфазных помех и шума. Вычитание среднего сигнала по всем трассам профиля из каждой отдельной трассы (удаление фона) только частично убрало эти помехи и шум с изображений из Марианы и Ксантена, потому что амплитуды таких полос часто сильно варьировали вдоль профилей и к тому же они иногда были не совсем горизонтальными. Это хорошо видно на рис. 4, б, а также на более глубоких горизонтальных срезах с низким соотношением «сигнал/шум» (см., например, рис. 6, a).

Авторы работы [1] напомнили, что Лекебуш (Leckebusch) в 2005 году при георадиолокационных исследованиях археологических площадок, чтобы удалить с временных срезов полосы, вызванные различиями между антеннами, комбинированно использовал «спектральное отбеливание» (полосовую фильтрацию и обратную фильтрацию, или деконволюцию) и медианную фильтрацию (energy matching – вычисление медиан энергии по каждому профилю для определенного момента времени или определенной глубины и замену медианными значениями аномальных точек и выбросов, то есть подавление некоррелированных и слабо коррелированных шумов и помех и малоразмерных деталей).

Есть и другие методы удаления полос от синфазных помех и шума, переотражений, «звона»:

• f-k-фильтрация (называемая также пространственно-временной фильтрацией или частотно-волновочисленной фильтрацией, когда на основе двумерного преобразования Фурье от временных функций переходят к циклической частоте f, а от координаты пространства х переходят к пространственной частоте (волновому числу) k, а затем выполняют фильтрацию в f-k-области);
• предсказывающая деконволюция (прогнозирующая обратная фильтрация);
• фильтрация собственных изображений с использованием сингулярной декомпозиции (разложения);
• дискретное вейвлет-преобразование;
• дискретное риджлет-преобразование;
• τ-p-фильтрация (фильтрация в τ-p-области на основе преобразования Радона, где τ – время прохождения сигнала, p – величина, обратная кажущейся скорости).

На рисунках 6–8 показаны результаты применения авторами работы [1] некоторых из этих методов к временным срезам трехмерных георадиолокационных записей из Марианы и Ксантена после удаления фона путем вычитания среднего по всем трассам.

 

Рис. 6. Сравнение результатов применения различных методов ослабления полос от синфазных помех и шума к энергетически-временному срезу (горизонтальному срезу примерно на уровне 27–28 нс) трехмерной георадиолокационной записи из Марианы: а – тот же срез, что и на рисунке 5, после удаления фона путем вычитания среднего по всем трассам (он послужил в качестве входных данных для последующей фильтрации); б – результат f-k-фильтрации с помощью простого фильтра высоких частот, как показано на рисунке 9, б (стрелка указывает на место, где не совсем ясен ход конструкции стены); в – изображение, полученное путем вычитания результата использования фильтра скользящего среднего в направлении x с применением окна в 100 отсчетов (5 м); г – результат f-k-фильтрации с использованием комбинации веерного фильтра высоких частот и круглого фильтра низких частот, как показано на рисунке 9, г

 

Рис. 7. Сравнение результатов применения различных методов ослабления полос от синфазных помех и шума к энергетически-временному срезу (горизонтальному срезу примерно на уровне 27–28 нс) трехмерной георадиолокационной записи из Марианы. Данные после удаления фона путем вычитания среднего по всем трассам (см. рис. 6, а) послужили в качестве входных для фильтрации с получением горизонтальных срезов (а, г), а также для фильтрации вертикальных профилей с последующим «извлечением» горизонтальных срезов (б, в): a – результат фильтрации с помощью 2D дискретного вейвлет-преобразования с использованием вейвлета Добеши 4 (db4), где коэффициенты детализации по горизонтали на уровнях разложения 1, 2, 3 (соответственно h1, h2 и h3) были обнулены; б – срез, «извлеченный» из вертикальных профилей, которые были подвергнуты τ-p-фильтрации с использованием параметров, перечисленных в таблице; в – cрез, «извлеченный» из вертикальных профилей, подвергнутых сингулярной декомпозиции, с удалением первых восьми собственных изображений; г – результат применения фильтра среднего по профилю (пространственного двумерного фильтра, сглаживания в окне – mean profile filter) с использованием фильтра низких частот с окном 15 х 1 в направлении x и фильтра высоких частот с окном 1 х 11 в направлении y. Примечание: фильтрация в пунктах б–г проводилась в комбинации с использованием круглого фильтра низких частот в волновочисленной области, как показано на рисунке 9, г

 

Рис. 8. Сравнение результатов применения различных методов ослабления полос от синфазных шума и помех к энергетически-временному срезу (горизонтальному срезу примерно на уровне 27–28 нс) фрагмента трехмерной георадиолокационной записи из Ксантена: a – временной срез после удаления фона путем вычитания среднего по всем трассам, который послужил в качестве входных данных для последующей фильтрации с получением горизонтальных срезов (б, в, д, е), а также для фильтрации вертикальных профилей с последующим «извлечением» горизонтального среза (г); б – результат f-k-фильтрации с применением простого фильтра высоких частот, подавляющего значения от минус 0,05 до плюс 0,05 м^-1, как показано на рисунке 9, б; в – фильтрация с помощью 2D дискретного вейвлет-преобразования с использованием вейвлета Добеши 4 (db4), где коэффициент детализации по горизонтали на уровне разложения 2 (то есть h2) был обнулен; г – cрез, «извлеченный» из вертикальных профилей, которые были подвергнуты τ-p-фильтрации с использованием параметров, перечисленных в таблице (см. также рис. 4, д, е); д – cрез, «извлеченный» из вертикальных профилей, подвергнутых сингулярной декомпозиции, с вычитанием первого собственного изображения; е – результат применения фильтра среднего по профилю (пространственного двумерного фильтра, сглаживания в окне – mean profile filter) вкупе с использованием фильтра низких частот с окном 63 х 1 в направлении x и фильтра высоких частот с окном 1 х 7 в направлении y

 

Таблица. Параметры для τ-p-фильтрации георадиолокационных данных, полученных в Мариане и Ксантене

 

3D частотно-волновочисленная фильтрация

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

При обработке сейсмических данных для снижения избыточной амплитудной разрешенности до настоящего времени широко использовалась и используется фильтрация в горизонтальной волновочисленной области k_x-k_y. Но это может затруднить интерпретацию данных и расчет сейсмических атрибутов. Шум и помехи, создаваемые геометрией сейсмической съемки, являются очень повторяющимися, поэтому в волновочисленной области они проявляются в виде четких всплесков. С другой стороны, схемы сбора данных по сейсмическим полям часто бывают сложными, с неортогональными линиями источника и приемника. Следовательно, большинство решений для ослабления избыточной амплитудной разрешенности в сейсморазведке приходится адаптировать к конкретным данным.

В случае же георадарных данных, полученных Вердонком с коллегами [1], напротив, картины шумов и помех, которые возникли из-за неодинаковых спектров антенн, были менее регулярны и не могли быть легко автоматически локализованы в волновочисленной области. Поэтому авторы работы [1] тот или иной фильтр предназначали для каждого конкретного случая. Для полученных ими георадиолокационных записей наиболее эффективной оказалась фильтрация в волновочисленной области, выполненная на частотных срезах. Линейные полосы в направлении x (см. рис. 6, a; 8, a) в горизонтальной волновочисленной области идут вдоль оси k_x = 0, как это видно по спектрам амплитуд на рисунке 9, а, б на записях из Марианы. Фильтрация высоких частот, подавляющая значения k_x от минус 0,05 до плюс 0,05 м^-1, выполнялась (определялась) визуально (рис. 9, б). Область около центра с радиусом 0,75 м^-1 была исключена из любой фильтрации. После фильтрации осталось несколько неопределенностей (например, в месте, на которое указывает стрелка на рисунке 6, б). Похожий эффект был получен при вычитании результата использования фильтра скользящего среднего в направлении x с окном в 100 отсчетов (что соответствует 5 м, см. рис. 6, в).

 

Рис. 9. Полосовые фильтры в частотно-волновочисленной области: а – горизонтальный волновочисленный спектр всего набора данных из Марианы для частоты 371 МГц (после удаления фона путем вычитания среднего по всем трассам в пространственно-временной области); б – простой вертикальный фильтр высоких частот для ослабления шума и помех вдоль оси k_x = 0; в – данные с частотой 283 МГц после удаления фона путем вычитания среднего по всем трассам показывают более широкую полосу шума вдоль оси k_x = 0; г – комбинация веерного фильтра и круглого фильтра низких частот для ослабления полос

 

Уточнить результаты для Марианы Вердонку с соавторами [1] позволили две модификации. Для частот ниже 350 МГц видна более широкая полоса шума и помех вдоль оси k_x = 0 (рис. 9, в). Удаление всей этой полосы привело бы к повреждению записей отражений от участков стен с ориентацией, близкой к направлению линий съемки. Поэтому вместо этого использовали фильтр, расширяющийся к концам оси k_x = 0 (а также для уменьшения частот) (рис. 9, г). Эта фильтрация была объединена с применением круглого фильтра низких частот с радиусом 6 м^-1 (см. рис. 9, г). После фильтрации стала видна протяженность едва различимых археологических сооружений (см. рис. 6, г).

Что касается данных из Ксантена, то использование простого направленного фильтра (см. рис. 9, б) не помогло удалить артефактные полосы с горизонтального среза (см. рис. 8, б). Результаты не улучшились и при увеличении ширины фильтра или применении его в виде веерного фильтра (см. рис. 9, г). Ответ на причины этого дало более внимательное изучение данных в пространственно-временной области. Поскольку в Ксантене поверхность земли была более неровной, чем в Мариане, контакт с ней антенн изменялся вдоль профилей. Поэтому после удаления фона путем вычитания среднего по всем трассам шум и помехи на профилях проявляются в виде сложных меняющихся по горизонтали картин, причем в основном между 20 и 40 нс (см. рис. 4, б). На срезах (см. рис. 8, а) это видно в виде серии точек, а не полос. В отличие от линейных структур округлые аномалии отображаются по всей плоскости Фурье (во всей частотной области).

Второй фактор, который усложнил f-k-фильтрацию данных для Ксантена, – это частота дискретизации в поперечном направлении, которая была в пять раз ниже, чем для Марианы. Следовательно, использование круглого фильтра низких частот, который для записей из Марианы способствовал эффективности f-k-фильтрации, не было подходящим для Ксантена, а наоборот, вызвало сильное ухудшение четкости в направлении у и потерю возможности выявления археологических деталей.

 

Дискретное вейвлет-преобразование и дискретное риджлет-преобразование

Авторы работы [1] напоминают, что при вейвлет-преобразовании (которое бывает непрерывным или дискретным) сигнал раскладывается по базисным функциям (вейвлетам), полученным из некоторой исходной базисной функции (материнского вейвлета) путем изменений масштаба (сжатий, растяжений) и сдвигов. Эту декомпозицию, которая приводит к аппроксимации (низкочастотным компонентам) и детализации (высокочастотным компонентам), можно итерировать, взяв результаты аппроксимации в качестве входных данных для еще одной фильтрации – многоуровневого вейвлет-преобразования.

Примерами использования вейвлет-преобразования в георадиолокационных исследованиях являются: выделение признаков при распознавании образов (feature extraction); объединение данных, собранных с помощью двух перпендикулярно ориентированных антенн; ослабление когерентных шума и помех на вертикальных георадарных профилях.

Для двумерного четырехуровневого дискретного вейвлет-преобразования данных из Марианы и Ксантена использовался вейвлет Добеши 4 (db4).

Для Марианы при получении вейвлет-реконструкций (обратных вейвлет-преобразований) H1, H2 и H3 для амплитудно-временного среза на 27,6 нс (рис. 10, а) использовали горизонтальные детализирующие коэффициенты h1, h2 и h3 только на уровнях разложения 1, 2 и 3 соответственно (рис. 10, б–г). Обнуление коэффициентов h1, h2 и h3 (то есть вычитание H1, H2 и H3 из исходного изображения) подавило полосы шума и помех (см. рис. 7, a) и дало результат, аналогичный результату комбинированного использования веерного f-k-фильтра и круглого фильтра низких частот (см. рис. 6, г).

 

Рис. 10. Двумерное четырехуровневое дискретное вейвлет-преобразование фрагмента записи из Марианы с использованием вейвлета Добеши 4 (db4): а – энергетически-временной срез на 27,6 нс после удаления фона путем вычитания среднего по всем трассам; б–г – вейвлет-реконструкции H1, H2 и H3, основанные на горизонтальных детализирующих коэффициентах h1, h2 и h3 соответственно (в этих реконструкциях концентрируется «полосчатость» в направлении x при уменьшении частоты в направлении y)

 

Аналогичная процедура была применена и к временным срезам данных из Ксантена, но обнуление было выполнено только для горизонтального детализирующего коэффициента h2. Подобно тому, что наблюдалось для f-k-фильтрации, из-за грубой поперечной дискретизации произошла некоторая потеря четкости для подповерхностных сооружений в направлении y (см. 8, в). Не лучше были и результаты при использовании других вейвлетов из семейства Добеши.

Хотя 2D дискретное вейвлет-преобразование дало хорошие результаты по ослаблению «полосчатости» на временных срезах данных из Марианы, обнуление целых горизонтальных полос (panels) (рис. 10, б–г) может привести к артефактам на вертикальных профилях. В связи с этим авторы статьи [1] напоминают, что другими исследователями ранее было разработано несколько методов восстановления полезного сигнала из зашумленных вейвлет-коэффициентов с помощью выбора порога. Широко используемая стратегия – это «мягкое» задание пороговой функции в пространстве вейвлет-коэффициентов на основе разброса (дисперсии) шумов и помех, минимизации несмещенной оценки риска Штейна или минимизации байесовского риска. Однако эти методы имеют дело только с белым гауссовским шумом. Более того, было показано, что вейвлеты неэффективно представляют анизотропные элементы изображения, такие как линии и кривые. Их можно лучше описать с помощью избирательно ориентированных преобразований, таких как риджлет- и курвлет-преобразования.

Далее Вердонк и коллеги [1] напоминают, что базовая функция (риджлет) риджлет-преобразования постоянна вдоль линий и представляет собой вейвлет, поперечный гребням колебаний (ridges). Непрерывное риджлет-преобразование может быть выполнено путем выполнения преобразования Радона R(θ,t) (где θ – угол линий проекций; t – перпендикулярное расстояние от линии проекции до начала координат) и применения одномерного вейвлет-преобразования к результату преобразования Радона по t.

На временных срезах данных из Марианы авторами работы [1] было проведено цифровое риджлет-преобразование (рис. 11), основанное на выполнении преобразования Радона, описанного в работе Авербуха и др. (Averbuch et al.) 2008 года. С помощью вейвлета Добеши 4 (db4) было выполнено одномерное восьмиуровневое дискретное вейвлет-преобразование, а уровни детализации (details) D1–D4 были обнулены для тех столбцов, которые представляют углы от минус 2,5 до плюс 7,5 град. (где 0 град. представляет направление x). Эта процедура эффективно удалила полосы и дала примерно такое же изображение, как и при большинстве других применяемых методов ослабления «полосчатости», хотя отношение «сигнал/шум» в описанном случае, кажется, несколько ниже, а элемент изображения, на который указывает стрелка на рисунке 11, является менее четким, чем при использовании других методов фильтрации.

 

Рис. 11. Энергетически-временной срез данных из Марианы на уровне около 27–28 нс после фильтрации с помощью цифрового риджлет-преобразования, состоящей из следующих этапов. Сначала к временным срезам применили дискретное преобразование Радона. На полученных изображениях (plots) столбцы, соответствующие углам от минус 2,5 до плюс 7,5 град., были отфильтрованы путем одномерного дискретного вейвлет-преобразования с использованием вейвлета Добеши 4 (db4) и обнуления уровней детализации (details) D1–D4. И наконец, было выполнено обратное преобразование Радона. Элемент изображения, на который указывает стрелка, получился несколько менее четким, чем на срезах, полученных при фильтрации другими методами

 

На вертикальных профилях из Ксантена авторы статьи [1] также исследовали эффективность 1D дискретного вейвлет-преобразования. На рисунке 4, в показан результат такого преобразования одной из трасс профиля, показанного на рисунке 4, б, для 37 нс. Из исходного сигнала S результат аппроксимации A8 и уровни детализации (details) D7 и D8 были вычтены путем обнуления их соответствующих коэффициентов. То же самое было сделано для каждой трассы георадарного профиля. Как показали ранее, например, Нуццо и Куарта (Nuzzo, Quarta, 2004), эта процедура может эффективно удалить горизонтальные полосы с записанных профилей (см. рис. 4, г). Однако для данных из Ксантена после «извлечения» временных срезов из вертикальных профилей Вердонку и коллегам [1] стало ясно, что «негоризонтальные остатки» шума и помех, которые различаются по амплитуде на разных поперечных разрезах (transects), делают невозможным удаление полос также и с горизонтальных срезов.

 

Линейное преобразование Радона

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

Подобно преобразованию Радона, используемому при обработке изображений и рассмотренному выше при обсуждении дискретного риджлет-преобразования, τ-p-преобразование (линейное преобразование Радона, или преобразование наклонного суммирования) при обработке сейсмических данных определяется как суммирование амплитуд в t-x-области вдоль или по касательной к линии с наклоном p, пересекающей ось времени в момент времени τ. Суммированные данные затем отображаются в точке S в τ-p-области:

 

где t – время, x – смещение, p – параметр луча (ray parameter).

 

Линейное преобразование Радона в τ-p-области часто применяется к сейсмограммам, полученным методом общей глубинной точки, но может быть расширено до сейсмических профилей с нулевым удалением источника от приемника или до георадарных профилей. При обработке сейсмических данных – это ослабление кратных отражений и анализ скоростей. Преобразование τ-p выполняется как форма дискретного преобразования Радона посредством обращения методом наименьших квадратов, выполняемого в частотной области.

При выполнении τ-p-преобразований данных из Марианы и Ксантена Вердонк и соавторы [1] сохранили интервал дискретизации для параметра луча (Δp) меньше 1/(x_rf_max), где x_r – длина георадарных профилей, f_max – максимальная частота. Максимальный диапазон для параметра луча (p_max) был выбран так, чтобы он не превышал 1/(2Δxf_max), где Δx – шаг сканирования (интервал дискретизации) вдоль линии профиля (см. таблицу). Наибольшую часть полезных сигналов содержала полоса частот 0–1 ГГц, поэтому именно она и была выбрана для преобразования. Дискретизация в направлении τ была такой же, как и для t. Был использован коэффициент затухания (damping factor), равный 0,001.

На рисунке 4, д показан георадиолокационный профиль, полученный авторами статьи [1] в результате τ-p-преобразования профиля из Ксантена, представленного на рисунке 4, б. Почти горизонтальные полосы на рисунке 4, б соответствуют узкой области вокруг p = 0 нс/м. Из τ-p-сечений были сохранены только данные в небольшом прямоугольном окне (см. таблицу) со значениями τ больше 20 нс, содержащие синфазные помехи и шум (linear noise). После обратного преобразования эти сечения были вычтены из исходных профилей. Итог (см. рис. 4, е) был аналогичен результату 1D вейвлет-преобразования (см. рис. 4, г), но и в этом случае остаточные негоризонтальные компоненты шума и помех не позволили полностью подавить полосы на временных срезах (см. рис. 8, г).

Профили из Марианы авторы работы [1] подвергли τ-p-фильтрации аналогичным образом. После создания горизонтальных срезов оптимальное ослабление полос от синфазных помех и шума (см. рис. 7, б) было достигнуто только после низкочастотной фильтрации данных в волновочисленной области с применением круглого фильтра, уже использованного в комбинации с веерным f-k-фильтром.

 

Сингулярная декомпозиция

Вердонк и соавторы [1] напоминают, что сингулярная декомпозиция (разложение по сингулярным, или особым, значениям) – это факторизация входной матрицы A:

 

где столбцы U – собственные векторы (eigenvectors) матрицы AA^T; столбцы V – собственные векторы матрицы A^TA;  Σ – прямоугольная диагональная матрица с сингулярными значениями матрицы A, расположенными в порядке убывания.

 

Это также можно записать следующим образом:

 

 

Здесь вклад собственного изображения (eigenimage) u_iv_i^T в матрицу A пропорционален сингулярному значению σ_i. Как было показано другими авторами ранее, первые собственные изображения содержат наиболее высококоррелированные части набора данных и их вычитание может быть мощным методом. При обработке георадарных записей он, например, использовался разными исследователями для оценки количества отражений волн, для отделения воздушных (air waves) и приземных (ground waves) волн от отражений, вызванных приповерхностными границами раздела, для ослабления аппаратного «звона» (собственных шума и помех от аппаратуры) на вертикальных георадарных профилях.

Сингулярная декомпозиция тесно связана с анализом с помощью метода главных компонент и с преобразованием Карунена – Лоэва. Применение метода главных компонент в георадиолокационных исследованиях до сих пор заключалось в удалении разными исследователями отражений от поверхности земли (ground-bounce) с записей, полученных с помощью георадара непрерывного излучения со ступенчато изменяемой частотой (SFCW GPR) воздушного базирования.

Авторы работы [1] выполнили сингулярную декомпозицию как на профилях, так и на горизонтальных срезах. Ее применение к временным срезам из Марианы (см. рис. 6, а) и вычитание первого собственного изображения ослабили полосы, хотя некоторые помехи и шум все еще присутствовали (рис. 12, a). И этот результат заметно не улучшился, когда было вычтено больше собственных изображений. Более того, на некоторых участках съемки стали видны линейные артефакты (рис. 12, б, № 1), вызванные высокоамплитудными аномалиями (например, на которую показывает стрелка № 2 на рисунке 12, б). На рисунках 12 в, г показано, как эта высокоамплитудная особенность отражается на левом и правом сингулярных векторах.

 

Рис. 12. Результаты обработки георадиолокационных данных из Марианы: а – энергетически-временной срез примерно на уровне 27–28 нс после сингулярной декомпозиции и обратного преобразования (реконструкции) без первого собственного изображения (видно, что на этом срезе все же присутствуют некоторые полосы); б – амплитудно-временной срез из Марианы на уровне 27,4 нс после реконструкции без первых трех собственных изображений (видно, что на этом срезе высокоамплитудная особенность, на которую указывет срелка № 2, отражается на левых сингулярных векторах (см. «в») и в большей степени на правых сингулярных векторах (см. «г»), в результате после вычитания собственных изображений на временном срезе появился линейный артефакт в направлении y (см. «б», № 1)); в – второй левый сингулярный вектор; г – второй правый сингулярный вектор

 

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

Авторы статьи [1] также выполнили сингулярную декомпозицию на данных из Ксантена. Как показано на рисунке 8, д, вычитание первого собственного изображения из горизонтальных срезов эффективно удалило шум и помехи, хотя в некоторых областях этот фильтр также частично подавил и полезные в археологическом отношении сигналы (рис. 13).

 

 

Использование фильтра среднего по профилю (пространственного двумерного фильтра, сглаживания в окне - mean profile filter)

Для записей из Ксантена авторы работы [1] получили относительно хорошие результаты с помощью фильтра среднего по профилю. Сначала вдоль линий съемки (survey lines) был использован одномерный фильтр низких частот (рис. 14, б), затем в поперечном направлении - одномерный фильтр высоких частот (рис. 14, в). Наконец, изображение, содержащее только шум и помехи (isolated noise), вычиталось из исходной картины (рис. 14, г). Такая фильтрация требовала применения различных коэффициентов и оценки результата, но адекватно отделяла шум и помехи от желаемых характеристик, как показано на рисунке 8, е.

 

Рис. 14. Использование фильтра среднего по профилю (пространственного двумерного фильтра, сглаживания в окне – mean profile filter) для фрагмента записи из Ксантена: а – изображение, полученное на частоте 410 МГц, после удаления фона путем вычитания среднего по всем трассам во временной области; б – результат применения фильтра низких частот с окном 63 х 1 (в направлении x); в – результат использования фильтра высоких частот с окном 1 х 7 (в направлении y); г – результат вычитания шума и помех (см. «в») из исходных данных (см. «а») только для характерных линий, изначально затронутых «полосчатостью»

 

Применение фильтра среднего по профилю для записи из Марианы, опять же в сочетании с круглым фильтром низких частот (см. рис. 9, г), дало итог, очень похожий на результаты использования других методов, испытанных Вердонком и соавторами [1] (см. рис. 7, г).

 

ОБСУЖДЕНИЕ РЕЗУЛЬТАТОВ И ИНТЕРПРЕТАЦИЯ ГЕОРАДАРНЫХ ДАННЫХ В АРХЕОЛОГИЧЕСКОМ ОТНОШЕНИИ

Итак, на временных срезах георадиолокационных данных, полученных авторами работы [1] в Мариане и Ксантене, даже после удаления фона путем вычитания среднего по всем трассам оставались неинформативные полосы шума и помех, которые к тому же отличались от «идеальных» артефактных полос с постоянными амплитудой и шириной на большей части направления x, в особенности для данных из Ксантена.

Дискретное вейвлет-преобразование и сингулярная декомпозиция, выполненные Вердонком и коллегами [1] при обработке данных, заняли меньше времени, чем f-k-фильтрация, использование фильтра среднего по профилю и τ-p-преобразование, где необходимо было выбирать соответствующие окна фильтров (filter windows). С одной стороны, медленная обработка менее привлекательна при больших наборах георадарных данных. Но с другой стороны, более быстрая сингулярная декомпозиция показалась авторам работы [1] наиболее склонной к удалению полезных аномалий и возникновению артефактов, особенно когда данные содержали высокоамплитудные элементы.

Однако в целом большинство испытанных авторами статьи [1] методов дали сопоставимые результаты: ни один из них не помог полностью удалить полосы шума и помех с временных срезов, даже если горизонтальные полосы были эффективно устранены с профилей. Исключением явилось 2D дискретное вейвлет-преобразование, хотя обнуление целых горизонтальных полос (panels) для фильтрации временных срезов может привести к появлению артефактов на профилях. И здесь полезной альтернативой, с точки зрения Вердонка и коллег [1], может быть дискретное риджлет-преобразование, поскольку оно более эффективно описывает линейные элементы изображения.

При этом окончательное качество фильтрации в значительной степени определялось методиками полевых измерений, выполнявшихся авторами работы [1].

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

Что касается записей из Ксантена, то это было невозможно, поскольку применение фильтра низких частот ухудшило разрешение в направлении y, учитывая более грубую дискретизацию в поперечном профилированию направлении. Хорошую основу для подходящего решения здесь обеспечило использование фильтра среднего по профилю (пространственного двумерного фильтра, сглаживания в окне – mean profile filter) с дополнительным уточнением результатов путем вычитания изолированных шума и помех (изображения, содержащего только шум и помехи, – isolated noise) только из линий, затронутых полосами (см. рис. 14, г; рис. 15, б). А применение других методов фильтрации только к зашумленным линиям не привело к заметному улучшению результатов.

 

 

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

На рисунке 16, а (для съемки в Мариане) показан результат моделирования шага между профилями 0,25 м путем уменьшения частоты дискретизации (decimation) данных, изначально полученных с полным разрешением (после f-k-фильтрации, как на рисунке 6, г). После этого были выполнены 2D-миграция и интерполяция на сетку 0,05 м х 0,05 м.

В отличие от смоделированной картины, изображение с полным разрешением (рис. 16, б) показало более резко очерченные края стен и более четкие дифракционные картины, вызванные отдельными камнями. В частности, по нему было легче интерпретировать особенности, параллельные или почти параллельные направлению профилирования. Например, элемент, на который указывает стрелка на рисунке 16, б, не может быть легко интерпретирован при пониженной частоте дискретизации (см. рис. 16, а), а на изображении с шагом между профилями 0,05 м можно выделить две отдельные параллельные аномалии – и эту картину, по мнению Вердонка с соавторами [1], можно было бы интерпретировать как водовод. Для примера на рисунке 16, в показан обнаруженный при раскопках в Мариане водовод с вертикально поставленными кирпичами на расстоянии около 0,2 м друг от друга в канаве с теми же размерами, что и упомянутая выше георадиолокационная аномалия.

 

 

Результат всей 3D-съемки, выполненной в Мариане, после удаления полос путем f-k-фильтрации (см. рис. 6, г; 9, г) и 3D-миграции показан на рисунке 17. Видны две параллельные улицы (рис. 17, № 1, 2) с той же ориентацией, что и в остальной части города. Пространство между ними в основном занимает большой дом, состоящий из ряда комнат вокруг внутреннего дворика (рис. 17, № 3). Одну из построек в западной части исследуемой территории можно с равным успехом интерпретировать как дом. Под входом со стороны улицы (№ 5) проходит дренажная или сточная канава. В некоторых местах есть признаки более поздних изменений (например, № 6-8).

 

 

Хотя в Ксантене не проводилась съемка с полным разрешением, полученные авторами работы [1] результаты все-таки добавили археологических знаний по территории исследования. Ее северо-западные, северо-восточные и юго-восточные границы на рисунке 15 примерно совпадают с краями улиц, окружающих квартал 31 в северной части древнеримского города. Этот блок зданий выполнял в основном жилую функцию, и некоторые из границ между прямоугольными участками видны на георадиолокационных изображениях (см. рис. 15, б, № 1-3). Часть границ расположена не под прямым углом к улицам (например, № 1; 3), что является напоминанием о более ранних этапах истории города. Также видны несколько построек 20 века (см. рис. 15, б, № 4).

 

ЗАКЛЮЧЕНИЕ

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

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

В статье [1] Вердонк и соавторы как раз и представили результаты 2D- и 3D-съемки с высоким разрешением с использованием модульной антенной системы при археологических исследованиях.

При выполнении 3D-съемки поперечное расстояние между средними точками соседних антенн было больше, чем требуемый шаг между профилями, поэтому авторам работы [1] пришлось выполнить дополнительные проходы для получения георадиолокационных записей между ранее полученными профилями. И здесь решающее значение имела точность позиционирования и записи. Оператор при профилировании должен был достаточно точно следовать вдоль теоретически выбранных линий сетки. Для съемок с использованием антенн с более высокой частотой (не менее 500 МГц) все же было сложно удерживать систему, буксируемую транспортным средством, строго на нужной линии сетки, чтобы соблюдался критерий выбранного шага дискретизации в четверть длины волны.

Даже если антенны были тщательно настроены производителем, на горизонтальных срезах профилей, полученных с помощью системы пар одиночных антенн, могут иметься неинформативные полосы от синфазных помех и шума в направлении линий съемки, что и произошло при работе Вердонка и коллег [1]. Они применили при обработке 2D и 3D георадарных записей несколько методов удаления этой «полосчатости», переотражений и «звона». Большинство испытанных фильтров оказалось достаточно эффективным в отношении ослабления или удаления горизонтальных полос, но на временных срезах остались проявления «негоризонтальных остатков» шума и помех с различными амплитудами в разных поперечных сечениях (transects) – остаточные линейные артефакты.

Оптимальное ослабление шума и помех на срезах оказалось возможным в основном при сочетании методов удаления полос с использованием фильтра низких частот в поперечном направлении. В 3D-данных это не вызвало заметного искажения полезной информации, но в 2D-данных это вызвало потерю деталей в археологическом отношении и поэтому не применялось. Для 2D-данных наиболее успешным оказалось использование фильтра среднего по профилю (пространственного двумерного фильтра, сглаживания в окне – mean profile filter), особенно когда изолированные шум и помехи (изображение, содержащее только шум и помехи, – isolated noise) вычитались только из линий, затронутых полосами.

---

Источник

 

1. Verdonck L., Vermeulen F., Docter R., Meyer C., Kniess R. 2D and 3D ground-penetrating radar surveys with a modular system: Data processing strategies and results from archaeological field tests // Near Surface Geophysics. European Association of Geoscientists & Engineers, 2013. Vol. 11 (2). P. 239-252. URL: researchgate.net/publication/261179966.

 

Список литературы, использованной авторами статьи [1]

 

Annan A.P. 2009. Electromagnetic principles of ground penetrating radar. In: Ground Penetrating Radar Theory and Applications (ed. H.M. Jol). Elsevier, Amsterdam-Oxford. P. 3-40.

Annan A.P., Davis J.L. 1992. Design and development of a digital ground penetrating radar system. In: Ground Penetrating Radar (ed. J.A. Pilon). Geological Survey of Canada. Paper 90-4. P. 15-23

Averbuch A., Coifman R.R., Donoho D.L., Israeli M., Shkolnisky Y., Sedelnikov I. 2008. A framework for discrete integral transformations II - The 2D discrete Radon transform. SIAM Journal on Scientific Computing 30, 785-803.

Beylkin G. 1987. Discrete Radon Transform. IEEE Transactions on Acoustics, Speech, and Signal Processing 35, 162-172.

Boeniger U., Tronicke J. 2010. On the potential of kinematic GPR surveying using a self-tracking total station: Evaluating system crosstalk and latency. IEEE Transactions on Geoscience and Remote Sensing 48, 3792-3798.

Cagnoli B., Ulrych T.J. 2001. Singular value decomposition and wavy reflections in ground-penetrating radar images of base surge deposits. Journal of Applied Geophysics 48, 175-182.

Chang S.G., Yu B., Vetterli M. 2000. Adaptive wavelet thresholding for image denoising and compression. IEEE Transactions on Image Processing 9, 1532-1546.

Di Renzo F. 2005. Nef sud: Etude des salles B-C-E. In: Rapport interm?diaire de fouilles programmees sur le site de Mariana (ed. P. Pergola). Service Regional de l'archeologie, Ajaccio. P. 36-58

Donoho D.L., Johnstone I.M. 1995. Adapting to unknown smoothness via wavelet shrinkage. Journal of the American Statistical Association 90, 1200-1224.

Ernenwein E.G., Kvamme K.L. 2008. Data processing issues in large-area GPR surveys: Correcting trace misalignments, edge discontinuities and striping. Archaeological Prospection 15, 133-149.

Fadili M.J., Starck J.L. 2009. Curvelets and ridgelets. In: Encyclopedia of Complexity and Systems Science (ed. R.A. Meyers). Springer, New York. P. 1718-1738.

Francese R.G., Finzi E., Morelli G. 2009. 3-D high-resolution multi-channel radar investigation of a Roman village in Northern Italy. Journal of Applied Geophysics 67, 44-51.

Freire S.L.M., Ulrych T.J. 1988. Application of singular value decomposition to vertical seismic profiling. Geophysics 53, 778-785.

Gerbrands J.J. 1981. On the relationships between SVD, KLT and PCA. Pattern Recognition 14, 375-381.

Gerlitz K., Knoll M.D., Cross G.M., Luzitano R.D., Knight R. 1993. Processing ground penetrating radar data to improve resolution of near-surface targets. Proceedings of the Sixth Symposium on the Application of Geophysics to Engineering and Environmental Problems. P. 561-574.

Grasmueck M., Viggiano D. 2007. Integration of ground-penetrating radar and laser positioning sensors for real-time 3-D data fusion. IEEE Transactions on Geoscience and Remote Sensing 45, 130-137.

Gulunay N., Benjamin N., Magesan M. 2006. Acquisition footprint suppression on 3D land surveys. First Break 24(2), 71-77.

Hill S., Shultz M., Brewer J. 1999. Acquisition footprint and fold-of-stack plots. The Leading Edge 18(6), 686-695.

Hugenschmidt J., Kalogeropoulos A. 2009. The inspection of retaining walls using GPR. Journal of Applied Geophysics 67, 335-344.

Karagul A., Crawford R., Sinden J., Ali S. 2004. Recent advances in 3D land processing: Examples from the Pakistan Badin area. First Break 22(9), 37-40.

Kim J., Cho S., Yi M. 2007. Removal of ringing noise in GPR data by signal processing. Geosciences Journal 11, 75-81.

Leckebusch J. 2005. Use of antenna arrays for GPR surveying in archaeology. Near Surface Geophysics 3(2), 111-115. doi: 10.3997/1873-0604.2005006

Leckebusch J. 2011. Comparison of a stepped-frequency continuous wave and a pulsed GPR system. Archaeological Prospection 18, 15-25.

Lines L.R., Treitel S. 1984. Tutorial: A review of least-squares inversion and its application to geophysical problems. Geophysical Prospecting 32, 159-186.

Linford N., Linford P., Martin L., Payne A. 2010. Stepped frequency ground-penetrating radar survey with a multi-element array antenna: Results from field application on archaeological sites. Archaeological Prospection 17, 187-198.

Linford N., Linford P., Payne A., David A., Martin L. 2011. Stonehenge: Recent results from a ground-penetrating radar survey of the monument. In: 9th International Conference on Archaeological Prospection, Izmir 19-24 September 2011 (eds M.G. Drahor and M.A. Berge). Archaeology and Art Publications, Istanbul. P. 86-89.

Marfurt K.J., Scheet R.M., Sharp J.A., Harper M.G. 1998. Suppression of the acquisition footprint for seismic sequence attribute mapping. Geophysics 63, 1024-1035.

Misiti M., Misiti Y., Oppenheim G., Poggi J. 2009. Matlab Wavelet Toolbox 4 User's Guide. The MathWorks Inc., Natick.

Muller M. 2008. Die stadtebauliche Entwicklung von der Coloniagr?ndung bis zur Spatantike. In: Colonia Ulpia Traiana. Xanten und sein Umland in romischer Zeit (eds M. M?ller, H-J. Schalles and N. Zieling). Verlag Philipp Von Zabern, Mainz am Rhein. P. 269-275.

Novo A., Lorenzo H., Rial F.I., Solla M. 2010. Three-dimensional ground-penetrating radar strategies over an indoor archaeological site: Convent of Santo Domingo (Lugo, Spain). Archaeological Prospection 17, 213-222.

Nuzzo L. 2003. Coherent noise attenuation in GPR data by linear and parabolic Radon transform techniques. Annals of Geophysics 46, 533-547.

Nuzzo L. 2005. Identification and removal of above-ground spurious signals in GPR archaeological prospecting. Archaeological Prospection 12, 93-103.

Nuzzo L., Quarta T. 2004. Improvement in GPR coherent noise attenuation using ?-p and wavelet transforms. Geophysics 69, 789-802.

Oimoen M.J. 2000. An effective filter for removal of production artefacts in U.S. Geological Survey 7.5-minute digital elevation models. 14th International Conference on Applied Geologic Remote Sensing, Las Vegas 6-8 November 2000, pp. 311-319.

Perrin S., Duflos E., Vanheeghe P., Bibaut A. 2004. Multisensor fusion in the frame of evidence theory for landmines detection. IEEE Transactions on Systems, Man and Cybernetics - Part C: Applications and Reviews 34, 485-498.

Radzevicius S.J., Guy E.D., Daniels J.J. 2000. Pitfalls in GPR data interpretation: Differentiating stratigraphy and buried objects from periodic antenna and target effects. Geophysical Research Letters 27, 3393-3396.

Sala J., Linford N. 2012. Processing stepped frequency continuous wave GPR systems to obtain maximum value from archaeological data sets. Near Surface Geophysics 10(1), 3-10. doi: 10.3997/1873-0604.2011046

Seren S., Loecker K., Hinterleitner A., Neubauer W., Ladst?tter S. 2011. Ephesos revisited: A decade of geophysical prospection on different field conditions in Ephesos, Turkey. In: 9th International Conference on Archaeological Prospection, Izmir 19-24 September 2011 (eds M.G. Drahor and M.A. Berge). Archaeology and Art Publications, Istanbul. P. 90-92.

Sensors & Software. 2010. Spidar user's guide. Sensors & Software Inc., Mississauga.

Soubaras R. 2002. Attenuation of acquisition footprint for non-orthogonal 3D geometries. EAGE 64th Conference & Exhibition, Florence, 27-30 May 2002, C-09.

Stempfhuber W., Schnaedelbach K., Maurer W. 2000. Genaue Positionierung von bewegten Objekten met zielverfolgenden Tachymetern. 13th International Course on Engineering Surveying. P. 144-153.

Strang G. 2006. Linear Algebra and its Applications. 4th Edition. Thomson Brooks/Cole, Pacific Grove.

Streich R., van der Kruk J., Green A.G. 2006. Three-dimensional multicomponent georadar imaging of sedimentary structures. Near Surface Geophysics 4(1), 39-48. doi: 10.3997/1873-0604.2005030

Tjora S., Eide E., Lundheim L. 2004. Evaluation of methods for ground bounce removal in GPR utility mapping. In: Proceedings of the Tenth International Conference on Ground Penetrating Radar, Delft, 21-24 June 2004 (eds E.C. Slob, A.G. Yarovoy and K.B. Rhebergen). Delft University of Technology. P. 379-382.

Trinks I., Johansson B., Gustafsson J., Emilsson J., Friborg J., Gustafsson C. et al. 2010. Efficient, large-scale archaeological prospection using a true three-dimensional ground-penetrating radar array system. Archaeological Prospection 17, 175-186.

Turner G. 1990. Aliasing in the tau-p transform and the removal of spatially coherent noise. Geophysics 55, 1496-1503.

---

Рецензии

В статье авторы предлагают различные способы обработки, а также варианты фильтрации синфазных помех, шумов для представления оптимального результата интерпретации трехмерной георадарной съемки.
Чаще всего для решения различных поисковых или археологических задач (поиск фундаментов исторических зданий и сооружений, поиск захоронений, склепов и т.д.) применяется георадариолокационное обследование территории по системе параллельных профилей. При этом расстояние между параллельными профилями выбирается оператором исходя из размера искомого локального объекта. Далее после обработки каждой отдельной радарограммы программа позволяет собрать все выполненные профили в квазитрехмерную модель. И результат интерпретации представляется обычно в виде амплитудных срезов, выполненных на определенной глубине (временные срезы).
Такая методика обработки и интерпретации георадарных данных безусловно достаточно трудоемка по времени. При этом требуется достаточно много времени как для выполнения полевых работ (каждый профиль необходимо пройти с антенным блоком), также уходит достаточно много времени на обработку и интерпретацию всех отснятых георадарных профилей. При этом, например, при выполнении полевых и камеральных работ учитывается именно количество погонных метров георадарной съемки.
В данной статье авторы рассматривают способ выполнения полевых работ для площадной съемки с помощью нескольких пар георадарных антенн импульсного типа, при этом перемещая всю модульную систему параллельно по площади исследования с помощью специального буксировочного транспортного средства.
Также авторами предлагаются различные способы и процедуры обработки полученных данных, которые применяются сразу ко всему объему георадарных данных, а именно к полученным в результате построения 3Д модели временным срезам.
Перечисленные в статье способы удаления различных помех и шумов достаточно давно известны и широко применяются при обработке отдельных радарограмм. Так, например, для подавления синфазных помех, а также аппаратного "звона" достаточно эффективно используется процедура "вычитания среднего" в окне по всем трассам. При этом пользователь, конечно, должен помнить, что при выполнении данной процедуры в заданном окне есть вероятность удалить горизонтальный участок реальной отражающей границы (особенно если это отражение получено от подошвы плиты). Для сглаживания часто применяется пространственная 2D-фильтрация в заданном окне. Для удаления многократных переотражений и повышения вертикальной разрешающей способности достаточно эффективно применяется предсказывающая деконволюция.
В статье авторы предлагают использовать перечисленные типы фильтраций, а также различные способы обработки, например вейвлет-преобразование, для обработки уже собранной 3D-модели, применяя указанные процедуры к временным срезам. Также проводится сравнение эффективности тех или иных процедур обработки.