CC BY
14АПОСТЕРИОРНЫЕ АЛГОРИТМЫ В РЕШЕНИИ ЗАДАЧ ГЕОФИЗИЧЕСКОГО МОНИТОРИНГА
Г. М. Воскобойникова, М.С. Хайретдинов*
Институт вычислительной математики и математической геофизики СО РАН,
630090, Новосибирск, Россия * Новосибирский государственный технический университет, 630073, Новосибирск, Россия
УДК 519.254
Рассматривается проблема повышения точности решения обратной задачи восстановления параметров сейсмических источников за счет применения более совершенного в сравнении с известными алгоритма обработки сейсмических данных. В связи с этим предложен и исследуется новый подход по определению параметров сейсмических волн, в рамках которого решение задач обнаружения и выделения волновых форм в шумах находится в едином процессе дискретной оптимизации. Эффективность подхода иллюстрируется на ряде численных экспериментов и примере решения модельной задачи мониторинга положения скважинного источника, вытекающей из проблемы локации в процессе нефтепромыслового бурения.
Ключевые слова: геофизический мониторинг, природные и техногенные события, апостериорные алгоритмы, обнаружение и выделение, численные эксперименты, сейсмическая локация, скважинный источник.
The problem of accuracy decision of the inverse problem for restoration of seismic sources parameters at the expense of application more perfect in comparison with known algorithm for processing of seismic data is considered. In this connection it is offered and the new approach by definition of seismic waves parameters in which frameworks the decision of detection and allocation problems of wave forms in noise are in uniform process of discrete optimization is investigated. Efficiency of the approach is illustrated on a number of numerical experiments and an example of a modelling problem the decision for monitoring position borehole a source. It is following from a problem in the course of oil-field drilling.
Key words: geophysical monitoring, natural and technogenic events, posteriori algorithms, detection and allocation, numerical experiments, seismic location, borehole source.
Введение. Проблемы геофизического мониторинга, прогнозирования и предупреждения чрезвычайных ситуаций природного и техногенного характера относятся к числу приоритетных современных проблем. Под мониторингом здесь подразумевается процесс регулярного слежения за возникновением природных и техногенных событий, а также за предшествующими им геодинамическими процессами в земле. К числу таких событий относятся землетрясения, извержения вулканов, лунно-солнечные приливы, оползни, падающие небесные тела, карьерные взрывы и др. Мониторинг базируется на ряде последовательных этапов, включающих в себя удаленную регистрацию откликов событий —
Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (гранты № 14-07-00518, № 15-07-10120), Междисциплинарного интеграционного проекта СО РАН № 14
сейсмических волн, измерение их времен прихода, либо же исходных волновых форм. На конечном этапе решаются обратные задачи восстановления параметров источника и вмещающей его среды, связанные либо с обращением времен пробега волн для обратного проектирования невязок времен пробега, либо же с обращением волновых форм для обратного проектирования невязок записей сейсмических волн [1]. При этом популярным методом решения обратных задач является метод наименьших квадратов, В то же время известна его критичность к грубым ошибкам измерений в исходных данных, что определяет ограниченность метода [1]. В связи с этим повышение точности оценивания параметров волн в шумах является одной из первостепенных задач.
Существуют два подхода к решению данной задачи: последовательный („on-line") и
апостериорный („off-line"). Последовательный подход ориентирован па получение наибо-
"
'Такой подход, в частности, лежит в основе семейства алгоритмов обнаружения моментов изменения свойств сигналов [2]. К их числу относится алгоритм авторегрессии проинтегрированного скользящего среднего (АРПСС) [3], нашедший свое применение для обнаружения землетрясений [3], промышленных взрывов [4, 5]. В последние годы при решении задач геофизического мониторинга, анализа сейсмических данных нашли широкое практическое применение алгоритмы вейвлет-фильтрации [6, 7, 8].
Главное преимущество апостериорного подхода перед последовательным состоит в том, что он ориентирован на получение оптимального (по всем накопленным данным) решения. Этот подход потенциально более точен, однако его алгоритмическая реализация сопряжена с решением трудоемких в вычислительном плане задач дискретной оптимизации. Это
обусловливает необходимость поэтапного решения задач с применением существующих
"
пия подзадач на каждом из этапов полученное в результате итоговое решение может не совпадать с оптимальным, т, к, решение, найденное по условным экстремумам, как известно, в общем случае не обязано совпадать с оптимальным, С учетом этого в данной работе предлагается и исследуется иной, елабоизученный применительно к геофизическому мониторингу подход, в рамках которого решение задачи находится в едином процессе дискретной оптимизации без разбиения задачи на этапы. При этом возможны две формы обнаружения: либо оценить непосредственно времена прихода волн, либо получить одновременно оценки времен прихода и формы волновых импульсов. Среди этого класса алгоритмов будем в дальнейшем рассматривать алгоритмы, предназначенные для обработки последовательностей, изменяющих свои свойства квазипериодически [9, 10]. Это означает, что временной интервал между двумя последовательными импульсами ограничен сверху и снизу заданными константами,
В работе приводятся результаты численных экспериментов по оцениванию точности и помехоустойчивости предложенных алгоритмов и их приложения к решению реальной задачи локации скважинного источника в виде маломощного взрыва, производимого на разных глубинах,
1. Апостериорные алгоритмы определения параметров волновых форм в шумах. В этом разделе приводится обоснование апостериорных алгоритмов, обеспечивающих решение задач обнаружения и оценивания волновых форм, представленных квазипериодической последовательностью и искаженных гауссовой помехой. Под волновыми формами подразумевается набор откликов среды, образуемых при многоканальной регистрации определенного события. Рассматриваются два варианта представления волновых
форм в квазипериодической последовательности — одинаковых и разных. При решении задач предлагается следующая модель анализируемых данных. Пусть компоненты вектора X = (я0,... ,хм-1) € образуют последовательность волновых форм в виде
M
Y,un-nm (m), n = 0,...,N - 1, (1)
m=1
где (n1,..., nN) E 0
M
Mm_____
M,
Qj
Mm in
0m = { (ni,..., nN) | 0 < ni < Tmax - q; N - Tmax < Пм < N - q; q < Tmin < nm - nm-i < Tmax, m = 2, . . . , M }
(2)
Дополнительно примем, что щ(т) = 0, если j = 0,... ,q- 1, при каждом m = 1,..., M, a Mmin и Mmax находятся из решения системы неравенств, входящих в определение (2), в котором q, Tmin и Tmax — натуральные числа, Положим, Um = (u0(m),... ,uq-1(m)),m = 1,..., M, Допустим, что 0 < ||Um||2 < = 1,..., M. В веде м w = (U1,..., UM ) и n = (n1,..., nN), Согласно введенным обозначениям, вектор X зависит от пары наборов n и w, содержащих одинаковое число M элементов, т.е. X = X(n,w). Пусть случайный вектор Y = (y0,...,yN-1) есть сумма двух независимых векторов Y = X(n, w) + E, гДе E = (e0,... ,eN-1) e Фх,а21< ж, Здесь через Фх,а21 обозначено нормальное распределение с параметрами (0,а2/),
С учетом этого задача обнаружения квазипериодических последовательностей волно-
n
с которым был порожден ненаблюдаемый вектор X(n, w), В этой модели компоненты векторов Y и X соответствуют наблюдаемому и ненаблюдаемому сигналу, а компоненты вектора Е — помехе. Номера компонент векторов ассоциируются с равномерным дискретным временем. Элементам набора (n1,...,nM) сопоставляются моменты времени вступления (начала) волновых форм т) = (П1,...,Пм)т; q-мерный набор Um,m = 1,...,M соответствует волновой форме. Значения Tmin и Tmax интерпретируются как максимальный и минимальный интервалы между двумя последовательными формами.
Для решения задач применяется принцип максимального правдоподобия. Авторами показано, что помехоустойчивое максимально правдоподобное обнаружение заданного числа неизвестных волновых форм моделируется следующей дискретной экстремальной задачей:
Задача 1, Дано: числовая последовательность Y=(y0,...,yN-1), натуральные числа q, M, Tmin и Tmax, Найти: набор n = (n1,..., nN) E 0M, такой, что
M q-1
F(n1,..., nM) = Y1 X^m +k ^ max.
m=1k=0
В случае, когда все волновые формы идентичны, т.е. Um = U = (u0,...,uq-1) для каждого m = 1,..., M, а их число M неизвестно, проблема обнаружения этих форм индуцирует следующую экстремальную задачу.
Задача 2, Дано: числовая последовательность Y = (y0,... ,Vn -1), вектор (u0,..., uq-1), натуральные числа Tmin и Tmax. Найти: набор n = (n1,..., nN) G Пм и его размерность, такие, что
M q-1
S(ni,.. .,пм) = ^ - 2yni+fc) ^ min. (3)
m=1k=0
Функции F и S аддитивны, поэтому обе задачи 1 и 2 решаются точно одним и тем же методом динамического программирования, но с использованием различных рекуррентных формул [11], Задача 1 решается за время O(MN2), а задача 2 — за время O(N2). Основную трудность представляет задача, которая индуцируется проблемой совместного обнаружения и оценивания повторяющейся формы при неизвестном числе повторов. Эта задача имеет следующую формулировку.
Задача 3, Дано: числовая последовательность Y=(y0,...,yN-1), натуральные числа q, Tmin и Tmax, Найти: набор n = (п1,..., nN) Е Пм и его размерность, такие, что
M M q-1
G(n1, . . . , nM) = Vnm+kVnj +k ^ max .
m=1j=1 k=0
Оптимальные значения компонент искомого набора U = (u0,... ,uq-1), соответствующего волновой форме, находятся по формуле:
uk = ¿ЕMM=1 ynm+k, k = 0,..., q — 1, где nm, m = 1,..., M и M — элементы оптимального решения задачи 3,
В общем случае эта задача NP-трудна [12], поэтому представляют интерес приближенные алгоритмы. Один из таких эвристических алгоритмов предложен в настоящей работе. Идея алгоритма состоит в следующем. Сначала находится решение задачи 1 при M = 1 для начального участка последовательноети Y, содержащего Tmax — q + 1 элементов, По найденному значению n1 находится набор (yni,... ,yni+q-1), Далее, используя этот набор, решаем задачу 2, положив U = (yni,... ,yni +q-1). Наконец, по найденному набору (n1,...,njM) вычисляем оценки компонент вектора U. С учетом обозначенного подхода предложен двухэтапный (локально-оптимальный) алгоритм нахождения оценки, а именно: на первом этапе определяется грубая оценка формы импульса, которая на следующем этапе уточняется в процессе решения задачи совместного оценивания формы импульсов и обнаружения моментов времени начала импульсов. При этом сущность первого этапа — решение задачи проверки гипотез, а сущность второго — решение задачи оценивания, которая сводится к минимизации аддитивного функционала (3),
Для проверки работоспособности и исследования точности работы предложенного алгоритма были выполнены численные эксперименты с моделированием различных волновых форм, одинаковых по форме и длительности, осложненных гауссовым шумом различного уровня. При этом задавались образцы реальных волновых форм, которые ранее были зарегистрированы от взрывных и вибрационных источников, и различные соотношения сигнал/шум. По сгенерированному набору (n1 ,...,nM) случайных номеров формировалась последовательность компонент вектора X, В соответствии с принятой моделью анализируемая последовательность компонент вектора Y синтезировалась как сумма вектора X и гауссовского вектора Е с параметрами распределения (0, а2/), В качестве примера на
Смоделированная последовательность импульсов в шумах
1.09 2.95 4.932 6.863 8.402 10.26 11.646 13.774 15.576 17.308 18.768 ""
Последовательность импульсов после применения алгоритма
1 096 2.948 5 026 6 88 8.412 10 258 11 646 13 778 15 576 17,4 18-774 сек
а)
Заданная и вычисленная формы импульса
сек
Ь)
Рис. 1. Соотношение сигнал/шум=1,25, Тт^п = 1,3 с, Ттах = 2,2 с, д = 1с; N = 20 с, М =11;
5и(М) = 6 • 10-2
рис. 1 приведены в графическом виде результаты совместного обнаружения и выделения волновых форм с помощью алгоритма решения задачи 2.
На рисунке изображены: а) сгенерированная модельная зашумлеппая последовательность и последовательность, найденная алгоритмом решения задачи 2; Ь) результаты численного оценивания погрешностей выделения одинаковых волновых форм в квазипериодической последовательности па фоне шума дня случая отношения еигпал/шум=1,25. Времена вступлений дня всех выделенных импульсов по отношению к обеим последовательностям проставлены па оси абцисс вначале каждого из импульсов. В серии численных экспериментов показано, что средняя абсолютная погрешность оценивания времени вступления волновой формы составляет 0,047 с, что в 3 раза меньше, чем у алгоритма вейвлет-фильтрации с пороговым обнаружителем, использованного для решения этой же задачи. Дня проверки качества алгоритма оценивания волновых форм использовалась мера среднеквадратического уклонения в виде ёц(М) = 1/д • Хд—о(ик — ик)2, где ик , ик, к = 0,..., ^ — 1 — заданные и вычисленные компоненты вол новой формы Ц Относительная средпеквадратическая погрешность оценивания волновой формы дня данных рис. 1 не превышает 6 %.
2. Фракталы в апостериорных алгоритмах. В задачах 13 и алгоритмах их решения параметры д,ТтгП и Ттах, соответствующие длительности волновой формы и границам снизу и сверху на интервал между двумя последовательными волновыми формами, являются входными данными. Однако в практических задачах эти параметры зачастую оказываются неизвестны заранее.
Для устранения этой априорной неопределенности авторами предложен подход для предварительного оценивания указанных параметров, основанный на фрактальном представлении волновых форм. Подход опирается на отображении волновых форм на двумерную плоскость „частота-время" с применением двумерного преобразования Фурье вида
1 N-1 , \
Р(к1 Л) = £ Р[П1,П2] ■ <П1<пз,ь,„ = ехр ( -¥ ) (4)
* 1 2 п1 =0 п2=0 ^ '
(4)
"
амплитудных значений будут соответствовать уровни яркостей. Получаемые таким образом изображения волновых форм служат для предварительного оценивания границ волновых импульсов. Последующее уточненное вычисление реализуется с помощью алгоритма дискретной оптимизации путем решения задачи 3,
Ниже приводятся результаты численного моделирования фрактального подхода к выделению границ волновых форм в шумах (рис, 2), Порядок моделирования состоял в следующем. Задавались образцы реальных волновых форм, взятых из экспериментов. Из набора выбирался образец, соответствующий конкретной задаче. Далее формировался кадр при различных значениях параметров N, М, Ттах, Ттгп и д в соответствии с (1), (2), На выбранные образцы накладывались шумы с гауссовым распределением с параметрами (0, ст). Выбором ст задавалось соотношение сигнал/шум.
Качественную картину сказанному дает рис, 2, В верхней части рис, 26 изображены за-шумленные последовательности волновых форм, подлежащие обработке. Запись содержит 8 волновых форм, вырезанных из реальных сейсмограмм, которые зарегистрированы от взрывных источников. На рис, 2а представлен результат двумерного Фурье-отображения записи согласно (4), Здесь уверенно выделяются из шумов начала и концы волновых импульсов, включая начало квазипериодической импульсной последовательности. Это способствует повышению надежности работы алгоритма дискретной оптимизации, В средней части рис, 26 представлены выделенные путем решения задачи 3 последовательности волновых форм без применения и с использованием фрактального подхода соответственно. Времена вступлений волновых форм представлены в начале каждого волнового импульса. На рис, 2с отображены в виде гистограмм погрешности вычисления времен вступлений для вариантов вычислений без использования фрактального представления — красные столбики — и с использованием его — синие столбики. Пунктирными рамками на рис, 26 выделены границы расположения волновых форм — исходных и вычисленных с применением фрактального подхода. Как следует из гистограмм, использование фрактального представления импульсной последовательности позволяет существенно повысить точность и надежность определения времен вступлений алгоритмом дискретной оптимизации, В ряде случаев наблюдается уменьшение погрешности на порядок [11].
3. Мониторинг положения скважинного источника сейсмических волн. Было рассмотрено приложение апостериорного алгоритма к решению одной из задач геофизического мониторинга, связанной с определением положения скважинного источника и
скоростей распространения сейсмических волн по глубине. Решение задачи имеет отношение к проблеме локации бурового инструмента в процессе нефтепромыслового бурения на основе определения времен вступлений сейсмических волн, порождаемых источником [13]. Ниже приводятся результаты модельного физического эксперимента, связанные с определением положения скважиппого источника (взрывов пороха массой 30 г), заглубляемого по вертикали в диапазоне 0-120 м.
Задача оценивания неизвестных параметров события сводится к решению нелинейной системы уравнений:
п = П(7,0) + e, (5)
где п = (П1,..., ПN)Т — вектор измеренных времен пробега сейсмических волн, n(Y, 0) = (Äi,..., ÄN)Т — N-мерный вектор вычисляемых времен пробега (теоретический годограф) или функция регрессии, е = (е1,..., eN)Т — вектор невязок, 0 = (x, y, z, v, t)T — ш-мерный вектор оцениваемых параметров, y = (y1, ... , yn) _ матрица координат датчиков, N — число датчиков, В качестве оцениваемых параметров выступают пространственные координаты источника — x, y, z, скоростная характеристика среды v и время в псточ-t
ei = fji(xi,0) — Пг(хг, 0). Будем в дальнейшем предполагать, что Eei = 0 Eeiej = а2^-, ai = a(xi), öij — символ Кронекера, i = 1,..., N.
0
решением служат оценки метода наименьших квадратов:
N
0 = arg min Q(0), Q(0) = ^ а"2 (Äi — n(Yi,0))2. (6)
öen i=i
Для отыскания минимума функционала Q(0) применяется итерационный метод Гаусса-Ньютона, основанный на линейной аппроксимации функции регрессии в окрестности точки 0 с применением на каждом шаге итерационного процесса методом псевдообращения (или обобщенного обращения), основанного на сингулярном разложении (SVD-разложении). Как известно, для вычисления SVD существует процедура SSVDC в библиотеке Linpaek [14], В [15] представлена стандартная процедура S FD-разложения на языке Фортран-IV, которая использовалась в данной работе.
Исходные данные для решения задачи в виде вектора n = (Ä1,..., ÄN)Т измеренных времен вступлений сейсмических волн определялись как результат решения задачи минимизации (3), В результате найдены оптимальный набор моментов времен вступлений
сейсмических волн и их количество: (n1,..., nм, M) = Arg min S1(n1,..., nM) и искомых
п
компонент волновой формы U : = m=1 ynm+fc, k = 0,... ,q — 1, Полученные времена были использованы в решении обратной задачи (5) с целью определения пространственных координат источника и скорости в среде. При определении скорости распространения волны по глубине используется метод „пристрелки" при лучевом трассировании в трехмерной слоисто-неоднородной среде [16], Результаты вычислений сведены в таблицу. Здесь последовательно по столбцам представлены глубины погружения источника, погрешности
xyz
волн. Данные таблицы иллюстрируют довольно высокую точность определения координат источника (погрешность по координате z на максимальных глубинах не превышает 1 %, отклонение по горизонтали не превышает 2 м).
с)
Рис.2. Результаты выделения волновых форм — с применением фрактального подхода (б) и без его использования (а), погрешности оценивания времен вступлений волн (с)
Таблица
Погрешности вычислений положения скважинного источника по глубине и вычисленные скоростные характеристики окружающей среды
Погрешность определения положения (м)
№ положения источника по глубине
X
у
Скорость сейсм. волн V (м/с)
Ьг1
Ьг5 ЬгЮ Ьг25 ЬгЮО
Ьг120
0,550 0,678 0,789 0,877 1,070 1,577
0,757 0,927 1,070 1,208 1,478 2,163
0,080 3,116 3,783 1,985 0,887 1,014
1950 1962 2740 2680 3770 3980
Заключение. 1. Для повышения точности решения задач обнаружения и выделения сейсмических волновых форм предложены и проанализированы апостериорные алгоритмы дискретной оптимизации, совмещающие в себе решение задач помехоустойчивого оценивания времен вступлений волн и выделения волновых форм в шумах. Высокая точность работы предложенных алгоритмов доказана на основе численных экспериментов, В частности, показано, что относительное ереднеквадратичеекое уклонение оценивания волновых форм не превышает 6 %, а относительные погрешности оценивания их времен вступления не хуже 0,1 %,
2, Рассмотрено применение предложенных алгоритмов к решению обратной задачи, связанной с определением положения скважинного источника и скоростей распространения сейсмических волн по глубине. Решение задачи имеет отношение к проблеме локации бурового инструмента в процессе нефтепромыслового бурения по сейсмическим данным, В модельном физическом эксперименте по определению положения скважинного источника (взрывов пороха массой 30 г), заглубляемого по вертикали в диапазоне 0-120 м, достигнута высокая точность определения координат источника (погрешность по координате г на максимальных глубинах не превышает 1 %, отклонение по горизонтали не превышает 2 м).
Оценены скорости распространения волн, нарастающие по глубине с 1950 до 3980 м/с. Список литературы
1. Сейсмическая томография. Под. ред. Г. Нолета. М.: Мир, 1990. С. 415.
2. Обнаружение изменения свойств сигналов и динамических систем. Под. ред. М. Бассвиль, А. Банвениста. М.: Мир, 1989. С. 278.
3. Никифоров И. В. Последовательное обнаружение изменения свойств временных рядов. М.: Наука, 1983.
4. Хайретдинов М. С., Омельченко О. К., Родионов Ю. И. Автоматизированная технология локации сейсмического источника / Тр. межд. конф. „Математические методы в геофизике". Новосибирск, 2003. Ч. II. С. 529-535.
5. Хайретдинов М. С., Воскобойникова Г. М., Седухина Г. Ф. Информационная технология сейсмолокации импульсных источников // Вестник национального ядерного центра Республики Казахстан. 2010. Вып. 3. С. 32-39.
6. Khairetdinov М. S., Avrorov S. A., Livenets A. A. Computing technology in seismic monitoring networks and systems // Bulletin of the Novosibirsk Computing Center. Series: Mathematical Modeling in Geophysics. 2010. I. 13. P. 51-69.
7. Любушин А. А. Вейвлет-агрегированный сигнал и синхронные всплески в задачах геофизического мониторинга и прогноза землетрясений // Физика Земли. 2000. № 3. С. 20-30.
8. Сагайдачная О. \!.. Дунаева К. А., Сальников А. С. Декомпозиция и анализ сейсмических полей на основе слоев вейвлет-разложения // Геофизика. 2010. № 5. С. 9-17.
9. Kel'manov А. V., Лкох В. A posteriori joint detection and discrimination of pulses in a quasiperiodic pulse train // IEEE Trans. Signal Processing. 2004. V. 52, N 3. P. 1-12.
10. Gruber P. and Todtli J. Estimation of Quasiperiodic Signal parameters by Means of Dynamic Signal Modes // IEEE Trans. Signal Processing. 1994. V. 42, N 3. P. 552-562.
11. Воскобойникова Г. M. Апостериорные вычислительные алгоритмы и программы в задачах геофизического мониторинга: дисс. канд. техн. наук / Сибирский государственный технический университет телекоммуникаций и информатики. Новосибирск, 2014. С. 119.
12. С11 mai.ii К. Кн.. Kel'manov A.V., Kel'Manova М.А., Khamidullin S. A. A Posteriori Detecting a Quasiperiodic Fragment in a Numerical Sequence // Pattern Recognition and Image Analysis. 2008. V. 18. N. 1. P. 30-42.
13. Illextмah Г. А. Определение параметров среды и траектории ствола скважины методом ВСП // Геофизика. 1996. № 5-6. С. 59-64.
14. Каханер Д., Моулер К., Нэш С. Численные методы и программное обеспечение. М.: Мир, 2001. С. 575.
15. Omelchenko О. К. Numerical implementation of wave mode of definition of bottom hole coordinates // Bulletin of the Novosibirsk Computing Center. Series: Mathematical Modeling in Geophysics. 1999. I. 5. P. 121-126.
16. Цецохо В. А., Виноградов С. П. О площадной пристрелке при лучевом трассировании в трехмерной слоисто-однородной среде / Препринт ИВМ и МГ СО РАН № 1093. 1997. С. 21.
Воскобойникова Гюльнара Маратовна, — науч. сотр. Института вычислительной математики и математической геофизики СО РАН,
тел. (383) 330-87-4-3; e-mail: gulya@opg.sscc.ru Хайретдинов Марат Саматович — д-р техн. наук, главн. науч. сотр. Института вычислительной математики и математической геофизики СО РАН, зав. кафедрой Новосибирского государственного технического университета;
тел. (383) 330-87-4-3; e-mail: marat@opg.sscc.ru
Дат,а, поступления — 24-06.2015
