А.А. Булычев, И.В. Лыгин, В.Р. Мелихов - Численные методы решения прямых задач грави- и магниторазведки, страница 21

Вычислим значения поля в расчетных точках поаналитической формуле и с помощью спектрального представления.Результат, на первый взгляд, оказывается неожиданным.Как и следовало ожидать, теоретическое поле в начале профиляблизко к нулевым значениям и возрастает к его окончанию. Поле,рассчитанное через спектральную функцию, имеет максимум при значенииx = 7 км, далее оно убывает и в конце профиля вновь возрастает,практически совпадая с теоретическим.Для того чтобы понять, с чем связано появление “ложного”максимума в начале профиля, следует рассмотреть формулу численногообратного преобразования Фурье:1401f ( x) =2π∞ˆ (ω )e iω x dω ≈ 1f∫2π−∞Nω∑ fˆ (ωj = − Nωj)eiω j x∆ω .2π, где N –N∆xчисло расчетных точек на профиле, ∆x – расстояние между расчетнымиточками (в нашем примере N = 128, ∆x =1).

Длина расчетного профиля –2π. Вычислим значение функцииX l = N∆x , значение частоты – ω j = jN∆xf ( x ) за пределами профиля в точке x = X l + x k = N∆x + x k . Тогда:Определим шаг по частоте следующим образом: ∆ω =1f ( x) =2πNω1iω xfˆ (ω j )e j ∆ω =∑2πj = − Nω1=N∆xNω∑ fˆ (ωj = − Nωj)eNω∑j = − Nω⎛ 2π ⎞i⎜ j⎟ ( N∆x + xk )⎝ N∆x ⎠fˆ (ω j )e1=N∆x⎛ 2π ⎞i⎜ j⎟x⎝ N∆x ⎠2π=N∆xNω∑ fˆ (ωj = − Nω1=N∆xNω∑j)e2π⎛⎞i ⎜ j 2π + jxk ⎟N∆x ⎠⎝fˆ (ω j )ej = − Nω⎛ 2π⎞i⎜ jxk ⎟⎝ N∆x ⎠== f ( xk ) .Из проведенных выкладок следует, что при таком задании расчетныхпараметров значение поля в точках x = x k + mX l (m = 0, ±1, ±2,…) будутповторяться, т.е.

рассчитанная таким образом функция будетпериодической.17. Приведенные примеры говорят о том, что вычисление аномальныхзначений поля силы тяжести или магнитного поля через их спектральныехарактеристики требует значительных как технических затрат, связанных свыбором шага дискретизации спектральных функций по частоте, выборомграничной частоты, выше которой значениями спектральной функцииможно пренебречь, так и временными затратами, поскольку для каждойрасчетной точки необходимо суммировать ряд Фурье, коэффициентыкоторого должны рассчитываться.18. В то же время, как мы уже отмечали, существует ряд прямых задач,где применение спектральных преобразований позволяет создатьэффективные алгоритмы по их решению.

Эти задачи описываютсяинтегралами типа свертки:f ( x) =∞∫ g( x )h( x − ξ )dξ = g( x ) ∗ h( x ) ,−∞141f ( x, y) =∞ ∞∫ ∫ g( x , y )h( x − ξ , y − η )dξdη = g( x , y ) ∗ h( x , y ) .− ∞− ∞К числу таких задач относится задача вычисления поля силы тяжестиот горизонтальной плоскости, расположенной на глубине ζ с заданной наней распределением плотности δ ( x ) :∞g z ( x , z ) = 2G ∫ δ (ξ )−∞(ζ − z )dξ ,(ξ − x ) 2 + (ζ − z ) 2задача определения поля силы тяжести, создаваемого горизонтальнымслоем с глубинами кромок ζ1 и ζ2 с известным в нем распределениемплотности:(ξ − x ) 2 + (ζ 2 − z ) 2dξ ,g z ( x , z ) = 2G ∫ δ (ξ ) ln(ξ − x ) 2 + (ζ 1 − z ) 2−∞∞и ряд других.Эффективность решения этих и аналогичных им задач основана натом, что если поле вычисляется также на горизонтальной плоскости(z=const), то для его вычисления необходимо вычислить спектры функцийплотности, ядра преобразования, их перемножить, и от полученногопроизведения вычислить обратное преобразование Фурье.19.

Рассмотрим этот алгоритм на примере решения прямой двухмернойзадачи гравиразведки от горизонтального слоя. Обычно распределениеплотности в таком слое задается в виде дискретных величин сравномерным шагом вдоль профиля. Для того, чтобы рассчитатьспектральнуюфункциюплотностислояприменяютпрямоепреобразование Фурье в “дискретном” аналоге (ДПФ), которое имеет вид:N −1δˆ (ω j ) = ∑ δ ( x k )e− iω j x k∆x , j = 0,..., N − 1 .k =0В этой формуле предполагается, что x k = k∆x , k = 0,..., N − 1 (N – число2π2π, частота – ω j = j.точек на профиле), шаг по частое – ∆ω =N∆xN∆xДля вычисления дискретного спектра функции можно положить∆x = 1, тогда формула прямого дискретного преобразования Фурьеприобретет вид:142N −1δˆ j = ∑ δ k e⎛ 2π ⎞jk ⎟−i⎜⎝ N⎠,k =0где δˆ j = δˆ (ω j ) , δ k = δ ( x k ) .

Соответственно, обратноепреобразование Фурье определяется следующим образом:1δk =NN −1∑ δˆ j e⎛ 2π ⎞i⎜jk ⎟⎝ N⎠дискретное.k =0Для вычисления как прямого, так и обратного преобразований Фурьеиспользуется алгоритм быстрого преобразования Фурье (БПФ), которыйпозволят во много раз сократить вычислительное время по сравнению свычислением этих преобразований непосредственно по представленнымформулам. Особенно ощутимым этот выигрыш во времени становится прибольших значениях N.20. Таким образом, решение прямой задачи можно было быорганизовать следующим образом: 1) рассчитать дискретный спектр отфункции плотности; 2) вычислить спектр ядра преобразования по заданнойаналитической формуле; 3) перемножить спектры; 4) осуществитьобратное дискретное преобразование Фурье.

В результате таких действийбудут рассчитаны значения поля с равномерным шагом по профилю в Nточках.Однако такой подход может привести и к ошибочному результату,связанному с тем, что рассчитанная спектральная функция ядрапреобразования может затухать медленно. В результате в решениипроявится эффект Гиббса, о котором нами уже упоминалось.С тем, чтобы избежать возникновения таких погрешностей, следуетиспользовать алгоритм, основанном на дискретной свертке двух функций –дискретной функции плотности и ядра преобразования:g z ( x j ) = ∑ δ (ξ k ) g 1 ( x j − ξ k , z ) = ∑ g 1 (ξ k , z )δ ( x j − ξ k ) .kkЕсли рассматривать эту формулу применительно к моделигоризонтального слоя, то δ ( x k ) соответствует значениям плотности,заданным с равным интервалом вдоль профиля, g 1 ( x k ) – поле силытяжести, создаваемое прямоугольной призмой (ячейкой), на которыеразбит слой, с постоянной единичной плотностью.

Для быстроговычисления такой свертки используется алгоритм “быстрой дискретнойсвертки”. Идея такой свертки состоит в том, что с помощью БПФвычисляются дискретные спектры ядра преобразования и дискретной143функции плотности, которые перемножаются, и затем от полученногопроизведения осуществляется обратное преобразование Фурье. Такойалгоритм вычисления дискретной свертки помогает во много раз сократитьвычислительное время по сравнению с ее прямым вычислением. Такаяэкономия времени особенно важна при решении трехмерных задач.Кроме того, экономия во времени может быть достигнута также засчет того, что, создав комплексную функцию, действительная частькоторой будет, например, соответствовать функции плотности, а мнимая –функции ядра преобразования, с помощью одного ДПФ можно вычислитьдискретные спектры сразу двух функций одновременно. Для того чтобывыделить эти спектры из полученного спектра комплексной функциинеобходимовоспользоватьсясвойствамисимметрии,которымудовлетворяют действительные и мнимые части спектральных функций.Дальнейший выигрыш во времени может быть достигнут за счеттого, что при расчете эффекта от системы горизонтальных слоев, нетнеобходимости для каждого отдельного слоя осуществлять обратноепреобразование Фурье.

Можно суммировать эффекты слоев в виде ихспектральных представлений и только на последней стадии осуществитьобратное ДПФ с целью вычисления пространственного распределенияполя.Такой подход к решению прямой задачи гравиразведки имагниторазведка оказывается особенно эффективным в случае расчетаполя от трехмерного куба с заданным в нем распределением свойств(плотности или намагниченности) среды.Литература.1. Алексидзе М.А. Приближенные методы решения прямых и обратныхзадач гравиметрии. – М: Наука. 1987.

336 с.2. Булычев А.А., Кривошея К.В., Мелихов В.Р., Зальцман Р.В. Вычислениеаномального гравитационного потенциала и его производных насфере. // Вестн. Моск. Ун-та. Сер.4. Геология. 1998. №2, С. 42 − 46.3. Геофизические поля и строение дна океанских котловин. / Под ред.Непрочнова Ю.П. – М.: Наука.

1990. 220 с.4. Голиздра Г.Я. Основные методы решения прямой задачигравиразведки на ЭВМ. – Обзор ОНТИ ВИЭМС. Сер. Регион., развед.и промысл. геофиз. М. 1977. 99с.5. Голиздра Г.Я. Вычисление магнитного поля масс переменногонамагничения. // Геофиз. журнал. 1981. Т. 3. №5. С 98-104.6.

Гравиразведка. Справочник геофизика. – М.: Недра. 1990. 607 с.7. Каппелини В., Константинидис А.Дж., Эмилиани П. Цифровыефильтры и их применение. – М. 1983.8. Мелихов В.Р., Булычев А.А., Састри Р.Г. Особенности решенияпрямой задачи гравиразведки с использованием быстрого144преобразования Фурье // Морские гравиметрические исследования.М., 1984. С. 73 − 80.9. Старостенко В.И., Манукян А.Г., Заворотько А.Н. Методы решенияпрямых задач гравиметрии и магнитометрии на шарообразныхпланетах. Киев.: Наукова думка. 1986. 112 с.10.

Страхов В.Н., Романюк Т.В., Фролова Н.К. Методы решения прямыхзадач гравиметрии, используемые при моделировании глобальных ирегиональных гравитационных аномалий // Новые методыинтерпретации гравитационных и магнитных аномалий. М: ИФЗ АНСССР.

1989. с. 118-235.11. Parker R.L. A new method for modeling marine gravity and magneticanomalies. // Jour. of Geophys. Research. 1974. V 79. N 14. P. 2014–2016.12. Talwany M., Ewing M. Rapid computation of gravitational attraction ofthree dimensional bodies of arbitrary shape. // Geophys. 1960. V.25. N1. P.203 – 225.Лекция 10. Решение прямой задачи магниторазведки для тел с высокоймагнитной восприимчивостью.В предыдущих лекциях нами рассматривались случаи решенияпрямой задачи от тел с заданной постоянной намагниченностью.