Анализ алгоритмов лучевой томографии для прогноза нарушенности выемочного столба Текст научной статьи по специальности «Математика»

-------------------------------------------- © В.Н. Захаров, А.П. Аверин,

С.А. Вартанов, 2010

УДК 550.832

В.Н. Захаров, А.П. Аверин, С.А. Вартанов

АНАЛИЗ АЛГОРИТМОВ ЛУЧЕВОЙ ТОМОГРАФИИ ДЛЯ ПРОГНОЗА НАРУШЕННОСТИ ВЫЕМОЧНОГО СТОЛБА*

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

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

Семинар № 3

ТЪ практике решения важнейшей

А# задачи шахтной сейсморазведки

- прогноза нарушенности угольных пластов - с середины 80-х годов 20 века начали активно применяться программноаппаратные комплексы на основе персональных ЭВМ для обработки данных полевых измерений, что позволило начать широкое использование одного из наиболее совершенных методов исследования внутренней структуры массива

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

Методы томографии можно разделить по следующим признакам:

- рабочая частота волн: ультразвуковая и сейсмоакустическая томография;

- лучевые траектории: прямолинейные и криволинейные;

- регистрируемые параметры волн: кинематические (скоростные) и динамические (амплитудные, частотные, энер-

гетические).

В угольном пласте отклонение упругих характеристик (модуля упругости, модуля сдвига) для сейсмоакустических волн обычно не превышает 20-30% от средних значений, поэтому в большинстве случаев обратная задача может быть поставлена в предположении прямолинейности лучевых траекторий. В тоже время, в зонах тектонических нарушений и замещений угольного пласта боковыми породами необходимо учитывать кривизну сейсмоакустических лучей для уточнения координат области с резким изменением скорости упругих волн.

В простейших, с точки зрения интерпретации результатов измерений по методу сейсмоакустической томографии, случаях достаточно использовать кинематические параметры, такие как скорость распространения волн различных типов (головных, типа Лява, Рэлея). Однако, как показывают результаты исследований [6, 7, 16], весьма перспективным направлением следует считать применение амплитудной сейсмотомографии, поэтому рассмотрим методы

*Работа выполнена при поддержке гранта РФФИ №»07-05-00718-а и №08-05-90437 - Укр-а

сейсмической томографии более подробно.

Используемые типы волн. Известно [2, 3, 4, 10, 22, 23, 24], что в угольном пласте распространяются волны следующих типов:

- объемные продольные, поперечные;

- головные;

- каналовые типа Лява и Рэлея.

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

В частности, для решения задачи прогноза нарушений в кровле угольного пласта применяются головные волны, а при анализе трещиноватости и смене литотипа пород в кровле - каналовые волны, обладающие дисперсией скорости от частоты. Объемные волны используются для оценки анизотропии свойств угольного пласта.

Сущность задачи сейсмоакустиче-ской томографии. Задача сейсмоаку-стической томографии сводится к решению системы уравнений для лучей, пересекающих N блоков массива под различными углами и различной длины. Рассмотрим простейшую систему уравнений:

=

Г—,і = 1,..., N Г C(?)

(1)

(1) используют их приближённые значения. Рассмотрим систему уравнений:

N Л?

іг =У-------Ц і = 1,..., N (2)

г 6 С(5у)’ ’ ’

Она получена из системы уравнений

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

Существуют два основных способа решения системы уравнений (1): решение системы интегральных уравнений и решение в первом приближении системы линейных алгебраических уравнений

(2).

Система уравнений (1) также может быть записана в виде:

или

Г^-,і = 1,...,М, і = 1,...,N (3) {С(*У

(4)

где і — номер источника, і — номер сейсмоприёмника, т( ?) — медленность

(величина, да(?) = 1

обратная скорости: ,).

где ь номер луча, ^ — полное время вдоль луча для первого вступления волны, С(а) — скорость распространения волны в массиве размерности [1 х N],

L — кривая, совпадающая с лучевой траекторией.

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

С {ау

Методы решения систем линейных алгебраических уравнений. Задачи сейсмоакустической томографии очень часто сводятся к системам линейных алгебраических уравнений (СЛАУ). Поэтому необходимо рассмотрение основных методов решения СЛАУ.

Методы исключения. Одним из наиболее ранних и известных методов последовательного исключения неизвестных является метод Г аусса [8].

Рассмотрим задачу: по известной

матрице

и вектору правой части

В

найти неизвестный вектор

||Х;.| ,т.е. АХ=В.

Идея метода Гаусса заключается в преобразовании матрицы ||Д^| к треугольному виду с элементами, расположенными выше ее диагонали, с помощью наибольшего ведущего делителя. При этом последнее уравнение имеет

элементарное решение Хы = Вы / Лтт .

Возвращаясь обратным ходом вверх, последовательно находим весь вектор

IIх, II.

Главное неудобство метода Гаусса заключается в невозможности решения СЛАУ при несоответствии числа уравнений числу неизвестных (то есть в том случае, когда матрица системы вырождена, т.е. det Л = 0), либо если в процессе преобразований получаем уравнение, в котором коэффициенты при всех неизвестных равны нулю, а свободный член не равен нулю.

Использование правила Крамера для решения такой системы крайне неэффективно и медленно (требуется найти N определителей матриц размерности N х N). Кроме того, определитель системы d должен быть отличен от нуля

( X: = —;, где

Л

А і =

а і—1

N, і—1

Ь1

то есть определитель матрицы, полученной из А заменой і-ого столбца на вектор правой части В, а А = det А).

Иной вывод правила Крамера связан с нахождением обратной матрицы А 1,

что дает решение X = Л 1В при det Л Ф 0 и вычисления идут быстрее.

Для решения СЛАУ с плохо обусловленной матрицей часто используется метод исключения типа ортогонализа-ции, предложенный А.А. Абрамовым [1]. При этом процесс исключения при фиксированной точности завершается преждевременно. Недостатком этого метода является невысокая скорость - она ниже скорости методов типа Гаусса и поэтому его следует применять в том случае, когда известно, что фактическое число шагов существенно меньше полного числа шагов.

Итерационные методы решения СЛАУ. При решении задач, в которых все величины, в том числе коэффициенты и свободные члены, известны лишь приближенно, методы Гаусса и обратной матрицы неприменимы, так как приводят к результату с плохой точностью. Однако многие задачи, в том числе обратные задачи сейсмоакустической томографии, являются некорректно поставленными задачами. Из трёх условий корректно поставленной задачи в них наиболее часто нарушается условие устойчивости, то есть соответствия малым ошибкам исходных данных малых ошибок решения. Поэтому для решения задач, сводящихся к СЛАУ, были разработаны итерационные методы, основанные на последовательном приближении неизвестных при решении некорректно поставленных задач [17-20].

Одним из таких итерационных методов является минимизация функционала Тихонова методом сопряженных градиентов [20]. В этом методе строится приближенное решение уравнения:

лг = и,

где Z и и определены в гильбертовом пространстве, т.е. с интегрируемым квадратом модуля, А- линейный либо нелинейный оператор. Кроме того, вво-

а

а

11

а

а

дится в рассмотрение сглаживающий функционал:

II II? II II?

ма[1 ] = |Л• 2-и5\в + а||г\\о, где

а > 0 - параметр регуляризации. Затем отыскивается минимум функционала методом сопряженных градиентов.

Рассмотрим функционал невязки

(р( z) = ІІА • z — и

|2

5 їїD

преобразуемый к

виду: р(г) = (г, Qz) + (d, г) + е , т.е. формируется матрицу Q = Л*Л и вектор d = 2 • Л*и8

Минимизация начинается с произвольной допустимой точки множества О и проводится в N - мерном пространстве Rn методом сопряженных градиентов [20]:

Вначале, для 1 итерации направление спуска вычисляется как

р(к) = —^айр(г(к)), для последующих итераций

||gradр( г(к

р(к} = —gradр(г(к}) + ^ У( )|[ 2 • р(к—».

||gradр( г(к ц)|

Определяется

1 (gradр( г(к)), р(к))

ак =

2

(Ор(к), р(к))

т.е. величина оптимального шага вдоль этого направления. Определяется атах -

величина максимально возможного ша-

~(к)

га вдоль найденного направления р , не выходящего за границу множества

а

Если ак < ат

то

г(к+1) = 7(к) + ак • р(к),

иначе

(к+1) (к) ,

1 + а„

• Р

(к)

и переход на метод проекции сопряженных градиентов. Останов итераци-

онного процесса происходит либо по невязке, либо по норме градиента.

Необходимо заметить, что в случае нелинейных задач обобщенный принцип невязки для выбора параметра регуляризации неприменим. Поэтому параметр регуляризации а выбирается левее точки разрыва согласно принципу сглаживающего функционала [20]. Способы выбора параметра регуляризации в нелинейных задачах отражены в работах А.С. Леонова [11-15], рассмотревшего случай с приближенно заданным оператором А.

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

Неявный итерационный метод для решения СЛАУ с комплексной матрицей [9]

Рассмотрим следующую СЛАУ:

(М + iQ)u = f, где и = иКе + т 1т,

1 = Ле +'/:т..

Комплексная матрица М + iQ имеет вещественную и мнимую компоненты, которые симметричны, и соответственно положительно и неотрицательно определены:

М = МТ > 0, Q = QT > 0.

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

ЛиКе = р, Л = ЛТ > 0.

Такая система уравнений решается быстро сходящимся методом с матрицей верхнего слоя М + а • Q с оптимальным параметром а . Норма начальной погрешности уменьшается в 1/ £ раз, 0 < £ < 1, за п = [0,57 • 1п(2/ £)] итераций. При этом число итераций огра-

ничено сверху равномерно по норме а = \Ммнимой компоненты мат-

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

(М + аQ) • ик+1 — ик + Л • ик = р,

^к+1

к = 0,1,...п —1, где и0 - начальное приближение.

Методика томографической интерпретации результатов измерений. Алгоритмы обработки сейсмограмм Обратное проецирование. Самый простой алгоритм, который широко применяется на практике. Алгебраическая запись его имеет вид:

Л = 1 •V Л • V

тп ^ у утп ’

Ь тп У

где Л.. — медленность или коэффициент

У

поглощения-рассеяния для у луча,

Ьтп =V Ьг]тп - сумма весовых коэф-

. 3

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

даннуЮ ячейку, . = Кг]тп / Яг] - весовой коэффициент для ячейки тп, Яупт — отрезок луча у , — длина луча

у.

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

Итерационное восстановление. Идея всех методов итерационного восстановления заключается в решении обратной задачи методом последовательных приближений. Последовательные приближения осуществляются по алгоритму А.Н.Тихонова [20]. Способ сведения интегрального уравнения к системе

линейных алгебраических уравнений описан выше. Наиболее распространённым из методов этого класса является метод одновременного итерационного восстановления. Рассмотрим данный метод подробнее.

Пусть в начале было получено решение методом обратного проецирования. Будем считать это нулевой итерацией:

Л'Чтп = Л'°тп = Л тп , где Лтп ПОЛДНЯ

на предыдущем шаге, а q = 0 . Рассчитаем для q-й итерации модельное среднее параметра по каждому лучу:

ЛЧ = V Л'Ч

У ~ о /—I № к1 ■

к ,1

После этого найдём расхождение между экспериментальным и расчётным лучевыми средними:

dЛЧ = Л.. — ЛЧ . Также введём среднеквадратическое расхождение для конкретной итерации:

Мч = N )2 , где N - количе-

ство лучей.

Если Мч устраивает нас, то итерационный процесс останавливается, и полученное значение Л'чтп является решением (на практике также часто используется правило останова, связанное с разностью расхождений на соседних итерациях: процесс продолжается до тех пор, пока Мч —Мч—1 не станет достаточно малой величиной).

В противном случае необходимо рассчитать сечение среднего расхождения:

dЛч = 1 У Лч • V

тп V /ши V .тп * тп

После этого можно найти следующее

приближение: Л'^ = Л'чтп +dA'Чmn, и

тп тп тп

перейти к следующей итерации.

Другие алгоритмы этого класса отличаются способами нахождения поправок dЛ'Чmn (например, с помощью метода Гаусса решения СЛАУ). Также отличие может заключаться в остановке процесса, например, существуют алгоритмы итерационного восстановления, основанные на минимизации функции

Мч .

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

Алгоритм суммирования разностей._Этот алгоритм представляет собой аналог свёртки обратной проекции (ОП) в алгебраическом представлении. Суть его заключается в коррекции сечения обратного проецирования. При этом предполагается уникальность (единственность) аномальных значений Л в

тп

каждой конкретной ячейке.

Рассмотрим, как и ранее, модельное лучевое среднее параметра:

Л у Ь у к1Л к1

к ,1

Определим для каждого луча лучевую сумму весовых коэффициентов по всем ячейкам за исключением ячейки тп:

У = V Ьгзк1 , и аналогичным образом

. , .

. Фт,

3 Ф п

рассчитаем лучевое среднее для всех ячеек, кроме тп:

л ~ = 1 V V л

утп ф. / ! ук1 к1

утп к ,1-к Ф т IФ п

Если предположить, что для каждого луча аномальные значения лучевого среднего связаны только с одной ячейкой, то можно записать величину этой аномалии как разность:

Л • V = Л • V — Л~ • V~

тп утп у у утп утп '

Если значительная часть рассматриваемого луча проходит вне рассматриваемого сечения (то есть V.. < 1), то эта

разность может оказаться достаточно большой из-за проходящей вне сечения части луча. В этом случае будем считать ^ =1.

Вычислим среднюю величину аномальных значений в ячейке тп по всем проходящим через неё лучам:

г]тп

Л" = 1 ~^(л Л~ • V V V

тп ^ \ у утп утп / ;

^ тп . 3

или, что эквивалентно,

Л'' = — е (Л.. - Л \.Б.. + V.. Л1 )£.. .

тп V-- У У Утп тп / утп

тп

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

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

Алгоритм коррекции ОП. Коррекция обратной проекции связана с обратной фильтрацией сечения [16, 21]. При этом само сечение обратной проекции рассматривается как результат фильтрации истинного сечения фильтром, к которому ищется обратный. После этого для восстановления истинного сечения производится обратная фильтрация найденного сечения. К недостатку метода относится зависимость от близости оценки сечения обратной проекции к действительному распределению свойств среды.

Заключение

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

1. Абрамов А.А. Об одном методе решения плохо обусловленных систем линейных алгебраических уравнений. //ЖВМ и МФ, 1991.Т.31, N 4.- С.483-491.

2. Азаров Н.Я. Шахтная сейсморазведка угольного пласта: Автореф. дис... докт. геол.-минер. наук / МГУ.- М., 1985.-36с.

3. Азаров Н.Я., Гильберштейн П.Г. Интерференционные волны, используемые при сейсмопросвечивании угольных пластов // Прикладная геофизика.-М.: Недра, 1978. -Вып.92. -С.42-57.

4. Аки К., Ричардс П. Количественная сейсмология. - М.: Мир,1983. -Т.1. -519 с.

5. Ефимова Е.А., Рудерман Е.Н. Возможности применения цифровой томографии для интерпретации геофизических данных. - М.: ВИЭМС, 1982, 55 с., ил.- (Обзор ВИЭМС).

6. Кабанихин С.И. Проекционноразностные методы определения коэффициентов гиперболических уравнений. -Новосибирск: Наука. Сиб. отд-ние, 1988. -С.131-143.

7. Колонин А.Г. Автоматическая обработка данных просвечивания // Передовой на-уч.-произв. опыт/ВИЭМС.-М., 1989.-Вып. 7.-С.3-6.

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

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

8. Курош А.Г. Курс высшей алгебры. -М: Наука. - 1975.-431с.

9. Кучеров А.Б. Быстрый итерационный метод решения в вещественной арифметике системы линейных уравнений с комплексной симметричной матрицей. //ДАН СССР, 1991. - Т.317, N 2.- С.294-296.

10. Левшин А.Л. Поверхностные и кана-ловые сейсмические волны . -М.: Наука, 1973.176с.

11. Леонов А.С. О построении устойчивых разностных схем решения нелинейных краевых задач // ДАН СССР.-1975.- Т.224, N 3.

- С.525-528.

12. Леонов А.С. Об алгоритмах приближенного решения нелинейных некорректных задач с приближенно заданным оператором // ДАН СССР. -1979.- Т.245, N2.- С.300-304.

13. Леонов А.С. О выборе параметра регуляризации для нелинейных некорректных задач с приближенно заданным оператором // ЖВМ и МФ. -1979.- Т.19, N 6. - С.1363-1376.

14. Леонов А.С. О связи метода обобщенной невязки и обобщенного принципа невязки для нелинейных некорректных задач // ЖВМ и МФ. - 1982. -Т.22, N 4. - С.783-790.

15. Леонов А.С. О некоторых алгоритмах решения некорректных экстремальных за-

дач //Мат. сб.- 1986. -Т. 129(171), N 2. - С.218-231.

16. Романов М.Е., Колонии А.Г. Криволи-

нейно-лучевая кинематическая и амплитудная сейсмотомография .-Новосибирск: ИМ СО

РАН, 1997.-40с.

17. Тихонов А.Н. О решении некорректно поставленных задач и методе регуляризации.// ДАН СССР, 1963.- Т.151, N 3.- С.501-504.

18. Тихонов А.Н. О регуляризации некорректно поставленных задач. // ДАН СССР, 1963.- Т.153, N 1.- С.545-548.

19. Тихонов А.Н., Арсенин В.Я. Методы решения некорректных задач. 2-е изд.- М.: Наука, 1979.- 288с.

20. Тихонов А.Н., Гончарский А.В., Степанов В.В., Ягола А.Г. Численные методы ре-

шения некорректных задач. - М.: Наука. Гл. ред. физ. - мат. лит., 1990.- 232 с.

21. Хермен Г. Восстановление изображений по проекциям: Основы реконструктивной томографии.- М.: Мир. 1983.-353 с.

22. Ямщиков В.С., Данилов В.Н., Вартанов А.З. О волнах Лява в тонком цилиндрическом волноводе в массиве // Изв. АН СССР. Физика Земли. - 1986.- N 2.- С.108-111.

23. Ямщиков В.С., Данилов В.Н., Вартанов А.З. Особенности характеристик канало-вых волн Лява в многослойном массиве //Изв. вузов. Горный журнал.-1986.- N 4. -С. 1-5.

24. Krey T. Channel waves as a tool of applied geophysics in coal mining // Geophysics, 1963, 28, N 6, p. 701-712. Ш

— Коротко об авторах ------------------------------------------------------

Захаров B.H. - доктор технических наук,

Аверин А.П. - кандидат технических наук,

Bартанов С.А. - УРАН Институт проблем комплексного освоения недр РАН, e-mail: averin.andrey@gmail.com

----------------------------------- ДИССЕРТАЦИИ

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

Автор Название работы Специальность Ученая степень

УЧРЕЖДЕНИЕ РОССИЙСКОЙ АКАДЕМИИ НАУК ГЕОФИЗИЧЕСКИЙ ЦЕНТР РАН

МИШИН Дмитрий Юрьевич Информационные геофизические модели и потоки данных в среде Г рид 25.GG.1G к.т.н.