Используется набор данных, который получил обозначение EHB98. Набор данных включает данные ISC (International Seismological Center) и NEIC (National Earthquake Information Center) бюллетеней, с переопределенными параметрами очагов землетрясений[Engdahl, E.R.,van der Hilst, R., Buland, R., 1998] в соответствии с моделью AK135 [Kennett, B. L. N., Engdahl, E. R., Buland R., 1995]. Из набора EHB98 авторы использовали почти 7 млн. отсчетов времени вступления Р - волн и 430 тыс. отсчетов времени вступления фаз рР – волн.
Для обнаружения вступлений коровых фаз применялся специальный согласующий фильтр, в котором цифровая сейсмограмма сворачивалась с функцией, задающей форму исходной волны (волноформой). Для вычисления отклонения времени вступления волны от стандартного годографа (модель АК135) применялся промежуточный параметр - "A differential travel time residual". При расчете этого разностного параметра используется одна сейсмограмма, на которой вычисляются временные сдвиги (которые и являются разностными параметрами) между вступлениями различных фаз, рефрагированных в ядре.
В работе приводится пример формулы расчета этого разностного параметра для двух фаз на одной сейсмограмме PKP_AB и PKP_DF. Физически, этот параметр является разностью отклонений каждой из этих фаз от стандартного годографа соответствующей фазы. Разностный параметр не содержит погрешности, с которой определяются такие параметры очага землетрясения, как время и координаты, так как они одинаковые для пары волн, выбранных для расчета разностного параметра. Так же разностный параметр свободен от влияния приповерхностных скоростных неоднородностей, в той части среды, где пути лучей пары волн совпадают с определенной точностью.
В работе не приводится информация о том, как этот разностный параметр пересчитывается в отклонение времени одной из пары волн, участвовавших в расчете параметра, то есть в или в . Технология обращения, которая используется в работе, описана в работах [Nolet, 1987, Spakman and Nolet, 1988]. Не учитывая регуляризацию, система нормальных уравнений, описывающая линеаризованную связь между параметрами модели (As) и данными (At) имеет вид:
LAs = At, где:L – матрица чувствительности, включающая производные Фреше, содержит длины лучевых сегментов в различных блоках и частные производные для источников; As – возмущения медленности относительно референтной модели ( АК135) и изменения в гипоцентре землетрясения (время в очаге, широта, долгота, глубина); At – отклонения времени.
Для решения этого уравнения применен итерационный LSQR алгоритм [Paige and Saunders, 1982, Nolet, 1985, van der Sluis and van der Vorst, 1987]. При обращении At определялись одновременно и структура среды, и уточненные значения положения источников. Минимизировалось значение дисперсии: ?2 = ?Т? , где ? = At – LAs – различия между экспериментальными данными и данными предсказанными на основе референтной модели. Было выполнено 200 итераций. Но примерно 98% сходимости было достигнуто после 25 итерационных циклов вычислений.
Количество исходных данных было значительно больше, чем параметров модели. Но большая часть модели была недоопределена из–за недостаточного покрытия блоков лучами. Так же данные были загружены шумом. Поэтому не представлялось возможным найти единственное решение без применения регуляризации. Применялось нормальное подавление, которое подавляет аномалии (разброс) в блоках с малым количеством данных и поэтому снижает амплитуду аномалий. Так же применялась градиентная (по первой производной) минимизация, которая обеспечивала сглаживание вариаций по широте и долготе. В радиальном направлении градиентная регуляризация не применялась.

Параметризация.
Модель (As) описывается линейной комбинацией ортогональных базисных функций, область определения которых охватывает изучаемое пространство. Многие подобные исследования в качестве базисных функций используют сферические гармоники [Su et al., 1994; Li and Romanowicz, 1996; Master et al., 1996]. В этом исследовании использовалось блоковое представление модели с непересекающимися блоками с постоянной скоростью, с заданными значениями широты, долготы и радиуса. [Inoue et al. 1990; Pulliam et al. 1993; Vasco et al. 1994, 1995a; Grand et al. 1997; van der Hilst et al, 1997; Vasco and Jonson, 1998; Bijwaard et al. 1998 ]. Локальный базис (блоковая параметризация) облегчает подавление артифактов в регионах с плохим покрытием и позволяет выполнять эффективную регионализацию. Но слишком большое количество параметров требует использование итеративных методов, которые, в свою очередь, требуют применения некоторых субъективных критериев сходимости и регуляризации и не позволяют рассчитать матрицу разрешения.
Использовалась регулярная трехмерная сеть, состоящая из 20 слоев, мощностью приблизительно 150 километров. Слои разделялись на 7200 ячеек размером 3⁰х3⁰ по широте и долготе.

Композитные лучи.
С целью уменьшения случайных погрешностей в определении времен вступления волн, содержащихся в бюллетенях МЦД, использовалось избыточное количество, почти 7х10₆ EHB98 P+ residuals и рассчитывались (конструировались) композитные лучи по избыточному набору данных [Spakman and Nolet, 1988].
Для построения композитных лучей землетрясения группировались в мантийные объемы размером 1⁰х1⁰х75км, а станции объединялись в приемные кластеры. Композиционные лучи соединяли мантийные объемы и приемные кластеры. Пути всех отдельных лучей, соединяющих конкретные источники из заданного мантийного объема и конкретные сейсмические станции, принадлежащие к заданному кластеру, добавлялись к одной и той же строке матрицы L. Из таких пучков лучей и складывались композитные лучи. Если в так называемых суммарных лучах объединяются лучи от конкретного землетрясения для всех станций приемного кластера, то есть с одной стороны, со стороны приемников, то в композитных лучах объединение производится с двух сторон, и со стороны источников, принадлежащих мантийному объему, и со стороны приемников из одного приемного кластера. Было показано, что применение этих различных подходов к объединению лучей дает немного различающиеся результаты [Widiyantoro, 1997].
Количество рядовых лучей, объединяемых в композитные лучи, доходило до тысячи и более. Но в среднем, эта цифра не превысила четырех. Для уменьшения эффекта разброса значений в пучке лучей, включаемых в композитный луч, использовалось медианное значение отклонений времени. Для того чтобы не исказить распределение значений, данные из выборки не исключались до определения медианного значения. Композитные лучи с медианными значениями, превышающими 5 секунд, исключались из дальнейшей обработки (таких данных было менее 1% от общего объема). Каждая строка матриц масштабировалась (scale) таким образом, чтобы общая длина пути равнялась корню квадратному из суммы всех длин путей всех лучей, вошедших в композитный луч. Соответствующим образом взвешивалось и отклонение (стоящее справа от знака равенства в линейном уравнении). Это грубый эквивалент использования индивидуальных лучей, но с повторением медиан отклонений в инверсии наименьших квадратов. Однако при превышении количества лучей в пучке 10 преимущества от добавления лучей уменьшаются по причине таких эффектов, как смещения положения источников, инструментальных эффектов, неоднородностей в районе приемников, несовершенством выбранной теоретической аппроксимации. величина этих эффектов превышает случайной ошибки определения отклонений времени. Поэтому мы применяем постоянный вес ко всем композитным лучам с 10 лучами и более. Несмотря на то, что только 10% композитных лучей имеют больший вес, Эти лучи включают в себя более 50% всех данных. Сужение диапазона изменчивости (vr) определяется по формуле:
VR = 1- Var(?)/Var(At)
Вычисляется в процентах и показывает насколько хорошо рассчитанная модель среды объясняет исходные данные. Наилучшему объяснению соответствует значение этого параметра, равное 100%. Расчеты показали, что параметр VR находится в диапазоне от 30% до 63%, в зависимости от количества лучей, включаемых в композитные лучи.

Зоны Френеля (Fresnel Zones) и чувствительные ядра (Sensitivity Kernels)
Пространство вокруг пути геометрического луча между источником и приемником, которое вносит основной вклад в формирование сейсмограммы, связано с зоной Френеля. В однородных средах максимальная ширина (?m) зоны Френеля может быть определена в лучевом приближении по формуле [Nolet, G.,1987]:
,
где: ? – длина волны и S – длина пути от источника до точки наблюдения. Зоны Френеля не получили широкого применения в кинематической томографии (travel time tomography). Как правило, вместо зон Френеля рассчитываются производные Фреше (то есть элементы матрицы L) путем простого трассирования лучей и применяется затухание или наложение некоторых пространственно корреляционных длин для получения гладких результатов.
Для Р – волн с частотой ~1 Hz, которые достигают границы с ядром, зона Френеля примерно равна 200 км. Это немного больше размеров ячеек, но меньше размеров структур, которые предполагается изучить.
Для интерпретации PKP_DF - P_diff данных теория лучевого приближения не дает удовлетворительных результатов. Это заключение следует из того обстоятельства, что записи с центральной частотой ~ 50 mHz имеют зону Френеля порядка 1000 км для волн, которые только коснулись границы ядра Земли и значительно большую зону Френеля для волн, дифрагируемых вокруг ядра. К тому же, P_diff не является фазой, описываемой в лучевом приближении. Это исчезающие, дисперсионные волны от границы мантии с ядром, такие, что чувствительность их годографов к неоднородной структуре частотно-зависимая. Для подобных волн рассчитываются частотно-зависимые чувствительные ядра [Zhao and Jordan, 1998; Zhao et al, 2000]
Сравнение результатов обращения только по P, pP и pwP данным, и обращения по этим же данным, но дополненным PKP_AB - PKP_DF , PKP_AB - PKP_BC и PKP_DF - P_diff данными, показало, что улучшилось разрешение для горизонтов, залегающих ниже 2200 километров от поверхности Земли.

Литература.
Nolet, G., Seismic wave propation and seismic tomography, in Seismic Tomography, With Application, in Global Seismology and Exploration Geophysics, edited by G. Nolet, pp.1-27, D. Reidel, Norwell, Mass., 1987.
Spakman, W., and Nolet G. Imaging algorithms, accuracy and resolution in delay time tomography, in Mathematical Geophysics: A survey of Recent Developments in Seismology and Geodynamics, edited by N. J. Vlaar, pp 155-188, D.Reidel, Norwell, Mass., 1988.
Paige, C. C., and Saunders, M. A., LSQR: An algorithm for sparse linear equations and sparse least squares, TOMS 8(1), 43-71 (1982).
Paige C. C., and Saunders M. A., Algorithm 583; LSQR: Sparse linear equations and least-squares problems, TOMS 8(2), 195-209 (1982).
Nolet, G., Solving or resolving inadequate and noisy tomographic systems, J. Comput. Phys, 61, 463-468, 1985.
van der Sluis, A. and van der Vorst H. A., Numerical solution of large sparse linear systems arising from tomographic problems, in Seismic Tomography, With Application in Global Seismology and Exploration Geophysics, edited by G. Nolet, pp. 53-87, D. Reidel, Norwell, Mass., 1987.

Используются сферические гармоники для параметризации модели:
Su, W.J., Woodward R., and Dziewonski A. M., Degree-12 model of shear velocity heterogeneity in the mantle, J. Geophys. Res., 99, 6945-6980, 1994.
Li, C. D., and Romanowicw, B., Global mantle shear velocity model developed using non-linear asymptotic coupling theory, J. Geophys. Res., 101, 22,245-22,272, 1996.
Masters, G., Laske, G. and Bolton, H., A shear-velocity model of the mantle, Philos. Trans. R. Soc. London, Ser. A, 354, 1385-1411, 1996.

Используются блоки для параметризации модели:
Inoue, H., Fukao Y., Tanabe K., and Ogata Y., Whole mantle P-wave travel time tomography, Phys. Earth Planet. Inter., 59, 294-328, 1990.
Vasco, D. W., Johnson, L. R., Pulliam, R. J. and Earl, P.S., Robust inversion of iasp91 travel time residuals for mantle P and S velocity structure, earthquake mislocations, and station corrections, J. Geophys. Res., 99, 13,727-13,755,1994.
Vasco, D. W., Johnson, L. R., Pulliam, R. J., Lateral variation in mantle velocity structure discontinuities determinated from P, PP, S, SS, and SS-SdS travel time residuals, J. Geophys. Res., 100, 24,037-24,059, 1995.
Grand, S.P., van der Hilst, R.D., and Widiyantoro S., Global seismic tomography: A snapshot of convection in the earth, GSA Today, 7, 1-7, 1997.
van der Hilst, R.D., Widiyantoro S. and Engdahl, E., Evidence for deep mantle circulation from global tomography, Nature, 386, 578-584, 1997.
Vasco, D. W., and Johnson, L. R., Whole earth structure estimated from seismic arrival times, J. Geophys. Res., 103, 2633-2671, 1998.
Bijwaard, H., Spakman, W., and Engdahl, E. R., Closing the gap between regional and global travel time tomography, J. Geophys. Res., 103, 30,055-30,079, 1998.

Композитные лучи.
Spakman, W., and Nolet G., Imaging algoritms, accuracy and resolution in delay time tomography, in Mathematical Geophysics: A survey of Recent Developments in Seismology and Geodynamics, edited by N. J. Vlaar, pp 155-188, D.Reidel, Norwell, Mass., 1988

Сопоставление "композитных" и "суммарных" лучей.
Widiyantoro, S., Studies of seismic tomography on regional and global scale, Ph.D thesis, 256 pp., Aust. Nat. Univ., Canberra, 1997.

Frequency-dependent sensitivity kernels.
Zhao, L., and Jordan T. H., Sensitivity of frequency-dependent traveltimes to laterally heterogeneous, anisotropic earth structure, Geophys. J. Int., 133, 683-704, 1998.
Zhao, L., Jordan T. H., and Chapman C. H., Three-dimensional Frechet differential kernels for the seismic delay times., Geophys. J. Int., 141, 558-576, 2000.

3D Earth's Models Player