CC BY
51
УДК 551.2+528.9:004+519.876.5 Б.Т. Мазуров, В.К. Панкрушин СГГ А, Новосибирск
ВЫЧИСЛИТЕЛЬНЫЙ ИДЕНТИФИКАЦИОННЫЙ ЭКСПЕРИМЕНТ: СОВМЕСТНАЯ МАТЕМАТИЧЕСКАЯ ОБРАБОТКА И ИНТЕРПРЕТАЦИЯ НИВЕЛИРНЫХ И ГРАВИМЕТРИЧЕСКИХ НАБЛЮДЕНИЙ ЗА ДИНАМИКОЙ ЗЕМНОЙ ПОВЕРХНОСТИ И ГРАВИТАЦИОННОГО ПОЛЯ В ВУЛКАНИЧЕСКОЙ ОБЛАСТИ
Настоящая работа является продолжением работ [1, 2] и описывает вычислительный идентификационный эксперимент по совместной математической обработке и интерпретации смоделированных в работе [2] пространственно-временных рядов комплексных наблюдений нивелирных превышений и абсолютных значений ускорения силы тяжести. Исследуемой геодинамической системой являются вертикальные движения земной поверхности и гравитационное поле в районе действующего вулкана [1].
Математическая обработка выполнялась с использованием алгоритма рекуррентного фильтра Калмана-Бьюси (ФКБ) [3] аналогично тому, как это описано в работе [4].
Для первой эпохи (до извержения) расширенный вектор параметров будет состоять только из отметок мобильных пунктов Им\, Нм 2, Нм з и массы конуса вулкана МК:
Хя( = 1) = (Нм 1, Нм2, Нм3,МК)Т.
Матрица Л(1) коэффициентов уравнений наблюдений формируется известным образом: каждое измерение у- (г) выражается как функция от параметров (г). Частные производные этих функций будут коэффициентами
матрицы Л(1). Покажем это на примере превышения /^(г = 1) . С учетом схемы сети (рис. 1 в [1]) это превышение выражаем как функцию от вектора параметров:
/1(г = 1) = Нм2(1 = 1) - На + (г = 1)мк(г = 1) + 8/, (1)
где Нм2 (г = 1) - отметка пункта М2 в эпоху 1; Н^\ - отметка исходного (стабильного) пункта С1; мК(г = 1) - масса маскона, моделирующего конус вулкана в эпоху 1; 8/^ - ошибка измерения превышения между пунктами С1 и
М2; кмк (г = 1) - коэффициент, учитывающий влияние этого маскона на "1
превышение /_(г = 1). Величина кмК (г = 1) вычисляется по известным
"1
координатам начального, конечного пунктов линии нивелирования и координатам центра шарового маскона Sk1.
Пример уравнения наблюдений абсолютного значения силы тяжести g -уравнение наблюдений g на пункте С1 в эпоху ? = 1:
§а(1 = 1) = Г+Д?”+С (г = 1)мк (г = 1)+8&_. (2)
В этом уравнении 8^ - ошибка измерения силы тяжести на пункте С1;
коэффициент км (г = 1) при мК (г = 1) вычисляем по известным координатам
пункта С1 линии нивелирования и координатам центра шарового маскона Sk1.
В эпохи 2, 3 и 4 в уравнениях наблюдений присутствуют слагаемые, учитывающие влияние не только увеличивающейся массы конуса вулкана (маскон МК) но и уменьшение массы в магматическом очаге (маскон МО), а также изменение высоты мобильных пунктов (через вертикальный градиент).
Координаты центра шарового маскона Sk1 достаточно знать приближенно (не грубее 50 метров). Для их определения используем результаты геодезических измерений размеров и формы конуса вулкана. Например, в работе [5] описана методика и результаты определения расхода магмы и размеров конуса при извержении Толбачинского вулкана.
Соотношение весов превышений и измерений силы тяжести устанавливалось с учетом точности измерений: ц/ = 0.5мм / км, = 5мкгал .
Веса превышений р/ = —, где - длина хода в км (табл. 1 в [2]). Вес
1 ^
измерения абсолютного значения силы тяжести для всех пунктов сети р& =001.
Через сутки после наблюдений первой эпохи началось извержение вулкана. При этом происходят одновременно вертикальные движения мобильных пунктов М1, М2, М3, уменьшение массы магматической камеры на величину МО и увеличение массы конуса мК на величину 8мК. Соответственно меняется и гравитационное поле.
Для создания модели вертикальных движений мобильных пунктов и изменения гравитационного поля на протяжении всего периода наблюдений, мы воспользуемся понятием ступенчатой функции [3] - функции, которая меняет свои значения только в дискретной последовательности точек разрыва. Это означает, что для каждого интервала времени между соседними эпохами принимается свое значение смещения или изменения массы, независимое от смещения или изменения массы в другие интервалы времени.
Расширенный вектор параметров для эпохи 2 будет:
ХЯ1 = (Нм 1 > Нм 2 > Нм 3 > им1 > им 2 > им 3 > мк 8мк >мо)Т,
где им1,им2,имз - вертикальные смещения мобильных пунктов между
эпохами 1 и 2.
Масса конуса увеличивается в результате излияния части вещества из магматической камеры. Определение характеристик процесса излияния является одной из важнейших задач при изучении вулканизма. Заметим, что адекватная математическая модель этого процесса будет работать на уточнение параметров, которые мы оцениваем по геодезическим и гравиметрическим наблюдениям. Однако, параметры экспоненциального закона можно определить только после нескольких эпох наблюдений, причем приближениями - модель должна настраиваться. Одновременно должна настраиваться еще одна
важнейшая величина - процент оставшегося в конусе вулкана и на его поверхности вещества, т.к. его часть выбрасывается в атмосферу.
Пока мы примем, что в конусе остается 70% от излившегося из магматического резервуара вещества. Неточность априорного задания этой величины отразится на поправках в результаты измерений после математической обработки - они резко возрастут.
В качестве начальных оценок нм1,нм2,нм3 и массы конуса вулкана мК и их ковариаций мы используем оценки, полученные при обработке наблюдений первой эпохи. Начальные оценки им1,им2,им3 и мО найдем по четырем необходимым измерениям второй эпохи /1, /2, /7, ЯС1. Используем при этом уравнения наблюдений, аналогичные (1) и (2).
Затем обрабатываем результаты наблюдений третьей и четвертой эпох. Заметим, что при обработке третьей эпохи вектор параметров расширится за счет включения параметра 8мО, который останется и в составе вектора параметров четвертой эпохи:
ХЯ1 = (Нм1 , Нм2 , Нм3 , им 1 , им2 , им3 , мк> 8мк,м0,8мО)Т.
Напомним, что мы моделируем изменения отметок и масс масконов с использованием ступенчатых функций - допускаем независимость этих величин в разные интервалы времени. Поэтому и в третью эпоху и в четвертую вычисляем вначале по необходимым измерениям начальные оценки им1,им2,имз, 8мК, 8м0. За исходные оценки параметров
Нм 1, Нм2, Нм3, мК, мО берутся оценки, полученные по наблюдениям предыдущей эпохи.
Полученные оценки отметок мобильных пунктов, их вертикальных смещений, масс масконов и их изменений от эпохи к эпохе позволяют выдвигать конкурирующие гипотезы закономерностей движений и изменений масс. Замечаем, что оценки смещений им1, им2, имз примерно равны в разные эпохи -отметки мобильных пунктов меняются с постоянной скоростью. Сложнее подобрать закон изменения во времени масс мК, мО. Увеличение массы конуса 8мК от эпохи к эпохе будет снижаться. Причиной этому является падение со временем избыточного давления в магматическом очаге, и как следствие, уменьшение излияния магмы, пропорционального этому давлению. Закон расхода магмы при извержении близок к экспоненциальному [5]. Для определения закона расхода магмы воспользуемся оценками параметров, вычисленных по модели со ступенчатыми функциями. Расчет количества магмы, извергнутой ко времени т с
начала извержения выполняется по формуле мизл (т) = ДВ(1 — е ^°т / Д ), где W0 - начальный расход магмы; ДВ - избыток магмы в очаге. Для трех эпох будем иметь три таких нелинейных уравнения с двумя неизвестными ДВ и Ж0. Мы определили эти неизвестные графическим способом. Избыток магмы в очаге ДВ = 1.12 * 107 тонн. Начальный расход магмы
Ж0 = 0.90 *107 тонн / сутки = 1042тонн / секунда.
Экспоненциальный закон излияния магмы с этими параметрами мы использовали для составления переходных матриц при обработке наблюдений с другими тремя моделями (произвольного смещения, движений жестких блоков и вращения жесткой плиты). Необходимо учесть, что в общем случае закон изменения параметров от эпохи к эпохе может быть различным. В нашем примере это касается масс масконов - увеличение этих масс от эпохи к эпохе неодинаково для разных пар соседних эпох. Сложно определить 8мК по оценке мО - необходимо знать процент оставшегося на конусе вещества в каждую эпоху. Эти величины могут быть определены по оценкам параметров модели 1. Но мы выполнили адаптацию (настройку) этих величин по критерию шт ггК^к (X, г) с использованием фильтра Калмана-Бьюси. В результате
установили, что в эпоху 2 и 3 доля оставшегося на конусе вещества составляла 70%, а в эпоху 4 - 90%.
Эта модель перераспределения масс масконов, аппроксимирующая динамику изменения гравитационного поля, сочеталась с тремя моделями вертикальных движений земной поверхности. Если принимать гипотезу произвольного смещения пунктов (модель 2), то расширенный вектор параметров будет:
ХЯ1 = (Нм1, Нм2, Нм3, им 1, им2, им3, мк, 8мК)Т,
где им1,им2,имз - вертикальные смещения мобильных пунктов между
соседними эпохами.
Третья модель соответствует ситуации, когда пункты С1 и С2 расположены на неподвижном жестком блоке, пункт М1 расположен на подвижном жестком блоке Л, а пункты М2 и М3 расположены на другом подвижном жестком блоке В. Расширенный вектор параметров будет:
ХЯ2 = (HM1, Нм2, Нм3 , иА, иВ, мК,8мК)Т,
где иА,иВ - вертикальные смещения блоков между соседними эпохами.
Согласно гипотезе 4 вертикальные движения мобильных пунктов М1, М2, М3 можно объяснить вращением жесткой плиты вокруг оси параллельной линии М2М3. Соотношение полученных по предыдущим моделям оценок смещений пункта М1 и пунктов М2, М3 позволяет принять при моделировании расположение оси вращения на расстоянии 0.50 км от М1 и 0.95 км от линии М2М3.
Учитывая малость угла поворота cоДt, связь отметок пунктов двух эпох обеспечивается введением угловой скорости вращения плиты о (в рад/год). Расширенный вектор параметров в этом случае будет:
Хя3 = (Нм1, Нм2, Нм3, о, мк, 8мк)т .
При обработке наблюдений по всем конкурирующим моделям получены оценки параметров состояния и их точностные характеристики. Приведем в табл. 1 те из них, которые отражают динамику земной поверхности и гравитационного поля - скорости вертикальных движений, массы масконов и их изменения для четырех эпох наблюдений по трем конкурирующим моделям.
Используя F-критерий аналогично тому, как это описано в [4], была выбрана адекватной модель 4. Согласно этой модели динамика земной поверхности описывается равномерным вращением блока с мобильными пунктами, а динамика гравитационного поля - изменениями масс поверхностного и глубинного масконов в соответствии с экспоненциальным законом излияния магмы и значениями доли остающегося на поверхности конуса излившегося вещества в эпохи 2 и 3 - 70%, в эпоху 4 - 90%.
По оценкам масс этих масконов и ковариационным матрицам векторов параметров в каждую эпоху для каждого пункта нами были вычислены оценки силы тяжести g, компонент уклонений отвесной линии £ 1 и аномалия высоты ^на пунктах геодезической сети с оценкой их точности (табл. 2).
XR t = 2 t = 3 t = 4
Модельные Оценки Точность Модельные Оценки Точность Модельные Оценки Точность
Модель произвольного смещения пунктов и экспоненциального закона излияния магмы
иМ 1 -16.23 мм -16.09 мм 0.89 мм -16.23 мм -16.09 мм 0.39 мм -16.23 мм -16.13 мм 0.89 мм
иМ 2 30.84 мм 30.86 мм 1.07 мм 30.84 мм 30.86 мм 0.47 мм 30.84 мм 30.82 мм 1.07 мм
иМ 3 30.84 мм 30.05 мм 1.05 мм 30.84 мм 30.05 мм 0.46 мм 30.84 мм 31.09 мм 1.05 мм
МК 4.39*109 тонн 4.38*109 тонн 4.90*106 тонн 6.20*109 тонн 6.19*109 тонн 3.55*106 тонн 4.39*109 тонн 7.74*109 тонн 2.94*106 тонн
8МК 3.01*109 тонн 3.02*109 тонн 4.90*106 тонн 1.81*109 тонн 1.81*109 тонн 1.33*106 тонн 3.01*109 тонн 1.55*109 тонн 0.72*106 тонн
Модель движения жестких блоков и экспоненциального закона излияния магмы
иМ 1 -16.23 мм -16.12 мм 0.85 мм -16.23 мм -15.66 мм 0.38 -16.23 мм -16.13 мм 0.23 мм
иМ 2 = иМ3 30.84 мм 30.43 мм 0.89 мм 30.84 мм 31.16 мм 0.39 30.84 мм 30.96 мм 0.24 мм
МК 4.39*109 тонн 4.38*109 тонн 4.71*106 тонн 6.20*109 тонн 6.19*109 тонн 3.47*106 тонн 4.39*109 тонн 7.74*109 тонн 2.92*106 тонн
5МК 3.01*109 тонн 3.02*109 тонн 4.71*106 тонн 1.81*109 тонн 1.81*109 тонн 1.30*106 тонн 3.01*109 тонн 1.55*109 тонн 0.71*106 тонн
Модель вращения жесткой плиты и экспоненциального закона излияния магмы
ш 1.95*10-5 рад/год 1.92*10-5 рад/год 3.05*10-7 рад/год 1.95*10-5 рад/год 1.94*10-5 рад/год 1.38*10-7 рад/год 1.95*10-5 рад/год 1.95*10-5 рад/год 0.82*10-7 рад/год
иМ 1 -16.23 мм -16.04 мм 0.25 мм -16.23 мм -16.16 мм 0.12 мм -16.23 мм -16.24 мм 0.07 мм
иМ 2 = иМ3 30.84 мм 30.48 мм 0.48 мм 30.84 мм 30.70 мм 0.22 мм 30.84 мм 30.86 мм 0.13 мм
МК 4.39*109 тонн 4.38*109 тонн 4.60*106 тонн 6.20*109 тонн 6.18*109 тонн 3.50*106 тонн 4.39*109 тонн 7.74*109 тонн 2.92 *106тонн
5МК 3.01*109 тонн 3.02*109 тонн 4.60*106 тонн 1.81*109 тонн 1.81*109 тонн 1.31*106 тонн 3.01*109 тонн 1.55*109 тонн 0.71*106 тонн
Таблица 2. Значения основных трансформант гравитационного поля, (модельные и полученные в результате обработки
по ФКБ)
Пункт gz (мкгал) Л С(см) gг (мкгал) л" С(см) gz (мкгал) л" С(см)
і = 2 і = 3 і = 4
Модель 980881768.0 5.3*10-2 0.0 0.92 980881531.5 5.6*10-2 0.0 1.38 980881354.7 7.2*10-2 0.0 1.72
Сі Оценка 980881767.4 5.2*10-2 0.0 0.92 980881529.8 5.6*10-2 0.0 1.38 980881355.4 7.2*10-2 0.0 1.72
Точность 0.37 3.7*10-4 0.0 1.2*10-3 0.70 3.3*10-4 0.0 1.2*10-3 0.72 3.0*10-4 0.0 1.1*10-3
Модель 980881821.5 3.6*10-2 -2.6*10-2 0.88 980881612.8 3.7*10-2 -2.6*10-2 1.31 980881466.8 4.7*10-2 -3.4*10-2 1.63
С 2 Оценка 980881821.0 3.6*10-2 -2.6*10-2 0.88 980881614.5 3.6*10-2 -2.6*10-2 1.31 980881463.2 4.7*10-2 -3.4*10-2 1.64
Точность 0.33 2.8*10-4 2.0*10-4 1.2*10-3 0.61 2.5*10-4 1.8*10-4 1.1*10-3 0.63 2.3*10-4 1.6*10-4 1.1*10-3
Модель 980880797.6 2.5*10-1 -1.2*10-1 1.46 980880027,0 3.2*10-1 -1.5*10-1 2.15 980879449.3 4.0*10-1 -1.9*10-1 2.67
М1 Оценка 980880796.4 2.5*10-1 -1.2*10-1 1.46 980880027.6 3.2*10-1 -1.5*10-1 2.15 980879453.7 4.0*10-1 -1.9*10-1 2.67
Точность 0.86 7.7*10-4 3.6*10-4 1.9*10-3 2.00 6.0*10-4 2.8*10-4 1.8*10-3 2.05 5.2*10-4 2.4*10-4 1.6*10-3
Модель 980878007.6 1.5 0.0 2.35 980875727,8 2.0 0.0 3.39 980873902.7 2.4 0.0 4.18
М2 Оценка 980878008.5 1.5 0.0 2.35 980875723.0 2.0 0.0 3.39 980873897.3 2.4 0.0 4.18
Точность 0.31 2.2*10-3 0.0 3.0*10-3 4.00 1.4*10-3 0.0 2.4*10-3 1.55 1.1*10-3 0.0 1.2*10-3
Модель 980880687.9 1.8*10-1 -3.0*10-1 1.53 980879836.5 2.4*10-1 -4.0*10-1 2.25 980879208.2 3.0*10-1 -5.0*10-1 2.80
М3 Оценка 980880686.6 1.8*10-1 -3.0*10-1 1.53 980879839.6 2.4*10-1 -4.0*10-1 2.25 980879201.0 3.0*10-1 -5.0*10-1 2.80
Точность 0.98 4.8*10-3 8.0*10-4 2.0*10-3 2.16 3.6*10-4 6.0*10-4 1.8*10-3 2.21 3.1*10-4 5.1*10-4 1.7*10-3
Заключение. Результатом совместной математической обработки и интерпретации смоделированных пространственно-временных рядов комплексных геодезических и геофизических наблюдений по алгоритму ФКБ явились оценки отметок мобильных пунктов, оценки вертикальных смещений этих пунктов, оценки аномалий силы тяжести, компонент уклонений отвесной линии и аномалия высоты на пунктах геодезической сети с оценкой их точности в 4-х эпохах. Совместная обработка нивелирных и гравиметрических наблюдений с использованием рекуррентного ФКБ позволила при адекватных моделях закономерностей движений земной поверхности и изменений масс поверхностного и глубинного масконов в соответствии с экспоненциальным законом и высокой точности наблюдений нивелирных превышений значительно повысить точность оценивания. Выполнен вычислительный эксперимент: параметрическая и структурная идентификация динамики земной поверхности и гравитационного поля в вулканической области.
БИБЛИОГРАФИЧЕСКИЙ СПИСОК
1. Мазуров Б.Т., Панкрушин В.К. Идентификационный эксперимент: построение физико-математической модели динамики земной поверхности и гравитационного поля в вулканической области. Сб. материалов Международной выставки и научного конгресса «ГЕО-СИБИРЬ-2005». - Новосибирск, 2005
2. Мазуров Б.Т., Панкрушин В.К. Идентификационный эксперимент:
моделирование системы наблюдений за динамикой земной поверхности и
гравитационного поля в вулканической области. Сб. материалов Международной
выставки и научного конгресса «ГЕО-СИБИРЬ-2005». - Новосибирск, 2005
3. Панкрушин В.К. Математическое моделирование и идентификация
геодинамических систем. - Новосибирск: СГГА, 2002.
4. Панкрушин В.К. Губонин П.Н., Дементьев Ю.В. Автоматизация математической обработки и интерпретации геодезических наблюдений за движениями и деформациями / Под ред. Панкрушина В.К.: Учебное пособие. - Новосибирск: НИИГАиК, 1989.
5. Большое трещинное Толбачинское извержение.(1975 - 1976 гг., Камчатка). Отв. редактор С.А.Федотов. - М.: Наука, 1984.
© Б.Т. Мазуров, В.К. Панкрушин, 2005
