Журнал вычислительной математики и математической физики, 2022, T. 62, № 9, стр. 1458-1472

1. ВВЕДЕНИЕ

При численном решении задач рассеяния монохроматических электромагнитных волн на идеально проводящих объектах широкое применение нашли методы, основанные на интегральном представлении электрического и магнитного полей через поверхностные токи [1], [2]. В точной постановке задача сводится к решению интегрального уравнения, записанного на поверхности облучаемых объектов, относительно поверхностных токов. При дискретизации такого уравнения возникает система линейных уравнений с плотной матрицей, причем размерность этой системы определяется числом ячеек разбиения поверхности, а число ячеек разбиения в свою очередь определяется как сложностью геометрии облучаемых объектов, так и требованием малости размера ячеек по отношению к длине волны. Поэтому при применении метода граничных интегральных уравнений вычислительная сложность алгоритма растет с увеличением частоты электромагнитного поля, что ограничивает применение метода для высоких частот.

При численном моделировании рассеяния волн высокой частоты широко применяются давно известные асимптотические методы. Одним из таких методов является метод физической оптики, основанный на приближении Кирхгофа–Макдональда [3]–[5]. В основе этого метода лежит гипотеза о том, что поле, индуцируемое источником излучения, взаимодействует с поверхностью идеально проводящего тела таким образом, что на затененной стороне тела поверхностные токи не возникают, а на каждом освещенном участке тела индуцируются поверхностные токи, такие же, как если бы инициирующее поле воздействовало на касательную плоскость. С математической точки зрения это, по сути, означает, что используется то же самое интегральное представление электромагнитного поля, что и при применении метода граничных интегральных уравнений, однако, само интегральное уравнение не решается, а его решение – поверхностные токи, определяется сразу на основании приближенной формулы, дающей связь между поверхностными токами на ячейках разбиения и инициирующим магнитным полем. В классическом методе физической оптики предполагается, что в качестве этого инициирующего поля выступает первичное поле, в качестве которого может выступать плоская волна, приходящая с бесконечности, или известное поле, возбуждаемое сторонними источниками излучения. Практика показывает, что такой подход применим только при моделировании облучения выпуклых тел.

Более сложная модель возникает при учете переотражений первичного внешнего поля. В этой модели предполагается, что для каждой ячейки, аппроксимирующей поверхность, в качестве инициирующих полей, возбуждающих поверхностный ток на этой ячейке, выступают как первичное внешнее поле, так и поля, индуцируемые поверхностными токами, сформировавшимися на других участках поверхности. В статье дается математическая формулировка такой модели и описывается численный итерационный алгоритм нахождения поверхностных токов. При этом на каждом шаге итераций определение поверхностных токов осуществляется умножением матрицы влияния, описывающей взаимодействие ячеек, на токи, найденные на предыдущей итерации. В результате вычислительная сложность алгоритма при прямом вычислении матрицы влияния становится порядка $\mathcal{O}({{n}^{2}})$, где $n$ – число ячеек разбиения.

Существенного ускорения работы такого алгоритма можно добиться за счет применения быстрых матричных алгоритмов. В настоящей работе к решению задачи осуществлено приложение метода мозаично-скелетонных аппроксимаций [6], [7]. Этот метод нашел широкое применение в различных прикладаных задачах [8]–[10]. В данной работе учитывается специфика аппроксимируемой матрицы, которая состоит в том, что ее элементы определяются аналитическим выражением для ядра интегрального оператора с учетом матрицы “видимости” ячеек разбиения (матрица “видимости” указывает, имеет ли отрезок, соединяющий центры двух ячеек, пересечение с облучаемой поверхностью в своих внутренних точках).

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

2. МЕТОД ФИЗИЧЕСКОЙ ОПТИКИ

Решается задача рассеяния заданного монохроматического первичного электромагнитного поля на системе идеально проводящих объектов (тел и экранов), суммарную поверхность которых обозначим через $\Sigma $. Поверхность $\Sigma $ может состоять из одной или нескольких замкнутых компонент, которые являются поверхностями идеально проводящих тел, и разомкнутых компонент с краем. Каждая такая компонента моделирует идеально проводящий экран.

Напряженности электрического и магнитного полей ищутся в виде ${\mathbf{E}}({\mathbf{x}}){{e}^{{ - i\omega t}}}$, ${\mathbf{H}}({\mathbf{x}}){{e}^{{ - i\omega t}}}$, где ${\mathbf{x}} = ({{x}_{1}},{{x}_{2}},{{x}_{3}})$ – точки пространства, $t$ – время, $\omega $ – круговая частота колебаний.

Для пространственных составляющих электрического и магнитного полей ${\mathbf{E}}$ и ${\mathbf{H}}$ в области вне облучаемых тел справедливы уравнения Максвелла [11, c. 17]:

где ${{\varepsilon }_{0}}$ и ${{\mu }_{0}}$ – электрическая и магнитная постоянные, $\varepsilon $ и $\mu $ – относительные диэлектрическая и магнитная проницаемости окружающей среды. Далее будем полагать, что среда вне тела имеет относительные диэлектрическую проницаемость $\varepsilon = 1$ и магнитную проницаемость $\mu = 1$.

Предполагается, что источником возбуждения является заданное первичное поле $\left( {{{{\mathbf{E}}}_{0}},{{{\mathbf{H}}}_{0}}} \right)$. В качестве первичного поля может рассматриваться, например, плоская волна, для которой

где ${\mathbf{k}}$ – заданный волновой вектор, определяющий направление движения волны, причем ${\text{|}}{\mathbf{k}}{\text{|}} = k = \omega \sqrt {{{\varepsilon }_{0}}{{\mu }_{0}}} $, ${\mathbf{r}}$ – радиус-вектор точки $x$, а ${{{\mathbf{E}}}_{i}}$ – заданный постоянный вектор, определяющий поляризацию волны, причем $({{{\mathbf{E}}}_{i}},{\mathbf{k}}) = 0$.

Задача состоит в отыскании полных электрического и магнитного полей (total fields)

где $\left( {{\mathbf{E}},{\mathbf{H}}} \right)$ – неизвестное вторичное поле.

При этом на идеально проводящей поверхности должно выполняться граничное условие для полного электрического поля [11, с. 17]:

где ${\mathbf{n}}$ – орт вектора нормали к поверхности, а вторичное поле должно удовлетворять условию излучения на бесконечности, которое мы запишем в форме Зоммерфельда [12, с. 110]:

Электрическое и магнитное поля в задаче (1)–(5) будем искать в виде [11, с. 39]:

где ${\mathbf{j}}$ – неизвестное касательное векторное поле на поверхности $\Sigma $. С физической точки зрения ${\mathbf{j}}$ есть поле поверхностных токов.

В теории дифракции большое значение имеют приближенные методы, позволяющие исследовать дифракцию коротких волн на различных телах. Одним из таких методов, применимых к широкому классу задач, является приближение Кирхгофа–Макдональда, известное в электродинамике как приближение физической оптики [3]–[5].

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

В первом приближении мы рассматриваем в качестве инициирующего поля $\left( {{{{\mathbf{E}}}_{{{\text{inc}}}}},{{{\mathbf{H}}}_{{{\text{inc}}}}}} \right)$ для каждого участка поверхности тела основное первичное поле $\left( {{{{\mathbf{E}}}_{0}},{{{\mathbf{H}}}_{0}}} \right)$ вида (2). В случае, если поле ${{{\mathbf{E}}}_{0}}$ есть плоская волна, определяемая формулой (2), мы будем называть точку $x \in \Sigma $ освещенной этим полем, если луч, проведенный из точки $x$ в направлении вектора $ - {\mathbf{k}}$, не имеет пересечений с поверхностью $\Sigma $, кроме самой точки $x$. В соответствии с теорией Кирхгоффа–Макдональда в каждой освещенной точке $x$ поверхности тела $\Sigma $ возбуждается поверхностный ток

Такая модель применима для решения задачи дифракции волны на выпуклом теле (когда поверхность $\Sigma $ ограничивает одну выпуклую область).

В более общем случае предполагаем, что для каждой точки $x \in \Sigma $ в качестве инициирующих полей, возбуждающих поверхностные токи в этой точке, выступают как первичное поле $\left( {{{{\mathbf{E}}}_{0}},{{{\mathbf{H}}}_{0}}} \right)$, так и поля, индуцируемые поверхностными токами, возникшими на других участках поверхности $\Sigma $, которые “видны” из точки $x$.

Если точка $x$ лежит на разомкнутой компоненте поверхности $\Sigma $, то будем предполагать, что в этой точке могут возникать токи ${{{\mathbf{j}}}^{ + }}(x)$ и ${{{\mathbf{j}}}^{ - }}(x)$ на положительной и отрицательной сторонах поверхности соответственно.

Если точка $x$ лежит на замкнутой компоненте поверхности $\Sigma $, то считаем, что в такой точке поверхностный ток возникает только на внешней (положительной) стороне поверхности. Этот ток обозначаем как ${{{\mathbf{j}}}^{ + }}(x)$ и при этом считаем, что ${{{\mathbf{j}}}^{ - }}(x) = 0$.

Предполагается, что в каждой точке $x \in \Sigma $ на каждой стороне поверхности $\Sigma $ возбуждается поверхностный ток в соответствии с формулой (7), где ${{{\mathbf{H}}}_{{{\text{inc}}}}}(x)$ есть сумма первичного поля ${{{\mathbf{H}}}_{0}}$, если точка $x$ освещена этим полем с соответствующей стороны, и магнитного поля, индуцируемого в соответствии с формулой (6) поверхностными токами, размещенными в тех точках поверхности $\Sigma $, которые видны из точки $x$ с рассматриваемой стороны (если условие видимости не выполнено, соответствующее слагаемое не прибавляется). Тогда возникает следующее уравнение для поверхностных токов:

где ${{s}^{ \pm }}(x)$ – параметры освещенности поверхности первичным полем, ${{a}^{ \pm }}(x,y)$ – функции видимости точки $y \in \Sigma $ наблюдателем, находящимся в точке $x \in \Sigma $ на положительной или отрицательной стороне поверхности соответственно. Мы полагаем, что:

• ${{s}^{ + }}(x) = 1$, если точка $x$ освещена первичным полем и при этом $({\mathbf{k}},{\mathbf{n}}(x)){\text{/|}}{\mathbf{k}}{\kern 1pt} {\text{|}} < - \delta $;

• ${{s}^{ - }}(x) = 1$, если точка $x$ освещена первичным полем и при этом $({\mathbf{k}},{\mathbf{n}}(x)){\text{/|}}{\mathbf{k}}{\kern 1pt} {\text{|}} > \delta $;

• ${{a}^{ + }}(x,y) = 1$, если отрезок $[x,y]$ не имеет пересечений с поверхностью $\Sigma $, кроме как в точках $x$ и $y$, и при этом $(y - x,{\mathbf{n}}(x)) > \delta {\text{|}}y - x{\kern 1pt} {\text{|}}$;

• ${{a}^{ - }}(x,y) = 1$, если отрезок $[x,y]$ не имеет пересечений с поверхностью $\Sigma $, кроме как в точках $x$ и $y$, и при этом $(y - x,{\mathbf{n}}(x)) < - \delta {\text{|}}y - x{\kern 1pt} {\text{|}}$;

• коэффициенты ${{s}^{ \pm }}(x)$ и ${{a}^{ \pm }}(x,y)$ равны 0 во всех остальных случаях,

в исходной постановке метода в данных условиях мы полагаем $\delta = 0$.

Соотношения (8) и (9) представляют собой систему интегральных уравнений, которая решается итерационно. На нулевом шаге полагаем ${\mathbf{j}}_{0}^{ + }(x) = {\mathbf{j}}_{0}^{ - }(x) = 0$ для всех $x \in \Sigma $. Далее для каждого $m = 1,2,...$ находим токи ${\mathbf{j}}_{m}^{ + }(x)$ и ${\mathbf{j}}_{m}^{ - }(x)$ по формуле:

Для численной реализации описанного итерационного метода (10), (11) поверхность $\Sigma $ аппроксимируется системой ячеек ${{\sigma }_{i}},$ $i = 1, \ldots ,n$, а для токов строится их кусочно-постоянная аппроксимация: предполагается, что на каждой ячейке ${{\sigma }_{i}}$ токи ${\mathbf{j}}_{m}^{ \pm }(x)$ принимают постоянные значения ${\mathbf{j}}_{{i,(m)}}^{ \pm },$ $i = 1, \ldots ,n$. При этом на каждой ячейке ${{\sigma }_{i}}$ выбирается точка коллокации ${{x}_{i}}$ и рассчитывается вектор ${{{\mathbf{n}}}_{i}}$ – приближение вектора нормали к ячейке. Далее на каждом шаге итерационного процесса дискретные значения токов рассчитываются по формуле:

где коэффициенты $s_{i}^{ \pm }$ и $a_{{ik}}^{ \pm }$ аппроксимируют значения ${{s}^{ \pm }}({{x}_{i}})$ и ${{a}^{ \pm }}({{x}_{i}},{{x}_{k}})$ соответственно. При этом мы учитываем, что поверхность $\Sigma $ заменяется при численном решении приближенной поверхностью $\tilde {\Sigma } = \bigcup\nolimits_{i = 1,...,n} {{\sigma }_{i}}$. Коэффициенты $s_{i}^{ \pm }$ и $a_{{ik}}^{ \pm }$ вычисляются как значения ${{s}^{ \pm }}({{x}_{i}})$ и ${{a}^{ \pm }}({{x}_{i}},{{x}_{k}})$ по описанному выше алгоритму, но с проверкой освещенности и видимости по отношению к поверхности $\tilde {\Sigma }$ и с использованием значения параметра $\delta {{ = 10}^{{ - 4}}}$.

В формуле (12) оператор $R$ вычисляется согласно выражению (6). В основе численного метода вычисления интегралов по ячейкам разбиения лежат квадратурные формулы, описанные в статье [13]. Заметим, что в силу определения функций ${{a}^{ \pm }}(x,y)$ справедливо равенство $a_{{ii}}^{ \pm } = 0$, поэтому в формуле (12) нет необходимости выделять особенность в операторе $R[{{\sigma }_{k}},{\mathbf{j}}_{{k,(m - 1)}}^{ \pm }]({{x}_{i}})$ для $i = k$.

Вычислительная сложность решаемой задачи определяется свойствами поверхности $\Sigma $ и волновым числом $k$. Для хорошей аппроксимации поверхностных токов необходимо иметь разбиение, при котором на длину волны приходится несколько ячеек (вычислительные эксперименты показывают, что в методах описываемого типа нужно брать не менее 6 ячеек на длину волны). Это приводит к необходимости увеличивать число ячеек при уменьшении длины волны, что в свою очередь приводит к увеличению вычислительных затрат. Нахождение значений токов ${\mathbf{j}}_{{i,(m)}}^{ \pm }$ по формулам (12) через значения этих же токов на предыдущем шаге итераций непосредственно на основании записанных формул (12) требует порядка $\mathcal{O}({{n}^{2}})$ вычислений значения оператора $R[{{\sigma }_{k}},{\mathbf{j}}_{{k,(m - 1)}}^{ + }]({{x}_{i}})$. При этом с вычислительной точки зрения выполнение такого шага итераций можно представить в виде умножения некоторой матрицы на вектор, причем эта матрица является плотной. Поэтому существенного повышения вычислительной эффективности метода можно добиться за счет использования быстрых матричных алгоритмов – это алгоритмы, в которых осуществляется аппроксимация матрицы, позволяющая приближенно выполнять умножение матрицы на вектор с вычислением небольшого числа элементов матрицы и имеющие вычислительную сложность меньше, чем $\mathcal{O}({{n}^{2}})$.

Более того, расчет матриц видимости ${{A}^{ \pm }} = (a_{{ik}}^{ \pm })$ по прямому алгоритму, при котором каждый элемент $a_{{ik}}^{ \pm }$ вычисляется путем проверки пересечения отрезка, соединяющего точки ${{x}_{i}}$ и ${{x}_{k}}$, c каждой из ячеек разбиения поверхности $\tilde {\Sigma }$, требует $\mathcal{O}({{n}^{3}})$ операций.

В настоящей работе для ускорения итерационного процесса, определяемого формулами (12), применяется метод мозаично-скелетных аппроксимаций. Существенная особенность данной задачи с точки зрения применения метода мозаично-скелетных аппроксимаций обусловлена дискретным характером функций видимости ${{a}^{ \pm }}(x,y)$, а так же необходимостью применения быстрых алгоритмов и к приближенному вычислению самих этих функции. Ниже сначала дадим общее описание метода мозаично-скелетных аппроксимаций матриц, а затем опишем особенности применения этого метода в рассматриваемой задаче.

3. МЕТОД МОЗАИЧНО-СКЕЛЕТНЫХ АППРОКСИМАЦИЙ

Алгоритм построения аппроксимации $\bar {A}$ для исходной матрицы $A$ на основе метода мозаично-скелетных аппроксимаций состоит из следующих этапов [7].

Шаг 1. Построение дерева кластеров.

Шаг 2. Построение списка блоков матрицы.

Шаг 3. Аппроксимации блоков матрицы.

На этапе построения дерева кластеров используется информация о сетке, на которой проведена дискретизация исходного оператора. Сетка рассматривается как набор точек коллокации ${{x}_{i}}$, $i = 1, \ldots ,n$, используемых в численной схеме (12) . Так как для токов на каждой ячейке рассматриваются постоянные значения ${\mathbf{j}}_{{j,(m)}}^{ \pm },$ $j = 1,...,n$, то каждому столбцу аппроксимируемой матрицы $A$ можно поставить в соответствие произвольную точку ${{y}_{j}}$ в ячейке ${{\sigma }_{j}}$. Мы в дальнейшем будем считать, что ${{y}_{j}} = {{x}_{j}}$ (это используется при учете видимости ячеек). Таким образом, квадратная матрица $A$ размера $n$ определяет некоторое взаимодействие “источников”, заданных на сетке $\{ {{y}_{j}}\} $, $j = 1, \ldots ,n$, с точками коллокации ${{x}_{i}}$, $i = 1, \ldots ,n$. Точки каждой из сеток (источников и точек коллокации) группируются по степени геометрической близости в кластеры. Кластеры в свою очередь формируют бинарное дерево, в корне которого лежит кластер, содержащий все узлы сетки, и при переходе на новый уровень дерева текущий кластер разбивается на два дочерних.

Имея два дерева кластеров, отвечающих паре сеток, на которых проведена дискретизация, мы разбиваем матрицу на иерархический набор блоков разного размера, каждый из которых отвечает взаимодействию некоторой группы источников $\{ {{y}_{j}}\} $ с некоторой группой точек коллокации $\{ {{x}_{i}}\} $. Блоки, отвечающие взаимодействию достаточно удаленных друг от друга кластеров точек, могут быть приближены матрицей малого ранга.

На фиг. 1 представлен пример разбиения матрицы на блоки в мозаично-скелетонном методе для частного случая, когда область разбиения является частью плоскости – квадратом. На этом рисунке представлено соответствие разбиения области, в которой находятся точки приемники (с источниками), дереву кластера и разбиению матрицы на блоки. Число уровней разбиения задается пользователем. Серым цветом показаны плотные блоки, остальные блоки малоранговые.

Используя алгоритм неполной крестовой аппроксимации [7], каждый малоранговый блок $B$ размера $p \times q$ с заданной точностью $\varepsilon $ представляется в виде скелетного разложения $B = U{{V}^{{\text{т}}}}$ (см. [6], [14]), где множители $U$ и $V$ имеют размеры $p \times r$ и $q \times r$ соответственно, а $r$ – ранг $B$.

Скелетное представление позволяет хранить матрицу в сжатом виде, так как вместо $p \cdot q$ элементов блока хранится только $(p + q)r$ элементов, что дает значительный выигрыш при $r \ll \min (p,q)$. Кроме того, малоранговый блок умножается на вектор за $(p + q)r$ операций вместо $p \cdot q$ операций для плотного блока.

4. БЫСТРЫЕ АЛГОРИТМЫ В ФИЗИЧЕСКОЙ ОПТИКЕ С ИСПОЛЬЗОВАНИЕМ ИЕРАРХИЧЕСКОГО ПРЕДСТАВЛЕНИЯ ДАННЫХ

Перепишем формулу (12) с учетом того, что токи ${\mathbf{j}}_{{i,(m)}}^{ \pm }$ постоянны на $i$-й ячейке:

где введены 4 матрицы, для которых необходимо применить мозаично-скелетный алгоритм,

В формуле (14) коэффициенты ${{R}_{{ \pm {\kern 1pt} \pm }}}[{{\sigma }_{k}}]({{x}_{i}})$ отличаются друг от друга только коэффициентами $a_{{ik}}^{ \pm }$ и $a_{{ki}}^{ \pm }$. Для удобства записи функцию видимости $a_{{ik}}^{ + }a_{{ki}}^{ + },\;a_{{ik}}^{ - }a_{{ki}}^{ + },\;a_{{ik}}^{ + }a_{{ki}}^{ - },\;a_{{ik}}^{ - }a_{{ki}}^{ - }$ соответствующего оператора ${{R}_{{ + {\kern 1pt} + }}},\;{{R}_{{ - {\kern 1pt} + }}},\;{{R}_{{ + {\kern 1pt} - }}},\;{{R}_{{ - {\kern 1pt} - }}}$ обозначим через ${{\hat {a}}_{{ik}}}$. Тогда формулы (14) запишем в виде одной:

где оператор $\hat {R}$ принимает одно из значений ${{R}_{{ + {\kern 1pt} + }}},\;{{R}_{{ - {\kern 1pt} + }}},\;{{R}_{{ + {\kern 1pt} - }}},\;{{R}_{{ - {\kern 1pt} - }}}$.

Заметим, что матрицы $(R[{{\sigma }_{j}}]({{x}_{i}})))$ являются векторными. Задача вычисления векторов ${\mathbf{H}}_{m}^{ \pm }({{x}_{i}})$, $i = 1,...,n$, по формулам (13) сводится к вычислению векторов

где ${{{\mathbf{j}}}_{k}}$ – один из наборов векторов ${\mathbf{j}}_{{k,(m - 1)}}^{ + }$ или ${\mathbf{j}}_{{k,(m - 1)}}^{ - }$. В свою очередь задача нахождения каждой из координат всех векторов $\widehat {{{{\mathbf{h}}}_{i}}}$ сводится к умножению матриц, составленных из координат векторов $\hat {R}[{{\sigma }_{k}}]({{x}_{i}})$, на соответствующие столбцы координат векторов ${{{\mathbf{j}}}_{k}}$ (нужно выполнить 6 умножений скалярных матриц на векторы). При этом нужно аппроксимировать 3 скалярные матрицы, составленные из координат векторов $\hat {R}[{{\sigma }_{k}}]({{x}_{i}})$.

Далее заметим, что можно аппроксимировать интегралы (15) так, чтобы задача свелась к выполнению умножений на одну скалярную квадратную матрицу размера $n$. Для этого запишем приближенное равенство:

(это выражение используем при $i \ne k$, при $i = k$ можем положить $\tilde {R}[{{\sigma }_{k}}]({{x}_{i}})) = 0$, так как $a_{{ik}}^{ \pm } = 0$ при $i = k$). Тогда где Задача нахождения векторов $\hat {R}[{{\sigma }_{k}}]({{x}_{i}})$ сводится к 6 умножениям на скалярную квадратную матрицу ${{\hat {R}}_{{ik}}}$ размера $n$, что видно из записи: Интегралы в выражениях (17) вычислялись по квадратурным формулам, описанным в статье [13].

Перейдем к применению метода мозаично-скелетных аппроксимаций. Для всех четырех матриц ${{\hat {R}}_{{ik}}}$ (14) строится одно и то же дерево кластеров. При построении блочной структуры матрицы $\hat {R} = ({{\hat {R}}_{{ij}}})$ необходимо учитывать дискретный характер функции видимости ${{\hat {a}}_{{ik}}}$. Если в малоранговый блок будут входить как видимые, так и невидимые пары ячеек, то в малоранговом блоке необходимо аппроксимировать нулевые подматрицы, что при заданной точности аппроксимации блока приводит к значительному росту ранга блока. Более того, если известно, что большая подматрица соответствует невидимым друг относительно друга областям, то умножение такой подматрицы на вектор в (13) будет давать нулевую добавку к вектору ${\mathbf{H}}_{m}^{ \pm }({{x}_{i}})$, поэтому не нужно аппроксимировать и хранить такую подматрицу.

При аппроксимации матрицы видимости полагаем, что точки-приемники и точки-источники совпадают (${{x}_{i}} = {{y}_{i}}$ для всех $i = 1,...,n$).

Согласно мозаично-скелетному методу при проходе по кластерным деревьям от корней к листьям каждой паре узлов дерева соответствуют некоторые области пространства ${{\Sigma }_{{\text{s}}}}$ и ${{\Sigma }_{{\text{r}}}}$, содержащие некоторую группу источников $\{ {{y}_{j}}\} $ и некоторую группу точек приемников (точек коллокации)$\{ {{x}_{i}}\} $. Этим группам точек соответствует некоторый блок матрицы. Введем следующие типы блоков с точки зрения выполнения условия видимости между их элементами.

1. Назовем блок “нулевым” , если для пары узлов деревьев, определяющих блок, соответствующая функция видимости ${{\hat {a}}_{{ij}}}$ равна нулю для любых $i$ и $j$.

2. Назовем блок “видимым”, если блок является малоранговым с точки зрения аппроксимации матрицы ${{R}_{{ik}}}$ и при этом для пары узлов деревьев, определяющих блок, соответствующая функция видимости ${{\hat {a}}_{{ij}}}$ равна $1$ для любых $i$ и $j$.

Сначала выполняется мозаичное биение матрицы на исходные малоранговые и плотные блоки с точки зрения аппроксимации матрицы ${{R}_{{ij}}}$ без учета видимости.

Далее все блоки проверяются на предмет отнесения к одному иp типов: “нулевой” или “видимый”. Если блок не принадлежит ни к одному из этих типов, то этот блок разбивается на подблоки следующего уровня и проверяется принадлежность этих подблоков к рассматриваемым типам. Если подблок не принадлежит ни к одному из указанных типов, то он опять разбивается на подблоки следующего уровня и т.д., до перехода на нижний используемый уровень. Блоки нижнего уровня, не относящиеся к рассматриваемым типам, объявляются плотными.

В результате возникло уточненное биение матрицы, которое содержит блоки: “нулевые”, “видимые” малоранговые и “плотные”.

“Нулевые” блоки не вычисляются, не хранятся и не используются при расчете векторов ${\mathbf{H}}_{m}^{ \pm }({{x}_{i}})$ из формулы (13) и токов ${\mathbf{j}}_{{i,(m)}}^{ \pm }$ из формулы (12). “Видимые” малоранговые блоки (далее будем называть их малоранговыми) аппроксимируются скелетонным представлением. Для плотных блоков вычисляются все элементы матрицы.

В расчетах проверка типа блока матрицы (“нулевой” или “видимый”) осуществляется приближенно. Вычислительные эксперименты показывают, что для определения свойства видимости на поверхности $\Sigma $ можно воспользоваться следующей процедурой. Для блока, заданного набором точек коллокации $\{ {{x}_{i}}\} ,$ $i = 1, \ldots ,p$, и группой источников $\{ {{y}_{j}}\} ,$ $j = 1, \ldots ,q$, рассмотрим подмножество точек коллокации $\{ x_{i}^{0}\} ,$ $i = 1, \ldots ,{{p}_{0}}$ (набор индексов обозначим через $\{ {{I}_{0}}\} $) и источников $\{ y_{j}^{0}\} ,$ $j = 1, \ldots ,{{q}_{0}}$ (набор индексов обозначим через $\{ {{J}_{0}}\} $), причем ${{p}_{0}} \ll p$ и ${{q}_{0}} \ll q$. Вычисляем функции видимости ${{\hat {a}}_{{ij}}}$ для $i = 1, \ldots p,$ $j \in {{J}_{0}},$ и для $i \in {{I}_{0}},$ $j = 1, \ldots q$. Если для всех выбранных индексов функция принимает одно и то же значение (видимости или невидимости), то предполагаем, что весь блок принимает то же самое значение. Наборы точек $\{ x_{i}^{0}\} ,$ $i = 1, \ldots ,{{p}_{0}}$, и $\{ y_{j}^{0}\} ,$ $j = 1, \ldots ,{{q}_{0}}$, для блока формировались следующим образом. Каждый блок определяет взаимодействие области ${{\Sigma }_{i}}$, в которой находятся точки коллокации $\{ {{x}_{i}}\} $, с областью ${{\Sigma }_{j}}$, в которой находятся источники $\{ {{y}_{j}}\} $. Область ${{\Sigma }_{i}}$ разбивалась на ${{p}_{0}}$ частей, в каждой из которой выбиралась точка $\{ x_{i}^{0}\} ,$ ближайшая к центральной точке подобласти. Аналогично выбиралась точка $\{ y_{j}^{0}\} $ из разбиения области ${{\Sigma }_{j}}$.

Такая проверка осуществляется только для блоков, начиная с некоторого назначаемого уровня ${{k}_{0}}$. Если блок имеет низший уровень, то такая приближенная проверка делается для всех его подблоков уровня ${{k}_{0}}$. Это актуально для тел со сложной геометрией, например, при наличии небольших отверстий и щелей.

В приводимых ниже примерах бралось значение параметра ${{k}_{0}} = 1$. При этом численные исследования показали, что для всех блоков достаточно выбрать параметры ${{p}_{0}} = {{q}_{0}} = 8.$

Теперь рассмотрим алгоритм построения функции видимости для двух точек.

При проверке видимости пары точек сразу проверятся условия:

• если $|{\kern 1pt} ({{y}_{k}} - {{x}_{i}},{\mathbf{n}}({{x}_{i}})){\kern 1pt} {\text{|}} \leqslant {{10}^{{ - 4}}}$ или $|{\kern 1pt} ({{y}_{k}} - {{x}_{i}},{\mathbf{n}}({{y}_{k}})){\kern 1pt} {\text{|}} \leqslant {{10}^{{ - 4}}}$, то $a_{{ik}}^{ \pm } = 0$;

• пусть $\Sigma $ ограничивает замкнутую (не обязательно выпуклую) область и эта поверхность аппроксимируется замкнутой же поверхностью $\tilde {\Sigma }$ и пусть точка ${{x}_{i}}$ лежит на поверхности $\tilde {\Sigma }$; так как на замкнутой поверхности поверхностные токи размещаются только на внешней стороне, следует положить $a_{{ik}}^{ - } = 0$, $k = 1,...,n$. Также для любой точки ${{y}_{k}}$, расположенной так, что $({{y}_{k}} - {{x}_{i}},{{{\mathbf{n}}}_{i}}) \leqslant {{10}^{{ - 4}}}$, сразу выполнено условие $a_{{ik}}^{ + } = 0$ (точка ${{x}_{i}}$ видна из точки ${{y}_{k}}$ изнутри). Для такой пары точек проверка видимости не проводится.

• если при проверке пересечения отрезка $[{{x}_{i}},{{y}_{k}}]$ с ячейкой ${{\sigma }_{j}}$ выяснилось, что точки ${{x}_{i}}$ и ${{y}_{k}}$ не видны, то нет необходимости проверять остальные ячейки поверхности $\tilde {\Sigma }$.

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

Для оставшихся пар точек воспользуемся иерархическими процедурами определения видимости.

Опять осуществим иерархическое разбиение области, содержащей поверхность $\tilde {\Sigma }$, на подобласти. Область, в которой находится поверхность $\tilde {\Sigma }$, помещается в куб. Далее строится восьмеричное дерево, согласно которому кубические области разбиваются на кубические подобласти с уменьшением размеров куба в 2 раза. Если подобласть не содержит ячейки поверхности $\tilde {\Sigma }$, то такая подобласть отмечается и дальше не разбивается. Пример иерархического разбиения области изображен на фиг. 2 (поверхность изображена жирной черной линией). Подготовительный этап выполняется один раз перед расчетом матриц видимости $a_{{ik}}^{ \pm }$.

Для проверки видимости точки ${{x}_{i}} \in \tilde {\Sigma }$ относительно точки ${{y}_{k}} \in \tilde {\Sigma }$ определяются кубические области $\{ {{\Omega }_{l}}\} $, которые соответствуют листьям восьмеричного дерева иерархичного разбиения области и которые могут содержать ячейки, пересекаемые отрезком $[{{x}_{i}},{{y}_{k}}]$. В список таких областей входят:

• области, содержащие точки ${{x}_{i}}$ и ${{y}_{k}}$;

• области, пересекаемые отрезком $[{{x}_{i}},{{y}_{k}}]$.

Набор областей получается при проходе дерева от корня к листьям. Если в результате прохода встретилась область, которую не пересекает отрезок $[{{x}_{i}},{{y}_{k}}]$, то далее подобласти такой области не проверяются. Если область не содержит ни одной ячейки поверхности $\tilde {\Sigma }$, то такая область в набор областей $\{ {{\Omega }_{l}}\} $ не включается.

После формирования списка областей $\{ {{\Omega }_{l}}\} $ проверяется пересечение отрезка $[{{x}_{i}},{{y}_{k}}]$ только с ячейками, лежащими в областях ${{\Omega }_{l}}$.

Решение задачи по формуле (13) с использованием малоранговых аппроксимаций показывает, что на аппроксимацию матрицы тратится гораздо больше времени, чем на умножение малоранговой матрицы на вектор. Поэтому перед итерациями один раз строится блочная малоранговая структура матриц ${{\hat {R}}_{{ij}}}$, а также малоранговая аппроксимация скалярной матрицы $({{R}_{{ij}}})$. Матрица в сжатом виде хранится в оперативной памяти, а на каждой итерации пересчет токов выполняется по формулам (12), (13).

Осуществим приближенную оценку сложности построенного вычислительного алгоритма. Сначала оценим сложность алгоритма выделения “нулевых” и “видимых” малоранговых блоков.

Число уровней восьмеричного дерева выбирается таким образом, чтобы с ростом числа ячеек $n$ среднее число ячеек на нижнем уровне оставалось примерно постоянным. Тогда число блоков в мозаичном биении матрицы растет как $\mathcal{O}(n\log n)$ (см. [15]). Значит, для выделения “нулевых” и “видимых” блоков требуется $\mathcal{O}(n\log n)$ вызовов функции видимости для пары расчетных точек ${{x}_{i}}$ и ${{y}_{j}}$. Сложность алгоритма выделения таких блоков можно оценить как $\mathcal{O}(n\log n)\Psi (n)$, где $\Psi (n)$ – среднее число операций при проверке видимости для пары проверяемых расчетных точек ${{x}_{i}}$ и ${{y}_{j}}$. Сразу можно сказать, что заведомо выполнено условие $\Psi (n) = \mathcal{O}(n)$.

Более сильную оценку для числа операций при проверке видимости для пары точек, лежащих на поверхности, можно получить, если предположить, что отрезок, соединяющий две любые такие точки, пересекает не более чем $L$ непустых областей одного уровня восьмеричного дерева (областей из списка $\{ {{\Omega }_{l}}\} $). Пусть $x$ и $y$ – две точки, лежащие на поверхности $\tilde {\Sigma }$, возникающей при аппроксимации поверхности $\Sigma $ системой из ${{n}_{1}}$ ячеек. Предположим, что для проверки функции видимости для этих точек потребовалось ${{\Psi }_{{x,y}}}({{n}_{1}})$ проверок пересечения отрезка $[x,y]$ с областями из списка $\{ \Omega _{l}^{1}\} $. Измельчим сетку на поверхности $\tilde {\Sigma }$, разбив каждую ячейку на 4 части. Тогда на поверхности $\tilde {\Sigma }$ возникнет новая поверхностная сетка с числом ячеек ${{n}_{2}} = 4{{n}_{1}}$. Увеличим число уровней восьмеричного дерева на 1, разбив каждую область нижнего уровня на 8 кубических областей следующего уровня, и опять выделим семейство областей $\{ \Sigma _{l}^{{(2)}}\} $, содержащих точки поверхности. При этом в каждой из областей семейства $\{ \Sigma _{l}^{{(2)}}\} $ опять окажется примерно $m$ ячеек. В этом случае для каждой области $\Omega $ из списка $\{ \Omega _{l}^{1}\} $, являющейся областью нижнего уровня для первого случая и имеющей пересечение с отрезком $[x,y]$, придется провести не более 8 проверок на пересечение этого же отрезка с областями нижнего уровня для второго дерева, лежащими в области $\Omega $. Поскольку число таких областей $\Omega $ не более $L$, получаем: ${{\Psi }_{{x,y}}}(4{{n}_{1}}) \leqslant {{\Psi }_{{x,y}}}({{n}_{1}}) + 8L$, откуда следует, что ${{\Psi }_{{x,y}}}(n) = \mathcal{O}(\log n)$.

Если гипотеза об ограниченности числа точек пересечения отрезков с кубами нижнего уровня верна, то сложность алгоритма выделения “нулевых” и “видимых” блоков оценивается в $\mathcal{O}(n\mathop {\log }\nolimits^2 n)$ операций.

После построения блочной структуры для матриц (14) с помощью алгоритма неполной крестовой аппроксимации необходимо построить скелетные разложения малоранговых блоков. При этом число операций при разложении малоранговых блоков не превосходит число операций при аппроксимации матрицы ${{R}_{{ij}}}$ (без множителей видимости), так как система малоранговых блоков матрицы ${{\hat {R}}_{{ij}}}$ по сути отличается от системы малоранговых блоков матрицы ${{R}_{{ij}}}$ выбрасыванием “нулевых” блоков и отнесением части малоранговых блоков к плотным. Поэтому для алгоритма построения скелетного разложения малоранговых блоков матрицы ${{\hat {R}}_{{ij}}}$ применима оценка $\mathcal{O}(rn{{\log }^{2}}n)$, где $r$ – мозаичный ранг [15]. Для числа операций при умножении малоранговой части (без учета плотных блоков) матрицы ${{\hat {R}}_{{ij}}}$ на вектор в сжатом формате применима та же оценка, что и для матрицы ${{R}_{{ij}}}$, и это число операций также оценивается как $\mathcal{O}(rn\log n)$ [15].

Таким образом, в случае если плотные блоки имеют малые размеры, можно надеяться на выполнение оценки $\mathcal{O}(rn{{\log }^{2}}n)$ как для числа операций по малоранговому разложению матриц ${{\hat {R}}_{{ij}}}$, так и для числа операций при выполнении каждой итерации по формулам (12).

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

5. ПРИМЕРЫ РАСЧЕТОВ И ОБСУЖДЕНИЕ РЕЗУЛЬТАТОВ

В результате решения дискретных уравнений (12) находятся приближенные значения поверхностных токов на ячейках разбиения. Далее по найденным поверхностным токам можно рассчитать распределения электрического и магнитного полей в окружающем пространстве, используя формулы (6).

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

Эффективной площадью рассеяния в направлении единичного вектора $\tau $ называется величина:

где ${{{\mathbf{E}}}_{0}}$ – напряженность падающего электрического поля, определяемого формулой (2). Величина ${{\sigma }_{0}}(\tau )$ приближенно выражается через поверхностные токи по формуле: где ${\text{|}}{{\sigma }_{l}}{\kern 1pt} {\text{|}}$ – площадь ячейки ${{\sigma }_{l}}$.

В приводимых ниже примерах положение облучаемых тел рассматривалось в декартовой системе координат $Oxyz$, направление облучения определялось волновым вектором ${\mathbf{k}}$, причем вектор ${\mathbf{k}}$ брался параллельным плоскости $Oxy$ в виде ${\mathbf{k}} = k( - \cos \alpha , - \sin \alpha ),0)$ ($\alpha $ – угол между вектором $ - {\mathbf{k}}$ и осью $Ox$), $k$ – волновое число (см. схему на фиг. 3а). Были получены диаграммы обратного рассеяния, показывающие зависимость эффективной площади рассеяния ${{\sigma }_{0}}(\tau )$ в направлении, противоположном направлению облучения ($\tau = - {\mathbf{k}}{\text{/}}k$) от угла $\alpha $ в децибеллах:

В первом примере рассматривалось рассеяние плоской волны на уголковом отражателе. Схема уголкового отражателя показана на фиг. 3а. На фиг. 3б приведены диаграммы обратного рассеяния для частоты $f = 10$ ГГц при $Oxy$ поляризации (в формуле (2) вектор ${{{\mathbf{E}}}_{i}}$ компланарен плоскости $Oxy$). В проведенном расчете использовалось разбиение на квадратные ячейки, длина волны составляла $\lambda = 3$ см, число ячеек – 12  800, что соответствовало стороне ячейки 0.25 см (длина волны определяется соотношением $\lambda = c{\text{/}}f$, $c = \sqrt {{{\varepsilon }_{0}}{{\mu }_{0}}} $ – скорость света в вакууме). На графиках пунктирной линией изображена диаграмма, полученная авторами методом граничных интегральных уравнений (использовалась численная схема, описанная в статьях [9], [16], [17]), а сплошной линией – график, полученный методом физической оптики, разработанным в данной статье. Наблюдаемое хорошее согласование приведенных результатов свидетельствует о применимости разработанной итерационной математической модели физической оптики с переотражениями при решении задач рассеяния на невыпуклых телах.

Следующий пример для тела более сложной геометрии иллюстрирует работу алгоритма аппроксимации матриц. Рассматривался двойной уголок, схема которого приведена на фиг. 4 (рисунок из статьи [18]). На фиг. 4, 5 показаны диаграммы обратного рассеяния для частоты 16.66 ГГц (длина волны – 1.8 см). На фиг. 4б показаны данные физического эксперимента, которые были приведены в статье [18]. На фиг. 5 приведены диаграммы обратного рассеяния, полученные авторами методом граничных интегральных уравнений и разработанным методом физической оптики.

В расчетах методом интегральных уравнений и методом физической оптики использовалась поверхностная сетка с числом ячеек 97  240, форма ячеек прямоугольная (близкая к квадратной), максимальная сторона ячейки 0.227 см. Видно хорошее согласование результатов, полученных указанными численными методами и в физическом эксперименте.

Вычислительные затраты при решении задачи дифракции на двойном уголке:

n = 1600,  6400,  25600,  102  400;

h = 0.866,  0.433,  0.217,  0.108;

T₁ = 75.5,  512.8,  2828.9,  23  459;

M₁ = 129,  650,  3673,  26  730;

T₂ = 7.3,  19.7,  75.9,  453.1;

M₂ = 89,  269,  983,  3322.

Здесь представлено:

$n$ – число ячеек;

$h$ – шаг сетки (в сантиметрах);

${{T}_{1}}$ – время (в секундах) решения задачи методом интегральных уравнений с применением малоранговых аппроксимаций;

${{M}_{1}}$ – память (в мегабайтах), необходимая для решения задачи методом интегральных уравнений с применением малоранговых аппроксимаций;

${{T}_{2}}$ – время (в секундах) решения задачи методом физической оптики с применением малоранговых аппроксимаций;

${{M}_{2}}$ – память (в мегабайтах), необходимая для решения задачи методом физической оптики с применением малоранговых аппроксимаций.

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

Для оценки вычислительной эффективности разработанного численного алгоритма оценки зависимости числа операций (по порядку) от числа ячеек разбиения была проведена серия расчетов для того же двойного уголока при частоте облучения 4 ГГц (длина волны $h = 7.5$ см). Были использованы 4 расчетные сетки с числом ячеек от 1600 (максимальная сторона ячейки $h = 0.886$ см) до 102  400 (максимальная сторона ячейки $0.108\;{\text{см}}$). При этом на самой грубой сетке выполнялось соотношение $h < 8\lambda $, на всех сетках были получены близкие диаграммы обратного рассеяния.

Выще приведены данные, характеризующие вычислительные затраты (время вычислений и затраты оперативной памяти) при применении разработанного алгоритма физической оптики, а также при применении метода граничных интегральных уравнений (использовался алгоритм авторов с применением мозаично-скелетонных аппроксимаций матриц [9], [16], [17]). Рассматривалась задача определения диаграммы обратного диаграммы рассеяния для 180 направлений облучения. Для каждого направления облучения решаются уравнения (13) с одними и теми же аппроксимированными заранее матрицами $\hat {R}[{{\sigma }_{k}}]({{x}_{i}})$ вида (16) и со своей функцией ${{{\mathbf{H}}}_{0}}(x)$. В случае применения метода интегральных уравнений многократно решается система линейных уравнений с одной и той же матрицей и различными правыми частями. Расчеты проводились на персональном компьютере с процессором Intel Core-I5 с оперативной памятью 32 Гб.

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

Для оценки зависимости числа операций (по порядку) от числа ячеек разбиения на фиг. 6 приведен график зависимости $\lg {{T}_{2}}$ от $\lg n$ на основе вышеприведенных данных (сплошная линия). На фиг. 6 штриховой линией показана линейная зависимость, построенная по крайним точкам – эта зависимость имеет вид $\lg {{T}_{2}} = A\lg n + B$, где $A \approx 0.992$. Тем самым наблюдается близкая к линейной зависимость времени вычислений от числа ячеек разбиения, что подтверждает гипотезу о теоретической оценке сложности алгоритма вида $\mathcal{O}(rn{{\log }^{2}}n)$.