Распространение радиоволн ВЧ/Рей-трейсинг

Модель

Графическая иллюстрация модели рейтрейсинга

Алгоритм

Идея последовательного приближения

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

$\left | F_i-F_{i-1} \right |<\Delta$, где

$\Delta$ - параметр моделирования, задаваемый пользователем.

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

$\left | F_i-F_{i-k} \right |<\Delta$, где

$k$ также будет задаваться пользователем.

$i$ здесь - это параметр цикла, стоящего над циклами основной программы, $i=\overline{0:N}$.

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

Шаги $\Delta_{\theta}\left(\theta\right)$ по азимуту и $\Delta_{\varphi}\left(\varphi\right)$ по зениту источника являются функциями от направления либо постоянными.

Угловой шаг дискретизации как функция ХН

Диаграмма направленности и испускаемые моделью лучи

Зависимость углового шага от азимутального угла

Большинство современных вещательных систем используют панельные антенны с ограниченными углами раствора диаграммы направленности (ДН) в горизонтальной плоскости (до 120°) и очень малыми углами в вертикальной плоскости (до 20°). Соответственно, в таких системах происходит серьезное перераспределение излучаемой энергии в пространстве. Типичные коэффициенты усиления: 16-18 dBi. Поэтому одним из решений задачи оптимизации является использование динамического углового шага дискретизации $\Delta_{\theta},~\Delta_{\varphi}$ как функции от характеристики направленности источника.

Изменение частоты дискретизации $\omega(\theta,\varphi, f)$ происходит по следующему закону:

$\omega(\theta,\varphi, f) = \omega_{min} + k F(\theta,\varphi, f)$ , где

$k = \omega_{max} - \omega_{min}$ - коэффициент девиации, где

$\omega_{max} = \frac{1}{\arcsin\left(\frac{\lambda}{R_{max}}\right)}\approx\frac{R_{max}}{\lambda}$ - максимальное и $\omega_{min} = \alpha\omega_{max}$ (для любого $0 \leq \alpha \leq 1$) минимальное значение частоты дискретизации;

$R_{max}$ - максимальное расстояние от источника до границ модели;

$F$ - функция ХН;

$\theta,~ \varphi$ - угол места и азимутальный угол;

$f$ - частота излучения;

$\lambda$ - длина излучаемой волны.

Тогда шаг дискретизации (угол) будет меняться по следующим образом:

$\Delta_{\theta},~\Delta_{\varphi} = \frac {1} {\omega(\theta,\varphi)}$

Общее кол-во лучей определяется выражением:

$N_{total} = \frac {1} {2 \pi^2} \int \limits_0^{\pi} \int \limits_0^{2 \pi} \omega(\theta,\varphi)~d \theta d \varphi$

С учетом, что максимум ДН находится в $\theta=0,~\varphi=0$, получим выражения:

$N_{\theta~total} = \frac {1} {\pi} \int \limits_0^{\pi} \omega(\theta,\varphi) \bigr|_{\varphi=0}~d \theta$,

$N_{\varphi~total}(\theta) = \frac {1} {2 \pi} \int \limits_0^{2 \pi} \omega(\theta,\varphi)~d \varphi$.

Инициализация геометрической модели

Перерасчет высот с учетом кривизны земли и рефракции радиоволн в тропосфере.

Входной параметр $h$ пересчитывается в соответствии с формулой:

$h'(x,y) = h(x,y)-\frac{-2 R_{eq}+\sqrt{(2 R_{eq})^{2}+4 r^{2}}}{2} \approx h(x,y)-\frac{r^{2}}{2 R_{eq}}$ , где

$R_{eq}=\frac{R_0}{1+R_0 \frac{dn}{dh}}$ - эквивалентный радиус Земли, где

$R_0=6371$ - радиус Земли (км),

$\frac{dn}{dh}=grad~n$ - изменение коэффициента преломления с высотой.

$r=\sqrt{x^{2}+y^{2}}$ - расстояние до точки с высотой $h(x,y)$.

Основной цикл программы | Вариант 1

Предусловия

• $G$ - входное описание среды распространения моделируемого поля, заданное геометрической моделью.
• $i$ - номер итерации моделирования, $i\geq 0$.

Основное течение

1. $S \leftarrow$ G.Множество первичных источников()
2. Если $S \colon ~\varnothing$ то ВЫХОД
3. $\forall ~ s_j \in S$
   1. $S \leftarrow S~ \backslash \left \{ s_j \right \}$
   2. $\forall ~ \omega_n \in$ $s_j$.Антенна().Тип антенны().Амплитудно-частотная характеристика()
      1. $\theta \leftarrow 0$
      2. $\zeta_{\theta} \leftarrow 0$
      3. Если $\theta < \pi \wedge (i=0 \vee \zeta_{\theta} \bmod 2 \neq 0)$
         1. $\varphi \leftarrow 0$
         2. $\zeta_{\varphi} \leftarrow 0$
         3. Если $\varphi<2\pi \wedge (i=0 \vee \zeta_{\varphi} \bmod 2 \neq 0)$
            1. Position $\leftarrow$ $s_j$.Антенна().Позиция()
            2. Ray $\leftarrow$ Луч().Создать(Position, $s_j$.Антенна().Мировая система координат(Вектор$(\theta,~\varphi, 1)$))
            3. Distance $\leftarrow \infty$
            4. $\forall ~ f_m\in$ G.Множество отражающих объектов()
               1. $\forall ~ t_{mn} \in f_m$.Множество отражающих поверхностей()
                  1. $P' \leftarrow$ Ray.Пересечение($t_{mn}$.Плоскость грани())
                  2. Если $t_{mn}$.Принадлежность($P'$)
                     1. Distance' $\leftarrow$ Расстояние(Position, $P'$)
                     2. Если Distance' $<$ Distance
                        1. Distance $\leftarrow$ Distance'
                        2. $t' \leftarrow t_{mn}$
                        3. $P \leftarrow P'$
            5. $\forall ~ \rho_k \in$ G.Множество контрольных точек()
               1. Если $\rho_k \in$ $s_j$.Область регистрации луча $(\theta,~\varphi,~\frac{\Delta_{\theta} (\theta)}{2^i},~\frac{\Delta_{\varphi} (\theta, \varphi)}{2^i},$ Distance$)$
                  1. $\rho_k$.Зарегистрировать($s_j$.Напряженность$(\omega_n,~\theta,~\varphi$, Расстояние(Position, $\rho_k$.Позиция()), G.Среда распространения()$)$)
            6. Если $s_j$.Напряженность$(\omega_n,~\theta,~\varphi$, Distance, G.Среда распространения()$)>E_{end}$
               1. Angle $\leftarrow$ Ray.Угол пересечения($t'$.Плоскость грани())
               2. $s' \leftarrow$ Вторичный источник при рейтрейсинге.Создать($s_j$.Напряженность$(\omega_n,~\theta,~\varphi$, Distance, G.Среда распространения()$),$ Angle, $P,~t'$)
               3. $S \leftarrow S \cup \{s'\}$
         4. $\varphi \leftarrow \varphi + \frac{\Delta_{\varphi} (\theta, \varphi)}{2^i}$
         5. $\zeta_{\varphi} \leftarrow \zeta_{\varphi}+1$
      4. $\theta \leftarrow \theta + \frac{\Delta_{\theta} (\theta)}{2^i}$
      5. $\zeta_{\theta} \leftarrow \zeta_{\theta}+1$
4. Переход на шаг 2

Основной цикл программы | Вариант 2

Геометрическая модель (рей-трейсинг 2). Наследует свойства Геометрической модели.

• Предикат завершения(). При $i=0$ возвращает истину, иначе возвращает значение функции SubApprox для результатов моделирования не предыдущем и на текущем этапах, а также для включаемого отклонения $\Delta$ и допустимого процента точек $k$.

Предусловия

• $G$ - входное описание среды распространения моделируемого поля, заданное геометрической моделью.
• $\Omega$ - список частот, на которых необходимо провести моделирование.
• $E_{end}$ - критическое значение напряженности, на котором дальнейшая трассировка луча прекращается.
• $\Delta$ - критерий разности между текущим и предыдущим результатом моделирования, при которой можно считать, что поле промоделированно точно.

Основное течение

1. $\forall ~ \omega_n \in \Omega$
   1. $i \leftarrow 0$
   2. Pre_SCP $\leftarrow$ G.Множество контрольных точек($\omega_n$) //Копирование состояния до изменения.
   3. Set_event $\leftarrow$ Preprocessing($G,~i,~\omega _{n}$)
   4. $\forall ~(s_{j},~\theta ,~\varphi )\in$ Set_event
      1. Ray $\leftarrow$ Луч().Создать($s_j$, $\theta ,~\varphi ,~\omega _{n}$)
      2. Trace(Ray, G, $E_{end}$)
   5. Postprocessing(G.Множество контрольных точек($\omega_n$))
   6. SAcheck $\leftarrow$ SucApprox(G.Множество контрольных точек($\omega_n$), Pre_SCP, $\Delta$)
   7. Если SAcheck
      1. $i \leftarrow i+1$
      2. Переход на 1.2
2. ВЫХОД

[Таблица событий] $\leftarrow$ Preprocessing(Геометрическая модель, Номер итерации, Частота)

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

1. $T\leftarrow$ Пустой односвязный список значений $(s_{j},~\theta ,~\varphi )$.
2. $G\leftarrow$ Геометрическая модель.
3. $i\leftarrow$ Номер итерации
4. $\omega \leftarrow$ Частота
5. $\forall ~ s_j \in$ G.Множество первичных источников()
   1. Если $s_j$.Антенна().Тип антенны().Амплитудно-частотная характеристика($\omega$) $\neq 0$
      1. $\theta \leftarrow 0$
      2. $\zeta_{\theta} \leftarrow 0$
      3. Пока $\theta < \pi$
         1. $\varphi \leftarrow 0$
         2. $\zeta_{\varphi} \leftarrow 0$
         3. Пока $\varphi<2\pi$
            1. Если $i=0\vee \zeta _{\theta }{\bmod {2}}\neq 0\vee \zeta _{\varphi }{\bmod {2}}\neq 0$
               1. $T$.Добавить$(s_{j},~\theta ,~\varphi )$
            2. $\varphi \leftarrow \varphi +{\frac {\Delta _{\varphi }(s_{j},~\theta ,~\varphi )}{2^{i}}}$
            3. $\zeta_{\varphi} \leftarrow \zeta_{\varphi}+1$
         4. $\theta \leftarrow \theta +{\frac {\Delta _{\theta }(s_{j},~\theta )}{2^{i}}}$
         5. $\zeta_{\theta} \leftarrow \zeta_{\theta}+1$
6. Вернуть Массив($T$)

Если $\Delta _{\varphi }(s_{j},~\theta ,~\varphi )={\rm {const}}$, и $\Delta _{\theta }(s_{j},~\theta )={\rm {const}}$, то есть углы распределены равномерно, таблица может быть реализована в виде массива значений $(s_{j},~\theta ,~\varphi )$, элементы которого могут заполняться параллельно, без синхронизации. Если же размер таблицы заранее не определен (общий случай), то таблица реализуется связанным списком, добавление элементов в который должно защищаться мьютексом на шаге 4.1.3.1.3.1.1, а алгоритм на заключительном шаге преобразует связанный список в массив с произвольным доступом. Второй путь при параллельном выполнении - создавать таблицы в локальной памяти потока, после выполнения которых - объединять таблицы в массив в основном потоке.

Параллельная форма

Параллельное выполнение в общем случае возможно только цикла, перебирающего первичные источники (шаг 4), так как распределение углов заранее не определено.

Пусть $P>0$ - число потоков, выполняющих функцию. Пусть $T'$ - массив изначально пустых связанных списков троек $(s_{j},~\theta ,~\varphi )$, так, что $T'_{i}$ - $i$-тый список массива, $0\leq i<P$.

1. $i\leftarrow$ Номер итерации
2. $\omega \leftarrow$ Частота
3. $T'\leftarrow$Создать массив.
4. $\forall ~ s_j \in$ G.Множество первичных источников()
   1. $t\leftarrow$ индекс потока.
   2. Если $s_j$.Антенна().Тип антенны().Амплитудно-частотная характеристика($\omega$) $\neq 0$
      1. $\theta \leftarrow 0$
      2. $\zeta_{\theta} \leftarrow 0$
      3. Пока $\theta < \pi$
         1. $\varphi \leftarrow 0$
         2. $\zeta_{\varphi} \leftarrow 0$
         3. Пока $\varphi<2\pi$
            1. Если $i=0\vee \zeta _{\theta }{\bmod {2}}\neq 0\vee \zeta _{\varphi }{\bmod {2}}\neq 0$
               1. $T'_{t}$.Добавить$(s_{j},~\theta ,~\varphi )$
            2. $\varphi \leftarrow \varphi +{\frac {\Delta _{\varphi }(s_{j},~\theta ,~\varphi )}{2^{i}}}$
            3. $\zeta_{\varphi} \leftarrow \zeta_{\varphi}+1$
         4. $\theta \leftarrow \theta +{\frac {\Delta _{\theta }(s_{j},~\theta )}{2^{i}}}$
         5. $\zeta_{\theta} \leftarrow \zeta_{\theta}+1$
5. Вернуть Массив$\left(\bigsqcup _{t=0}^{P-1}T'_{t}\right)$

Trace

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

1. $R\leftarrow \infty$
2. $T\leftarrow$ NaN
3. $\forall ~ f_m\in$ G.Множество отражающих объектов()
   1. $\forall ~ t_{mn} \in f_m$.Множество отражающих поверхностей()
      1. $p \leftarrow$ Ray.Пересечение($t_{mn}$.Плоскость())
      2. Если $t_{mn}$.Принадлежность($p$)
         1. $r\leftarrow$ Ray.Пробег($p$)
         2. Если $r < R$
            1. $R \leftarrow r$
            2. $T \leftarrow t_{mn}$
            3. $P \leftarrow p$
4. $\forall ~ O_k \in$ G.Множество плоскостей вывода результатов()
   1. $p \leftarrow$ Ray.Пересечение($O_k$.Плоскость())
   2. Если $O_k$.Принадлежность($p$)
      1. $r\leftarrow$ Ray.Пробег($p$)
      2. Если $r < R$
         1. $\rho \leftarrow$ $O_k$.Ближайшая контрольная точка($p$)
         2. $\rho$.Зарегистрировать(Ray, $r$, G.Среда распространения())
5. Если $T \not\equiv$ NaN
   1. $\gamma \leftarrow$ Ray.Угол пересечения($T$.Плоскость())
   2. Ray.Отразить($T,~P,~R,~\gamma$, G.Среда распространения())
   3. Если Ray.Напряженность $>E_{end}$
      1. Trace(Ray, G, $E_{end}$)

Postprocessing

Функция "проявляет" поле в контрольных точках по накопленным значениям напряженности от всех пришедших лучей.

SucApprox

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

Свойства алгоритма

Возможности распараллеливания