Сейсмическая томография земли (стр. 3 )

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

Возможны и другие способы дискретизации. В работе Michelena, Harris [1991] предложено разбиение области с помощью ячеек, растянутых вдоль лучевых траекторий. Такой способ дискретизации был назван natural pixels. В работе Vasco et al. [1995] этот подход развивается с использованием объёмов Френеля, при этом учитывается физическая природа волны: зоны Френеля зависят от длины волны. Преимуществом подхода является его быстрота. В работе Michelena, Harris [1991] на модельном примере показано, что для восстановления модельного объекта с одинаковым качеством подход, использующий natural pixels, требует существенно меньше времени, чем стандартный, основанный на использовании прямоугольной сетки. Дискретизация с помощью natural pixels, или объёмов Френеля, не подходит для анизотропных сред: это связано с тем, что при её использовании необходимо учитывать влияние ячеек друг на друга в зоне пересечений лучей. В изотропной среде это влияние легко учитывается. Использование прямоугольных сеток является более универсальным и простым при программной реализации алгоритмов.

В работе T. Bodin, M. Sambridge [2009] предложен алгоритм решения задачи томографии, на каждой итерации которого модифицируется не только скоростная модель, но и сетка дискретизации. При этом используется нерегулярная сетка на основе многоугольников Вороного. Такое разбиение определяется только количеством ячеек и их «центральными» точками. Этот подход является оригинальным и интересным, но очень сложен в реализации и требует большого количества вычислений на каждом шаге итераций.

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

В настоящей работе будет рассматриваться стандартная дискретизация целевой области на основе прямоугольной сетки.

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

Можно ли применять подобную аппроксимацию, насколько она точна и как связаны между собой томографическая матрица и исходный оператор, какой должна быть дискретизация? Эти вопросы в большинстве работ исследуются методом обращения синтетических данных для различных моделей. Понятно, что такой подход не является достаточно строгим и универсальным. Существует ряд теоретических результатов, связывающих дискретизацию и разрешение [Наттерер, 1990], вывод которых основан на преобразовании Фурье. Эти теоремы и оценки носят общий характер, они не учитывают компактности томографического оператора, не предполагают наличие ошибок в данных и ошибок аппроксимации оператора.

Оператор томографии является компактным, поэтому обратный оператор не является ограниченным. Получающаяся в результате дискретизации томографическая система наследует свойства исходного бесконечномерного оператора. Томографическая матрица плохо обусловлена. Известно, что оптимальным выбором, позволяющим получать устойчивые решения плохо обусловленной задачи, является базис в пространстве решений, образованный правыми сингулярными векторами исходного линейного оператора [Чеверда, Костин, 2010]. Устойчиво определяются проекции на подпространства, соответствующие старшим правым сингулярным векторам. В научной литературе эти проекции называются r-решениями [Годунов и др., 1992; Чеверда, Костин, 2010].

Существует множество программных пакетов, позволяющих проводить различные операции линейной алгебры, в частности находить сингулярное разложение матриц. Общепринятым стандартом является библиотека LAPACK. Эта библиотека находится в свободном доступе (netlib. org) и написана на фортране. Она входит в широко распространённые математические пакеты, например, в MATLAB и математическую библиотеку компании INTEL MKL (Math Kernel Library). Существуют различные модификации LAPACK: SCALAPACK-LAPACK, переписанный на языке C, PROPACK - параллельная версия и другие. LAPACK позволяет вычислять сингулярное разложение для комплексных матриц двойной точности с использованием процедуры ZGESVD, комплексных матриц одинарной точности CGESVD, вещественных матриц двойной точности DGESVD и вещественных матриц одинарной точности с помощью процедуры SGESVD. В MATLAB этим процедурам соответствует команда svd.

Во все перечисленные процедуры исходные матрицы необходимо передавать в "полном" виде (двумерные массивы соответствующих типов). Процедуры вычисляют сразу весь спектр. Указанные особенности алгоритмов накладывают существенные ограничение на размеры матриц и, следовательно, на размеры сетки дискретизации и количество пар источников-приёмников. Вместе с исходной матрицей необходимо также хранить вспомогательные массивы, которые служат «рабочим пространством» для процедуры SGESVD, а также получаемые сингулярные векторы и числа. Кроме этого, вычисление сингулярного разложения с применением процедур LAPACK требует большого количества времени, так как вычисляется сразу весь спектр.

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

LAPACK не позволяется вычислять сингулярное разложение для спарсовых матриц. Нет такой возможности и в MKL. Для вычисления сингулярного спектра спарсовых матриц применяется пакет ARPACK. Пакет ARPACK (http://www. caam. rice. edu/software/ARPACK) предназначен для вычисления части собственных чисел и векторов квадратных матриц. В пакете есть различные процедуры для вычисления собственных чисел и векторов соответственно для симметричных и несимметричных комплексных и вещественных матриц разной точности. Доступна и параллельная версия ARPAСK. Задача поиска части сингулярного разложения матрицы B может быть сведена к поиску собственных чисел симметричной квадратной матрицы BB*. В пакете ARPACK есть готовые подпрограммы ssvd, dsvd для вычисления первых сингулярных чисел и векторов соответствующих матриц ("s","d" вещественные одинарной и двойной точности). Пакет ARPACK используется в MATLAB для вычисления части сингулярного разложения спарсовых матриц (команда ssvd). Существенным недостатком MATLAB является невозможность использования спарсовых матриц одинарной точности. Если работать с библиотекой ARPACK, непосредственно требуется задавать лишь действие матрицы на вектор и действие сопряжённой матрицы в виде процедур. Таким образом, хранить матрицу вообще не требуется. Тем не менее в ходе выполнения сингулярного разложения действие матрицы и сопряжённой матрицы считается большое число раз, поэтому на практике всё же приходится вычислить элементы матрицы один раз, а затем хранить их в спарсовом формате для действия прямой и сопряжённых матриц.

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

Классическим подходом является регуляризация [Тихонов, Арсенин, 1986]. При использовании этого типа регуляризации возникает сложность правильного выбора регуляризующего параметра.

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

Использование многомасштабных базисов [Woodward et al., 2008] также является одним из вариантов регуляризации. Решение задачи отыскивается итерационно с использованием последовательности измельчающихся прямоугольных сеток дискретизации целевой области. Решение, полученное на крупной сетке, используется в качестве начального приближения для решения задачи на более мелкой. При переходе с крупной сетки на мелкую используется сглаживание. Этот метод регуляризации по свойствам схож с предыдущим.