CC BY
60
Построение скоростных моделей по данным вертикального сейсмического профилирования
1 2 Э.А. Бляс , A.C. Рожков
1 Судоводительский факультет МГТУ, кафедра высшей математики НИИ Моргеофизика
Аннотация. Рассматривается задача определения по первым вступлениям волн, регистрируемых в скважине. Ее решение разбивается на два шага. Сначала автоматически определяются первые вступления, затем полученные времена аппроксимируются непрерывной ломаной по критерию минимума среднеквадратичной погрешности аппроксимации.
Abstract. The problem of the layered model velocity identification from borehole seismic data is considered in the paper. For solving this problem two steps have been suggested. First, the travel-time arrivals have been defined, then a continuous curve has been constructed with the help of the least squared method.
1. Введение
Метод вертикального сейсмического профилирования (ВСП), основанный на регистрации сейсмических колебаний во внутренних точках среды, применяется для решения широкого круга задач при проведении геологоразведочных работ на нефть и газ. Одной из основных задач метода ВСП является задача определения скоростной модели среды по годографу прямой продольной волны (Murat et al., 1992). Эта задача разбивается на две независимые задачи. Первая состоит в снятии первых вступлений прямой волны, то есть в определении времени вступления сигнала на фоне помех. Вторая задача сводится к построению пластовой модели разреза, то есть к разбиению годографа на участки линейности и определению пластовых скоростей в выделенных пластах.
Задача автоматического определения первых вступлений и построения скоростной модели разреза рассматривалась в работе (Яновский, Боголюбский, 1976). Времена вступлений в дальнейшем использовались для вычисления пластовых скоростей. Здесь необходимо отметить, что прибавление постоянной величины ко всем временам вступлений не меняет пластовых скоростей, поэтому для построения скоростной модели достаточно знать времена вступлений с точностью до постоянного слагаемого.
В cTaHflapTHofi o6pa6oTKe пластовые cKopocTH определяются по первым вступлениям прямой волны. Для более точного определения пластовых скоростей разрабатываются различные способы снятия первых вступлений (Шопин, Литвиненко, 1976; Шехтман, 1971;
Чаплыгин и др., 1968), но их точность, очевидно, ограничена снизу. Фактически (npn дискретизации записей с шагом 1 мс) мощности слоев, в которых удается найти cKopocTH с точностью 3-4 %, составляют не менее 50-70 м. Оценим, с какой точностью необходимо oпpeдeлять время npoöera волн в слое мощностью 10 м, чтобы eKopoeTb в нем определялась с точностью до 3-5 %. Если принять, что eKopoeTb составляет величину пopядкa 2,5-4,5 м/мс, то время распространения пpямoй волны в таком слое составляет пopядкa 3 мс, поэтому погрешность его нахождения не должна превышать 0,1-0,15 мс. Таким образом, необходимо определять время вступления прямой волны с точностью до 0,1 шага дискретизации, который обычно равен 1мс.
Если находить времена первых вступлений по срыву трассы, то есть по амплитуде, превышающей некоторый порог (зависящий от уровня помех), то погрешности такого определения могут достигать величины порядка 1 мс. Вблизи первого вступления амплитуды сигнала ненамного превосходят амплитуды помех, в частности, первый отсчет сигнала сравним с амплитудой помехи. Это приводит к существенным погрешностям определения первого вступления (времени первого отсчета сигнала) по превышению амплитудой трассы уровня помех.
В то же время, как отмечалось, для построения скоростной модели среды можно определять вступления волны с точностью до постоянного слагаемого. Это означает, что можно находить характерные точки сигнала при больших амплитудах. Отметим, что использование максимальной амплитуды сигнала (то есть точки, в которой отношение сигнал/помеха максимально) не дает максимальную точность, так как на этих временах вместе с падающей волной проявляются отраженные
волны. Отсюда следует, что характерные точки сейсмических трасс необходимо искать в начальной части вступления прямой волны.
Вторая задача сводится к автоматическому определению точек излома годографа по некоторому критерию. Решение этой задачи дано в работе (Рожков, 1999). Она фактически разбивается на две подзадачи: выбор критерия оптимальности и нахождение разбиения годографа, минимизирующего целевую функцию.
В настоящей работе рассматривается решение обеих задач. Для решения первой задачи предлагается алгоритм, основанный на нахождении характерных точек кривой вблизи вступления прямой волны с последующей их аппроксимацией полиномом второй или третьей степени и занулении первой или второй производной полинома. Для второй задачи разработан алгоритм, основанный на кусочно-линейной аппроксимации годографа и предположении о нормальном распределении помех. По данным алгоритмам разработаны программы, применяющиеся при обработке данных ВСП, полученных на шельфе Баренцева моря.
2. Определение характерных точек сейсмических трасс
Рассмотрим задачу автоматического определения точки минимума и точки перегиба на участке трассы, ближайшем к первому вступлению. Реальные трассы имеют форму, показанную на рис. 1. Обычно время от первого вступления трассы до ближайшей точки минимума не превосходит 8-10 мс. Вступление прямой волны характеризуется сильным повышением амплитуд и регулярным видом записи. Для предварительного определения первых вступлений можно использовать увеличение амплитуд после вступления прямой волны, а также увеличение коэффициентов корреляции соседних трасс.
Рис. 1. Реальные трассы
Для каждой точки трассы справа и слева от этой точки отложим равные отрезки некоторой заданной длины А/. Далее возьмем среднее арифметическое абсолютных величин амплитуд всех точек первого и второго отрезков, обозначим их и •2 и найдем отношение В соответствии с
вышесказанным, вступление прямой волны характеризуется максимальным отношением •2/а,1, поэтому можно считать, что первое вступление соответствует той точке трассы, для которой это отношение максимально. При наличии на отдельных трассах сильных помех максимальное отношение может появляться до вступления прямой волны, поэтому в программу заложены дополнительные критерии. Они основаны на том, что после вступления прямой волны амплитуды имеют ненулевые значения одного знака на некотором промежутке, длина которого задается пользователем, исходя из вида трасс.
По данному алгоритму написана программа, которая прошла опробование на большом числе реальных данных. Ее тестирование показало, что погрешность для абсолютного большинства трасс не превосходит 1-2 отсчета.
После предварительного определения времен первых вступлений находим ближайший справа минимум (или максимум) трассы и между ними точку перегиба. Сначала данные точки находятся с
точностью до дискреты, а затем уточняются с использованием аппроксимации трассы полиномом 3-й степени.
Теперь рассмотрим решение второй задачи: разбиение полученных годографов на участки, соответствующие постоянной скорости. Другими словами, имея заданную в дискретных точках (с постоянным шагом) последовательность , где , - номер трассы (приемника), характеризующий глубину, ^ - время характерной точки, необходимо выбрать среди этих точек точки излома, в которых первая производная терпит разрыв. Эту же задачу можно сформулировать несколько по-иному. Требуется найти кусочно-линейную аппроксимацию функции с неизвестными точками начала и конца отрезков этой аппроксимации. Так как наблюдаемые времена известны с некоторыми погрешностями, то для решения этой задачи применим метод максимального правдоподобия. Можно считать (экспериментальные данные подтверждают это), что ошибки измерения вызваны аддитивной случайной помехой, имеющей нормальное распределение. В этом случае метод максимального правдоподобия сводится к минимизации суммы квадратов отклонений реальных данных от аппроксимируемой кривой.
3. Автоматическое определение точек излома
Так как полученные таким образом значения функции имеют погрешности измерения, то для нахождения кусочно-линейной аппроксимации поставим следующую задачу. Из бесконечного числа ломаных линий на области определения (обозначим эту область 3: , = 1,...,^ где N - число трасс) найти такую, среднеквадратичное отклонение которой от ^ минимально. В рассматриваемом случае, кроме параметров отрезков (участков ломаной, аппроксимирующей реальные данные), неизвестными являются также точки излома.
Итак, требуется найти минимум функции
п т (/' )
ф = ЕЕ - а (х, - х,) - ь ]2. (1)
/ •
Здесь п - число отрезков, на которые разбивается аппроксимирующая функция (ломаная), т() -число точек на/-ом участке, t = а/ (х - X/) + Ьj - уравнение отрезка ломаной на /-ом участке, X/ = (1+1 + I) /2 - центр /-го участка, 11,12, ..., 1п - точки разбиения ломаной. Уравнение отрезка на /'-ом участке берется в виде, содержащем середину этого участка, так как в дальнейшем это позволит упростить систему уравнений.
Длина ломаной равна N-1, где N - число трасс, начало первого отрезка имеет координату , = 1, конец последнего , = N. Пусть, например, взята некоторая ломаная, состоящая из п отрезков, удовлетворяющая требованиям, определенным выше. Таким образом, функция Ф, определенная равенством (1), зависит не только от неизвестных коэффициентов а/, Ь/, но также от точек 11, 12, ..., 1п разбиения ломаной на отрезки и числа п этих отрезков. Ясно, что без дополнительных ограничений данная задача имеет очевидное решение, когда между любыми двумя точками проводится отрезок, соединяющий эти точки, то есть все точки являются точками излома. Для практики такое решение не представляет интереса, поэтому необходимо задать дополнительные ограничения на искомые величины. С точки зрения математической статистики, получаемые после минимизации функции Ф значения коэффициентов а/, Ь/ можно рассматривать, как их оценки. В этом случае эти оценки (когда каждая точка кривой является точкой излома) имеют бесконечную дисперсию и поэтому представляют малый интерес. Другими словами, погрешности таких оценок настолько велики, что ими нельзя пользоваться. Фактически это означает, что скорости между соседними приемниками находятся с большими погрешностями. Поэтому в качестве дополнительной информации можно задавать минимальную длину отрезка ломаной.
Исходя из этого, задача сводится к нахождению минимума функции Ф, зависящей от точек разбиения ломаной и коэффициентов отрезков при условии, что задано минимальное расстояние между точками разбиения. Данную задачу будем решать в два этапа. Сначала для заранее заданного множества точек разбиения рассмотрим задачу нахождения ломаной, при которой Ф минимально. После этого рассмотрим задачу минимизации функции Ф при всех возможных точках разбиения.
Первая задача сводится к нахождению минимума функции Ф по переменным а/, Ь/. Рассмотрим ее решение.
4. Нахождение ломаной, при которой Ф минимальна, для заданного множества точек разбиения
Пусть даны опорные точки ломаной l1, l2, ..., ln+i. На каждом из участков [lm, lm+1] отрезок ломаной описывается уравнением
t = üm(x - Xm ) + bm , m = 1,2,...,П.
X1,...,Xm - координаты середин отрезков /1,...,/m на оси Ox (рис. 2), то есть
Xi = (l2 + ll) / 2, X2 = (l3 + l2) / 2, ..., Xm = (lm+1 + lm) / 2.
(2)
h x
12 "Л "4
Рис. 2. График ломаной с фиксированными опорными точками Из условия непрерывности ломаной получаем уравнения (уравнения связи):
а(2 - X;) + Ь = - Х2) + Ъ2, а( /з - Х2) + Ь2 = аз( /3 - Х3) + Ь3,
аш-1( 1т - Хш-1) + Ьш-1 = ат( 1т ~ Хт) + Ьт.
(3)
Так как l2 - X1 = l2 - (l2 + l1)/2 = (l2 - l1)/2, l2 - X2 = l2 - (l3 + l2)/2 = (l2 - l3)/2, и так далее, то уравнения (2) перепишутся в виде:
a1{l2 - 1/2 + b1 = a2(l2 - 1з)/2 + b2,
a(l3 - 2/2 + b2 = a3(l3 - l4)/2 + Ьз,
am-1 (lm - lm-1 )/2 + bm -1 = am(lm " L+i)/2 + bm.
(4)
Теперь потребуем, чтобы ломаная, заданная уравнениями (2), была такова, что ее среднеквадратичное отклонение от ^ было минимальным, то есть
ф = Ев - а;(Ху - X;) - Ь;]2+ - а*% ~ Х2) ~ Ь2 ]2 +...+ - ат (% - Хт) - Ьт]2 ^ Ш1П.
ieL
i^L
Для нахождения минимума применим метод Лагранжа. С учетом уравнений связи (3) функция Лагранжа записывается в виде
L = Ф + Л1[a1(/2 - X!) + Ь1 - а2 (/2 -X2) - Ь2] + Л2[а2(/3 -X2) + Ь2 - а3(/3 -X,,) - Ь3] +...
+ Лт-1[ат-1 (/т ~ Хт-1) + Ьт-1 ~ ат (/т ~ Xт) ~ Ьт],
где Л-1, Л-2, ..., - множители Лагранжа. Согласно методу Лагранжа, в точке минимума выполняются равенства д^да] = 0, 5L/5Ьj = 0, которые с учетом последнего равенства принимают вид
дЬ/дц = 21a-(Xi -Xj) + bj - ti ](Xi -Xj) - ÄJ.1(lj -Xj) + lj(lj+1 -Xj) = 0,
(5)
izL
дЫдЬ/ = 2Е[а/(х, ~X/) + Ь/- и] + Л) = 0, (6)
1&1/
из (5) и (6) и того, что 1/+1 -X, = -(/,- - X/), получаем:
2а, 1(х, - X/ )2 + 2Ь, !(х, - X/) = 2Ъ, (х, - X) - (¿1 + Л. )(1, +1 - X), (7)
1е1, ¡е!/ ¡е!/
2а/ Е(х, - X) + 2Ь/ Ц = 2^,- + 4^ - 4. (8)
/е!/ /еТ/
Для упрощения расчетов, для каждого /-го отрезка возьмем систему координат с началом в центре этого отрезка, т.е. в точке X/. Тогда в уравнениях (7) и (8)
I (х, -X/) = 0.
Так как I/ -X/ = -(// +1 - I) /2, то из (7) и (8) получим:
2а/ 1(х, - X / )2 = 2Ъ, (х, -X /) - Ц.-1 + 4 )(/+1 - / )/2
'Ц/ Щ
2Ь/11 = + 4-1 - 4
Положим
А/ = 2Е(х, - X )2, В/ = 2ЕА- (х, - X), С/ = Е1, = 2Е4 1&!/ 1&!/ /е^ 1&!/
Получим при / ^ 1, / ^ т:
а/ = В//А/ - (4-1 + 4)(//+1 - /)/(4А/),
Ь/ = Д /С/ + (4-1 )/(2 С/), (9)
при/ = 1 а1 = В1/А1 - 4(/2 - 11)/(4А 1), Ь1 = А/С1 - 4/(2СД
при/ = т ат = Вт /Ат ~ ^т-1 (1т+1 ~ Iт )/(4Ат ), Ьт = Дт / Ст + Кп-1/(2Ст ).
Подставив выражения а/, Ь/, / = 1,...,т в систему (3), получим систему да-1 уравнений с да-1 неизвестными Л1,...,Лт-1., причем эта система трехдиагональная, поэтому для ее решения применим метод прогонки. Найдя этим методом неизвестные Л1,...,Л1^ и подставив найденные значения в уравнения (9), мы найдем, тем самым, искомые коэффициенты а/, Ь/, которые полностью определяют искомую ломаную при данных опорных точках 11,...,1п.
4. Минимизация функции Ф при всех возможных точках разбиения
Рассмотрим теперь вторую задачу. Перебрав все возможные случаи таких ломаных, координаты соединения любых двух последовательных отрезков которых - целые числа между 1 и п (очевидно, что число таких возможных наборов координат ограничено), и для каждого из этих случаев найдя такую ломаную, среднеквадратичное отклонение которой от ^ минимально, и затем выбрав из полученных ломаных такую, у которой среднеквадратичное отклонение минимально по сравнению с другими ломаными, мы и найдем такую ломаную, которая из всех возможных ломаных (образно говоря) наиболее плотно прилегает к и. Тем самым мы и решим поставленную задачу.
Рассмотрим какую-либо ломаную, определенную на отрезке (1, п). Пусть ее начало имеет координату 11, конец - координату 1т+1, т - число отрезков ломаной, 12 - координата соединения 1-го и 2-го отрезка, 13- координата соединения 2-го и 3-го отрезка,..., 1т - координата соединения (да-1)-го и да-го отрезка.
Получим l1,...,lm+1 опорных точек ломаной. Из этих точек l1 и lm+1 фиксированы, a l2,..,lm могут варьироваться. Теперь достаточно перебрать все возможные варианты расположения точек l2,...,lm. Для простоты расчетов, не умаляя общности, можно поступить следующим образом.
1) Вводим минимальную длину отрезков ломаной mindl.
2) Будем разбивать не всю область точек i = 1,...,n. Вместо этого поступим следующим образом. Возьмем отрезок разбиения равный шести минимальным отрезкам 6xmindl=otrraz. Выберем точку l3, отстоящую от начальной координаты l1 на шесть минимальных отрезков, то есть на otrraz, и рассмотрим все допустимые разбиения участка [l1, l3] на отрезки, с учетом того, что минимально допустимый отрезок имеет длину mindl.
Сначала рассмотрим все возможные разбиения на два отрезка. Если у нас два отрезка, то ломаная имеет три опорные точки. Из них две точки будут фиксированными l1 и l3, a l2 может принимать лишь значения от k1 = l1 + mindl, до k2 = l3 - mindl. (Это в случае, если k2 > ki или k2 = k1, иначе у нас имеется лишь один отрезок, и для него существует тривиальное решение.)
Итак, перебирая значения l2 от k1 до k2, мы переберем все разбиения отрезка длиной otrraz с началом в точке l1 на два отрезка. Для каждого из этих разбиений отыщем ломаную с этими опорными точками, которая наиболее плотно к реальной функции ti на этом же отрезке. Из всех этих ломаных выберем, наконец, ту, которая плотнее всего прилегает к t. Запомним ее точки разбиения l1, l2, l3. Затем способом, аналогичным описанному выше, будем разбивать отрезок разбиения всеми возможными способами на три отрезка, и так же из всех таких разбиений выберем то, для которого существует ломаная, наиболее плотно прилегающая к t. Затем сравним полученную ломаную с той, которую выбрали при разбиениях на два отрезка и выберем из этих двух ломаных ту, среднеквадратичное отклонение которой от ti минимально. Запомним ее точки разбиения. Аналогично будем действовать и далее, деля отрезок разбиения на 4, 5, ив конце концов на шесть отрезков. Такое разбиение единственно, так как длина отрезка равна 6xmindl. Таким образом мы выберем из всевозможных ломаных на этом отрезке разбиения такую, которая на нем плотнее всего прилегает к ti.
Далее сдвинем отрезок разбиения на длину первого из отрезков найденной ломаной. То есть будем считать, что мы нашли вторую опорную точку l2 искомой ломаной. Зафиксируем найденное l2. Затем отложим от найденного l2 отрезок длиной otrraz. Получим следующий отрезок разбиения с началом в точке l2 и концом в точке l2 + otrraz. Поступая аналогично тому, что описывалось выше, получим в конце концов наиболее плотно прилегающую ломаную к ti на этом отрезке разбиения. Запомним длину 1-го отрезка этой ломаной и отложим его от найденного l2. Найдем таким образом l3 и будем считать, что нашли третью опорную точку искомой ломаной. Продолжая этот процесс дальше, найдем таким образом все опорные точки. Затем уже для найденных опорных точек, из всех ломаных, имеющих эти опорные точки, выберем наиболее плотно прилегающую к ti на отрезке (1, n). То есть найдем искомую функцию y = f(i), которая удовлетворяет всем требуемым условиям. Она является ломаной, наиболее плотно прилегающей к ti на отрезке (1, n).
Литература
Murat M.E. and Rudman A.J. Automated first arrival picking: a neural network approach. Geophysical
prospecting, v.40, 1 6, p.587-605, 1992. Волков E.A. Численные методы. M., Наука, с.248, 1987.
Рожков A.C. Способ автоматического определения точек излома вертикального годографа ВСП.
ВИНИТИ РАН, № 2268-В99, 8 е., 1999. Чаплыгин Э.Н., Клушин C.B., Федотова Г.А. Способы выделения сейсмических границ по годографам волн, регистрируемых при вертикальном сейсмическом профилировании. Материалы 2-й научной конференции молодых геологов Белоруссии, Минск, с.184-186, 1968. Шехтман Г.А. Использование кажущихся скоростей при интерпретации результатов ВСП.
М., Разведочная геофизика, вып. 44. с.18-24, 1971. Шопин Ю.Г., Литвиненко Е.Р. Определение положения отражающих границ в околоскважинном пространстве по данным вертикального сейсмического профилирования. Л., Записки Ленингр. горного ин-та, № 2, с.11-16, 1976. Ямпольская Е.Е. Определение скоростной характеристики среды по данным скважинных и наземных наблюдений. Новые геофизические методы изучения геоструктуры и физических свойств осадочного чехла Прикаспийской впадины. Саратов, с.61-71, 1986.
Яновский А.К., Боголюбский А.Д. Способ автоматической аппроксимации вертикального годографа, основанный на последовательном выделении пластов. Прикладная геофизика, вып. 82, с.95-100, 1976.
