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

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

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

Аннотация

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

Введение

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

Для оценки скорости разработано несколько методов, хорошо зарекомендовавших себя на практике. Например, корреляционный метод [2] на протяжении ряда лет используется при выполнении акустического каротажа фирмой 8сЫитЪе^ег. Для обработки данных с высокой зашумленностью был разработан метод предсказания [3], и его модификация [4]. Однако этот метод не позволяет проводить дисперсионный анализ и, таким образом, может с достаточной надежностью применяться в анализе только продольных и поперечных волн, исключая дисперсные волны Стоунли.

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

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

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

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

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

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

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

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

Расчет трубной волны

Необходимо рассчитать акустическую волну, распространяющуюся в скважине (рис.1). Известны следующие параметры: Я - радиус скважины, а -радиус зонда. Зонд предполагается идеально упругим. Ур и V, - продольная и поперечная скорости акустической волны в породе, р- плотность породы. V - скорость акустической волны в жидкости, р^ -плотность жидкости.

R

передатчик

A

приемник 1

приемник 2

Рис. 1. Вид скважины и зонда с одним передатчиком и двумя приемниками

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

р(а, г, /) =- 1

(П)

да да

: J J [[(a)A(a,к,а)]ехр{ikzz}exp{-imt}dkdm

(1)

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

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

N (а, к ,а)

A(a, к ,a) = -

(2)

D(a, к ,a) где

N(a, к, я) =

= K ( fa )/„( fR) -Ko( fK)U fa) - g[Ki(JR)I0 ( fa) + K0 ( fa)Il(fR)] (3)

fKi(fa) '

D(a,к,a) - дисперсионное уравнение:

D(a,к,a) = Io(fR) + ^^K0(fR) -Ki( fa)

ii( fR) - Kf Ki( fR)

Ki( fa)

где

g =

fP f 2^/ ps

pfp I к2c2

1 + 2V? Ko(sR)

sR c K,(sR)

Г 2Jl -V

Ko( pR)

Ki( pR)

(5)

p=iк-V2,s = *-w

f = 4к2 -

Vf2

'к-

(6)

Здесь 1п и Кп (п=0,1) - модифицированные

функции Бесселя первого и второго рода, соответственно. Затухание волны в жидкости и в породе моделировалось путем введения комплексных скоростей [14]

vp ^

Vf ^

V, ^

(7)

где Qp, Qs и Qf - параметры, характеризующие затухание продольных, поперечных волн и волн в жидкости.

Интеграл (1), как правило, вычисляют методом RAI (real-axis integration) [12]. В этом случае, в рассчитанном сигнале будут присутствовать все типы волн. Чтобы оценить качество выделения итеративным алгоритмом одного из типов волн, мы вычислили отдельно трубную волну. Для этого, в соответствии с теоремой о вычетах [13], был рассчитан вклад в общий сигнал полюсов выражения (2), которые соответствуют корням дисперсионного уравнения (4).

Тогда вклад трубной волны в общий сигнал записывается в виде:

ч 1

p(a, z, t ) = — > 2п

: J X(a){2niR(a,a,z)}exp{-iat}da

(8)

где вычет, соответствующий корню дисперсионного уравнения для частоты а , имеет вид:

R(a,a, z ) =

1 N (a, к,a)

2п dD(a,к,ю).

-exp faz}

/дк

(9)

к=кр (a)

2

2

2

Здесь

= ^ fЯ) + /Я) - K1( /Я) - я'

дk / ^ / К,2 (/а) / К/) ,и

Л( /Я) - ^^ К,( /Я)

К,( /а)

Кк_

/

10(/а)-/ + ^^х( К0(/Я) + К^ 1-/Я) /Я К,(/а) ( ^ ' /а ] ЯК2(/а) ^

/а

д = 1 10(/а) -1 К,(/а) + ( К0(/а) + 11,(/а)

/а

Р/

о+/

2К

, р 2 р, "¡Икр +к, "к"

1_ + К,(,Я)

с2 К^К)

2У.

— 1

Як

1 + -

т= - к - V Як

с2 ,

1-

К^Я)! К0(,Я) +

К.(,Я)

К2(,Я)

Корни дисперсионного уравнения D(a, к,ю) находились методом Ньютона по формуле:

к1+1 = к1

D(a, к ,ю)

дD (а, к ,ю) дк

Здесь I - номер итерации метода Ньютона.

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

,!( ю) = А(ю)ехр(р( ю)). (10)

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

,2(ю) = ю)ехр(/'к(ю)/), (11)

где I - расстояние между приемниками, к(ю) -комплексное волновое число. Запишем к (ю) в виде:

к(ю) = кг (ю) + а(ю), (12)

где кг (ю) - описывает скорость волны, а(ю) - затухание волны. Тогда (4) можно переписать в виде: ,2(ю) = А( ю)ехр(-а(ю)1)х х ехр(/'(р(ю) + кг (ю)1)) Необходимо оценить параметры кг (ю) и а(ю) по измеренным сигналам на приемниках во временной плоскости.

Вычислим натуральный логарифм ,2(ю) , и выделим действительную и мнимую части:

(13)

К0(рЯ) I К0(рЯ) +

К.( рЯ) рЯ

К2( рЯ)

ЭТ [1п ,2(ю)] = 1п А(ю) -а(ю)1. 3[1п ,2(ю)] = Развернутое значение [р(ю) - кг (ю)1 ]

(14)

(15)

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

Адаптивный итеративный алгоритм восстановления искаженного фрагмента трубной волны

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

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

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

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

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

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

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

1. Полагаем отклики на первом и втором приемниках эталонами и начальными приближениями итеративного процесса:

1 ^(а) = Ф{в($)}, 2Setl(а) = Ф{2)} ;

1 ^(О = ! S(t), 2 So(t) = 2 S(t).

2. Находим Фурье-спектры волн на приемниках:

1 ^(а) = Ф{1 ^(t)} = !Ап(а)ехр[/(1%(а))],

25(а) = Ф{2^)} = 2Ап(а)ехр[/(2^и(а))],

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

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

О = {а : Ап (а) > Ь};

(16)

шт

, 2 Ап (а)

1п 2 п +а_а- сп

1 Ап (а)

здесь Ь задает минимальный уровень спектральных компонент, при котором они учитываются в выражении (16).

Выполняется процедура адаптации, т.е. высокочастотные компоненты (а г О) сигналов принимаются равными, соответственно, ! 8е„ (а) и 2Беа(а), а затем умножаются на параметр адаптации ё (0 < ё < 1) .

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

А+ДО = Ф-111 А(а)ехр

2 ^(О = Ф-112 А(а)ехр

2 Р

(1+ап).

2Рап

а+сп

ехр(;^(а))|

(1+ап)

а+сп

ехр<2^(а))|

где р = [1А(а) + 2 А(аУ2.

5. Если сходимость итеративного алгоритма замедлилась, то выполняется процедура установки эталона:

1 ^(а) = ф{,^(О} , 2(а) = Ф{2^+1«} .

Далее выполняем замены

!spn+l(t) = 1 s(t), 2spn+l(t) = 2s(t), t еП ;

здесь П - определяет «достоверную» область сигнала, в которой, как полагается, отсутствуют наложения других типов волн.

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

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

Пусть трубная волна регистрируется двумя приемниками, отстоящими друг от друга на 1,22 м. Расстояние между излучателем и ближайшим к нему приемником составляет 4,88 м. Параметры, задающие затухание в породе и жидкости в соответствии с формулами (7), полагаются равными Qp = 100 , Qs = 50 и Qf = 20 . В качестве импульса

излучателя выбрана вторая производная функции Блэкмана-Харриса [14]. Для рассчитанных трубных волн вычисляется скорость, используя какой-либо из известных методов [2-7]. В соответствии с вычисленной скоростью две трубные волны совмещаются, как показано на рис. 2.

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

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

Рис. 2. Совмещенные по времени трубные волны на двух приемниках

26.0- IV

-1-1-1-г—-' -340.0 -630.0 -420.0 -21 0.0 Т I | | | - 21 0.0 420.0 630.0 840.0

-26.0-

-52.0-

-78.0- 1

Рис. 3. Совмещенные трубные волны после обнуления их значений во временном интервале, соответствующем зоне наложения других типов акустических волн Модуль спектра трубной волны на 1-м приемнике представлен на рис. 4. Интервал дискретизации в частотной плоскости уменьшен в 16 раз путем добавления соответствующего количества нулей к трубной волне во временной плоскости перед выполнением Фурье-преобразования.

Рис. 4. Модуль спектра трубной волны на 1-м приемнике

Действительная часть комплексного волнового числа, рассчитанная гомоморфным методом, приведена на рис. 5. Заметны отклонения функции от линейности, что приводит к погрешности начальной оценки величины затухания трубной волны. Линейная аппроксимация в соответствии с формулой (16) дает значение оценки 0,00095.

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

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

-840.0 -630.0 1 1 -420.0 -21 0.0 21 0.0 420.0 I 1 630.0 840.0

-0.1 р \

/0.24 - \

/ -0.36- \

/ -0.48- \

/ -0.6- \

Рис. 5. Действительная часть комплексного волнового числа, рассчитанная гомоморфным методом

В процессе итераций произошло уточнение начальной оценки величины затухания. Действительная часть комплексного волнового числа, рассчитанная адаптивным итеративным методом за 25 итераций, показана на рис. 6.

-840.0 -630.0

1 1 (20.0-210.0 { I 1 1 21 0.0 420.0

-ол/-

/22-

/ -0.33-

/ -0.44-

^ -0.55-

630.0 840.0

Рис. 6. Действительная часть комплексного волнового числа, рассчитанная адаптивным итеративным методом

Первые 10 итераций выполнялись с коэффициентом адаптивности равным нулю. Итеративный процесс сошелся к значению оценки величины затухания 0,00103. Затем была произведена установка эталона во временной плоскости и коэффициент адаптивности в спектральной плоскости положен равным единице. После этого были выполнены еще 15 итераций. Итеративный процесс сошелся к точному значению величины затухания 0,00113. Таким образом, предложенный в данной работе адаптивный итеративный метод позволил повысить точность оценки величины затухания трубной волны на 20% по сравнению с гомоморфным методом.

Модули восстановленных спектров трубных волн на 1-м и 2-м приемниках показаны на рис. 8.

Среднеквадратичное отклонение волн, восстановленных адаптивным итеративным методом и изображенных на рис. 9, от эталонных волн, изображенных на рис. 2, не превышает 1-2%.

Рис. 7. Зависимость величины затухания от количества итераций

й A

i 01П l\

// 0.08- -1-Г^-1-1-

Ц 0.04-

-840.0 -630.0 -420.0 -21 0.0

21 0.0 420.0 630.0 840.0

Рис. 8. Модули спектров трубных волн на 1-м и 2-м приемниках, рассчитанные адаптивным итеративным методом

Рис. 9. Трубные волны на 1-м и 2-м приемниках, восстановленные адаптивным итеративным методом

Выводы

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

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

Литература

1.

2.

3.

4.

5.

6.

7.

C.H. Cheng, J. Zhang, D.R. Burns Effects of in-situ permeability on the propagation of Stoneley (tube) waves in a borehole // Geophysics, 1987, v.52. C. V. Kimball, and T.L. Marzetta Semblance processing of borehole acoustic array data // Geophysics, 1986. V.49.

X.M. Tang Predictive processing of array acoustic waveform data // Geophysics, Soc. of Expl. Geo-phys., 1997. V.62.

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

R. Prony, Essai experimental et analytique // L'ecole Polytech., 1975. V.1.

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

K.J. Ellefsen, C.H. Chengand, K.M. Tubman Estimating phase velocity and attenuation of guided waves in acoustic logging data // Geophysics, 1989, v.54. K. J. Ellefsen, D.R. Burns, C.H. Cheng Homomorphic processing of the tube wave generated during acoustic logging // Geophysics, 1993. V. 58. Chang S.K. et al. Method and apparatus for discrete-frequency tube-wave logging of boreholes // US Patent 5,331,604, 1994.

10. Казанский Н.Л., Серафимович П.Г., Харитонов С. И. Итеративный алгоритм расчета скорости и затухания трубных волн по данным акустического каротажа // Известия Самарского научного центра Российской академии наук, 2001. Т. 3, №1.

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

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

13. X.M. Tang, C.H. Cheng Effects of a logging tool on the Stoneley waves in elastic and porous boreholes // The Log Analyst, 1993, September-October.

14. C.H. Cheng, M.N. Toksoz, M.E. Willis // Determination of in situ attenuation from full-waveform acoustic logs // Journal of Geophysical Research, 1982. V.87.

9.