авторефераты диссертаций БЕСПЛАТНАЯ БИБЛИОТЕКА РОССИИ

КОНФЕРЕНЦИИ, КНИГИ, ПОСОБИЯ, НАУЧНЫЕ ИЗДАНИЯ

<< ГЛАВНАЯ
АГРОИНЖЕНЕРИЯ
АСТРОНОМИЯ
БЕЗОПАСНОСТЬ
БИОЛОГИЯ
ЗЕМЛЯ
ИНФОРМАТИКА
ИСКУССТВОВЕДЕНИЕ
ИСТОРИЯ
КУЛЬТУРОЛОГИЯ
МАШИНОСТРОЕНИЕ
МЕДИЦИНА
МЕТАЛЛУРГИЯ
МЕХАНИКА
ПЕДАГОГИКА
ПОЛИТИКА
ПРИБОРОСТРОЕНИЕ
ПРОДОВОЛЬСТВИЕ
ПСИХОЛОГИЯ
РАДИОТЕХНИКА
СЕЛЬСКОЕ ХОЗЯЙСТВО
СОЦИОЛОГИЯ
СТРОИТЕЛЬСТВО
ТЕХНИЧЕСКИЕ НАУКИ
ТРАНСПОРТ
ФАРМАЦЕВТИКА
ФИЗИКА
ФИЗИОЛОГИЯ
ФИЛОЛОГИЯ
ФИЛОСОФИЯ
ХИМИЯ
ЭКОНОМИКА
ЭЛЕКТРОТЕХНИКА
ЭНЕРГЕТИКА
ЮРИСПРУДЕНЦИЯ
ЯЗЫКОЗНАНИЕ
РАЗНОЕ
КОНТАКТЫ


Pages:     | 1 || 3 |

«МИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИ РОССИЙСКОЙ ФЕДЕРАЦИИ РОССИЙСКИЙ ГОСУДАРСТВЕННЫЙ ГЕОЛОГОРАЗВЕДОЧНЫЙ УНИВЕРСИТЕТ ...»

-- [ Страница 2 ] --

2 x 2x В частности, при h0 = 0, z 0, k0 = 0, µ = µ0, компоненты вектора E в произвольной точке M(x,y,z) в земле описываются формулами [Юдин, Киселев, 1985]:

I 1 k A I K + A I K + A I K + A I K + Ex = 4 1 00 0 0 10 1 0 01 0 1 11 1 (2.1.2.20) 2 3 + 3k1R + k1 R 2 1 + k1R k1R + z e, 2 R k R 1 где ( ), ( ) k1zy2, z 3 y 2 R2 ( R + z ) 3 y2 R A= A= 00 01 k R R4 R ( ) + 2 x2 z, ( ) k1zy2 ;

z 3 y2 + R2 ( R z ) 3 y 2 R A= A= 10 11 k R R4 R 2r 2 R I xy 1 k B I K + B I K + B I K + B I K ].

Ey = 4 R5 1 00 0 0 10 1 0 01 0 1 11 1 Здесь B = 3( z R) + k 2 zR2, B = 3k zR, 00 1 10 B = 3( z + R) k 2 zR2, B = 3k zR + 3k zR3 / r 2.

01 1 11 1 Формула для расчета вертикальной компоненты вектора E приведена в работе [Ваньян, 1965]:

I 3 xz ( ) 1 + k R + k R 2 / 3 ek1R.

Ez = 1 (2.1.2.21) 2 R5 1 Расчетные формулы для компонент вектор-потенциала и полей E(–i, r, z) и H(–i, r, z) существенно упрощаются, если электрический диполь и точки наблюдения совпадают с поверхностью земли (z = 0, h0 = 0). Компоненты векторов E и H вычисляются по следующим формулам:

I x 2 1 3 2 + ek1r 1 + k r, ) ( Ex ( i, r,0 ) = (2.1.2.22) 2r 3 r 2 I 1 xy, E ( i, r,0 ) = 0, E y ( i, r,0 ) = (2.1.2.23) z 2r ( I0K1 I1K0 ), Iµ xy 8 I K k r H x ( i, r,0 ) = 4 11 4r I µ 1 4 y 2 / r 2 k ) ( x 1 I K I K H y ( i, r,0 ) = IK, 2 11 r 2 01 r2 ( ).

3I µ x kr 1 + k r + k 2 r 2 / H z ( i, r,0 ) = 1 e 1 2r 2 r ) ) ( ( Здесь I = I k r / 2, K = K k r / 2.

1 Б. Верхнее полупространство (z 0). Если же источником поля является магнитный диполь, то исследование поведения поля в рассматриваемой области представляет практический интерес для аэроэлектроразведки.

Изучение полей в воздухе, возбуждаемых электрическим диполем, не представляют практического интереса. Полезно получить формулы, когда и источник поля и приемник расположены в проводящей среде (земле). В этом случае нужно принять, что индекс 1 соответствует воздуху, а индекс 0 – земле.

Формулы для расчета компонент вектор-потенциала и скалярного потенциала в воздухе для горизонтально-однородной слоистой среды имеют вид:

I µ e p0|h0 + z| p p / R* p h p z 0 0 J (r ) d, C = e 0 0, A= +C e x0 4 0 p p p + p / R* 0 0 0 0 1 Iµ pz A=0 Z (0) e 0 J ( r ) d = z0 4 x 0 0 Iµ 2 p z =0 V0 (0) X 0 (0) / e 0 J 0 ( r ) d, 4 x 0 k I µ i p z V (0) 0 X (0) e 0 J ( r )d.

U (r, z ) = 4 k 2 x 0 0 2 0 0 В однородном полупространстве R* =1 и коэффициент C принимает вид p p 0 1 e p0h0.

C= ) ( 0 p p +p Рассмотрим вычисление интегралов, через которые выражаются потенциалы.

1) Компонента A в области z 0.

x Преобразуем часть выражения, входящего в C, ) ( 2 + k p p p p 01 1 = 01= = p0 2 p1 + ( ) ( ) p p 2 p 2 k0 k1 p p p +p 001 00 p 2 + (k 2 k 2 ) 2( p p ) 1 0 = p 2p + 0 1 0 1.

2 k2 0 1 2 k p p k k 0 1 0 Примем = 0 k = 0 p =.

1 1 При сделанных предположениях 2( p ) p p ( ) 2 0 1 = p0 p0 + k =.

) ( k2 k2 k2 k p p p p +p 0 001 01 После таких преобразований интегралы, посредством которых представляется компонента A, могут быть выражены через элементарные функции и x модифицированные функции Бесселя. Действительно, 2 + k 2 e p0 (h0 z ) J ( r ) d = ( ) pz 0 J (r ) d = C0e p0 p0 2 p 0 0 k0 k 0 2 3 2 S (, k0 ) (, k0 ) 2 (, k0 ) S (, k0 ), := h0 z.

= + k 2 2 k 0 (2.1.2.24) Итак, компонента A в области z 0 может быть рассчитана по формуле:

x Iµ { A = 0 S (| h + z |, k ) S (, k ) + x0 4 0 0 2 2 S ( 0, k0 ) ( 0, k0 ) 2 ( 0, k0 ), 0 = h0 z.

+ + k k 2 2 3 0 0 0 Если источник в земле на глубине h0 в полупространстве с волновым числом k1, то для получения аналога последней формулы нужно сделать замену k0 на k1 и z на –z. В результате получим:

Iµ { A = 0 S (| h z |, k ) S (, k ) + x1 4 0 1 2 2 S (1, k1) (1, k1) 2 (1, k1),1 := h0 + z.

+ + k 2 2 k 1 1 1 Компонента Ах вектор-потенциала на постоянном токе равна Iµ, R = x 2 + y 2 + ( z h)2.

Ax (r, z, h) = 4 R В. Компонента A в области z 0. Полагая = 0 и p =, x0 1 k p V (0) + 1 X (0) = V (0) = 0, 2 X (0) p h p h ) ( 2 1 = e 0 0 = p e 0 0, p + k2 2 0 Iµ 1 p z A = 0 X 0 (0) e 0 J 0 ( r ) d = 4 x z Iµ p (h z ) 0 p e 0 0 J ( r ) d.

= 2 x 0 0 2 k Формулы Зоммерфельда и Фока дают I µ 3(, k ) 2 S (, k ) 0 0 0, := h z.

0 0+ A= (2.1.2.25) x z0 2 k 2 x 2 0 Если источник в земле на глубине h в полупространстве с волновым числом k1, то нужно поменять местами k и k1 заменить на. В результате 0 0 получим:

I µ 3(, k ) 2 S (, k ) 0 1 1, := h + z.

11+ A = x z1 2 k 2 x 2 1 1 Компонента Аz вектор-потенциала на постоянном токе равна Iµ x z + h 2 2 Az (r, z, h) = 1, R+ = x + y + ( z + h).

4 r 2 R+ Г. Скалярный потенциал в области z 0. Источник в нижнем полупространстве.

В этой модели 1 Ax Az { S (| h z |, k ) S (, k ) + U = divA = + = q µ x µ µ x 1 0 1 z 3 2 (1, k1) S (1, k1) 2 (1, k1) + + k k 2 3 1 1 1 3 2 (1, k1) S (1, k1) +.

x z k 2 x 2 1 1 После преобразования получим I S (| h z |, k ) S (, k ) 2(, k ) 1 1 1.

0 1 1 1 U ( z) = (2.1.2.26) 4 x 1 x x Найдем интегральное представление для потенциала погруженного в нижнее полупространство источника. Найдем скалярный потенциал U :

Ax Az q ( X + Z ) J 0 ( r )d.

U = divA = + = µ µ x µ x 1 z p p p1 z h 1 0 e p1 z +h, Z = p Z = p p z +h e X1 = e +.

) ( 1 1 1 1p +p p p p +p 1 11 ( ) p1 p0 2p 1 = 2.

+ ) ( p +p p p p +p 10 11 Таким образом, q p (h + z ) S (| h z |, k ) + 2 e 1 J ( r )d = U = µ x 1 0 1 0 p 1 2(, k ) I S (| h0 z |, k1) S (, k ) 1 1 1.

1 1 = 4 x x x На постоянном токе получим x 1 1 q, U = + =q + µ1 x R R+ µ1 R R 1 + U 2 1 = q 1 + 1 = q 1 1 3x + 1 1 3x, µ x2 R R+ µ R R R3 R x 3 1 1 + + U 1 = q 1 + 1 = q 3xy 1 + 1.

µ1 x 2 R R+ µ1 R R x + Полагая = k = 0 и p =, получим [Ваньян, 1965] 11 p z k 2 p0|h0 + z| k I µ i V (0) 0 C e 0 0 e 0 J0 ( r ) d.

U= 0 4 k 2 x 0 0 2 0 2 p 0 4. Электрическое поле с дипольным источником в нижнем полупространстве.

В соответствии с формулой E = i A gradU, имеем U U U E ( x, y, z ) = i A 1 ;

E ( x, y, z ) = 1 ;

E ( x, y, z ) = i A 1.

x1 x1 x y1 z1 z1 z y Компоненты вектора Е вычисляются посредством формул (2.1.2.24)- (2.1.2.26).

I µ i E ( x, y, z ) = 0 S (| h z |, k ) S (, k ) 4 1 x1 0 I 2 S (, k ) 3(, k ) 2 (, k ) 1 1 1 + 11+ 1 1 k (2.1.27) 2 2 3 1 I 2 S (| h0 z |, k1) 2 S (, k ) 3(, k ) + 1 1 1.

2 2x2 2x2 x Отметим I µ i 0 = k1 I 4 2 I 2 S (| h z |, k ) 2 S (h + z, k ) 3(, k ) E ( x, y, z ) = 1 1 1, (2.1.28) 0 1 0 xy y1 2xy 2xy I 3(, k ) 2 S (, k ) E ( x, y, z ) = 1 1 1 + 11+ xz z1 xz I 2 S (| h z |, k ) 2 S (, k ) 3(, k ) + 1 1 1.

0 1 2 2xz 2xz xz 2 После преобразований получим I 2 S (| h z |, k ) 2 S (, k ) E ( x, y, z ) = 1 1 1.

0 1+ (2.1.29) xz z1 xz Найдем производные, необходимые для вычисления электрических полей 2 S ( z, k ) 3 + kR + k 2 R 2 1 + kR = x2, R x 2 R 2 S ( z, k ) 3 + kR + k 2 R = xy.

xy R Далее k z z = 1 1 I K + 1 + I K = 1 0 0 z x 2 x R R k xz ) ) ( ( z z 1 I K I K + 1 x I1K0 + 1 + R x I0 K1.

3 1 0 0 2 R R С учетом правил вычисления производных от модифицированных бесселевых функций I = I, I ( x) = I ( x) I ( x) / x;

K = K, K ( x) = K ( x) K ( x) / x 0 11 0 1 0 11 0 найдем ) ) ( ( kx x IK = I K I K IK, 1 0 2 R 0 0 1 1 R( R z ) 1 x ) ) ( ( kx x I K = I K I K IK.

01 0 0 1 1 R( R + z ) 0 x 2R Подстановка последних соотношений и простые преобразования дают k xz kx z x = 1 I K I K + 1 I K I K I K + I K, z 2 R3 1 0 0 1 R 1 1 R 0 0 R 2 1 0 0 x k1x z k1 xz ) ( z = 1 I K + 1+ I K + I K I K.

x z 2 R 2 R 1 0 R 0 1 2 R2 0 0 1 Сопоставляя последнюю формулу с производной / z, получаем x k1 xz ) ( = + I K I K.

x z R 2 z 2 R 2 0 0 1 Так как х и у входят в функцию / z одинаково, то сразу можем записать y k1 yz ) ( = + I K I K.

y z R2 z 2 R2 0 0 1 Если использовать интегральное представление производной и воспользоваться вариантами интеграла Фока, то можно получить 2 1 2 x 2 2 x 2 2 = 1 + k, 2 z r 2 rz r 2 z z x r 2 xy 2 xy 2 = 2 + k.

xy z r 3 r z r 2 z z Приведенные формулы позволяют рассчитать электрические поля, создаваемые горизонтальным электрическим диполем, погруженным в нижнее полупространство на глубину h0 от поверхности земли.

Отметим, что при h0 = 0 легко убедиться, что расчеты по общим формулам (2.1.27) – (2.1.29) совпадают с вычислениями по более простым формулам (2.1.2.20) – (2.1.2.21) 5. Поле погруженного электрического диполя с постоянным током в нижнем полупространстве При = 0 приведенные выше формулы существенно упрощаются.

Компоненты вектор-потенциала в полупространстве с источником имеют вид [Кауфман, 1997]:

x z+h A ( x, y, z ) = q, A ( x, y, z ) = q 1.

x1 R z1 r R+ Скалярный потенциал равен:

x1 1.

U ( x, y, z ) = divA = q + µ µ R R + Компоненты электрического поля вычисляются по формулам 1 3 x 2 1 3x 2 U ( x, y, z ) + 1 E ( x, y, z ) = =q 1, µ R R R 2 R x1 x 2 + + 3xy 1 U ( x, y, z ), E ( x, y, z ) = = q + µ y1 R5 R y + U ( x, y, z ) z h + h+ z.

3x E ( x, y, z ) = = q µ z1 R5 z R+ Пример. Приведем графики изменения с глубиной действительных и мнимых частей компонент напряженности электрического поля для следующих параметров.

Частоты. Расчеты выполнены на постоянном токе и для двух различных групп частот. Каждая группа состоит из 4 частот. Первая группа начинается с частоты 10 Кгц (Т = 0.0001 с), вторая – с частоты 100 гц (Т = 0.01 с). Далее частоты уменьшаются в геометрической прогрессии со знаменателем 0.5.

Таким образом, имеем частоты в герцах:

• I группа частот: f1 = 10000, f2 =5000, f3 =2500, f4 =1250, f0 = 0;

• II группа частот: f1 = 100, f2 =50, f3 =25, f4 =12.5, f0 = 0.

Индекс кривых на графиках соответствует номеру частоты в группе.

2. Положение источника и приемников. Диполь расположен на глубине 100 м. Точки наблюдения имеют координаты х = 10 м, у = 20 м и переменные значения по оси z, изменяющиеся от 0 до 200 м. Расчеты выполнены на сетке с шагом дискретизации 2 м, количество дискретных значений равно 101. По оси абсцисс отложены номера узлов дискретизации. По оси ординат – значения компонент электрического поля.

3. Удельное сопротивление среды равно 1 омм.

Результаты расчетов изображены на рис. 2.4–2.7.

6. Поле электрического диполя с постоянным током в слоистой среде Устремляя частоту к нулю, получают решение задачи для диполя постоянного тока [Ваньян, 1965]. Ввиду того, что E = i A gradU, при 0 для компоненты Ex получим I x U = 1 RJ1 ( r ) d, Ex = x 2 x r где = cth[ h + arcth 2 cth( h +... + arcth n )].

R := lim 1 n 1 0 R* Рис. 2.4. Поле погруженного диполя. Графики изменения Re Ex с глубиной частотах f1=10 Кгц, f2=5 Кгц, f3=2.5 Кгц, f4=1.25 Кгц, f0=0 гц.

Шифр кривых соответствует номеру частоты в группе.

Меняя порядок интегрирования и дифференцирования, получим U I 1 1 x 2 1 J (r ) d.

Ex = = R + (2.1.2.30) x 2 0 r r r r В согласии с (2.1.2.20) при = 0 для однородного полупространства получаем I x 1 3 1.

E x := 2 r 3 r Относительное кажущееся сопротивление k /1 обычно определяется посредством деления напряженности электрического поля на поверхности горизонтально-слоистой среды на напряженность электрического поля на поверхности однородного полупространства Рис. 2.5. Поле погруженного диполя. Графики изменения Re Ex с глубиной на второй группе частот: f1= 100гц, f2=50 гц, f3 =25 гц, f4 =12.5 гц, f0 = 0 гц.

Шифр кривых соответствует номеру частоты в группе.

Рис. 2.6. Поле погруженного диполя. Графики изменения Im Ex с глубиной на частотах f1=10 Кгц, f2=5 Кгц, f3=2.5 Кгц, f4=1.25 Кгц, f0=0 гц.

Шифр кривых соответствует номеру частоты в группе.

Рис. 2.7. Поле погруженного диполя. Графики изменения Im Ex с глубиной на частотах f1= 100 гц, f2=50 гц, f3 =25 гц, f4 =12.5 гц, f0 = 0 гц.

Шифр кривых соответствует номеру частоты в группе.

k = Ex = r 3 R 1 + x 1 J ( r ) d.

1 Ex 0 r r r r В экваторе диполя при х = 0, у = r I 1 RJ ( r ) d.

Ex = 2 r 0 В методе вертикальных электрических зондирований (ВЭЗ) обычно используются симметричные питающие линии конечной длины (установка Шлюмберже). Напряженность электрического поля E в экваторе такой установки (т.е. при x = 0) находится посредством интегрирования поля диполя, определяемого формулой (2.1.2.24), по длине линии. Если линия имеет длину 2r и момент диполя I = Jdx, то в результате простых преобразований получим:

J +r = Ex dx = 1 R J ( r ) d.

E 0 AB r Следовательно, в методе ВЭЗ кажущееся сопротивление может быть рассчитано путем вычисления интеграла k = r 2 R J ( r ) d.

1 Применяя обратное преобразование Ханкеля, найдем k R= J r r dr.

1( ) 0 1r По сравнению с кажущимся сопротивлением функция R просто связана с параметрами горизонтально-слоистого разреза и используется для интерпретации данных ВЭЗ.

В случае двухслойной среды можно получить выражение для кажущегося сопротивления в виде ряда [Заборовский, 1963;

Ваньян, 1965]. В этом частном случае 2 mh R = cth( h + atcth 2 ) = 1 + k m e 1, k := 2 1.

1 12 + 1 m=1 12 Выполняя почленное интегрирование ряда, получим k = r 2 R J ( r ) d = r 2 1 + 2 k m e2 mh J ( r ) d = 1 1 m= 0 0 3/ 3 k m r 2 + mh () = 1 + 2r.

m=1 12 Здесь использован интеграл Вебера-Липшица z J ( r ) d = e z J ( r ) d = r =, z 0.

e 1 r 0 r r 2 + z 2 3/ ( ) 0 r 2 + z Основная трудность построения численных квадратур вычисления интегралов, содержащих функции Бесселя, состоит в осциллирующем характере подынтегральной функции. Подынтегральное выражение представляет собой произведение двух функций, одна из которых меняется медленно (это функции вида X, Z и их производных), второй сомножитель – это функция Бесселя первого рода нулевого или первого порядка.

7. Поле погруженной линии АВ конечной длины, заземленной на концах.

Будем считать, что электрический ток в кабеле изменяется гармонически, синфазно во всех точках линии АВ. Общий метод решения этой задачи предложен В.А.Фоком [Vock, 1933]. Решение задачи по методу В.А.Фока приведено в монографии [Заборовский, 1960].

Приведем основные формулы этого решения, необходимые для выполнения расчетов.

Пусть кабель имеет вид гладкой дуги и расположен на поверхности земли (на границе двух однородных полупространств, рис. 2.9). В соответствии с правилом суперпозиции токовую линию можно рассматривать как последовательность диполей, расположенных вдоль кабеля. Рассмотрим один из диполей с моментом Jds, находящийся в точке с координатами (,,0) на кабеле. Если поле рассчитывается в точке М с координатами (x,y,z), то расстояние R между выделенным на кабеле элементом и точкой М равно R = r 2 + z 2, r 2 = ( x )2 + ( y )2.

Рис. 2.8. К задаче о поле заземленной линии АВ Обозначим Jds = Jd i + Jd j. d :=| ds | cos(ds, i), d :=| ds | cos(ds, j).

Таким образом, можно рассматривать поле, создаваемое системой из двух взаимно перпендикулярных диполей с моментами Jd, Jd соответственно.

Будем полагать, что магнитная проницаемость среды всюду постоянна и равно магнитной проницаемости воздуха.

В рамках принятой модели условия сопряжения для компонент dAx и dAy вектор-потенциала для системы диполей на границе раздела земля-воздух (z = 0) независимы, поэтому отыскание решения для них можно вести автономно и отдельно от поисков решения для компоненты dAz.

Примем 2 p |z| e J ( r ) d, = 0,1.

P (r, z ) := 0 p0 + p Тогда решение для тангенциальных компонент вектор-потенциала дипольного элемента имеет вид [Заборовский, 1960, с. 140-141] Jµ Jµ d P (r, z ), dAy = d P (r, z ), = 0,1.

dAx = 4 Поле, создаваемое линией АВ, получается в результате интегрирования по длине кабеля B B Ax = q P (r, z )d, Ay = q P (r, z )d, = 0,1.

A A Компонента Az вычисляется по формуле [Заборовский, 1960, с. 143, 145] Az = q Q (r ) Q (r ), = 0,1, B A где rB, rA - расстояния от точки наблюдения до точек В и А соответственно, ( ) 0 ( k 2 p + k 22p )( p p |z| Q (r, z ) := k 2 k 0( r ) d.

e J ) +p 01 1 0 0 Приведем явные выражения для функций P (r, z ) и Q (r, z ).

Вычисляя последний интеграл, находят 2 1 + k0 R 2 3 + 3k0 R + k0 R k0R P (r, z ) = z + e 2 k 2 R 0 R k 1 (2.1.2.31) z k 2 k 2 k r 2 k 2 z 1 + k0 R + 1 + k0r e 0 1 e 1.

+k 1 0 R3 r3 Формула для функции P (r, z ) получается из P (r, z ), если обменять местами 1 величины k и k. Если принять k = 0, то расчет функции P (r, z ) сведется к 0 1 0 вычислению интеграла 2 p z e 1 J ( r ) d, z 0, P (r, z ) := 1 0 p0 + p который был ранее вычислен (см. (1.24)). Поэтому имеем 2 2 S 3 P (r, z ) = + k. (2.1.2.32) 1 2 z 2 z3 1 z k Функция Q (r, z ) при k0 = 0 соответственно равны 2 k1R + z 1 + k1R k1(r + z ) Q (r, z ) +, k =0, e (2.1.2.33) 0 k 2 R3 k r3 1 1 1 2 S p z e 1 J ( r )d = + Q (r, z ) = 2. (2.1.2.34) + p 1 0 2 z z k 0 1 Компоненты электрического поля в нижнем полупространстве.

Согласно общей формулы E = i A gradU имеем U E ( x, y, z ) = i A 1 = i A + divA, x1 µ x x1 x1 x U E ( x, y, z ) = 1 + divA, y µ y y1 E ( x, y, z ) = i A + divA.

z1 µ z z1 С учетом найденных выше решений для компонент вектор-потенциала divA1 принимает вид:

A A A x1 + y1 + z1 = divA = 1 x y z x y d + Q (rB, z ) Q (r, z ).

d + = q P (r, z ) P (r, z ) 1 A AB r 1 1 z AB r r r Составим выражение для компоненты E ( x, y, z ) x 1 B P (r, z ) x B P (r, z ) y B E ( x, y, z ) = i q P (r, z )d 1 d + d + x1 2 x A r x A r r r k A 2 Q (rB, z ) Q (r, z ).

1 A rz 1 Запишем последнее выражение в более удобном для вычислений виде B E ( x, y, z ) = i q P (r, z )d x1 A 2 B P ( r, z ) r 2 ( x ) 1 B P (r, z ) x 1 d + 1 d + 2 A r 2 r A r r k 1 (2.1.2.35) 2 P (r, z ) B B P (r, z) ( x )( y ) ( y )( x ) 1 d 1 d + + r 2 2 A r r r A x ( ) Q (r, z ) Q (r, z ).

+ r zr 1 B 1A Аналогично B E ( x, y, z ) = i q P (r, z )d y1 A 1 B P (r, z ) ( y )( x ) B P (r, z ) ( x )( y ) 1 d 1 d + A r k 2 A r 2 r2 r 1 (2.1.2.36) B 2 P ( r, z ) y 2 B P (r, z ) r 2 ( y ) 1 d + 1 d + 2 r r A r r A y Q (rB, z ) Q (r, z ).

+ r zr 1 1 A Если кабель длиной 2L прямолинейный и ориентирован в направлении оси х и начало координат совпадает с его центром (рис. 2.10), то координаты концов линии АВ равны:

A( L,0,0), B( L,0,0).

В точках этой линии = 0, d = 0, а изменяется в от – L до + L. Кроме того, 2 r = ( x + L ) + y2, r = ( x L ) + y 2.

B A Рис. 2.9. Частный случай расположения заземленной линии АВ.

Формулы для расчета компонент напряженности электрического поля существенно упрощаются:

1 L P (r, z ) x L E ( x, y, z ) = i q P (r, z )d d + x1 1 k 2 L r 2 r L 1 (2.1.2.34) L P (r, z ) 1 x ) ( Q (rB, z ) Q (r, z ).

+ y2 1 d + + r zr 1 1A L r r 3 1 L P (r, z ) y ( x ) L P (r, z ) ( x ) y 1 d E ( x, y, z ) = i q d + y1 r k 2 L r 2 r2 r L y ( ) Q ( r, z ) Q (r, z ).

+ (2.1.2.35) r zr 1 B 1A dE, dE, dE поля диполя происходит Если вычисление компонент x1 y1 x достаточно быстро, то проще найти все компоненты полей посредством численного интегрирования полей по длине питающей линии L E ( x, y, z ) = dE ( x, y, z;

,0,0)d. (2.1.2.36) 1 L Пример. Приведем графики компонент электрического поля, создаваемого линией конечной длины АВ для той же модели среды, для которой были приведены выше графики дипольного источника на двух группах частот и на постоянном токе.

1. Частоты. Расчеты выполнены на частотах в герцах:

• I группа: f1 = 10000, f2 =5000, f3 =2500, f4 =1250, f0 = 0;

• II группа: f1 = 100, f2 =50, f3 =25, f4 =12.5, f0 = 0.

Индекс кривых на графиках соответствует номеру частоты в группе.

2. Положение источника и приемников. Длина питающей линии равна 100 м и погружена на глубину 100 м. Точки наблюдения имеют координаты х = 10 м, у = 20 м и переменные значения по оси z, изменяющиеся от 0 до 200 м.

Расчеты выполнены на сетке с шагом дискретизации 2 м, количество дискретных значений равно 101. По оси абсцисс отложены номера узлов дискретизации. По оси ординат – значения компонент электрического поля.

3. Удельное сопротивление среды равно 1 омм.

Результаты расчетов приведены на рис. 2.10–2.13.

Рис. 2.10. Поле погруженной линии АВ. Графики изменения Re Ex с глубиной на частотах f1=10 Кгц, f2=5 Кгц, f3=2.5 Кгц, f4=1.25 Кгц, f0=0 гц.

Шифр кривых соответствует номеру частоты в группе.

Рис. 2.11. Поле погруженной линии АВ. Графики изменения Re Ex с глубиной на частотах f1= 100гц, f2=50 гц, f3 =25 гц, f4 =12.5 гц, f0 = 0 гц.

Шифр кривых соответствует номеру частоты в группе.

2.1.3. Поле вертикального электрического диполя Поле вертикального электрического диполя (ВЭД) полезно рассмотреть в связи с скважинно-наземной электроразведкой. В рамках этой модификации Рис. 2.12. Поле погруженной линии АВ. Графики изменения Im Ex с глубиной на частотах f1=10 Кгц, f2=5 Кгц, f3=2.5 Кгц, f4=1.25 Кгц, f0=0 гц.

Шифр кривых соответствует номеру частоты в группе.

источник поля помещают в скважину, а измерения выполняют на поверхности земли. Эти решения также могут быть использованы при построении тензора Грина при решении задач геоэлектрики по методу интегральных уравнений.

Рис. 2.13. Поле погруженной линии АВ. Графики изменения Im Ex с глубиной на частотах f1= 100гц, f2=50 гц, f3 =25 гц, f4 =12.5 гц, f0 = 0 гц.

Шифр кривых соответствует номеру частоты в группе.

Будем полагать, что ось диполя коллинеарна оси z. Влияние самой скважины на электромагнитное поле учитывать не будем. В этом случае электромагнитные поля полностью могут быть вычислены, если известно решение для одной компоненты Az (r, z ) векторного потенциала A=(0, 0, Az).

Напомним, что электрическое и магнитное поле связано с вектором потенциалом посредством соотношений B = rotA, E = i (A grad divA), k поэтому 1 2 Az 1 2 Az 1 2 Az, E = i Az Ex =,E =, (2.1.3.1) µ yz z µ xz y µ z 1 Az 1 Az Hx =, Hy =.

µ y µ x Непрерывность Ех, Еу, Нх и Ну на границах разрыва свойств среды будет обеспечена, если на этих границах принять условия сопряжения в следующем виде 1 Az 1 Az = 0, =0.

µ µ z Положим µ = const, тогда условия сопряжения несколько упрощаются:

1 Az Az = 0, =0.

z Очевидно, в однородном пространстве с волновым числом k компонента Az равна kR I µ e p |z h| o Iµ e J (r ) d, Az = = 4 R 4 0 p (2.1.3.2) R := r 2 + ( z h)2, p = 2 + k 2.

1) Вертикальный электрический диполь в нижнем однородном полупространстве (земле).

Пусть волновое число воздуха (в области z 0) равно k, а земли (при z0) равно k (см. рис.2.14).

Тогда решения задачи в воздухе имеет вид pz A = e 0 J ( r ) d, z 0;

z0 0 0 а в земле p z A = A0 + e 1 J ( r ) d, z 0.

Рис. 2.14. Модель среды и источника.

z1 z1 1 Здесь k R p |z h| 0 = q e 1 = q e 1 J (r ) d, A z1 R p 0 Iµ R := r 2 + ( z h)2, q := 1.

Коэффициенты и найдем из условий сопряжения решений на границе 0 раздела земля-воздух. Условия сопряжения на этой границе дают систему p h e 0 1 = q p, p h p 0 + p1 = q e 1, 0 1 0 1 решая которую найдем p h 0e 1 q p1 0 p01 p1h q p h k e 1. (2.1.3.3) 0 = 2q, = = e p + p 1 p p +p p 10 01 1 10 01 Здесь использовано обозначение p n+1 pn p n pn n := n+1 n+1 = n+1.

k nn+1 p n + pn p pn n+1 + n+1 n+ n n+ Следовательно, p1|h z| p (h + z ) +k e 1 J (r ) d.

A =q e (2.1.3.4) z1 01 p 0 1 При 0 = 0 получим p |h z| p (h+ z ) 1 e 0( r ) d A =q [e ]J (2.1.3.5) z1 p 0 и 2 p h e 1 chp z J ( r ) d, z h, 2 A I 0 q z1 = E = (2.1.3.6) 2 I 2 p z µ1 r z r 1 1 chp h J ( r ) d, z h.

2 e В частности, на поверхности земли получим I p h E = 1 2e 1 J ( r ) d, z = 0.

r1 2 0 Вычисления дают I 1 rh 3 + 3kR + k R ekR.

E= r1 2 R На постоянном токе I 3rhI 3 + 3kR + k 2 R 2 kR E = 1 rh 1.

= e (2.1.3.7) r1 2 5 2 R R k = 2) Вертикальный электрический диполь в двухслойной среде А. Диполь в первом слое (рис. 2.15).

Решение в воздухе, первом и втором слое дают формулы (2.1.3.8) Коэффициенты, и найдем из условий сопряжения решений на 0 1 z = 0 и z = H. При z = 0 имеем границах pz A = e 0 J ( r ) d, z 0;

z0 0 0 0 p z +p z A = Az + ( e 1 + e 1 ) J ( r ) d, z1 1 1 0 z H;

p z A = e 2 J (r ) d, z H.

z2 0 2 Рис. 2.15. Модель среды и (2.1.3.8) источника p h e 0 = q p + +, p h p 0 = q e 1 p1 + p1, 0 1 0 1 1 Аналогично, при z = H получаем p ( H h ) p H p H p H e + e 1 + e 1 = e 2, q 1 1 p p ( H h ) p p H p p H p pH e 1 e 1 + 1 e 1 = 2 e 2, q 1 1 1 1 1 Коэффициент 1 равен p h 2 p ( H h ) k e 1 (1 k e 1 ) q 01 =.

2 p H 1p 1) 1 (1 + k k e 01 При 0 = 0 k = –1, поэтому p h 2 p ( H h) e 1 (1 k e 1 ) q =.

2 p H 1 p 1) 1 (1 k e Когда 1 = 2, k12 = 0, выражение для 1 должно совпадать с соответствующим коэффициентом однородного полупространства:

2 p ( H h) p 2H 1 +k e e p h p h 1, = q k e q.

= k e 2 p H 1 p 1 p 1 1+ k k e 01 При 0 = 0 k = –1, поэтому 2 p H 2 p ( H h) p h e 1 e q = k e 1.

2 p H 1 p 1 1 k e Когда 1 = 2, k12 = 0 и 1 = 0. Можно убедиться, что = 0. Подстановка и 0 в (2.1.3.8) дает p h 2 p ( H h) p |h z| e 1 (1 k e 1 ) p z 1 e 1+ A = q [e 2 p H z 0 p1 1 k e 2 p H 2 p ( H h ) p h e 1 e 1 pz +k e 1 e 1 ] J ( r ) d.

2 p H 12 1 k e Постоянный ток, двухслойная палетка Диполь в первом слое 2 = 2**(-10)-2**10, q=1,4141, n= 1 = 1E+ 1E+ 1/ 1/ Er/ Er 1/ 0. 1/ 1/ 0. 1/ 0. 0. 0.1 1 10 1E+002 1E+003 1E+004 1E+ r/ H Рис. 2.16. Постоянный ток. Двухслойная палетка кажущегося сопротивления. Шифр кривых соответствует 2 / 1.

Диполь в середине первого пласта.

При z h I p |h z| E (2) (r, z ) = 1 2[e 1 + 4 r p ( h+ z ) 2 p ( H h) p (h z ) 2 p H 2 p ( H h) e1 1 )+k e 1 1 e (1 k e (e ) 12 12 ]J ( r ) d.

+ 2 p H 1) (1 k e Полагая z = 0, получим 2 p ( H h) I 2 p h 1 k e (2) (r,0) = 1 e 1 12 J (r ) d.

E (2.1.3.9) 2 p H 2 r1 1 k e На постоянном токе 2 ( H h) (2) (r,0) = I 1 2e h 1 k12e J (r ) d.

E (2.1.3.10) 1 k e2 H 2 r1 Примем E (2) (r, z ) E (2) (r,0) k := r1 = r1. (2.1.3.11) 1 E (r, z ) E (r,0) r1 r z = Палетка кажущихся сопротивлений приведена на рис. 2.16. На палетке видим, что в ее средней части при 2 / 1 1 кривые сильно дифференцированы и имеют глубокие фиктивные минимумы, не отражающие распределение сопротивлений в двухслойной модели среды.

Б. Диполь во втором слое (рис.2.17).

Очевидно, в этом случае можем записать pz A = e 0 J ( r ) d, z 0;

z0 0 0 p z +p z A = ( e 1 + e 1 ) J ( r ) d, z1 1 1 0 z H.

p z A = A0 + e 2 J ( r ) d, zH.

Рис. 2.17. Модель среды и источника z2 z2 0 2 Здесь k R p |z h| 0 = Iµ e 2 = Iµ e 2 J (r ) d, A z 2 4 R 4 0 p R := r 2 + ( z h)2, p = 2 + k 2.

2 При z = 0 из условий сопряжения следуют уравнения = +, 0 p p p 0 0 = 1 1 + 1 1.

0 1 Аналогично, при z = H имеем p (h H ) p H +p H p H 1 + e 1 = q e + e 2, e 1 1 p p (h H ) p p H p p H p pH e 1 e 1 + 1 e 1 = q 2 e 2.

1 1 1 1 2 Из системы найдем p H 1k01e p (h H ) = 2q e 2, 2 p H 1) ( p + p )(1 + k k e 12 12 10 p H 1(1 + k01)e p (h H ) = 2q e 2, 2 p H 1) ( p + p )(1 + k k e 12 12 10 p (h H ) p H +p H +p H e = ( e 1 + e 1 q )e 2.

2 1 1 p На основании тождеств 2 p 2 p 1 2, 1 k = 1+ k = 12 p + p 12 p + p 12 12 12 имеем 2 p H p (h H ) p H k12 + k01e q e = e.

2 p H 2 p 2 1+ k k e 10 Выражения для коэффициентов могут быть записаны в несколько ином виде p (h H ) (1 + k )e p1H q e 2 =, 2 p H 1 p 2 1+ k k e 10 p (h H ) (1 + k )k e p1H q e 12 =, (2.1.3.12) 2 p H 1 p 2 1+ k k e 10 p (h H ) (1 + k )(1 + k )e p1H q e 2 12 =.

2 p H 0 p 2 1+ k k e 10 При 1 = p h p h p h q e 2 q e 1 q e = k, = k,=. (2.1.3.13) 2 01 1 01 p p p 2 1 Вектор-потенциал в первом слое равен p z +p z p H I e p2 (h H ) (1 + k )(k e 1 + e 1 )e 12 A = 1 J ( r ) d, (2.1.3.14) 2 p H z1 4 p 0 1) 2 (1 + k k e 10 0 z H.

а во втором 2 p H I p |h z| k12 + k01e p ( z +h2 H ) A = 1 2 e2 J 0 ( r ) d.

+ e 2 p H z 2 4 0 p 2 1+ k k e 10 z H.

Электрическое поле в первом слое равно 2 p z p ( H z ) 1 )e I 2 p (h H ) p (1 + k )(1 k e 12 1 e 2 1 J ( r ) d.

E= 2 p H r1 4 0 p 2 1+ k k e 10 При H z h 2 p H I 2 p (h z ) k12 + k01e p ( z + h2 H ) E = 1 e 2 e2 J 0 ( r ) d.

2 p H r 2 4 0 1+ k k e 10 При z h 2 p H I 2 p ( z h) k12 + k01e p ( z + h 2 H ) E = 1 e 2 e2 J 0 ( r ) d.

2 p H 4 r2 1+ k k e 10 В частности, при z = 0 и 0 = p H p (1 + k )e I p (h H ) E = 1 2e 2 1 2 p H 1 ( ) J r d.

r1 2 0 p 2 1 k e 2.1.4. Поле вертикального магнитного диполя Поле вертикального магнитного диполя, приподнятого на высоту h над землей, описывается одной компонентой Az* вектор-потенциала A *, которая в области с постоянными свойствами удовлетворяет уравнению A* = k 2 A*. (2.1.4.1) z z где k := iµ, Re k 0.

Условия сопряжения на границах пластов для Az* имеют вид * * ] = 0, 1 Az = 0. (2.1.4.2) [ Az µ z Примем Mµ * A* (r, z, ) = a (, z, ) J 0 ( r )d, (2.1.4.3) z 4 0 z относительно a * получают задачу z d 2a* z = 2a* (2.1.4.4) iz dz где m = 2 + km, km = iµ m, -пространственная частота, с условиями сопряжения на границах пластов da* 1 z [a* ] = 0, =0 (2.1.4.5) z µ dz и краевым условиям и условию на бесконечности a* * z z =0 = a0, az 0, z + (2.1.4.6) Здесь величина a должна учитывать присутствие источника в верхнем полупространстве.

Частные случаи.

Пусть диполь находится в начале цилиндрической системы координат.

Запишем расчетные формулы для различных геоэлектрических моделей, опираясь на известные решения и полученные здесь формулы.

1. Однородное пространство Решение задачи имеет вид (Л.Л.Ваньян,1965, А. И. Заборовский,1960) ekR Mµ A* = p, R := r 2 + z 2, p :=. (2.1.4.7) z0 R Интегральное представление (интеграл Зоммерфельда) e |z| k R * = pe, := k 2 + = p A (2.1.4.8) z0 R 0 или |z| * =e a.

z 2. Однородное полупространство.

Пусть диполь находится на высоте h над поверхностью земли (Заборовский, 1960). Тогда в воздухе (z 0) будем иметь 2 h z a* ( z ) = 1 + 0 1 e 0 e0. (2.1.4.9) z + 01 Коэффициент a0 вычисляется по формуле h e 0, a= (2.1.4.10) 0 + поэтому в земле (z0) ( z h) a* ( z ) = a e 1. (2.1.4.12) z 3. Горизонтально-однородная слоистая модель среды.

В пределах m-го пласта конечной мощности с постоянными свойствами решением уравнения (2.1.4.4) является a* ( z ) = a q ( z ) + am q ( z), z := z z, (2.1.4.13) z m1 1m m 2m где a m 1, a m - значения функции a * ( z ) соответственно в кровле подошве пласта, z а sh[m (hm z )] shm z q ( z) =, q ( z) =. (2.1.4.14) shm hm shm hm 1m 2m В подстилающем N-м пласте бесконечной мощности с учетом поведения поля на бесконечности ( zH ) q ( z ) = e N =1, q ( z) = 0. (2.1.4.15) 1m 2m Выбор функций (2.1.4.14) в качестве частных решений уравнения (2.1.4.4) и общего решения в виде (2.1.4.15) автоматически обеспечивает непрерывность 1 da* функции a* ( z ). Требование непрерывности z на границах пластов z µ dz приводит к системе алгебраических уравнений относительно величин a m, m=1,…,N В соответствии с формулой (2.1.4.13) на границе первого и второго пластов должно выполняться равенство 1 (a q (h ) + a q (h )) = (a q (0) + a q (0)). (2.1.4.16) µ1 0 1,1 1 1 2,1 1 µ2 1 1,2 2 2, Если a0 известно, то 1 1 1 a ( q (h ) q (0)) + a ( (0)) = a q q (h ).

1 µ 2,1 1 µ 1,2 2µ 0 µ 1,1 2, 1 2 2 Введем обозначения 1 q (0) = m cthm hm, cm := q (hm ) = µm 2,m µm 1,m µm (2.1.4.17) m 1 bm := q (hm ) = (0) = q.

µm 1,m µm 2,m µm shm hm В этих обозначениях уравнение (2.1.4.16) примет вид a (c + c ) + a b = a b.

11 2 22 Система уравнений для вычисления величин am, m=1,...,N- a (c + c ) + a b = a b, 11 2 22.................

a m1bm + am (cm + cm+1) + am+1bm+1 = 0, m = 2,..., N 2, (2.1.4.18).................

+a + c ) = 0.

a b (c N 2 N 1 N 1 N 1 N Коэффициент cN получается из формул (2.1.4.17) при hN. Система имеет трехдиагональную матрицу коэффициентов и решается прогонкой. В согласии с работой [Ваньян,1965, 1997] в системе (2.1.4.18) коэффициент a0 равен R h h 0 e 0 = e 0.

a = + (2.1.4.19) 0 0 + 1 / R 0 + 0 R Поле в воздухе описывает соотношение / R 2 h z a* ( z ) = 1 + 0 1 e 0 e0. (2.1.4.20) z + / R 01 В формулах (2.1.4.18), (2.1.4.19) функция R выражается через гиперболические функции и зависит от свойств горизонтально-слоистой модели среды (Ваньян, 1965).

Если известен вектор-потенциал A*, то электрическое и магнитное поля вычисляются по формулам (Ваньян, 1965) E = i rotA*, (2.1.4.21) B = k 2 A* + grad divA*.

С учетом того, что A* = (0,0, A* ) получаем z A* E = i z, (2.1.4.22) r 2 A* * * = k 2 A* + z = 1 (r Az ), Bz (2.1.4.23) z r r r z 2* * Az.

Br = (2.1.4.24) zr Компоненты электрического и магнитного поля в m-том пласте равны B* = q (a ( z ))(m k 2 ) J ( r )d = q ( z ) + am q z m1 1,m 2,m 0 (2.1.4.25) = q 2 (a ( z )) J ( r )d, q ( z ) + am q m1 1,m 2,m * Br = q (a q ( z ) + am q ( z )) J ( r )d, (2.1.4.26) m1 1,m 2,m * E = q (a q ( z ) + am q ( z )) J ( r )d. (2.1.4.27) m1 1,m 2,m 2.1.5. Поле токовой линии (кабеля) Электромагнитное поле токовой линии найдем посредством интегрирования поля диполя в направлении оси X в бесконечных пределах.

При этом нужно учесть, что компонента вектор-потенциала A x W Az = x = x r r является нечетной функцией по переменной x, то + V. p. Az dx = 0.

Выражение для компоненты Ax получается в результате преобразования интеграла:

J µ + Jµ + 0 Ax dx = 4 X ( z, ) J 0 ( r ) d dx = 4 X ( z, ) d J 0 ( r ) dx, 0 где внутренний интеграл приводится к табличному:

rJ ( r ) dr 2cos y 22 J 0 ( r ) dx = 2 J 0 ( r )d r y = 2 =.

y r 2 y Окончательно получим формулу для единственной компоненты вектор потенциала при z = 0:

J µ e p0h Ax = 0 cos y d. (2.1.5.1) µp p + 0 µR В произвольном m-том слое нижнего полупространства значения вектора A можно найти численно посредством вычисления интеграла Jµ Ax = 0 X m ( z, ) cos yd, 2 Xm (z,) в котором функция определяется выражением (2.1.1.4).

Коэффициенты, содержащиеся в этой функции, находятся путем решения системы (2.1.1.6) или по формулам (2.1.2.10)-(2.1.2.11) с указанными выше заменами.

Расчет поля кабеля в горизонтально-слоистой среде в общем случае производится численно посредством вычисления несобственных интегралов вида (2.1.5.1).

В случае однородной земли при k0 = 0,-h z 0 формула (2.1.5.1) принимает вид h + z J µ e 0 1 p / h0 + z Jµ 0 I 2 I + 2I, Ax = 0 1e cos yd = + 2 1 k 2 2 0 + p1 где I1, I2, I3 - табличные интегралы h cos y 2 zh I = sh ze d = arth, 1 2 + y 2 + z h 0 ( ) h z2 y h0 z I = e cos yd =, ( ) h z2 + y ( ) ( ) S i S i p h0 z cos yd = 1,1 + 1, I = 1e.

3 0 k2 2 i 2 i ) ( В последней формуле S1,1 – функция Ломмеля, = µ h z + iy.

0 2.2. Аналитические решения для нестационарных полей в горизонтально-слоистой среде 2.2.0. Введение При решении нестационарных задач токами смещения будем пренебрегать. Будем предполагать, что изменение электрического тока в источнике описывается произведением силы тока J на функцию Хевисайда 1(t):

1, t 0, 1(t ) = 0, t 0.

Если форма импульса тока в источнике имеет более сложную форму, то функции, описывающие переходной процесс, можно преобразовать к переходным процессам, соответствующим импульсу Хевисайда [Соколов и др, 196?].

Аналитические решения для нестационарных электромагнитных полей обычно получаются путем применения обратного преобразования Фурье или Лапласа к решениям, полученным в частотной области. Далеко не для всех решений удается получить оригиналы во временной области. Как правило, лишь для простейших моделей это удается сделать. Для нахождения оригиналов полезно пользоваться справочниками по операционному исчислению [Диткин, Прудников,1965].В частности, при использовании справочников по преобразованию Лапласа или Лапласа-Карсона нужно в формулах, являющихся функцией от круговой частоты, сделать замену:

-i на p, если p – параметр преобразования Лапласа F( p) = f(t )e pt dt.

Символом будем обозначать соответствие между изображением и оригиналом:

F(p) f(t) В большинстве случаев приходится прибегать к численному обращению преобразования Фурье.

2.2.1. Одномерные задачи 1. Поле плоской волны. Для функции R при z = 0 для двухслойной модели среды получен оригинал следующим образом [Ваньян, 1965]. Эту функцию можно представить в виде ряда 2 (k h + arcth 2 ) = 1 + Q n e2nk1h = 1 + 4 nQ n e2nk1h, 1 R(i ) = cth 1 n=1 n= где 2 Q=.

+ 2 Затем по таблицам [Диткин, Прудников, 1965] находим:

p erfc, e (2.2.1.1) 2t где erfc(z) - функция ошибок.

Следовательно, 4 n R 2 ( i ) 1+4 nQn erfc, := 10 12 t. (2.2.1.2) / h n=1 Здесь - аналог длины волны во временной области. Этот результат обобщен и на случай произвольного числа слоев [Румшиская, 1975], однако из-за громоздкости выражений мы его не приводим. В этой работе. Приведены палетки кривые кажущегося сопротивления для двух трехслойных моделей:

типа H с параметрами 1 = 3 = 1, 2 = 1/16 и типа K с сопротивлениями 1 = = 1, 2 = 64. Цифры у кривых соответствуют мощности промежуточного пласта.

Сплошной линией изображены частотные графики амплитуды кажущегося сопротивления T МТЗ, пунктиром - графики /1. Отметим, что функция (2.2.1.2) является аналогом относительного кажущегося сопротивления / становления поля электрического диполя в дальней зоне источника (при r ) [Ваньян, 1965].

2. Начально-краевые задачи для полуплоскости и полупространства.

Для случаев E-поляризации и H-поляризации таблица 2.1.1.1 во временной области внешне остается без изменения, если вместо функции G(k,y,z) подставить соответствующий ей оригинал:

2 µ y + z exp K k y2 + z2 4t kz 1 z.

G(k, y, z ) = y2 + z2 y2 + z Трехмерные задачи для аномальных полей при u0. = 1. Для трехмерного случая соответствующие решения для тангенциальных составляющих электрического поля, приведенные в таблице 2.1.1.2, изменяются с учетом соответствия:

µ R R µ R µ (1+kR)ekR erfc + exp.

t 4t 2t 2.2.2. Поле горизонтального электрического диполя.

Вектор-потенциал в однородной среде.

В согласии с формулами (2.1.2.0) и ( 2.2.1.1) получаем I µ ekr r µ Iµ Ax = erfc. (2.2.2.1) 4 r 4 r 2t Электромагнитное поле электрического диполя в проводящем полупространстве (земле) Пусть электрический диполь находится на границе раздела земля-воздух.

Волновое число воздуха и земли равны соответственно k0 и k1.

Диполь расположен в начале декартовой прямоугольной системы координат и его ориентация совпадает с осью ОХ. Ось z направлена вниз.

Электромагнитное поле электрического диполя в присутствии однородного проводящего полупространства выражается через интегралы Зоммерфельда 2 +k 2 |z| J ( r )d = e kR, R := r 2 + z 2, P(i, R) := e 0 R 0 2 + k и Фостера e +k |z| J ( r )d = I ( R z) K ( R + z ).

1 k k Q (i, R ) := 0 0 2 0 0 2 + k2 Для получения оригиналов P, Q функций P и Q воспользуемся таблицами [Диткин, Прудников, 1965]. По ним находим:

R µ P( - i,R) P(p,R) = P(t,R) = erfc, 2t R t1 µ 2 2 µ 2 Q(t,r,z) = exp 8 ( R + z ) I0 8 ( R z ) d.

Q(p r, z) 0 2 Согласно [Заборовский,1972], при k0 = 0 в нижнем полупространстве компоненты Ax и Az вектор-потенциала могут быть записаны через P и Q в следующем виде:

I µ 3Q 2 Q 2 P Ax = +k +, 2 z3 1 z z 2 k I µ 3Q 2P Az = +.

2 x z 2 x z 2 k Обозначим A x (t, R), Az (p,R) A z (t, R), Ax (p,R) тогда I µ t 3Q 2Q 2 P + µ Ax = + dt, (2.2.2.2) 2µ 0 z3 z t z I µ t 3Q 2P Az = + dt. (2.2.2.3) 2µ 0 x z 2 x z На основании последних формул, пользуясь теоремами операционного исчисления, легко находятся соотношения для компонент векторов электрического и магнитного (нестационарных) полей в однородном полупространстве. Так как в частотной области в однородной среде E = i ( A + grad div A), k H= rotA, µ то соответствующие оригиналы E, H равны A E =- + grad div A, (2.2.2.4) t µ H= rotA. (2.2.2.5) µ Отметим, что 2Q div A = 2.

x z В частности, нестационарное поле на поверхности земли при условии, что диполь также находится на границе раздела земля-воздух (т.е. ho = 0) на основе формул (2.2.2.2) - (2.2.2.5) получим I x 2 u u 2) + erfc( ) + u exp( ), ( E x (t,r,0) = 3 r2 2 r I 2 r xy, E z (t,r,0) =0;

u := E y (t,r,0) =, 2 r Iµ xy u u2 u2 u H x (t,r,0) 4exp( ) I ( ) + exp( )I ( ), =- t 4 t r 4 r 4 14 4 0 H y (t,r,0) u Iµ 1 4 y 2 / r 2 u2 u 2 x2 u )I ( ), exp( ) I ( ) exp( =- t 4 t r3 r 4 14 4 0 r2 u2 u H z (t,r,0) 3I u 1 ecfc( ) exp( )(1 + ).

= t 2r 5 2 2 Электромагнитное поле электрического диполя в горизонтально-слоистой среде В общем случае нестационарное электромагнитное поле электрического диполя в горизонтально-слоистой среде находится численно на основе известных решений в области Фурье-изображений. Основой для тестирования этих программ и оценки точности расчетов по ним могут служить приведенные в этом разделе аналитические выражения для компонент поля в пространстве оригиналов. Алгоритмы вычисления интегралов обсуждаются в разделе 2.4.

2.2.3. Поле вертикального электрического диполя Однородное полупространство.

Электрическое поле описывается соотношениями I p h Er = 1 2e 1 J ( r ) d, z = 0.

2 0 и I 3 + 3kR + k 2 R2 kR E = 1 rh e, R = r 2 + h2.

r1 2 R Согласно [Диткин, Прудников, 1965, c.253, формула 23.119] p h e 1 = exp( h2 µ ( p + )) (, t ), µ 1 µ µ 1(, t ) = e herfc(h ) + e herfc(h 1 t t 1 1 + ), (2.2.3.1) µ1 µ 2 t t где p := i. При t или t 2 µ µ1e 1(, t ) sh h.

t С учетом соотношений [Диткин, Прудников, 1965] R µ ekR = exp( R µ p ) erfc( 1 ), c.247,23.89, 1 2t µ 1 exp( R µ1 ), c.248,23.96, k R k e 1 = µ p exp( R µ p ) 1 1 1 t 4t R µ 1 exp( R µ1 ), c.248,23. 2ek1R = µ p exp( R µ p ) µ k 1 1 1 1 2t t 4t нестационарное электрическое поле над однородным полупространством принимает вид R µ R µ I 1 exp( R µ1 )(1 + µ1R ). (2.2.3.2) 1 rh 3 erfc( 1)+ Er (r,0, t ) = 2 R5 t 4t 6t 2t Определение. Назовем функцию R µ R µ 1 exp( R µ1 )(1 + µ1R ) 1 (r,0, t ) := rh 3 erfc( 1)+ Er (2.2.3.3) R5 t 4t 6t 2t основным решением рассматриваемой задачи для полупространства.

Тогда радиальная компонента электрического поля на поверхности земли определяется соотношением:

I Er (r,0, t ) = 1 E1(r,0, t ), z = 0.

2 t При t получим I lim Er (r,0, t ) = 1 rh. (2.2.3.4) t R Как известно, При t 3 R µ R µ1 1 R µ1 R µ 1 ) 1 2 +1 1, erfc( 3 2 t 10 2 t 2 t 2t 2 2 µ 2 µ R2 µ R 2 µ R R 1 + 1 1 1 1) exp( 4t 6 4t, 4t 4t 2 поэтому () 8 R5 µ µ I 1 + O(t 2 )).

1 rh 3 (1 Er (r,0, t ) = 2 R5 15t 2 t I 3 5 32 2 R, := 107 2 t.

Er (r,0, t ) 1 rh 1 (2.2.3.5) 2 15 1 1 R При t получим I lim Er (r,0, t ) = 1 rh. (2.2.3.6) t R Двухслойная среда.

Электрическое поле на поверхности земли описывается соотношением:

2 p ( H h) I 2 p h 1 k e 1 e 1 12 J ( r ) d.

Er (r,0) = (2.2.3.7) 2 p H 2 0 1 k e 2a. В основании лежит изолятор. В этом случае 2 = 0, k12 = 1, поэтому I 2 p h 1 e2 p1( H h) Er (r,0) = 1 e 1 J ( r ) d. (2.2.3.8) 2 p H 2 0 1 e С учетом того, что 2 p1nH =e.

2 p H n= 1 e Подынтегральную функцию представим в следующем виде:

2 p ( H h) p h 1 e p (2nH +h) p1(2(n+1) H h) e1 =e 1 e. (2.2.3.9) 2 p H n= 1 e Обратное преобразование Лапласа последнего выражения имеет вид:

2 p ( H h) p h 1 e 1 (, t ), e 2 p H 1 e µ (, t ) = e (2nH +h)erfc((2nH + h) 2 2 n=0 t µ ) + e (2nH +h)erfc((2nH + h) t t 1 + ) µ µ t 1 (2.2.3.10) µ e (2(n+1) H h)erfc((2(n + 1) H h) t 1 + ) µ t µ e (2(n+1) H h)erfc((2(n + 1) H h) t 1 + ).

µ t Функция erfc( z ) при z имеет асимптотическое разложение e z (2m + 1)!!

1 + (1) m erfc( z ). (2.2.3.11) z m=1 2 )m (2 z В разложении (4) основной вклад в сумму ряда вносят несколько первых его членов, поэтому аргумент функции erfc() может неограниченно возрастать лишь за счет времени t. На основании такого предположения µ t t 1 + erfc + erfc (2nH +h) µ t µ t 1 и t 2 µ µ1e 3 µ t erfc(+ ) 1. (2.2.3.12) µ1 2 t + t Подставляя ((2.2.3.12) в (2.2.3.10), получим t 2 µ µ1e 3 µ (, t ) sh( (2nH + h)) sh( (2(n + 1) H h)) 2 t 2 n= t или t 2 µ µ e 3 µ (, t ) {sh( h) + [ sh( (2 H + h)) sh( (2 H h))] +...}.

2 t t Так как sh( (2 H + h)) sh( (2 H h)) = 2sh ( h ) ch ( 2 H ), то t 2 µ µ1e 3 µ (, t ) {sh( h)[1 + ch( 2 H ) +...]}.

2 t t При h 0 и h H при (, t ) 0.

Двухслойная среда с изолятором в основании.

В этом случае 2 = 0, k12 = 1, поэтому с учетом соответствия 1 E (r, h) Et (r, h) из (2.4) получаем I E r ( r, 0, t ) = 1 Et ( r, 2 nH + h ) E1 ( r, 2( n + 1) H h ).

t 2 n = 0 Обозначим R := r 2 + (2nH + h)2, R := r 2 + (2(n + 1) H h)2.

n1 n При t () 8R5 µ µ 3rh Etq (r, h) 1 1), ( 5 2 t R 15t тогда 3I r Er (r,0, t ) () 5 2 µ 2(n + 1) H h 8R5 µ N 2nH + h 8Rn1(µ1) µ 1 1 n2 1 1 = n=0 R 5 R 15t 2 t 15t 2 t n1 n ) ( µ N 2(n + 1) H h 2nH + h 4 I r µ 1 1 = Er (r,0,0) +.

5 t 2 t R2 R n= n2 n Очевидно, для рассматриваемого частного случая характер асимптотики не изменяется по отношению к асимптотике однородного полупространства I r ( H h)µ µ R2 4nH 2 (n + 1) 1 Er (r,0, t ) Er (r,0,0) + 1 1, R = r + h.

2R 2 t t n=0 R n1 n При t R2µ R2 µ R µ rhµ µ 4t 6t µ 1 ) 2 te exp 1.

1, E (r, h) erfc( t t µ R 2 2tR R µ 4t 2t 1 I r µ µ Er (r,0, t ) = 1 4 t t R 2 µ (2.2.3.13) R 2 µ (2nH + h) (2(n + 1) H h) exp n2 1.

exp n1 4t n=0 R 2 R 4t n1 n На постоянном токе 3I (2nH + h) / r (2(n + 1) H h) / r Er (r,0,0) = (2.2.3.14) 5 3 n= 2 r 2nH + h 2 2 2(n + 1) H h 2 1 + ( 1 + ( ) ) r r Двухслойная среда с идеальным проводником в основании.

В этом случае поле описывает соотношение I Er (r,0, t ) = 1 (1)n ( E1(r,2nH + h) + E1(r,2(n + 1) H h)). (2.2.3.15) t t 2 n= 2.2.4. Поле токовой линии (кабеля) Однородное полупространство.

Нестационарное электромагнитное поле токовой линии над однородным полупространством найдем посредством таблиц обратного преобразования Лапласа. Будем полагать магнитную проницаемость постоянной и равной магнитной проницаемости воздуха, k0 = 0.

1. Поле в земле ( z 0) [Юдин, 1970].

С учетом сделанных упрощающих предположений формула (2.1.4.1) для компоненты Ax вектор-потенциала принимает вид:

h z 2 +k e 0e Jµ cos yd.

Ax = (2.2.4.1) 0 + 2 + k Запишем в более удобном виде часть подынтегральной функции, зависящую от частоты, и по таблицам [Диткин, Прудников, 1965] найдем соответствующий ей оригинал:

z 2+k 2 z2µ ( p+ 2 / µ ) e e = p 2 + k 2 + µ ( p + 2 / µ ) + 2t z 2 µ z µ t exp( z )erfc exp +.

2t µ µ t 4t µ Подставим последнее выражение в (2.2.4.1) и запишем A x (t, y, z ) в виде суммы двух интегралов. После простых преобразований получим:

J µ t e[ z /(2a)] 1 t S dt, a := A x (t, y, z ) = S, µ 1 µ µ 0 t где 2t h cos y d, S = exp µ 1 0 z µ ( ) t = exp (h z ) erfc cos y d.

+ S 2t 2 0 µ 0 Введем обозначения := µ [(h z ) iy], := µ (h iy).

0 Интегралы S1 и S2 табличные [Диткин, Прудников,1965]:

µ 2 * * erfc, exp + exp erfc S= 1 t 2 t 4t 4t 2 t f (,h ) + f ( *,h ).

S= 2 2a h 0 Здесь 2 1 z /2a z µ u erfc exp erfc(u )du.

f(, h ) := exp 1 exp 2a t 0 2 t a 4t t 2. Поле в воздухе (z 0) Нестационарное поле в воздухе найдем на основе формул (2.1.2.29)–(2.1.2.30). В этих формулах интегралы I1 и I2 не зависят от частоты, поэтому при переходе к нестационарным полям они остаются без изменения. Для вычисления третьего интеграла воспользуемся соответствием [Диткин, Прудников, 1965]:

erfc + 1 := (,t ).

p S ( p) exp 1,1 t 4t 2 t В результате получим:

t (h0 z ) y 2 zh Jµ Ax = arth 2 + µ (h z )2 + y 2 + y2 + z h 0 (2.2.4.2) t 1 *, t ) + 2Re dt.

(, t ) + ( + 2 t * 0 Электромагнитное поле кабеля в горизонтально-слоистой среде В общем случае электромагнитной поле в произвольной горизонтально слоистой среде находится численно.

2.3. Некоторые свойства классов функций, связанных с решением задач в горизонтально-слоистой среде При решении прямых задач геоэлектрики нередко применяются различного рода математические операции над векторными и скалярными полями, правомерность которых требует знание тех или иных функциональных свойств Обсуждение рассмотренных в этом разделе вопросов приведено в работе [Юдин, 1970].

Функцию X (, p) := X (, p) = (, p) ( p = i, pm = + km ) + p /R можно рассматривать в качестве функции, заданной в области изображений по Лапласу (или Лапласу–Карсону) (аргумент p ) и Ханкелю (аргумент ). Для вычисления оригиналов по этой функции применяются либо численные методы, либо приемы, основанные на использовании асимптотических свойств интегральных преобразований и функции X (, p). И в том и в другом случае необходимо знать свойства основных классов функций в области изображений A ( r, t ) { } X := X (, p), класса T := x в области оригиналов и класса функций t, являющегося преобразованием Лапласа функций из T = L[T ].

2.3.1. Свойства класса функций X и связанных с ним подклассов Функцию двух переменных X (, p) будем обозначать символом X p ( ), когда будем говорить о ней как о функции переменной с параметром р и X ( p) - когда будем рассматривать ее как функцию от р с параметром.

Наряду с классом X X := { X (, p)} (2.3.1.1) будем рассматривать частные классы { } { } X := X p ( ), X p := X ( p). (2.3.1.2) Напомним (формула (2.1.2.9)), что X (, p) =.

+ p / R(, p) (, p) (формула R Отметим, что функция (2.1.2.8)) бесконечно дифференцируема и имеет не равные нулю асимптотические значения при, p 0 и, p.

Рассмотрим частные классы функций.

А. Подкласс X.

При X p ( ) = O(1), так как в этом случае p / R (, p ).

При 0 X p ( ) = O( ), так как в этом случае p / R (, p) k / R( p ).

1 Таким образом, X p ( ) X C ( 0, ) V ( 0, ), (2.3.1.3) где V ( 0, ) – пространство функций с ограниченной вариацией на интервале (0,).


Б. Подкласс X. К этому классу отнесем функции переменной X (, p ):= X (, p) X (, p) = X p ( ) X ( ).

p Функция X ( ) имеет следующие асимптотические оценки:

p • при X ( ) = O( 2 ), p • при 0 X ( ) = O(1).

p Таким образом, X p ( ) X C ( 0, ) V ( 0, ) L ( 0, ) L ( 0, ). (2.3.1.4) 1 В. Подкласс X. К этому классу отнесем функции переменной, равные разности функций для слоистой среды и полупространства:

X (, p) := X (, p) X (, p) = X p ( ) X p ( ).

+ p Очевидно ( ) p R X p ( ) = + p = R + p + p.

( )( ) + p /R 1 1 Можно показать, что 1 2 h 2 h 1 = O (e 1), • при X p ( ) e R • при 0 X p ( ) = O( ).

k Отсюда следует X p ( ) X C ( 0, ) V ( 0, ) L ( 0, ) L ( 0, ). (2.3.1.5) 1 t В. Подкласс X. К этому классу отнесем функции X (, t ) X (, t ) = C 1 pX (, p ) = C 1 p, 1 t 2 + k + где C 1 – оператор обратного преобразования Лапласа-Карсона.

Путем тождественных преобразований получим a a= 1, =p p µ + 2 + k12 a + a2 + p и в таблицах [Диткин, Прудников, 1965, формула 22.48] найдем 21 2 a X (, t ) = C 1 pX (, p ) = aea t aea t erfc. (2.3.1.6) 1 t t Очевидно, что µ • при X (, t ) = O e, • при 0 X (, t ) = O ( ).

Поэтому X (, t ) Xt C ( 0, ) V ( 0, ) L ( 0, ) L ( 0, ). 2.3.1.7) 1 Таким образом, большинство классов функций, являющихся решениями прямых задач для горизонтально-слоистых моделей среды, в области изображений по Ханкелю принадлежат пересечению пространств ( 0, ) V ( 0, ) L ( 0, ) L ( 0, ) непосредственно, либо путем несложных C 1 преобразований могут быть сведены к классу функций, принадлежащих этому пересечению.

2.3.2. Свойства классов функций T и Основным физическим параметром, с которым имеют дело в структурной электроразведке, является кажущееся сопротивление ( t ), которое пропорционально первой производной по времени от компоненты Ax ( r, t ) при регистрации вертикальной составляющей магнитной индукции Bz ( r, t ). В этом случае Ax ( r, t ) / t при t имеет оценку Ax ( r, t ) = O t 3/2. (2.3.1.8) t Класс функций, являющихся решением телеграфного уравнения (с учетом токов смещения) и имеющих на бесконечности оценку (2.3.1.8), изучался П. П.

Макагоновым [1966]. Воспроизведем далее некоторые результаты, приведенные в цитированной работе.

1. Класс T имеет следующие функциональные свойства ( ) T L ( 0, ) L ( 0, ) B ( 0, ) C r µ,, 1 где B - пространство ограниченных функций, - диэлектрическая проницаемость.

2. Класс функций = L[ T ], получаемый применением прямого преобразования Лапласа к классу T, характеризуется следующими свойствами.

Любая F • аналитическая в правой полуплоскости и непрерывная в замкнутой полуплоскости, • суммируемая в L по прямым, параллельным мнимой оси, • суммируемая в L, ограниченная и равномерно непрерывная на мнимой оси, • стремится равномерно к нулю при p в области Re p 0 и имеет там () оценку O p 1, • L -норма b+i F ( s + i ) d F= bi не только ограниченная, но и убывающая функция от s при Re p = s 0.

3. По приближенному решению F ( p), полученному в области изображений и близкому в метрике L2 ( p p ) к точному решению, применяя обратное преобразование Лапласа, получим приближенное решение f (t ) T, близкое в метрике L (t t ) к точному решению f (t ).

Рассмотренные в этом разделе основные классы функций и их свойства будут полезными и при решении более сложных двумерных и трехмерных задач, связанных с исследований локальных неоднородностей в горизонтально слоистой среде.

2.4. О вычислении интегралов, содержащих функции Бесселя Расчет поля диполя в горизонтально-слоистой среде в общем случае производится посредством вычисления несобственных интегралов, содержащих функции Бесселя.

Существует несколько алгоритмов вычисления подобного рода интегралов. Основные из них базируются на • применении математических фильтров [Андерсен, 1975];

• интерполяции полиномом не всей подынтегральной функции, а лишь медленно меняющегося сомножителя с последующим вычислением получившихся интегралов в аналитическом виде (метод Филона) [Ваньян, 1965];

• сведении несобственного интеграла к ряду, членами которого являются интегралы по интервалам, расположенным между корнями функции Бесселя, с последующим суммированием этого ряда по методу Эйлера;

• аналитическом продолжении подынтегральной функции вещественной переменной интегрирования в комплексную плоскость и выбором в ней оптимального в вычислительном отношении пути интегрирования [Заборовский, 1963;

Табаровский, 1975].

Здесь мы рассмотрим два алгоритма приближенного полуаналитического вычисления интегралов, содержащих функции Бесселя. Эти интегралы будут вычислены таким образом, что в результате будут получены формулы, допускающие дальнейшее вычисление производных по переменной r, соответствующей расстоянию до оси в цилиндрической системе координат.

Первый алгоритм основан на разложении подынтегральной функции в ряд Фурье по системе собственных функций интегрального оператора с ядром в виде функции Бесселя первого рода.

Второй алгоритм базируется на экспоненциальной аппроксимации части подынтегральной функции с последующим аналитически вычислением интегралов от произведения экспоненты на функцию Бесселя.

Каждый алгоритм имеет свой диапазон действия по переменной r.

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

2.4.1. Вычисление преобразования Фурье-Бесселя на основе разложения по собственным функциям интегрального оператора Согласно [Бэйтмен, Эрдейи, 1970,формула 8.9(3)] () () 2 H m me /2 L(n ) 2 := me /2 L(n ) 2 J m ( r )d = m m 0 (2.4.1) () 2 /2 ( m ) = ( 1)n r mer Ln r, Re m 1, m где L(n ) ( z ) – полином Лагерра.

Отметим, что m m L( ) ( z ) = 1, L( ) ( z ) = z + m + 1, 0 ( m ) z = z + 2n + m 1 L( m ) z n + m 1 L( m ) z, n = 2,3,… nLn ( ) ( ) n1 ( ) ( ) n2 ( ) Корни многочленов Лагерра [Справочник… под ред. Абрамовица, 1979, с.594].

Все нули являются действительными, простыми и лежат в интервале (0,).

Примем:

m • z (n) - k-тый нуль полинома L(n ) ( z ) z (n) z (n)... znn), ( 1 k •j k-тый положительный нуль функции Бесселя m,k J m ( z ), 0 j j...

m,1 m, Тогда j (n) m,k z 4 n k (n ) k 2 + 1/ 4 m2 = r + m + 1.

2 + zk r n k k ( ) j2 2 1 + j 2 2m m,k m,k + O(n5 ) 4 1+ 48n n { }n=0 (m) Согласно (2.4.1), система функций dn dn ) = me /2 L(n ) ( 2 ) (m m (2.4.2) представлена собственными функциями интегрального оператора Фурье Бесселя. Эта система функций ортогональна на [0, ) [Янке и др., 1968], причем ( m + 1) C n, m z ( m) m = n, m z e Ln ( z )L( ) ( z ) dz = n+m k m n.

0, 0 Согласно [Сегё, 1962, стр.116, теорема 5.7.1], она (при m -1) замкнута в {} (m L2(0, ). Эти свойства системы dn ) позволяют любую функцию из L2(0, ) приблизить в метрике этого пространства с любой точностью:

() () (m X 2 = ndn ) 2, (2.4.3) n= где коэффициенты Фурье равны скалярному произведению функции X(2) на dn ) (m () n = X, d n = X ( 2 )d n 2 d 2. (2.4.4) {}(m dn ) Система функций приводится к ортонормированной системе n= {}(m dn ), если положить:

n= (m dn ) dn = m.

dn ) ( L2 0, В геоэлектрике чаще всего m = 0 или m = 1, поэтому далее будем полагать, что (m величина m равна именно этим значениям. В этом случае d n ) = n + 1, следовательно (m (m d n ) = d n ) / n + 1, m = 0,1.

{}(m dn ) Образом ортонормированной системы тоже является n= {} (m ортонормированная система ln ), с помощью которой любую функцию n= f(r2)L2(0, ) можно в свою очередь разложить в сходящийся в L2- норме ряд f (r 2 ) = nlnm) (r 2 ), ( (2.4.5) n= где () () (m m ln ) r 2 = r mer /2 L(n ) r 2 / n + 1. (2.4.6) Приведем несколько графиков этих функций и оператор на языке MathCAD для вычисления функций базиса (таблица 2.4.1).

() Покажем, что коэффициенты разложений в ряд по системам { d n ) 2 } и (m () (m { ln ) r 2 } соответственно функций X(2) и f(r2) с точностью до знака равны.

Это утверждение следует из цепочки равенств () n (m n (m (m (m n = X, dn ) = Hm f r 2, dn ) = f,H1 dn ) = f, ( 1) ln ) = ( 1) n.

m В процессе преобразований учтено, что прямое H m и обратное H m преобразования Ханкеля совпадают и подынтегральные функции удовлетворяют теореме Фубини. Итак, получаем:

() () () () (m (m f r 2 = nln ) r 2 = H m X 2 = H m n dn ) 2 = n=1 n=1 () () n (m (m = n H m d n ) 2 = ( 1 ) nln ) r 2.

n= n= Таблица 2.4. Оператор на языке MathCAD Графики функций d n ( ) d (n, ) :

для n = 0, 1,2 (слева) и n = 0, 4, 7 (справа).

Таким образом, если известны коэффициенты Фурье n разложения в ряд (2.4.3) функции X ( 2 ), то преобразование Фурье-Бесселя f ( r 2 ) этой функции можно записать в виде ряда:

() () 1 2 (m) n ( 1 ) n r mer /2 Ln r 2.

F (r ) = f r2 = (2.4.7) n + 1 n= Аналитические решения осесимметричных задач геоэлектрики выражаются через интегралы, содержащие функции Бесселя J0(r) и J1(r). Им соответствуют собственные функции () () () () () 2 d 2 = e 2 /2 L 2, L 2 L 0 2, dn n n n n d n ( 2 ) = e /2 L(n ) ( 2 ) 2 m с собственными числами n = (-1)n.

( m ) ( 2 ) присутствует множитель exp( 2 / 2) = d 0 ( ), В функциях d n который при вычислении коэффициентов Фурье обеспечивает достаточно быстрое убывание подынтегральной функции. Однако при больших значениях индекса n многочлен Лагерра имеет осцилляции, затрудняющие вычисление интегралов. Поэтому алгоритм будет давать лучшие результаты, если ряд (2.4.3) достаточно быстро сходится.

Замечание. В результате вычисления интеграла получаем функциональный ряд, коэффициенты Фурье которого зависят от геометрических и электромагнитных параметров модели. Эти коэффициенты умножаются на собственные функции, зависящие только от величины r (разноса). Такое «разделение переменных» по одним и тем же коэффициентам позволяет:


• вычислять интеграл для всех требуемых разносов, • находить производные по пространственным переменным в аналитическом виде, что бывает необходимым при вычислении различных компонент электромагнитного поля.

Подынтегральные функции X ( z, ), Z ( z, ) не относятся к классу интегрируемых с квадратом функций, поэтому для применения рассмотренного здесь подхода нужно их модифицировать. С этой целью рассмотрим функции вида X a z, = X z, X (0) z,, Z a z, = Z z, Z (0) z,, ( ) ( ) ( ) ( ) ( ) ( ) в которых X (0) и Z (0) соответствуют X и Z для однородного проводящего полупространства (z 0). Легко показать, что X a, Z a L2 (0, ).

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

На простых примерах покажем возможности рассматриваемого подхода для вычисления известных интегралов. Вычисления будут выполняться средствами MathCAD, которые позволяют написать математически наглядные операторы, реализующие вычисления, но весьма ограниченные возможности в управлении вычислениями.

Пример 1. Пусть нужно вычислить интеграл f (r ) = e z J ( r )d.

0 Его точное значение дает формула, следующая из интеграла Вебера-Липшица, f (r ) = e z J ( r )d = e z J ( r )d = 1 z =.

0 z 0 z r 2 + z 2 3/ ( ) 0 r 2 + z Найдем коэффициенты Фурье () n = e z, dn = e z dn 2 d () 2 = e z e /2 Ln 2 d 2 = e z x e x /2 Ln ( x ) dx, 0 () () (m X 2 = ndn ) 2.

n= Окончательный результат дает формула f (r ) = n dn ( r 2 ) = ( 1) n dn ( r 2 ).

n n=1 n= Приведем операторы на языке MATHCAD, реализующие приближенное вычисление интеграла. Комментарии к операторам приведены в левом столбце таблицы 2.4.2, а операторы – в правом столбце.

Таблица 2.4. Комментарии Операторы языка MATHCAD Вычисление базисных функций Данные Коэффициенты Фурье n Реконструкция функции Коэффициенты n Функция для приближенного вычисления интеграла Функция FF(r, z), дающая точное значение интеграла Вычисление несобственного интеграла средствами MathCAD’а Тестовые расчеты. Вычисление коэффициентов Фурье и реконструкцию подынтегральной функции посредством отрезка ряда Фурье. Результаты расчетов посредством MATHCAD приведены на рис. 2.4.2 – 2.4.4. Рис. 2.4. иллюстрирует характер изменения коэффициентов n, n. Коэффициенты которые равны по модулю, но n знакопостоянны, а n – знакочередующиеся.

На рис. 2.4.3 сопоставляются графики функции X ( 2 ) и ее аппроксимации N отрезками рада Фурье X (, N ) = n d n( m ) ( 2 ) при N = 10 и 20. Визуально n = графики совпадают. Рис.2.4.4 показывает сравнение результатов приближенного вычисления интеграла с точными значениями. Здесь также отличия не видны.

Рис. 2.4.2. Коэффициенты Фурье n Рис. 2.4.3. Результат реконструкции экспоненты по 10 и 20 коэффициентам (положительные значения) и n в сравнении с точными значениями n (знакочередующиеся).

(графики совпадают).

Приближенное вычисление интеграла (функции F(r, N)) Рис. 2.4.4. Результаты расчета тестового интеграла по различным алгоритмам при z = 0.5 и различным разносам r (в интервале от 0.01 до 10).

Сравнения точности вычислений на некотором наборе разносов приведены на рис. 2.4.5.

Как показали тестовые расчеты, на относительно небольших разносах точность вычисления интеграла Фурье-Бесселя посредством разложения части подынтегральной функции в ряд по полиномам Лагерра и последующего вычисления интегралов от членов ряда аналитически позволяет получить достаточную высокую точность вычисления интеграла.

Рис. 2.4.5. Сравнение результатов расчета посредством ряда по собственным функциям (F(rri)) и точным формулам (FF(rri), V(rri)) для ряда разносов rr, изменяющихся от 0.01 до 5.12 со знаменателем геометрической прогрессии 2.

Пример 2. Вычислим интеграл Зоммерфельда для вещественного волнового числа. Пусть нужно вычислить интеграл e z + k J ( r ) d.

2 f (r ) = +k 2 Здесь e z + k.

2 X ( 2 ) = +k 2 Точное значение интеграла Зоммерфельда дает формула ek r + z.

2 f (r ) = r 2 + z Найдем коэффициенты Фурье e z + k, d = e z + k d ( 2 ) d 2 = 2 2 2 n = 2 + k2 n n 2 + k2 e z + k e 2 /2 L ( 2 )d 2 = e z x + k e x /2 L ( x ) dx.

2 2 = n n +k x+k 2 2 0 Восстановление функции Х дает ряд Фурье X ( 2 ) = n dn ( 2 ).

n = Окончательно находим f (r ) = (1)n n d n ( r 2 ).

n = Приведем вычисления интеграла Зоммерфельда посредством операторов языка MATHCAD.

Определим исходные данные Данные Вычислим коэффициенты Фурье и реконструируем подынтегральную функцию посредством отрезка ряда Фурье. Результаты вычислений коэффициентов Фурье и аппроксимации подынтегральной функции посредством операторов на языке MATHCAD сведены в таблицу 2.4.3.

Таблица 2.4. Комментарий Операторы на языке MATHCAD и результаты их работы Коэффициенты Фурье Графики коэффициентов Фурье n, n Реконструкция части подынтегральной функции Таблица 2.4.3 (продолжение) Сравнение результатов реконструкции для N = 10 и N = Вычислим интеграл посредством реконструкции его отрезком ряда Фурье. Результаты вычислений интеграла Фурье-Бесселя по разным алгоритмам посредством операторов языка MATHCAD сведены в таблицу 2.4.4.

Таблица 2.4. Комментарий Операторы на языке MATHCAD и результаты их работы Операторы вычисления интеграла Аналитическое Приближенное Оператор вычисления интеграла «в лоб»

Разносы Сравнение результаты расчетов по разным алгоритмам для нескольких «небольших»

разносов rr.

Таблица 2.4.4(продолжение) Разносы Сравнение результаты расчетов по разным алгоритмам для нескольких относительно «больших»

разносов rr.

Таблица 2.4. Комментарий Операторы на языке MATHCAD и результаты их работы Данные Оператор вычисления части подынтегральной функции, ее действительной и мнимой части Вычисление Фурье коэффициентов действительной и мнимой частей Графики Фурье коэффициентов подынтегральной функции (в области изображений) Пример 3. Вычислим интеграл Зоммерфельда для комплексного волнового числа Алгоритм вычислений совпадает с тем, который описан в Примере 2.

1. Вычисление коэффициентов Фурье и реконструкцию подынтегральной функции посредством отрезка ряда Фурье.

Результаты вычислений коэффициентов Фурье и аппроксимации подынтегральной функции посредством операторов на языке MATHCAD сведены в таблицу 2.4.5.

2. Вычисление интеграла посредством реконструкции его отрезком ряда Фурье.

Результаты вычислений интеграла Фурье-Бесселя по разным алгоритмам посредством операторов языка MATHCAD сведены в таблицу 2.4.6.

Таблица 2.4. Комментарий Операторы на языке MATHCAD и результаты их работы Реконструкция действительной и мнимой части.

Оператор вычисления интеграла по формуле Приближенное вычисление интеграла Графическая иллюстрация результатов вычисления интеграла Зоммерфельда по разным алгоритмам 2.4.2. Вычисление преобразования Фурье-Бесселя на основе экспоненциальной аппроксимации Аппроксимация табличных данных двумя экспонентами рассмотрена в книге [Хэмминг, 1964].

Опишем способ аппроксимации экспериментальных данных или функций экспоненциальными полиномами и изложим методику расчета входящих в него коэффициентов, приведенный в дипломной работе В.Сапожникова. На каждое слагаемое линейной комбинации экспоненциальных функций нужно иметь по крайней мере пару экспериментальных точек. Если число экспериментальных данных имеет значительно больший объем, то полезно подобрать аппроксимирующий экспоненциальный многочлен методом наименьших квадратов. Полученный полином дает наилучшее приближение, при котором среднеквадратическая погрешность по всем точкам минимальна. Требование о равенстве аппроксимирующей функции значениям вектора данных в узловых точках снимается.

Пусть задана таблица значений функции и независимого переменного х1 х1 х1 … хn У1 У2 У3 … Уn Необходимо аппроксимировать вектор данных Y = (Y1,..., Yn ) полиномом вида m ax y ( x) = Ai e i. (2.4.8) i = Многочлен содержит m членов. Нужно, чтобы объем выборки n был больше 2m. По вектору данных составим систему уравнений m Y = Ai exp(ai xi ), k = 1,..., n. (2.4.9a) k 1= Пусть данные заданы на равномерной сетке аргументов x = x + (k 1)h, k = 1,2,..., n. (2.4.9b) k где h – шаг независимой переменной.

Введем обозначения ci := Ai exp(ai x ), zi := exp(ai h).

После подстановки принятых обозначений в выражение (2.4.9) получим m Y = ci zik 1, k = 1,..., n. (2.4.10) k i= Решение системы (2.4.10) найдем в два этапа. На первом этапе найдем коэффициенты aк. На втором этапе, используя известные величины ak, найдем коэффициенты Аk.

Первый этап. Для решения системы (2.4.10) предположим, что z1, z2,..., zm являются корнями некоторого полинома Pm ( z ) степени m вида m z m(m1) + Sm = zi Smi = 0, S = 1.

Pm ( z ) = z m + S z m1 + S z m2 +... + S m 1 2 i = (2.4.11) Используя систему (2.4.10), сначала найдем коэффициенты S k (k = 1,..., m) полинома (2.4.11). Затем по известным величинам S k найдем корни этого многочлена.

Вычисление величин S k. Умножим первое уравнение системы (2.4.10) на S m, второе уравнение на Sm 1, …, m-е уравнение на S1, m+1–е уравнение на 1 и сложим первые m+1 уравнений:

mm m c zi S = Y S.

i=0 i +1 mi i=0 j =1 j j mi Меняя порядок суммирования и учитывая тот факт, что z j являются корнями уравнения (2.4.11), получим mm m m m c zi S = c j zij Smi = c j Pm ( z j ) =0.

i=0 j=1 j j mi j =1 i=0 j = Итак, m Yi +1Smi = 0, i= поэтому, имеем уравнение m = Y YS.

i=1 i mi +1 m+ Эту процедуру проделаем далее для всех уравнений системы (2.4.10), начиная со второго. Последний раз это будет возможно сделать, начиная преобразования с уравнения с номером n-m. В результате получим систему из n-m линейных уравнений относительно величин S k :

m Yi Smi +1 = Ym+1, i= m Yi +1Smi +1 = Ym+ 2, (2.4.12) i=...............................

m Yi +nm1Smi +1 = Yn.

i= Пусть B – матрица системы (2.4.12), имеющая размер (n-m)xm, а вектор b R n m – ее правая часть, S R m – искомый вектор. Обозначим Y Y Ym m S 1 m, b =, B :=, S := Yn Ynm Y S n1 тогда систему (2.4.12) можно записать в более компактном виде BS = b.

S (k = 1,..., m) можно найти Система переопределенная. Коэффициенты k методом наименьших квадратов. Как известно, система нормальных уравнений может быть записана в следующем виде BT BS = BT b, где BT - транспонированная матрица B.

Приведем ее в более подробной записи nm nm n m nm Sm Yi2 + S Yi Yi+1 +... + S1 Yi Yi+ m1 = Yi Yi+m, m1 i= i =1 i=1 i= nm n m 2 n m n m Sm Yi +1Yi + Sm1 Yi +1 +... + S1 Yi +1Yi +m1 = Yi+1Yi +m, i =1 i=1 i=1 i=...............................

nm nm n m 2 n m Sm Yi +m1Yi + Sm1 Yi +m1Yi +1 +... + S1 Yi +m1 = Yi+m1Yi +m.

i =1 i =1 i=1 i= (2.4.13) Отыскание корней полинома (2.4.11) и вычисление коэффициентов ak.

Для вычисления корней полинома с вещественными коэффициентами S k (k = 1,..., m) существуют библиотечные программы. Их можно найти, например, в библиотеке IMSL Visual Digital Fortran 6.X. Пусть найдены корни многочлена (2.4.10) и они равны z, z,..., zm. Корни могут оказаться комплексными. Запишем их в показательной форме i z = z e k.

k k Так как z := exp(a h), то величина a равна k k k ( ) a = ln | z | +i / h. (2.4.14) k k Вычисление a, k = 1, 2,..., m завершает первый этап расчетов.

k Второй этап. Вычисление величин Ak. Для определения величин Ak снова воспользуемся методом наименьших квадратов. Экспериментальные данные в узлах сетки будут приближаться функцией следующего вида:

m Y ( x) = A j f j ( x), f j ( x) = exp(a j x).

j = Здесь будем использовать коэффициенты ai, найденные на первом этапе вычислений. Посредством последней формулы на сетке аргументов (2.4.9b) получаем систему уравнений FA = Y, где F = fi, j := f j ( xi ), A := ( A,..., Am ), Y := (Y,..., Yn ).

1 i=1,...,n, j =1,...,m Следовательно, нормальная система уравнений в векторно-матричной записи примет вид FT FA = FT Y. (2.4.15) Если принять Ф:= FT F, g := FT Y и Ф = i, j, g = ( gi ), то i =1,...,n i =1,...,n;

j =1,...,m n n = f ( xi ) f j ( xi ), g j = f j ( xi )Yi, j, k = 1,2,..., m.

kj i =1 k i= Таким образом, решение систем (2.4.13) и (2.4.15) и формула (2.4.14) позволяют вычислить коэффициенты, участвующие в экспоненциальной аппроксимации функции.

Тестирование алгоритма. Вычислим интеграл Фока e z 2 + k 2 k k J ( r ) d = I r 2 + z2 z K r 2 + z2 + z.

F (r, k, z ) = 0 0 2 0 2 + k2 посредством аппроксимации функции e z + k f (, k, z ) = 2 + k линейной комбинацией экспонент a (k, z ) M f (, k, z ) f (, k, z ) = Ai (k, z )e i.

i = Напомним формулу Вебера-Липшица z J ( r ) d = 1/ z 2 + r 2.

e С помощью этой формулы получим m a F (r, k, z) F (r, k, z ) := Ai e i J ( r ) d = 0 i = M a m = Ai e i J ( r ) d = Ai / ai2 + r 2, Re ai 0.

i =1 0 i= 1. Будем полагать k = 0.1, z = 10, m = 3. В таблице 2.4.7 приведены величины коэффициентов Ai и ai.

Таблица 2.4. a(1) (-12.7106,0.0) A(1) (6.05972, 6.648E-015) a(2) (-34.1283, 19.2381) A(2) (-1.1886, 0.1302) a(3) (-34.1283, -19.2381) A(3) (-1.1886, -0.1302) Re ai Видим, что для всех коэффициентов. В таблице 2.3 пары коэффициентов a(2) и a(3), A(2) и A(3) комплексно сопряженные величины, поэтому в сумме они дают вещественные значения:

( A + iB ) e(aib) + ( A iB ) e(a+ib) = 2 Aea cos(b ) 2Bea sin(b ).

В таблице 2.4.8 приведены значения функции f (, k, z ) и ее аппроксимации f (, k, z ) тремя экспонентами, а их графики на рис. 2.1.16.

Таблица 2.4. 0.01 0.02 0.03 0.04 0.05 0.06 0. f (, k, z ) 3.64232 3.5366 3.3718 3.1624 2.924 2.6715 2. 6 18 64 34 078 33 3.64230 3.5366 3.3718 3.1623 2.923 2.6716 2. f (, k, z ) 3 26 81 07 974 35 Рис. 2.1.16. Графики функций f (, k, z ) и f (, k, z ) Согласно рис. 2.1.16 и Таблицы 2.4 три экспоненты удовлетворительно аппроксимируют часть подынтегральной функции f (, k, z ).

В таблице 2.4.9 приведены значения интеграла Фока F (r, k, z ) для нескольких разносов и результат вычисления интеграла на основе использования трех экспонент, а их графики на рис. 2.1. Таблица 2.4. r 1 2 3 4 5 6 F (r, k, z ) 0.414616 0.410349 0.403493 0.394391 0.383461 0.371144 0. F (r, k, z ) 0.419527 0.415124 0.408067 0.398727 0.387547 0.374986 0. Рис. 2.1.17. Графики функций F (r, k, z ) и F (r, k, z ) 2. Приведем результаты расчетов для одиннадцати экспонент (m = 11, z = 10, k = 0.1).

Согласно таблицы 2.4.10, для аппроксимации части подынтегральной функции можно ограничиться линейной комбинацией первых восьми экспонент. Экспоненты с большими номерами не вносят заметного вклада в результат: амплитуды их малы, а действительная часть показателя степени отрицательна и большая по абсолютной величине.

Таблица 2.4. a(1) (-10.64100,0.0) A(1) (1.593879770,3.6868E-9) a(2) (-13.40297,0.0) A(2) (3.156464840,-1.3834E-8) a(3) (-18.47613,0.0) A(3) (3.026429778, 2.4514E-8) a(4) (-36.43866,0.0) A(4) (-4.57653336, -1.8466E- a(5) (-51.88975,0.0) A(5) (-4.607062366,1.2175E-6) a(6) (-62.63337,0.0) A(6) (5.883841654,-1.8568E-6) a(7) (-96.09251, 26.67572) A(7) (-0.406060257,-1.9946E-2) a(8) (-96.09251,- 26.67572 A(8) (-0.40606227, 1.9947E-2, a(9) (-132.83724, 66.85405) A(9) (7.0911569E-3,-4.361E-3) a(10) (-132.83724,-66.85405) A(10) (7.090781E-3,4.3619E-3) a(11) (-357.94843, 0.0) A(11) (-4.393606E-4,9.7535E-7) В таблице 2.4.11. приведены значения функции f (, k, z ) для различных и функции f (, k, z ) – результат аппроксимации линейной комбинацией экспонент Таблица 2.4. 0.01 0.02 0.03 0.04 0.05 0.06 0. f (, k, z ) 3.642326 3.536618 3.371864 3.162434 2.924078 2.671533 2. f (, k, z ) 3.642326 3.536618 3.371864 3.162434 2.924078 2.671533 2. В таблице 2.4.12 приведено сравнение результатов вычисления интеграла по точной F (r, k, z ) и приближенной F (r, k, z ) формулам на основе использования 11 экспонент.

Таблица 2.4. r 1 2 3 4 5 6 F (r, k, z ) 0.419107 0.414705 0.407648 0.398309 0.387131 0.374571 0. F (r, k, z ) 0.419527 0.415124 0.408067 0.398727 0.387547 0.374986 0. Ввиду того, что графики функций f (, k, z ) и f (, k, z ), F (r, k, z ) и F (r, k, z ) при m = 11 визуально совпадают, мы их не приводим.

2.4.3. Алгоритм Андерсена.

Рассмотрим несобственный интеграл, содержащий функцию Бесселя первого рода F (r ) = f ( ) J n ( r )d, n = 0,1. (2.4.16) Сделаем в нем замену переменных r = e x, = e y и умножим обе части равенства на e x. Интеграл (2.4.16) примет вид ( ) ( ) e x F ( e x ) = f e y e x y J n e x y dy.

Обозначая ( x) := e x F (e x ), g ( y) := f (e y ), n ( ) := e J n (e ), сведем интеграл (2.4.16) к интегралу свертки ( x ) = ( g n ) ( x ) = g ( y) n ( x y) dy. (2.4.17) Примем ( ) = F [ ], g ( ) = F[ g ], n ( ) = F[ n ].

Тогда по теореме о Фурье-преобразовании свертки имеет место равенство ( ) = g ( ) n ( ). (2.4.17) Для построения функции n ( ) (ее называют откликом фильтра) используют интегралы r 2 / 4a a / ( 2a ), e J ( r )d = e (2.4.18) r 2 / 4a 2 a () / 4a 2.

e J ( r )d = re (2.4.19) Рассмотрим более подробно построение фильтра для случая n = 0. В этом случае, согласно (2.4.18) имеем r 2 / 4a, ( ) = e J e.

, g ( y) = ea / ( 2a ) ( x) = re =e y 0 r =e x Из равенства (2.4.17) найдем спектр фильтра ( ) ( ) = ( ) / g ( ).

Опасность деления на нуль в последней формуле автор удачно обходит.

Применяя обратное преобразование Фурье, находят ( x) = F 1 ( ) = F 1 ( ) / g ( ).

0 0 При построении ( x) используется алгоритм БПФ, поэтому фильтр получают в дискретном наборе точек, равномерно расположенных на оси переменой х.

{} n Дискретный аналог фильтра ( x) обозначим. Вычисление интеграла 0 0i i =n (2.4.16) сводится к дискретной свертке n 1 f e yi ln r, n,..., n ;

n 1, n 283.

F (r ) = r (2.4.20) i=n 0i 1 21 На основе интеграла (2.4.19) аналогичным образом строится фильтр, 1i с помощью которого выполняется дискретное преобразование Ханкеля при n=1.

Коэффициенты фильтров и приведены в работе [Андерсен, 1975].

0i 1i Глава 3. Электромагнитные поля в цилиндрически слоистой среде 3.0. Введение.

Мы будем обсуждать задачи, направленные на изучение полей точечного источника, электрического диполя и полубесконечной питающей линии, расположенных в скважине, на электромагнитное поле внутри и вне скважины.

Эти задачи возникают при оценке влияния обсаженной скважины на результаты измерений электромагнитного поля во вмещающей среде. К ним примыкают также проблемы изучения искажающего влияния трубопровода на результаты скважинно-наземных измерений полей на нефтегазовых месторождениях. В скважинно-наземной электроразведке нужно изучить влияние системы «скважина-обсадная труба» на результаты измерения электромагнитного поля источников, расположенных в скважине или вне ее (например, линия конечной длины на поверхности земли), на относительно больших расстояниях от вертикальной металлической трубы. В этой ситуации ее влияние можно учесть приближенно, пользуясь аналогиями между обсаженной скважиной и линиями электропередач. Обсуждение проблемы с этих позиций проводилось многими исследователями [Кауфман, 1997;

Могилатов, Гендельман, 1983;

Горюнов и др., 1991, 1993].

Один из важных запросов нефтегазовой промышленности - измерение удельного сопротивления среды вне стальной обсадной трубы. Удельное сопротивление среды – важный индикатор воды, газового и нефтяного насыщения, истощения продуктивной зоны, затопления и т.д.

Большинство работ по изучению влияния обсадной трубы выполнено не с точки зрения скважинно-наземной электроразведки, а с позиций исследования моделей каротажа скважин. Дадим краткий обзор этих работ.



Pages:     | 1 || 3 |
 





 
© 2013 www.libed.ru - «Бесплатная библиотека научно-практических конференций»

Материалы этого сайта размещены для ознакомления, все права принадлежат их авторам.
Если Вы не согласны с тем, что Ваш материал размещён на этом сайте, пожалуйста, напишите нам, мы в течении 1-2 рабочих дней удалим его.