CC BY
32
2007
НАУЧНЫЙ ВЕСТНИК МГТУ ГА серия Математика и физика
№ 114
УДК 528.8.04
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ОТРАЖЕНИЯ ПОЛЯРИЗОВАННОЙ ВОЛНЫ В СЛУЧАЕ СГЛАЖИВАЮЩЕГО ПЕРЕХОДНОГО СЛОЯ
А. В. САМОХИН
Для ситуации, при которой комплексная диэлектрическая проницаемость меняется плавно и даже гладко, без разрывов первой производной диэлектрической проницаемости по глубине (например, свежий снег, лежащий на увлажненном фирне, переходящем в лед) интерпретация результатов дистанционного зондирования в рамках стандартных моделей будет приводить к большим погрешностям. Работа посвящена численному моделированию отражения поляризованной волны в случае переходного слоя.
1. Постановка задачи
Для плоской волны и = elk^‘PxX+PvV+1PzZ+^lt\ где ^ _ (Ех, Еу, Ez) вектор поляризации и где 6 — показатель преломления, волновое уравнение сводится к уравнению Гельмгольца
V2 + k2e(z)~lt + grad (Ё, grad In e(z)) = 0.
Имеется в виду, что разделяющий слой расположен в области 0 ^ z ^ /г; V2 — это оператор Лапласа по ж, у и г. Ориентируя ось х вдоль проекции вектора (Px,Py,Pz) на плоскость z — 0, получаем ру = 0, и падающая волна (в полупространстве г < 0) может быть записана в виде
иПад — Еeik^,xX+PtZ\ При этомр2+р2 = 1. Компонента (0, Еу, 0) соответствует горизонтальной поляризации, a (Ex,0,Ez) — вертикальной.
Сам показатель преломления устроен, как и выше, следующим образом: e(z) = 1 при z < 0, e(z) = вк — const > 1 при z > h, а в промежутке 0 ^ z h это гладкая функция, с условиями непрерывности первой производной на границе разделяющего слоя: e'(z)|2=o =
е'(-г)\z=h = 0.
В этих условиях в [3] были оценены коэффициенты отражения R = \Еотр(9)/Епа^(в)\ (здесь Е ■— напряженность электрического поля, в — угол падения. Решение преобразованной волны в полупространстве z > 0 найдено в виде и = elkpxX(F(z), U(z),H(z)). Функции F, U, G связаны условиями
С & + к'л(е(2) ~ pI)f + ine(2)) =о
| %$ + k2(e(z)-pl)U =0 (1)
[ 0 + **(ф)-й)Я+£(Я£Inе(г)) =0.
Второе и третье уравнения этой системы независимы, и их можно решать по отдельности.
Перейдем сначала к безразмерным переменным. Обозначим £ = z/h и и — kh.
Тогда второе уравнение системы (1) перейдет в
+ *2(е«) - = 0, (2)
где и (С) = U(h() = U(z) и £г(0 = e(h() = e(z).
Для третьего уравнения
т.е. в безразмерных переменных
д"а + ^(е-рі)д + д'Л + д
.ее" - є'2
0.
Чтобы избавиться от первой производной, сделаем замену
У
(4)
(5)
В результате получим
У" +
2єє" - Зе 4е2
/2
+ Ає-рі)
У = 0.
Наконец, для первого уравнения в безразмерных переменных имеем
Т" + и2(е - р\)Т = -гирхд^,
т.е.
Т" + и2(е - р2х)Т = -гирхе'е г У. Решение этого последнего уравнения можно искать в виде
г = р + <3,
(6)
где Р и <3 — общее решение однородного уравнения и частное решение неоднородного соответственно.
Итак, в безразмерных переменных, в которых и будет осуществляться численное решение, имеем систему обыкновенных дифференциальных уравнений на комплекснозначные функции Т, и, У
Т" + и2(е - р2х)Т = -іирхє'є 2 У;
+ ^МС) - ЙМС) = 0;
(П
У" +
2££"-3£ы , т .2 „г',
4^2------Ы' (є -
У = 0.
Здесь ^-(0 = ^(Х) = ^(г), У = и д(С) = <3(АО = <ЭД.
Вводя обозначения для действительных и мнимых частей входящих в систему функций;
Т = к + И;
и = / + гд\
У = ш -I- гп;
е = а + гЬ,
получим следующую систему из шести уравнений;
' Ы'(г) + и2[(а(г) — р2х)К(г) — Ь(г)1(г)] = Ат 4- Вп\
1"{г) + и2\Ь(г)Н{г) + (а(г) - р2)1(г)] = Ст + Бп;
/"(г) + г/2[(а(г) - р2х)/(г) - Ь(г)д{г)] = 0; д"{г) + и2[Ъ(г)/(г) + (ф) - РІ)д(г)} = 0; т"(г) + [и2(а(г) - р2х) + 5]т(г) - [и2Ь(г) 4- Т]п(г) = 0; , п"(г) 4- [и2Ь(г) + Т]тп(г) + [и2(а(г) - р2х) + 5]п(г) = 0,
где
4 — ирх
(а2 + Ь2) 1
В = "Рх (а2 + Ъ2)
С = VРх
(а2 + Ь2)і
в = »Рх
(а2 + Ь2)%
5 = аа" + ЬЬ"
2 (а2 + Ь2)
т - ой” - Ьа"
2 (а2 + Ь2)
\/ V а2 +И2 — а(аа' + ЬЬ') — л/ а2 Л-~№ + а(аЬ' — а'Ь)
Ь2 — а(аЬ' — а'Ь) + у у/а2 + Ь2 + а(аа' + ЬЬ')
і2 + Ь2 + а(аа' + ЬЬ') + у у/а2 + Ь2 — а(аЬ' — а'Ь)
— у у/а2 + Ъ2 + а{аЬ' — а'Ь) + у V а2 + Ь2 — а(аа' + ЬЬ') + ЬЬ" 3 (аа' + ЬЬ')2 — (аЬ' — а'Ь)2
(а2 + Ь2)2
(9)
(10)
(П)
(12)
(13)
(14)
2. Матрица рассеивания
Коротко изложим основные результаты [4,7].
Коэффициент отражения Яу для горизонтально поляризованной компоненты равен
5є' (у) - 4є"(у)(є(у) - РІ) 2ІІ>$ у/е(т)-рІ<іт
------е 0
іб[ф) - р1\
<іу.
(15)
Для волны, падающей вертикально (т.е. для случая рх = 0), коэффициент отражения Иу равен
Ъе'2(у) — 4е"(у)е(у) / У/Ят) ёт
1{у
16є(у)
-е 0
<1у.
(16)
Отметим, что в этом случае Ех = 0 и уравнение на Ех такое же, как и на соответственно
совпадают и коэффициенты отражения для этих компонент: Яу — Кх.
В правой части последнего равенства стоит быстроосциллирующий интеграл, значение которого было асимптотически оценено по методу, близкому к методу стационарной фазы. Эта оценка дает
16и2
е"(0) -
е"(1)
2І" / \Л(Т) <*т
(17)
ІД
Для наклонно падающей волны (15)
е"{1)
1Ьи2
(г(1) - )-гё)
2ги / л/є{т) йт
е"(0)
(18)
Перейдем теперь к третьему уравнению
32Я , и2Г.,., 2\ тг , д, „э
Коэффициент отражения для ^-компоненты равен
о
2е£" + Ъе'2 4 (г — р1)е" — 5е/2 1 ^
4е2(е -р2х)*
Ще-р1)3
2IV / Рг (т) <1т
е 0 йу.
(20)
Быстроосциллирующий интеграл (20) оценивается тем же методом, что и (15). Получается следующая приближенная формула
1Д*
г"(0)(1 ~ 2^(0)Яр,(0)) е"(1)(е(1) - 2^(1)«р,(1)))
2рЗ(0)3?рг(0) 2е(1)р1(Щр2(1))
где р2(0) = у/1 -Р2Х и ря( 1) = у/е(1)-р%.
ж-компонента F есть решение неоднородного уравнения (1)
Я2 тр ,/
— + /с2(е(г) -р2)Р = -гкрхН
или в безразмерных переменных
Т" + г/2(е - р2).?7 = -гирхд—,
т.е.
(21)
.Т7" + /^2(ег - р2)^- = —%ирх£*£ 2 у. Решение этого последнего уравнения представим в виде
Т = Р-+ <5,
(22)
где Р и (3 — общее решение однородного уравнения и частное решение неоднородного соответственно. Поскольку однородное уравнение для (6) такое же, как уравнение на компоненту Еу, то выкладки, повторяющие таковые для (15), приводят к формулам Р = Р0 + где
р _ Угу/р1 ™1Уе(г)-Р2*<1т
V е(0 ~
Рх
И
А(С) =
«Рг\/Р1
2у
I
( у <^1рх(т)<1т } е 0 5е'2(у)—4е" [у)р\ (у) ги^п(С-у)/рг(т)<1т | С ь /
\/рг\у) 16Рг(2/) Л/Рг(0у/Рг(у)
(1у.
В частности, в полупространстве ( < у, где £ = 1, имеем отраженную волну
р1г>( о =
5£'2(!/)-4£»(у)р^!/)
1бр| (у)
2гу / рг (т) (1т
е 0 йу
(23)
и выражение в фигурных скобках дает коэффициент отражения собственно самой компоненты Ех, которое обозначим 11х->х. Как видно из (23), Ях->х = Яу. Однако частное решение неоднородного уравнения вносит дополнительную добавку за счет отражения компоненты
Ег.
Коэффициент отражения компоненты Ег в Ех
1
е"(у) «*Ыт)*г
£{у)Рг {у)ЩРг(у))
е"(0) е"(1) „2“'[Р*{т)Лт
ыЩЩрТШ ~ ?ИШ1не
(24>
1
Отметим, наконец, что отражение в х-компоненте — суммарное; главная часть отраженного сигнала имеет в выбранных осях компоненты, рассчитываемые по формуле
Н-х-лх 0 1 / м 1
0 11у 0 ЕУ (25)
0 0 Я, ) ' V Ел (
Все слагаемые имеют второй порядок малости по и 1, как это видно из формул (24), (17) и (21).
3. Численное моделирование формирования собственного излучения внутри переходного слоя по глубине
Специфика формирования собственного микроволнового излучения для слоев с плавным изменением показателя преломления и отсутствием резких скачков на границе состоит в том, что при дистанционном зондировании с борта воздушного судна виден сигнал, созданный всем слоем в целом (а не сигнал с границы, разделяющей среды с различным показателем преломления, как это обстоит в случае резких скачков). В задаче дифференциации подстилающих поверхностей с целью обеспечения безопасной посадки чрезвычайно важным является понимание процесса формирования сигнала по всей глубине слоя.
Выше предлагались асимптотические (при V —¥ оо) оценки. Коротковолновая асимптотика для коэффициента отражения была получена для слоя, разделяющего два полупространства с постоянными значениями коэффициента преломления и непрерывной производной коэффициента преломления внутри слоя. Удалось получить эффективную оценку матрицы отражения. Как показало рассмотрение примеров, определяющим в формулах для коэффициентов отражения оказывается скачок второй производной на границе воздух <—¥ слой и среднее значение показателя преломления. Однако для получения этих оценок дважды применялся асимптотический переход: при использовании квази-функции Грина и для вычисления быстроосциллирующего интеграла. Естественно, встает вопрос, в какой степени можно доверять полученным таким образом оценкам. Ниже приводятся результаты численного моделирования, позволяющего проверить эти оценки на решениях уравнения Гельмгольца.
Аналитическое решение уравнения Гельмгольца, описывающего движение электромагнитной волны внутри слоя, как правило, невозможно, поэтому в этом разделе предложены численные методы, позволяющие детально прослеживать движение электромагнитной волны внутри слоя и, как следствие, оценить качество оценок. Для численного решения уравнения Гельмгольца был использован алгоритм Рунге-Кутта-Фельберга порядка 4-5, реализованный при помощи пакета прикладных математических программ Мар1е.
В рассматриваемой модели, кроме падающей из воздушной (е = 1, г < 0) среды, в этом же полупространстве присутствует отраженная волна, в то время как за переходным слоем (г > 1) есть только прошедшая и экспоненциально убывающая волна. При пг = РхИ Ф
0 в общем решении при г > 1 для однородной среды
и = Сг ехр(шу/ек -р2х С) + С2 ехр(-г!/у/е* - рЮ
присутствуют как убывающая, так и растущая экспоненты. Поэтому начальные данные задавались при г > 1 в форме Бора-Зоммерфельда с тем, чтобы исключить рост прошедшей волны. При С > 1 имеем С2 = 0 и у/ёь = Щ + т2,
Второе уравнение системы (8) может решаться отдельно, и на его примере покажем, как задавать физически осмысленные начальные условия. Фактически, комплексное уравнение Гельмгольца на Еу было представлено в виде системы двух действительных уравнений второго порядка на действительную / и мнимую д части решения:
Приведенные графики иллюстрируют расчеты при следующих значениях параметров: и = 10, £к = 3 + г рх = 0, т.е. у/ёк = щ 4- гп2 = 2 + 1. Показатель преломления представляет собой одну полуволну косинуса:
ференция прямой и отраженной волны перед переходным слоем. За ним — экспоненциальное убывание вследствие выбора начальных условий по Бору-Зоммерфельду. Заметим, что интерференция отражается в модуле следующим образом. Если егкх +Ле~1к(х+в') — сумма прямой и отраженной волн в полупространстве, где е = 1, то \ё1кх + = |1 + Ле~%2Нх+Ф)|.
Поэтому вокруг 1 происходят колебания с частотой , удвоенной по отношению к частоте падающей волны. Амплитуда колебаний модуля совпадает с коэффициентом отражения. На рис. 4 интерференционная часть дана в увеличенном виде; удвоение частоты видно при сравнении с 2, а амплитуда может быть оценена как примерно равная 0.0089. Асимптотическая оценка для этого е = 3 + г равна 0.007. Сравнение коэффициента отражения, полученного
(26)
где щ = Щу/е - р1).
с (26)):
/(Со) + гд((о) = Сх ехр(й/^/е* -р^Со),
/'(Со) + *^(Со) = ги^Сг ехр(шу/ек ~Р2Со) = й'л/^Жо) + *Жо)) = у(-П2 + »П1)/(Со) + гд((о),
откуда
/Ы = /о, д{*о) = 0о, о) = -^(п2/0 + щд0), д'{г0) = 1/(щ/о - п2д0).
В частности, можно взять
/(-г о) = /о, = ~п2и/0, д(г 0) = 0, д'(г0) = пги/0.
его график для = 3 + г представлен на рис. 1.
Начальные данные следующие
/(1.509) = 0.000013, /'(1.509) = -0.00013, $(1.509) = 0, $'(1.509) - 0.00026.
(Начальная точка взята за переходным слоем, а данные Коши взяты с тем расчетом, чтобы было у(0) = 1, у'(0) = 0).
Рис. 2 иллюстрирует действительную часть решения Еу
На рис. 3 представлен график модуля решения, \Ы(г)\ = у/Р + д2. На нем видна интер-
Рис. 1. Действительная и мнимая часть е
при помощи асимптотики и при помощи прямого численного моделирования в соответствии с изложенной схемой для различных переходных слоев, собрано в таблице (таблица содержит также данные об эффективной безразмерной глубине максимума излучения, полученного при помощи методики, изложенной в [1]).
Как видно из таблицы, различие становится исчезающе малым при больших частотах, в то время как при и = 5 оно достаточно заметно. Это хорошо согласуется с характером асимптотики, взятым при V -» оо. Следует отметить, что в случае больших значений |е| (случай воды в таблице) численный метод становится достаточно чувствительным по отношению к начальным условиям и не очень устойчивым. Тем не менее, при больших частотах и в этом случае наблюдается удовлетворительное согласование оценок.
Рис. 3. График модуля Е(х). Видна интерференция прямой и отраженных волн в полупространстве, где е — 1
Рис. 4. Часть графика модуля. Интерференция прямой и отраженных волн. Видна длина волны и амплитуда, равная коэффициенту отражения
С точки зрения обеспечения безопасной посадки важна не безразмерная глубина, на которой формируется сигнал, а абсолютная (на эту величину нужно вносить поправку в показания высотомера). Ее нетрудно определить, исходя из того, что г/Н = С,, кН = и и, значит, V ■ С = {^/К) = кК. Сомножители в левой части этого равенства берутся из таблицы, и абсолютная глубина, на которой формируется сигнал при выбранном к, дается формулой
* = (*/-С)А»=^рЛ>
где А - длина волны,' £ — безразмерная глубина.
Таблица
Сравнение асимптотической оценки коэффициента отражения с оценкой, получаемой при численном решении при различных сочетаниях и и е
V 4 £ —> Пресный лед, 3 + г Сухой снег, 2 + г•10~3 Мокрый снег, 4 + 0.5г Вода, 60 + 40г
5 Численная модель Асимптотика 0.038 0.03 0.015 0.018 0.0476 0.0365 0.88 0.341
Глубина макс. излучения 0.56 1 0.55 0.18
10 Численная модель Асимптотика 0.0089 0.007 0.0042 0.0037 0.012 0.0095 0.20 0.22
Глубина макс. излучения 0.44 1 0.44 0.15
20 Численная модель Асимптотика 0.0018 0.0017 7.9 ■ 10~4 8.1 ■ 10~4 0.00234 0.00253 0.49 0.55
Глубина макс. излучения 0.36 1 0.35 0.11
40 Численная модель Асимптотика 4.2 ■ 10~4 4.3 • 10-4 1.42 • 10“4 1.42-10-4 5.8- 10-4 5.86 ■ 10-4 0.18 0.14
Глубина макс. излучения 0.28 1 0.32 0.08
80 Численная модель Асимптотика 1.06 • 10“4 1.08-10“4 6 ■ ю-5 6•10~5 1.38 • 10~4 1.46 • 10-4 0.00397 0.0034
Глубина макс. излучения 0.2 1 0.28 0.06
Очевидно, что для выбора поправки бортовой радиолокатор должен быть оснащен компьютером, в который должен быть заложен соответствующий банк данных возможных наблюдаемых показателей типа таблицы с целью оперативного выбора величины поправки.
Вернемся к проблеме задания начальных данных для системы (8).
Третье уравнение системы 8 также может решаться отдельно, и, поскольку при С > 1 выполняется е' — 0, для прошедшей волны оно будет иметь вид
гч Ие-й)]у = о-
Чтобы избавиться от первой производной, сделана замена
в = ^ (28) Vе
Поэтому в области С > 1 исходная амплитуда Я и У отличаются на постоянный (е|^>1 = £*) множитель, который можно спрятать в константы общего решения. Таким образом, для неизвестных функций тип для прошедшей волны начальные условия Бора-Зоммерфельда пишутся ровно в той же форме, что и для пары /, д.
Для первого уравнения в безразмерных переменных
Т" + и2(е - р1)Т = -шрхе'£ 2 У.
(29)
Его правая часть в области ( > 1 обращается в ноль, т.к. е' = 0. Решение уравнения имеет вид Т = P + Q, где Р и Q — общее решение однородного уравнения и частное решение неоднородного соответственно; однако Q|c>i = 0 т.к. s' = 0. Поэтому и в этом случае начальные условия для пары h,l неизвестных функций можно записать в форме Бора-Зоммерфельда точно также, как и в случае уравнения (2).
ЛИТЕРАТУРА
1. Самохин А. В. Отражение от сглаживающего переходного слоя // Электромагнитные волны и системы. Т. 6, № 2-3, 2001.
2. Samokhin А. V. Symmetries of linear and linearizable systems of differential equations // Acta Applicandae Mathematicae, Kluwer Academic Publishers, Dordrecht, Boston, London. V. 56, № 2&3, 1999.
3. Самохин А. В. Отражение поляризованной волны в случае сглаживающего переходного слоя // Научный Вестник МГТУ ГА, серия Математика и физика, № 91, 2005.
4. Самохин А. В. Коротковолновая асимптотика для уравнения Гельмгольца с потенциалом сглаживающего типа // Научный Вестник МГТУ ГА, серия Радиофизика и радиотехника, № 8, 1999.
NUMERIC SIMULATION OF A POLARIZED ELECTROMAGNETIC WAVE PASSING THROUGH A SMOOTHING TRANSITIONAL LAYER
Samokhin A. V.
The interpretation of a distant radio probing using the standard assumption of a permittivity jump at the boundary between different layers results in significant errors in such a cases as an ice with its composition changing with depth or in the cases of dry and wet snow, wet sands and soil. The situation when the permittivity changes continuously and even smoothly up to a first derivative is common enough. This paper deals with the numeric simulation of a polarized probing wave passing through a smoothing layer.
Сведения об авторе
Самохин Алексей Васильевич, 1947 г.р. окончил МГУ им. М.В. Ломоносова (1971), доктор технических наук, профессор, действительный член Московского и Американского математических обществ, Соросовский доцент, заведующий кафедрой высшей математики МГТУ ГА, автор 32 научных работ, область научных интересов — дифференциальные уравнения математической физики, симметрии, законы сохранения.
