CC BY
36
УДК 551.2+528.9:004+519.876.5 Б.Т. Мазуров, В.К. Панкрушин СГГ А, Новосибирск
ИДЕНТИФИКАЦИОННЫЙ ЭКСПЕРИМЕНТ: ПОСТРОЕНИЕ ФИЗИКОМАТЕМАТИЧЕСКОЙ МОДЕЛИ ДИНАМИКИ ЗЕМНОЙ ПОВЕРХНОСТИ И ГРАВИТАЦИОННОГО ПОЛЯ В ВУЛКАНИЧЕСКОЙ ОБЛАСТИ
В книге [1] в аспекте проблемы идентификации геодинамических систем (ГДС) дана теория метода и математическое обеспечение решения задачи совместной математической обработки и интерпретации пространственно -временных рядов комплексных нивелирных и гравиметрических наблюдений; при этом для сокращения размерности задачи в определяемые параметры вектора состояний вводятся масконы (mass concentration) -_аномалиеобразующие тела (АОТ), как функции пространственных координат и времени. С помощью таких масконов выполняется аппроксимация гравитационного поля в пространстве и времени.
Для разработки методики идентификационного эксперимента, алгоритмического обеспечения и компьютерной технологии решения задачи совместной математической обработки и пространственно-временной интерпретации результатов разнородных комплексных геодезических и гравиметрических наблюдений нами была построена модель вертикальных движений и временных вариаций гравитационного поля, аномалии которого аппроксимируются в пространстве и времени двумя масконами - глубинным и поверхностным.
В качестве физической модели геодинамической системы (ГДС) был выбран локальный участок на земной поверхности, расположенный в области действующего вулкана (рисунок). Система наблюдений за геодинамическим объектом включает нивелирные и гравиметрические наблюдения в пространстве и времени.
Рис. Физическая модель геодинамического объекта и системы наблюдений
Нивелирная сеть состоит из двух условностабильных пунктов (реперов) С1 и С2 (их взаимное положение со временем не меняется) и трех мобильных пунктов М1, М2 и М3 [2].
При моделировании ГДС «Движения физической поверхности и изменения гравитационного поля Земли» из множеств параметров состояний выделим следующие множества (векторы параметров состояний), которые необходимо будет определить в процессе идентификации системы [1]:
X, {р (1.)} - вектор пространственных координат пунктов р (?) сети наблюдений в дискретные моменты времени (эпохи) t = 1,2,... 'I = 1,2,...,^;
Ж,{р(X ,г)} - вектор внешнего потенциала силы тяжести в тех же пунктах или, при задании нормального потенциала и,{р (X,t)}, вектор возмущающего потенциала Т, {р (X, t)};
Х^ (X, t) - вектор параметров движений в широком смысле (изменений как геометрических, так и физических параметров).
Из перечисленных выше векторов формируется так называемый расширенный вектор определяемых параметров состояний ГДС:
XR(X,I) = [XT{р.(I)},ЖТ{р(X,0},XI(X,t)]Т е X^ (1)
или, при задании нормального потенциала и, {Р. (X, t)}, -
XR (X, t)=X {р ц )},ТТ {р (X, t)}, XI (X, 1)]т , (2)
где Т {р (X, t)} - возмущающий потенциал; верхний индекс Т — символ транспонирования.
Входящий в (1), (2) вектор Xв (X, t) формируется в следующем виде:
XВ (X,t) = ^ОфПЗ), XD(ВГПЗ) (X,t)} , (3)
где XD(ФПЗ)(X, t) - вектор параметров движений и деформаций
физической поверхности Земли или, в общем случае, вектор параметров напряженно-деформированного состояния земной коры (горных пород);
Xв(вгпз) - вектор параметров изменений (вариаций) во времени
внешнего гравитационного поля Земли.
Для построения конечномерных моделей с пространством состояний переходим от возмущающего потенциала Т{р(X,t)} к его основным трансформантам в пространстве и времени:
щ, {р, (X, t)} = [АЙТ {р (X, t )}л! {р (X, t )},ц,Т {р (X, t )},СТ {р (X, t )}]Т,
(4)
где А§, £, ц, С - соответственно, аномалии силы тяжести, компоненты уклонения отвесной линии, аномалии высот.
Пусть в состав разнородных комплексных наблюдений У(X,Ж,t) за геодинамической системой (объектом, процессом, явлением) входят измерения превышений (X ,Ж, t) методом геометрического нивелирования
в моменты времени (эпохи) t = 1,2,... между пунктами (реперами) сети
наблюдений на ФПЗ р(Xj,Жу,t) е ФПЗ и р^(Xj,Жу,t) е ФПЗ,
,, у = 1,2,...,N.
Наблюдаемые превышения Иу (X ,Ж, t) могут быть представлены
пространственно-временными геодезическими функционалами высот (геодезических, нормальных, ортометрических) и соответствующих основных трансформант щ{р(Т,X,!)} возмущающего потенциала Т{р(X,t)}.
Линейную модель пространственно-временных рядов результатов измерений превышений в гравитационном поле Земли представим в следующем виде [1]:
Иу (X, Ж, t) = Н у (X, t) - Н , (X, t) + [& (X, t) + А^т (X, t)]А% (X, t) +
+ ц* (X, t )Ауу (X, t) + еу (X, Ж, t), /, у = 1,2,...,
(5)
где еу (X ,Ж, t) - ошибки наблюдений превышений;
Аху (X, t) = ху (!) - х, (t); Ау у (X, t) = у у (t) - у, (!) - приращения
координат.
При моделировании абсолютных наблюдений силы тяжести, принимается, что в вектор параметров состояния ГДС, включены ускорения силы тяжести (XІ,WІ,t)} на пунктах сети наблюдений
р ^Ж, t), где X = (х, у, Н),, = 1,2,...,N, t = 1,2,...
Модель пространственно-временных рядов абсолютных наблюдений силы тяжести строится в следующем виде:
{р, (X, ж, 0} = ^ р (X/7, Ж,г, t - )}-р (X[, иг, t - р (X?, t - ))+4;,в {рг (X,, Ж, г, t - »,„■ +е8 (X, Ж, t).
дН
(6)
Г Г л Г л Г
В (6) {р, (X , ,Ж , ,t )} - счислимая одношаговая прогнозная
оценка фона наблюдаемой силы тяжести на эпоху t по всем наблюдениям до
д г
эпохи t - 1 включительно; ^ [рг (XF, иг, t-)} - вертикальный градиент
дН
нормальной силы тяжести; (X ,Ж, t) - ошибки наблюдений силы тяжести.
Для составления общей системы уравнений разнородных комплексных геодезических и геофизических наблюдений, которые были бы связаны общими параметрами, используется подход, который заключается в использовании понятий аномалиеобразующих тел (АОТ), масконов и точечных масс при аппроксимации аномального гравитационного поля.
Сущность его состоит в следующем: аномальное гравитационное поле в пространстве и времени, представленное возмущающим потенциалом Т[р^ , t)} или его основными трансформантами Щр^Ж,!)} = щ^^Т, t)} и их вариациями Зд{р(X,T, t)}, аппроксимируется комбинациями простых
распределений масконов (АОТ) и их вариаций Мк {Ок (X, 1)},к = 1,2,...,S, тк {Ок (X, t)} [1].
Для составления общей системы разнородных наблюдений Y(X,W,t) и аппроксимации ГПЗ совокупностью аномальных масс (масконов) М(ХД) переходят к представлению расширенного вектора параметров ГДС XR (X, t) в следующем виде:
XR (X, t) = [X, {р а)}, мТ {Ок (X, ()}, XI (X, ()]Т. (7)
Вектор координированных масконов, входящий в определяемые параметры расширенного вектора ГДС, позволяет существенно сократить размерность задачи, так как размерность вектора параметров Мк, к = 1,2,.. £ существенно меньше, чем, если бы в определяемые параметры входил вектор основных трансформант возмущающего потенциала силы тяжести щ , , = 1,2,...,N, у которого размерность намного больше.
В нашем примере изменение гравитационного поля на локальном участке местности характеризуется деятельностью ближайшего вулкана, иначе говоря, перемещением аномальных масс вдоль подводящего канала от зоны образования магмы до магматической камеры, заполнением камеры и выходом магмы на поверхность при извержении вулкана (формирование его конуса). Значения составляющих уклонений отвесных линий можно выразить через аномальные массы и координаты центров масс.
С учетом этого математическую модель состояний геодинамической системы (Рисунок) без учета детерминированных и стохастических воздействий внешней среды [1] представим уравнением XR()= Ф(|,t - 1^( - 1), (8)
где X R () - расширенный вектор определяемых параметров состояния
геодинамической системы, в который входят: геодезические отметки
мобильных пунктов НМ1, Нмг, НМ3, параметры вертикальных смещений мобильных пунктов им 1, им 2, иМ 3, масса МК маскона, аппроксимирующего аномальное гравитационное поле, и масса <МК, аппроксимирующая изменение (вариацию) аномального гравитационного поля. Масса маскона МК аппроксимирует влияние конуса вулкана, а через массу ЗМК может быть определена масса маскона МО, аппроксимирущего влияние пустого шарового маскона в верхней части магматической камеры,
Ф(1 , t -1) - переходная матрица, показывающая закономерность изменений значений параметров от эпохи t -1 к эпохе ! .
Интервал времени между эпохами А| постоянный и равен А1 = 2 месяца (61 сутки). Выполняются четыре эпохи нивелирных и гравиметрических наблюдений в течение шести месяцев. В соответствии с работой [4] нами была смоделирована следующая хронология событий эксперимента. Первая эпоха наблюдений была сделана 7 мая. Ровно через сутки - 8 мая началось извержение вулкана. Вторая, третья и четвертая эпохи наблюдений выполнялись соответственно 7 июля, 6 сентября и 6 ноября.
На протяжении всего периода наблюдений (с 1-й по 4-ю эпоху) происходят вертикальные движения физической поверхности вблизи вулкана. Модель вертикальных движений задана нами как равномерное вращение одного жесткого блока (плиты), на котором находятся пункты М1, М2, М3 [2]. Ось вращения параллельна линии М2-М3 и находится на расстоянии 0.95 км от нее и на расстоянии 0.50 км от пункта М1. Угловая скорость вращения
плиты О является постоянной и равна о = 1.95 *10-5 рад/год. При таком расположении оси и значении угловой скорости модельные смещения пунктов М1, М2, М3 за шесть месяцев составили им1 =-48.70мм, иМ2 = 92.53мм, им з = 92.53мм.
Модельные значения отметок пунктов в четыре эпохи представлены в таблице:
Таблица. Модельные значения отметок пунктов М1, М2 и М3
Пункт II +-> я о2 п +-> Я п К ^4 п К
M1 1151060.00 1151044.00 1151028.00 1151011.00
M2 1365190.00 1365221.00 1365252.00 1365283.00
M3 1280610.00 1280641.00 1280672.00 1280703.00
Гравитационное поле в первую эпоху наблюдений характеризуется значением силы тяжести, вызванной притяжением бесконечной горизонтальной плиты толщиной плиты dH = 8 км и шара,
аппроксимирующего конус вулкана. Плотность масс 8=2.63г/см3.
Для вычисления влияния конуса вулкана в 1-ю эпоху необходимо знать массу шара МК (г = 1), аппроксимирующего этот конус. Для радиуса шара
-5
ЯК (г = 1) = 500 м и плотности 8 = 2.63 г/см , получим массу:
МК (Г = 1) = у л[ЯК(Г = 1)]38 = 1.37706 * 1015 г.
Масса шара, аппроксимирующего магматическую камеру
(магматический очаг) вулкана:
МО(г = 1) = 4/ъ лЯ38 = 1.72133 *1017 г,
где Я = 2500 м - радиус шара, аппроксимирующего магматический очаг вулкана в 1-ю эпоху.
Аномалия силы тяжести для 2-й, 3-й и 4-й эпох будет вызвана излиянием вещества из магматической камеры и увеличившейся массой конуса вулкана. Для построения модельных значений соответствующих этим массам шаровых масконов необходимо задать закон расхода магмы. Моделирование расхода магмы Ж (г) в моменты времени г выполнялось в
соответствии с экспоненциальным законом [3]
Ж (г) = Ж0е ~Ж°11АБ, (9)
где Ж - начальный расход магмы; АБ - избыток магмы в очаге; ^ -время от начала извержения.
Были выбраны следующие параметры для закона расхода магмы. Начальный расход магмы Щ = 1042т / с (1042 тонны в секунду). Избыток
7
магмы в очаге АБ = 1.12 * 10 тонн.
Следствием закона (9) будет формула расчета количества магмы,
извергнутой ко времени г с начала извержения Мизл (г) = АБ(1 - е/ АБ ). Время с начала извержения составляет 60 суток для второй эпохи, 121 и 182 суток для третьей и четвертой эпох соответственно. Получаем с учетом заданных нами параметров и их размерностей следующие значения массы излившейся из резервуара магмы:
Мизл (г = 2) = 1.12 *107 * (1 - е-1000*60*2^*60*60/1 16*107 ) = 4.28 *109 тонн;
Мизл. (г = 3) = 1. 12 *107 * (1 -е^^^^Ч21*24*60*60 /1 16*10') = 6.96*109 тонн;
Мизл. (г = 4) = 1.12*107 *(1 -е-«м*182*24*60«./1. 16*107) = 8 . 61 * 109 тонн.
Этим значениям с учетом плотности 8 = 2.63 г/см3 соответствуют радиусы пустых шаровых масконов в верхней части магматического очага для 2-й, 3-й и 4-й эпох. Применяя формулу связи радиуса и объема шара, получаем:
ЯО(г = 2) = 724 м; ЯО(г = 3) = 855 м; ЯО(г = 4) = 921 м.
При извержении вулкана не все вещество остается в конусе вулкана и на его поверхности, т.к. его часть выбрасывается в атмосферу. Поэтому примем, что в конусе остается 70% от излившегося из магматического резервуара вещества в эпохи 2 и 3, и 90% в эпоху 4. Таким образом, получаем массу конуса вулкана в эпохи 2, 3 и 4:
МК(г = 2) = МК(г = 1) + Мизл (г = 2) * 0.7 = 4.305 *109 тонн; мк (г = 3) = мк (г = 1)+Мшл (г = 3) * 0.7 = 6.200 * 109 тонн;
МК (г = 4) = МК(г = 1) + Мизл (г = 4) * 0.9 = 7.746 * 109 тонн.
Этим значениям соответствуют радиусы шаровых масконов, аппроксимирующих конус для 2-й, 3-й и 4-й эпох:
ЯК(г = 2) = 731 м; ЯК(г = 3) = 826 м; ЯК(г = 4) = 889 м.
Смоделированные массы масконов, как функции пространства и времени, позволили нам описать дискретно гравитационное поле в каждую эпоху - через вычисленные значения компонент уклонений отвесной линии £ т] и силы тяжести g для всех пяти пунктов геодезической сети (значения здесь не приводятся ввиду ограничений по объему статьи).
Таким образом, нами построена физико-математическая модель вертикальных движений земной поверхности и изменений гравитационного поля в вулканической области, для корректной идентификации которой необходима совместная математическая обработка пространственновременных рядов разнородных комплексных нивелирных и гравиметрических наблюдений.
1. Панкрушин В.К. Математическое моделирование и идентификация геодинамических систем. - Новосибирск: СГГА, 2002.
2. Панкрушин В.К. Губонин П.Н., Дементьев Ю.В. Автоматизация математической обработки и интерпретации геодезических наблюдений за движениями и деформациями / Под ред. Панкрушина В.К.: Учебное пособие. - Новосибирск: НИИГАиК, 1989.
3. Большое трещинное Толбачинское извержение. (1975 - 1976 гг., Камчатка). Отв. редактор С.А.Федотов. - М.: Наука, 1984.
© Б.Т. Мазуров, В.К. Панкрушин, 2005
