CC BY
94
ГЕОДЕЗИЯ И МАРКШЕЙДЕРИЯ
УДК 622.2:528.061.4
АППРОКСИМАЦИЯ ГРАВИТАЦИОННОГО ВЛИЯНИЯ ЛОКАЛЬНОГО РЕЛЬЕФА С ИСПОЛЬЗОВАНИЕМ НЕКОТОРЫХ АНАЛИТИЧЕСКИХ МОДЕЛЕЙ И МЕТОДА КОНЕЧНЫХ ЭЛЕМЕНТОВ
Борис Тимофеевич Мазуров
Сибирский государственный университет геосистем и технологий, 630108, Россия, г. Новосибирск, ул. Плахотного, 10, доктор технических наук, профессор кафедры физической геодезии и дистанционного зондирования, тел. (383)343-29-11, e-mail: [email protected]
Определение поля силы тяжести, его трансформант не является тривиальной задачей. Для ее решения возможно сочетание аналитических описаний некоторых элементарных пространственных тел с конечно-элементным разбиением рельефа. Здесь описаны некоторые возможные к практическому использованию аналитические модели. А именно, кроме традиционных точечных представлений масс, качество аппроксимации может быть улучшено использованием в качестве конечных элементов параллелепипедов и цифровых моделей рельефа. Данный способ позволяет более детально описывать гравитационное влияние локального рельефа.
Ключевые слова: сила тяжести, аналитические модели, конечные элементы.
APPROXIMATION GRAVITATIONAL INFLUENCE LOCAL RELIEF USING SOME ANALYTICAL MODELS AND FINITE ELEMENT METHOD
Boris T. Mazurov
Siberian State University of Geosystems and Technologies, 630108, Russia, Novosibirsk, 10 Plakhotnogo St., Ph. D., Prof., Department of Physical Geodesy and Remote Sensing, tel. (383)343-29-11, e-mail: [email protected]
Gravity field definition, it is not trivial to transformant. To solve the combination of analytical descriptions of some elementary spatial solids with finite-element partitioning. Here are some possible practical use of analytical models. Namely, in addition to traditional bitmap representations of mass quality of approximation can be improved using a finite-element boxes and digital elevation models. This method allows more detail to describe the gravitational influence of the local topography.
Key words: gravity, analytical models, finite element.
5
Вестник СГУГиТ, вып. 3 (31), 2015
Для обеспечения желательной корректной интерпретации геодезических наблюдений необходимо отталкиваться от фундаментальной научной основы, как это представлено в [1]. Желание качественного и более точного количественного понимания меняющихся со временем окружающих Землю физических полей - это многие века главная составляющая часть парадигмы всех наук о Земле. Главная преследуемая исследователями цель - прогнозирование такой сложной природно-технической системы, какой является земной приповерхностный слой в сочетании с находящимися на поверхности и в недрах объектами, имеющими техногенное происхождение. Самым логичным методом прогнозирования по количественным данным является математическое моделирование. Корректно выполненная процедура моделирования с учетом статистических свойств комплексных геодезических и геофизических измерений и предположений о физической природе изменений геодинамических объектов и их напряженно-деформированного состояния делает более обоснованным прогноз природных и техногенных катастроф [2].
Поиски эффективных методов моделирования - объективная особенность человеческого познания окружающего мира, имеющая историю в несколько тысячелетий. Галилею принадлежит мысль, что книга Природы написана на языке математики. Современные темпы развития вычислительной техники объективно требуют математической и алгоритмической поддержки востребованных программных продуктов. Для удовлетворения этой технологической потребности ведется глобальная работа. В сфере геоинформатики некоторые математические приложения, алгоритмы предложены в [3-6]. Вопросы математического моделирования и визуализации динамики природных и техногенных объектов решались в работах [7, 8].
Многие науки о Земле используют для решения своих задач данные о поле силы тяжести в глобальном, региональном и локальном масштабах. Достоверные выводы о состоянии поверхностного слоя земной коры гравитационных полей и их динамике необходимы в геодинамике, для поиска предвестников землетрясений [1, 2, 8], вулканизма [9-14] и т. д. Хорошие технологические возможности имеет сочетание гравиметрических исследований с ГНСС-технологиями, имеющими массу положительных свойств. Очень важно выполнять адекватную математическую обработку измерений, которые часто представляют собой длинные пространственно-временные ряды. Для успешного решения задачи мониторинга нужны корректные математические методы обработки данных, алгоритмы. Такими могут быть метод наименьших квадратов, регрессионные модели и нейронные сети. Для достоверного результата исследований необходимо также привлечение аналитических моделей гравитационного влияния тел простой геометрической формы (сфера, цилиндр, призма и т. д.). Добавление метода конечных элементов позволяет получить достаточно точное описание гравитационного поля.
Очень важным объектом исследований как в научном аспекте, так и в практическом является меняющееся гравитационное поле. Приложения полученных
6
Геодезия и маркшейдерия
данных о силе тяжести и их интерпретации разнообразны. Например, геодинамические исследования. Значимые вертикальные смещения точек земной поверхности и смещения уровенных поверхностей во времени могут возникать, в частности, при разработке крупных месторождений полезных ископаемых [15]. Перераспределение больших объемов пород и руды вызывает значительное изменение поля силы тяжести, а недоучет влияния перемещаемых масс в результаты нивелирования может быть причиной неверного представления о картине вертикальных движений и оказывает влияние на результаты геодезических измерений и определяемые по ним деформации земной поверхности. Понимание локального гравитационного поля и его учет важен при мониторинге вулканизма и корректной интерпретации геодезических наблюдений различного вида [9-14]. Зависимость поправок за плоский и сферический слои в неполной топографической редукции от их толщины и радиуса учитываемой зоны описана в [16, 17].
Определение поля силы тяжести, его трансформант не является тривиальной задачей. Поэтому часто вынужденно для достижения цели делается сочетание аналитических описаний некоторых элементарных пространственных тел с последующим конечно-элементным разбиением сложного элемента рельефа. Многие тела простой формы и постоянной плотности имеют влияние на силу тяготения, выражаемое аналитически в замкнутой форме (в квадратурах). Для шара радиуса R с постоянной плотностью или состоящего из концентрических слоев известна формула
где Ар - разность плотности возмущающего тела и плотности окружающей его среды; G - гравитационная постоянная; х, y, z - координаты центра шара. Часто такие гравитирующие тела называют точечными массами, масконами.
Если гравитирующее тело представлять как бесконечно протяженное в горизонтальной плоскости, то формулы упрощаются. Есть технологическое допущение, что это возможно, когда горизонтальные размеры тела в два раза больше глубины его залегания. Это позволяет использовать аналитически выражаемую модель «бесконечного» плоского слоя. Есть еще примеры аналитического подхода к анализу гравитирующего влияния тел простой формы, например, вертикального цилиндра. Но формула работает только для точек, расположенных на его оси. В случае конуса это также ограничено расположением точки пространства и гравитирующего конуса.
Изучение такой сложной самоорганизующейся природной системы, которой является меняющаяся со временем земная поверхность, - очень актуальная задача геодинамических исследований. Ее успешное решение предполагает комплексность наблюдений различных видов с соответствующей последующей математической обработкой [1, 2, 12]. При решении некоторых задач геодина-
Ag = -4 kGR3 Ар
(1)
7
Вестник СГУГиТ, вып. 3 (31), 2015
мики, когда интерпретация результатов наблюдений за движениями земной поверхности должна учитывать изменения поля силы тяжести, встает вопрос аппроксимации гравитирующих тел в исследуемом районе. Для достижения цели моделирования выбор сложности модели должен основываться на точностных характеристиках экспериментальных данных. Например, какой сложности (подробности) должна быть модель, аппроксимирующая гравитационное влияние рельефа, чтобы высокоточные геодезические и гравиметрические измерения при их комплексной математической обработке соответствовали выбранной модели гравитирующего тела?
В работах [12, 14] даются примеры решения обратных задач, когда результатом совместной математической обработки многократных геодезических и гравиметрических наблюдений на земной поверхности являются оценки не только координат пунктов, их смещений, но и масс гравитирующих тел, а также изменений этих масс. Эти примеры даны в отношении вулканического извержения и подготовки к нему. Меняющимися гравитирующими телами являлись шарообразный глубинный магматический очаг вулкана и его поверхностный конус.
Достаточно типичной элементарной формой тела для рельефа земной поверхности является конус. Причем, более общий случай - усеченный конус. В естественной среде, кроме вулканов, это могут быть отдельные составляющие горных хребтов. В техногенной сфере мы имеем не только конусообразные возвышения, но и конусообразные выемки [15]. При открытых горных разработках пространство выработанной породы является конусообразным. Конусообразными являются, например, кимберлитовые месторождения. Разработка трубки Мир (Якутия) привела к образованию конусообразного карьера глубиной 520 метров. Отвалы также часто представляют собой усеченные конусы.
В работах [12, 14] конусообразные гравитирующие тела аппроксимировались шаром (одной точечной массой). Будет вполне разумным уточнить аппроксимирующую модель конуса путем увеличения числа точечных масс. Но при этом желательным является, чтобы число оцениваемых параметров осталось минимально необходимым. Тем самым будет выполняться требование большей избыточности измерений, которая необходима при математической обработке. Например, пусть оцениваемым параметром будет оставаться общая масса конуса, но рассредоточенная определенным образом в пространстве (в пяти точках) (рис. 1).
Таким образом, имеем следующую задачу. Известны размеры усеченного конуса: R1 - радиус нижнего основания, r - радиус верхнего основания, h1 - высота конуса (рис. 2).
Необходимо разбить конус на пять равнообъемных частей и для каждой части найти координаты центра тяжести. Вариантов разбиения на пять равнообъемных частей может быть много. Здесь выбран следующий (рис. 2). Горизонтальной плоскостью отсекается верхний усеченный конус объемом одну пятую от объема всего конуса. Нижний усеченный конус делим на четыре сектора двумя вертикальными взаимноперпендикулярными плоскостями.
8
Геодезия и маркшейдерия
а)
б)
в)
Рис. 1. Усеченный конус и модели аппроксимации его гравитационного влияния: а) усеченный конус; б) одноточечная модель; в) пятиточечная модель
Рис. 2. Усеченный конус и его разбиение на пять равнообъемных частей
Найдем радиус R2 горизонтального сечения и высоту верхнего конуса. Воспользуемся формулами вычисления объема усеченного конуса для исходного V1 и для верхнего V2, полученного после отсечения горизонтальной плоскостью
jt ^hi , 2 2 t-j ч Vi = 31(Ri +r + Rir); (2)
Tr ^2 „2 2 л \ V2 = 32(R2 + r + R2r) . (3)
9
Вестник СГУГиТ, вып. 3 (31), 2015
Мы выбрали вариант разбиения конуса на пять пространственных тел равного объема. Поэтому
^=5 и.
(4)
и если учесть изображение сечения конуса вертикальной плоскостью (рис. 3), мы получаем из отношения сторон треугольника
h2 = hl Rl
R1 - r
(5)
r
\ hl —i i i i i i i Ri \
/Ri-r, i hi i i \
j R\-r i i i i i i Ri \
Рис. 3. Сечение конуса вертикальной плоскостью
Объемы усеченных конусов могут быть также представлены через разность кубов радиусов. Делаем подстановку формулы (5) в формулу (3) и после этого формула (4) будет следующая:
nh (R2 - r) 2 2 „ ч 1 nhi /n2 2 n 4
---f (R2 + r + R2 r) = 1—I(R1 + r + Rir);
3 (R1 - r) 5 3
(R2 - r)(R2 + r2 + R2r) = 1(R1 - r)(R12 + r2 + R1r).
(6)
(7)
И теперь заменяем формулу (7), используя формулу разности кубов. Для наших пространственных тел (см. рис. 2) зависимость будет такая:
D3 3 1/ D3 3Ч
R2 - r = 5(R1 - r )■-
(8)
10
Геодезия и маркшейдерия
и, следовательно,
R = 35( R
3Ч 3
r ) + r
(9)
Теперь стало возможным аналитическое определение координат центра тяжести верхнего усеченного конуса. Для этого будем использовать классический подход из теории сопротивления материалов к таким вычислениям через суммы статических моментов элементарных геометрических фигур (прямоугольник, треугольник, сектор круга).
Вертикальное сечение усеченного конуса представляет собой трапецию, которую можно разделить на две симметричные прямоугольные трапеции (рис. 4). Прямоугольную трапецию можно представить сочетанием прямоугольника и треугольника.
Рис. 4. Усеченный конус:
а) после разбиения на пять равнообъемных частей; б) вертикальное сечение верхнего усеченного конуса; в) вертикальное сечение нижнего объемного сектора; г) горизонтальное сечение нижнего объемного сектора
Обращаясь к теории сопротивления материалов, находим статические моменты прямоугольника и треугольника. Это произведение площади этих элементарных фигур на расстояние от принятых осей координат. Для общего сечения в виде равнобедренной трапеции с осью z, проходящей через ось симметрии, координата x центра тяжести будет равна нулю. Координата z центра -расстояние от основания трапеции - вычисляемая величина. Согласно правилу теории сопротивления материалов, координаты центра тяжести всей трапеции
11
Вестник СГУГиТ, вып. 3 (31), 2015
будут находиться как результат деления суммы статических моментов этих фигур на сумму их площадей. Таким образом, для сечения верхнего объемного сектора центр тяжести будет иметь координаты
x = 0;
z =
rh2 + (R2 - r) h22 2__________6
rh2
+
(R2 - r )h2 2
(10)
С четырьмя нижними объемными секторами поступим следующим образом. Для 90-градусного сектора плоскость симметрии будет проходить под углом 45°. Сечение этой плоскостью будет представлять собой прямоугольную трапецию, которую мы представляем сочетанием прямоугольника и треугольника.
Вычислим площади фигур: S1 - прямоугольника, S2 - треугольника
S = R2(h -*2); S2 = ±(h -a2)(Ri -R2). (11)
Координаты центра тяжести трапеции в целом
x
S1+ S2 1 2 2
^ R1 + 2 R2^
у .
S1 + S2
(12)
S1
h - h2 „ h - h2
—----2 + S 2 —-2
z = •
S1 + S2
(13)
Определяем радиус горизонтального сечения конуса на высоте, равной z :
R3 = R1 (R1 R2)
h1 h2
(14)
Найдем координаты центра тяжести С четверти этого сечения. Это четверть круга радиусом R3. С учетом уравнения окружности, координаты будут вычисляться
X.
Ус
R3
Idx I
0 0
2 2
- X
ydy
1 ^R32 4 3
4R3
3n
(15)
12
Геодезия и маркшейдерия
Итоги определения координат пяти точечных масс представлены в таблице. Координаты соответствуют системе координат (см. рис. 1, а).
Таблица
Координаты центров тяжести точечных масс
Координаты Масса 1 Масса 2 Масса 3 Масса 4 Масса 5
X 0 Xc - Xc - Xc xc
Y 0 Vc Vc - Vc - Vc
Z h1 - h2 + z z' z' z' z'
Аппроксимация гравитационного влияния конуса, в частности, с использованием в качестве конечного элемента прямоугольной пластины с нулевой толщиной, но обладающей массой, предложена Ю. В. Дементьевым [16, 17]. Иное пространственное тело в качестве конечного элемента использовалось в [18, 19]. Весь объем конуса представлялся совокупностью кольцевых секторов (криволинейные параллелепипеды). Недостаточно технологически еще освоена формула гравитационного влияния параллелепипеда (рис. 5), которая представлена выраженным аналитически замкнутым интегралом [20].
Рис. 5. Гравитирующее тело в виде параллелепипеда
После решения интеграла получается формула вычисления аномальной силы тяжести в точке поверхности P
x У
Ag = G-Ар||| - x• ln(y + r) -y• ln(x + r) + z• arctg——
z r
X2 У2
X1 У1
z2
z, ’
(16)
где r = л[х2 + y2 + z2.
13
Вестник СГУГиТ, вып. 3 (31), 2015
В работе [23] предложена идея сочетания аналитической модели гравитационного влияния параллелепипеда (16) и цифровых моделей рельефа. Авторами выполнены вычислительные эксперименты для некоторых реально существующих локальных участков - вулкана Св. Елена (США, штат Вашингтон) и прибрежной территории полуострова Камчатка. Во втором случае были учтены рельеф суши и дна, а также различная плотность морской воды и горных пород.
Самым простым вариантом описания гравитационного влияния локальных элементов рельефа является использование точечной модели. Для некоторых прикладных задач это может быть достаточным. Но существуют ситуации, которые требуют более точной аппроксимации. Например, геодинамические исследования, прогноз землетрясений и вулканизма, вопросы горной механики, обеспечения безопасности при разработке месторождений полезных ископаемых, включая нефтяные и газоносные объекты. Правильное сочетание аналитических моделей и метода конечных элементов, как в [21], позволяет повысить точность описания локального гравитационного поля по сравнению с существующими методиками примерно на 3-4 % для локальных объектов, имеющих размеры территории около 20-100 км с горным и предгорным рельефом. Основой повышения точности является сочетание аналитических моделей, метода конечных элементов и цифровых моделей рельефа.
БИБЛИОГРАФИЧЕСКИЙ СПИСОК
1. Идентификация движений и напряженно-деформированного состояния самоорганизующихся геодинамических систем по комплексным геодезическим и геофизическим наблюдениям: монография / В. А. Середович, В. К. Панкрушин, Ю. И. Кузнецов, Б. Т. Мазуров,
B. Ф. Ловягин; Под общ. ред. В. К. Панкрушина. - Новосибирск: СГГА, 2004. - 356 с.
2. Мазуров Б. Т., Панкрушин В. К., Середович В. А. Математическое моделирование и идентификация напряженно-деформированного состояния геодинамических систем в аспекте прогноза природных и техногенных катастроф // Вестник СГГА. - 2004. - Вып. 9. -
C. 30-35.
3. Вовк И. Г. Геометрическое моделирование линейных объектов в прикладной геоинформатике // Вестник СГГА. - 2013. - Вып. 4 (24). - С. 57-62.
4. Вовк И. Г. Моделирование формы и оценка размеров систем в прикладной геоинформатике // Вестник СГГА. - 2013. - Вып. 2 (22). - С. 17-25.
5. Вовк И. Г. Определение геометрических инвариантов пространственной кривой в прикладной геоинформатике // Вестник СГГА. - 2012. - Вып. 3 (19). - С. 51-62.
6. Вовк И. Г. Определение геометрических инвариантов поверхности в прикладной геоинформатике // Вестник СГГА. - 2012. - Вып. 4 (20). - С. 59-69.
7. Дорогова И. Е. Изучение движений и деформаций земной коры на геодинамическом полигоне Таштагольского железорудного месторождения // Вестник СГГА. - 2010. -Вып. 2 (13). - С. 9-13.
8. Мазуров Б. Т., Дорогова И. Е., Дербенев К. В. Горизонтальные движения земной коры вращательного характера, наблюдаемые на геодинамических полигонах // Интерэкспо ГЕО-Сибирь-2012. VIII Междунар. науч. конгр. : Междунар. науч. конф. «Геодезия, геоинформатика, картография, маркшейдерия» : сб. материалов в 3 т. (Новосибирск, 10-20 апреля 2012 г.). - Новосибирск: СГГА, 2012. Т. 1. - С. 232-236.
14
Геодезия и маркшейдерия
9. Мазуров Б. Т. Модель вертикальных движений земной поверхности и изменений гравитационного поля в районе действующего вулкана // Изв. вузов. Геодезия и аэрофотосъемка. - 2007. - № 2. - С. 97-106.
10. Мазуров Б. Т. Модель системы наблюдений за вертикальными движениями земной поверхности и изменениями гравитационного поля в районе действующего вулкана // Изв. вузов. Геодезия и аэрофотосъемка. - 2007. - № 3. - C. 93-101.
11. Мазуров Б. Т. Совместная математическая обработка разнородных комплексных геодезических и геофизических наблюдений за движениями земной поверхности и изменениями аномальных масс // Изв. вузов. Горный журнал. - 2007. - № 6. - C. 30-38.
12. Мазуров Б. Т. Совместная математическая обработка и интерпретация нивелирных и гравиметрических наблюдений за вертикальными движениями земной поверхности и изменениями гравитационного поля в районе действующего вулкана // Изв. вузов. Геодезия и аэрофотосъемка. - 2007. - № 4. - С. 11-20.
13. Мазуров Б. Т. Моделирование и идентификация геодинамического объекта в вулканической области по комплексным нивелирным и гравиметрическим наблюдениям // Вестник СГГА. - 2006. - Вып. 11. - С. 84-94.
14. Мазуров Б. Т. Идентификация напряженно-деформированного состояния вулканической области по результатам математической обработки разнородных геодезических и геофизических наблюдений // Изв. вузов. Горный журнал. - 2007. - № 7. - С. 58-62.
15. Мазуров Б. Т. Математическая обработка нивелирных и гравиметрических наблюдений в условиях извлечения и перемещения больших объемов руды и пород // Изв. вузов. Горный журнал. - 2006. - № 4. - С. 99-104.
16. Дементьев Ю. В. Зависимость поправок за плоский и сферический слои в неполной топографической редукции от их толщины и радиуса учитываемой зоны // Вестник СГГА. -2010. - Вып. 2 (13). - С. 13-17.
17. Дементьев Ю. В. О возможности и необходимости определения аномалий силы тяжести в полной топографической редукции // Вестник СГГА. - 2011. - Вып. 3 (16). - С. 3-14.
18. Мазуров Б. Т., Некрасова О. И. Гравитирующее влияние конусообразных форм рельефа на результаты геодезических измерений // Геодезия и картография. - 2013. - № 5. -С. 2-6.
19. Мазуров Б. Т., Некрасова О. И. Конечно-элементная модель конусообразных форм рельефа для учета их гравитирующего влияния на результаты геодезических измерений // Геодезия и картография. - 2013. - № 6. - С. 42-45.
20. Nagy В. The gravitational attraction of a right angular prism / Geophysics 31. - 1966. -P. 362-371.
21. Мазуров Б. Т., Некрасова О. И. Аппроксимация гравитационного влияния локального рельефа по его цифровым моделям // Геодезия и картография. - 2014. - № 7. - С. 2-4.
Получено 12.08.2015
© Б. Т. Мазуров, 2015
15
