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

Отсюда следует вывод с одной стороны о том, чточрезмерно детальное представление модели может и не дать точногорезультата, а с другой – о том, что по возможности стоит использоватьпредставление чисел с удвоенной точностью.10. Рассмотрим теперь вопрос, связанный с временными затратами навычисление прямой задачи. Будем рассматривать туже модельпрямоугольной призмы, которую представим в виде наборапрямоугольных ячеек (призм) количеством Nx=1024, Nz=512.

Числорасчетных точек на профиле будет равно 100.Организуем вычислительный процесс следующим образом. Будемрассчитывать влияние каждой ячейки в расчетную точку по формуле,описывающей поле от прямоугольника в действительной области. Время,затраченное на такой расчет, оказалось равным 53,28 с. Теперь будемаппроксимировать каждую ячейку горизонтальным цилиндром.

Формулавычисления поля от горизонтального цилиндра значительно проще посравнению с формулой вычисления поля от прямоугольника, и она требуетменьших временных затрат. Время, затраченное ЭВМ при таком подходе крешению прямой задачи, составило 4,39 с.Ясно, что увеличение количества элементарных тел, с помощьюкоторых аппроксимируется реальная геологическая среда может привестик неоправданным временным затратам. Это особенно становитсяактуальным при решении трехмерных задач.134Изменим алгоритм расчета.

Для этого рассмотрим какой эффектбудетсоздаватьодингоризонтальныйрядпредставленныйпрямоугольными ячейками одинакового размера, при этом каждая ячейкабудет иметь свою постоянную плотность. Тогда для эффекта такого слояможно записать:∞ ζ2ζ −zdξζ =g z ( x , z ) = G ∫ ∫ δ (ξ , ζ )22()()xzξ−+ζ−− ∞ζ 1ξi −= G∑i∆ξ2 ζ2∫ξζξi −∆2∫ δ (ξ i )1ζ −zdξdζ =(ξ i − x ) 2 + (ζ − z ) 2= ∑ δ (ξ i ) g 1 ( x − ξ i , z ) = ∑ g 1 (ξ i , z )δ ( x − ξ i ) ,iiгде g 1 (ξ i , z ) – эффект призмы с единичной плотностью, горизонтальными∆ξ∆ξ ⎞⎛координатами ⎜ ξ i −,ξ i +⎟ и глубинами (ζ 1 , ζ 2 ) в точку с22 ⎠⎝координатами (0,z), или, что тоже самое, это – эффект призмы с⎛ ∆ξ ∆ξ ⎞,координатами ⎜ −⎟ в точку с координатами (ξi,z).

Если поле⎝ 2 2 ⎠рассчитывается на постоянном уровне z=const, и шаг между расчетнымиточками ∆x=∆ξ, то представленную формулу можно записать в следующемвиде:g z ( x j ) = ∑ δ (ξ k ) g 1 ( x j − ξ k , z ) = ∑ g 1 (ξ k , z )δ ( x j − ξ k ) ,kkт.е. в виде дискретной свертки двух функций – g 1 ( x k ) , δ ( x j ) . Тем самыммы пришли к идее палеток. Нам достаточно рассчитать эффект однойпризмы, расположенной в начале координат, в точки с координатами xj, идалее для каждой расчетной точки умножить эти значения насоответствующие плотности и складывать их между собой.

Такой подходзначительно сокращает вычислительное время по сравнению с расчетомэффекта от каждой ячейки в расчетную точку. Однако здесь надо обратитьвнимание на то, что число коэффициентов такой палетки всегда будетограничено, и результаты будут только там, где палетка полностьюложится на ячейки, которыми представляется модель нашего слоя.Для того, чтобы рассчитать эффект от нашей модели призмы,заданной в виде ячеек как по горизонтали, так и по вертикале, необходимонашу модель с краев дополнить пустыми ячейками, т.е. ячейками с135нулевой плотностью. При этом, поскольку модель призмы по горизонталипредставлена числом Nx=1024, то число ячеек, которыми следуетдополнить слои, должны также составлять 1024 как слева от модели, так исправа.

Не смотря на то, что число ячеек при таком подходе возросло, темне менее, выигрыш во времени оказывается значительным, поскольку нетнеобходимостирассчитыватьэффектоткаждойячейки,аппроксимирующей исходную модель, в каждую расчетную точку напрофиле, а достаточно рассчитать эффект от одной призмы каждогогоризонтального слоя.Выигрыш во времени может быть увеличен еще в большей степени,если для вычисления свертки использовать так называемый алгоритм“быстрой дискретной свертки”, который мы рассмотрим несколькопозднее.11. Рассмотрим некоторые особенности решения прямых задач наоснове спектральных представлений элементов гравитационного имагнитного полей.Напомним, что для абсолютно интегрируемой функции f(x),удовлетворяющейусловиямДирихле,существуетследующеепреобразование, называемое прямым преобразованием Фурье:fˆ (ω ) =∞∫ f ( x )e− iω xdx .−∞Функция fˆ (ω ) носит название спектра функции f(x), а параметрпреобразования ω обычно сопоставляют с частотой (временной илипространственной).По спектральной функции fˆ (ω ) можно определить исходнуюфункцию f(x) с помощью обратного преобразования Фурье:1f ( x) =2π∞∫ fˆ (ω )eiω xdω ,−∞которое в этом разделе и будет представлять для нас интерес.Как уже отмечалось в предыдущей лекции, это связано с тем, что вразные годы были предприняты попытки создания алгоритмов решенияпрямых задач по расчету элементов аномального гравитационного имагнитного полей через их спектральные функции, аналитическиепредставления которых могут быть легко получены для широкого классаобъектов.

Однако, большинство из предложенных алгоритмов не далижелаемого результата.13612. Отметим основные причины, которые не позволили создатьэффективных алгоритмов по решению прямых задач при таком подходе.Для этого используем модель бесконечной горизонтальной материальнойлинии и, создаваемую ею, вертикальную компоненту поля силыпритяжения:g z ( x , z ) = Vz ( x , z ) = 2Gδ лζ −z,(ξ − x ) 2 + (ζ − z ) 2где (ξ,ζ) – координаты положения материальной линии, (x,z) – координатырасчетной точки.

Будем предполагать, что ось oZ направлена вниз, и полебудем рассчитывать на оси абсцисс (z=0), тогда соответствующаяспектральная функция будет иметь вид:− ω ζ − iωξgˆ z (ω ) = 2πGδ л ee.13. Поместим нашу модель в точку с координатами (0, ζ). Выражениедля спектра поля силы тяжести приобретет вид:−ωζgˆ z (ω ) = 2πGδ л e.Из этого выражения следует, что спектральная функция такой моделиимеет только действительную часть, а ее мнимая часть равна нулю. Изэтого же выражения видно, что абсолютные значения спектральнойфункции убывают по экспоненциальному закону с увеличением частоты,значения которой меняются в пределах от −∞ до +∞.

Чем глубжерасположен источник, тем быстрее происходит затухание спектральнойфункции.13714. Рассчитаем по полученным спектральным характеристикам значенияфункции g z ( x ) в точках, расположенных на оси oX, с шагом 0.1 км. Длятого, чтобы по спектральной функции gˆ z (ω ) вычислить значения функцииg z ( x ) необходимо от интегральной формулы обратного преобразованияФурье перейти к ее численному аналогу.

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

Сравним результаты для случаев, когда для расчета1024поля будем использовать значения спектра до частоты, равной π (число Nω= 512), и до частоты 2π (Nω = 1024).Полученный результат дает наглядное представление об эффектеГиббса, т.е. ограничение числа суммируемых гармоник приводит куклонению решения от точного. Этот эффект будет присутствовать всегда,но может проявляться в большей или меньшей степени. Так на следующемрисунке представлены результаты расчетов для модели, расположенной наглубине 0.5 км.138Как видно полученного результата совпадение полей с увеличением числагармоник и глубины положения модели улучшается.15.

Пусть наш источник располагается в точке с координатами (ξ,ζ). Вэтом случае выражение для спектральной функции поля g z ( x ) приобрететвид:− ω ζ − iωξgˆ z (ω ) = 2πGδ л ee.Появление сомножителя e − iωξ приведет к тому, что спектральнаяфункция gˆ z (ω ) будет иметь как действительную, так и мнимую части,которые будут представлять собой синусоиды, затухающие по амплитуде сувеличением частоты по экспоненциальному закону. На следующемрисунке представлен характер поведения элементов функции gˆ z (ω ) дляслучая расположения источника в точке с координатами (10., 0.5).139Частота осцилляций будет тем выше, чем дальше будетрасполагаться источник от начала координат.

В свою очередь это приведетк тому, что необходимо будет увеличивать дискретизацию представленияспектральной функции по частоте с тем, чтобы правильно описатьхарактер ее поведения.16. К чему может привести неверное задание шага дискретизации почастоте, иллюстрирует следующий пример. Пусть шаг дискретизации по2π, и пусть расчетные точки располагаются начастоте равен ∆ω =128расстоянии ∆x = 1 друг от друга. Число точек на профиле будет равно 128.расположим наш источник в точке с координатами (135., 5.), т.е. запределами профиля.