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

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

УДК 550.83+535.12

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

© 2001 Н.Л. Казанский, П.Г. Серафимович, С.И. Харитонов Институт систем обработки изображений РАН, г. Самара

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

Введение

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

Распространенным методом обработки данных акустического каротажа являются метод Прони [3] и его модификации [4, 5]. Недостатками метода Прони являются большое количество выбросов в значениях оцениваемых параметров и большая погрешность в оценке коэффициента затухаемости. Чтобы преодолеть эти недостатки был разработан гомоморфный метод [1], который будет подробно рассмотрен далее. Однако при малом количестве приемников и высокой зашумленности данных, как будет показано ниже, гомоморфный метод приводит к большой погрешности в оценке затухаемости волны. В качестве шумов могут присутствовать как измерительная погрешность аппаратуры, так и искажения вследствие наложения различных типов волн.

Для обработки данных с высокой зашумленностью был разработан метод предсказания [6, 7]. Этот метод устойчив к шумам данных, но позволяет оценивать только скорость волны. Кроме того, он не позволяет проводить дисперсионный анализ и, таким образом, может с достаточной надежностью применяться в анализе только продольных и поперечных волн, исключая дисперсные волны Стоунли.

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

В методе предсказания [6] расчеты выполняются только во временной плоскости. Это исключает возможность проведения дисперсионного анализа.

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

Идея использования и временной и частотной плоскостей в итеративном подходе к оценке требуемых параметров не нова. Она применялась в оптике, астрономии [8, 9], при расчете дифракционных оптических элементов [10, 11].

Результаты работы предложенного итеративного алгоритма сравниваются с гомо-

морфным методом. Приведем описание гомоморфного метода.

Гомоморфный метод определения скорости и затухания трубной волны

В общем случае для нахождения требуемых параметров волны необходимо решить волновое уравнение [12]:

1

s( z’ *) = ТТ^ х (2р)

от от

х I da exp(-i at) | dkX(k, a) exp(-ikz)

(1)

где 5( z, ?) - отклик, записанный на приемнике, г - расстояние между излучателем и приемником, t - время, а - частота, к - волновое число. Выражение X(к,а) включает в себя информацию о геометрии модели, излучателе, граничных условиях и т.д. Для реализации и гомоморфного и итеративного алгоритмов не требуется знать точный вид выражения X(к,а).

Записывая (1) в частотной плоскости и выполняя интегрирование по контуру в плоскости комплексного волнового числа, получим выражение, состоящее из суммы вычетов, охваченных контуром, и двух ветвящих-

ся интегралов

1 [13]:

Jcut L J

Sj( a) = A(a) exp(ij(a)).

(3)

Спектр волны на втором приемнике описывается в виде:

s2 (а) = $( а) ехр(гк (а )1), (4)

где I - расстояние между приемниками, к (а)

(6)

1 от

S(z, а) = — I dkX(к, а) exp(-ikz) =

2p J

—от

= /VЯ, exp(-ikz) + I (2)

■* Jcut

j

В некоторых частотных диапазонах основной вклад в (2) вносит лишь один вычет. Экспериментально установлено [1], что так происходит для трубной волны в диапазоне от 0 до 5 kHz.

Пусть спектр волны на первом приемнике имеет вид:

- комплексное волновое число.

Запишем к(а) в виде:

к (а) = кг (а) + га (а), (5)

где кг (а) - описывает скорость волны, а (а)

- затухание волны. Тогда (4) можно переписать в виде:

52 (а) = А( а) ехр(-а(а)/) х х ехр(г( ^ (а) + кг (а)/))

Необходимо оценить параметры кг (а)

и а (а) по измеренным сигналам на приемниках во временной плоскости.

Вычислим натуральный логарифм

s2(а) и выделим действительную и мнимую части:

Я[1п s2(а )] = 1п А(а) - а (а )1, (7)

3[1п 5 2 (а)] =

= Развернутое значение \(р(а) - кг (а)/]. (8)

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

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

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

енты прямых - кг (со) и а(со), но не заканчиваем на этом расчеты, а переходим опять во временную плоскость. Здесь, получившуюся фазу мы оставляем неизменной, а результирующую амплитуду заменяем на амплитуду отклика на приемниках. На этом единичная итерация заканчивается.

Кроме того, в процессе итераций мы

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

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

Пусть 1 s(t) и 2 s(t) - отклики на первом и втором приемниках, соответственно.

15 (№>) = Ф{ ^0]=1 А(«)ехр[г(1ф(«))] - фу-рье-спектр отклика на первом приемнике.

1. Полагаем отклик на втором приемнике начальным приближением 2 s(t)=2 s0(t) .

2. Находим Фурье-спектр

2 5п (^) = ф{2 sn ^)}= 2 Ап (ехр[г(2 Фп («))] ,

где п - номер итерации.

3. Находим значения ап и сп из условия

Ш1П і

I

ln

2 An (W)

+ a nW - Cn

(9)

іA(w)

находим значения kn из условия

min ш 1(р(со) + кпсо-1(рп(со)^\ (10)

кп Y со J

4. Выполняем обратное Фурье-преобра-зованне

2 Sn +1

(t) = Ф

-1 Г1 A(w) exp(-a nW + Cn) х]

х ехр(! р(а) + к„а) 5. Выполняем замену

2 ^ )

, SPn+1 (t) :

(t)

[2 S„+1(t)]

Переходим к шагу 2 или заканчиваем итерации, если изменение 2 sn (t) меньше заданной погрешности.

Численные результаты

В качестве тестового сигнала на первом приемнике использовался (рис.1)

/20^ N

і s(ti) =cos(—(i - —)) exp

(i - N )2

2

А2

. (11)

где N = 128 - количество отсчетов сигнала,

А = 20, і = 1...128 .

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

2 ^(А-) = Ф"1 {ехР(«і)Ф[і s(ti ХІ^ (12)

На рис.2 показаны амплитуды сигнала на первом приемнике и, получившегося в результате операции (12), сигнала на втором приемнике.

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

Рис. 2. Сигнал на первом приемнике и сигнал на втором приемнике с коэффициентом затухания

а = 0.003

Рис. 3. Фаза сигнала на втором приемнике

при наложении друг на друга поперечной волны и волны Стоунли. Выражения (7) и (8) для модели (6) не учитывают возможности такого искажения фазы. Таким образом, цель предлагаемого итеративного процесса - точное восстановление сигнала (в том числе -фазы на втором приемнике, рис.3, включая оценку значений запаздывания и затухания сигнала, рис.4).

На рис.4 показаны амплитуды сигнала на первом приемнике и, получившегося в результате операции (12) и сдвига по времени для к = 0.1, сигнала на втором приемнике.

Эффективность работы предложенного итеративного алгоритма проверялась на зашумленных сигналах, показанных на рис.5. На этом рисунке представлены сигналы на первом и втором приемниках с коэффициентом затухания а = 0.003 и запаздыванием к=0.1, искаженные гауссовым шумом с дисперсией й = 0.01.

Для генерации гауссова шума из равномерной случайной выборки использовался алгоритм Бокса-Мюллера. Генерация равномерной случайной выборки выполнялась с помощью алгоритма Ранеку [14].

Рис. 4. Сигнал на первом приемнике и сигнал на втором приемнике с коэффициентом затухания

а = 0.003 и запаздыванием к = 0.1

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

Действительно, зашумленный сигнал на приемнике можно записать в виде:

s(а) = (А( а) ехр(-а(а)/) + еа) х х ехр(г(р (а) + кг (а)/ + е г)) ,

где еа и ег - шумовые компоненты амплитуды и фазы, соответственно. После логарифмирования вариация фазы записывается в виде:

А(3[1п ^(а)]} = А( ).

Таким образом, преобразование не меняет вариацию фазы, и, следовательно, оценку скорости.

Используя разложение в ряд, вариация амплитуды записывается в виде:

А(1пэ(а)]»А(еа )ехр<а(:)1).

А(а)

Т.е. вариация амплитуды растет экспоненциально с ростом расстояния / между приемниками, что затрудняет оценку величины затухания.

Для примера, представленного на рис.5, при оценке величины запаздывания алгоритм сошелся к правильному значению за две итерации, причем после первой итерации (что соответствует оценке гомоморфным методом)

втором приемнике с коэффициентом затухания а = 0.003 и запаздыванием к = 0.1, искаженные гауссовым шумом с дисперсией ё = 0.01

Рис.6. Зависимость оценки коэффициента затухаемости от числа итераций

ошибка составляла всего 2%. Для оценки коэффициента затухаемости потребовалось 10 итераций, чтобы получить ошибку менее 1 % (рис.6).

Выводы

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

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

1. Ellefsen K.J., Burns D.R.,Cheng C.H. Homomorphic processing of the tube wave generated durng acoustic logging // Geophysics. V.58. 1993.

2. Cheng C.H., Zhang J., Burns D.R. Effects of in-situ permeability on the propagation of

Stoneley (tube) waves in a borehole // Geophysics. V.52. 1987.

3. Prony R. Essai experimental et analytique // L’ ecole Polytech. V.1. 1975.

4. Lang S.W., Kurkjian A.L., McClellan J.H., Morris C.F., Parks T. W. Estimating slowness dispersion from arrays of sonic logging waveforms // Geophysics. V.52. 1987.

5. Ellefsen K.J., Cheng C.H., Tubman K.M. Estimating phase velocity and attenuation of guided waves in acoustic logging data // Geophysics. V.54. 1989.

6. Tang X.M. Predictive processing of array acoustic waveform data // Geophysics, Soc. of Expl. Geophys. V.62. 1997.

7. Chunduru R.K., Tang X.M. Automated velocity analysis of array acoustic waveform data // SEG Expanded Abstracts. 1998.

8. Gerchberg R.W., Saxton W.O. A practical algorithm for the determination ofthe phase from image and diffraction plane pictures // Optik. V.35. 1972. №237.

9. Fienup J.R. Phase retrieval algorithms: a comparison // Applied Optics. V.21. 1982. №15.

10. Методы компьютерной оптики // Под ред. В.А. Сойфера. М.: Физматлит, 2000.

11. Kotlyar V.V., SeraphimovichP.G., Soifer V.A. An iterative weight-based method for calculating kinoforms // Proc. SPIE "Image Processing and Computer Optics". V.2363. 1994.

12. Cheng C.H., Toksoz M.N. Elastic wave propagation in a fluid-filled borehole and synthetic acoustic logs // Geophysics. V.46. 1981.

13.Kurkjian A.L. Numerical computation of individual far-field arrivals excited by an acoustic source in a borehole // Geophysics. V.50. 1985.

14.L'EcuyerP. // Commun. ACM. 1988.

ITERATIVE ALGORITHM FOR CALCULATING VELOCITY AND ATTENUATION OF THE TUBE WAVE MEASURED DURING ACOUSTIC LOGGING

© 2001 N.L. Kazanskiy, P.G. Seraphimovich, S.I. Kharitonov

Image Processing Systems Institute of the Russian Academy of Sciences, Samara

We develop a new iterative algorithm to obtain estimates for the velocity and attenuation of the tube wave, which is measured during acoustic logging by a few receivers. A distinguishing feature of the algorithm is that optimization of the estimated parameters is performed in both the space plane and the frequency plane. On the one hand, it provides a high degree of tolerance of our algorithm of noise and on the other permits one to conduct dispersion analysis of wave modes, for example, Stoneley waves. A comparison of the new algorithm and the homomorphic algorithm is made.