Algorithm for filtering coordinates of a ground source unauthorized radio emission in a satellite communication system with direct retransmission

Full Text

Введение

Вопросы теории и практики оценки координат объектов находят широкое применение в различных отраслях науки, техники и промышленности (в радионавигации, радиолокации, радиоуправлении, радиоастрономии, радиогеодезии, в военной сфере). Для решения подобного класса задач учеными (К.М. Антонович, Дж. Вэббер, Д.А. Кошаев, С.П. Панько, В.В. Савельев, А.А. Савин, А.Г. Сайбель, Ю.Г. Сосулин, В.С. Черняк, Дж. Эффланд, М.С. Ярлыков и др.) отечественной и зарубежных научных школ статистической теории радиолокации и радионавигации разработан ряд адекватных математических моделей фильтрации координат подвижного объекта.

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

Известные алгоритмы и программно-аппаратные системы координатометрии наземных ИНР в ССС с прямой ретрансляцией, которые разработаны такими учеными, как Р.В. Волков, А.Х. Кельян, В.И. Могучев, В.Н. Саяпин, В.В. Сухотин, В.И. Тисленко, А.О. Чемарев, Wolfgang Koch, Darko Musicki, Hugo Seuté, Ali Khenchaf и др. [1–7], предполагают двухэтапное разделение. На первом решается задача оценки координат и вектора скорости возмущенного орбитального движения основного и смежного спутника-ретранслятора (СР) на различных интервалах времени наблюдения [8]. На втором этапе для определенных оценок кинематических параметров СР выполняется непосредственное решение задачи координатометрии ИНР. Подобное разделение несколько снижает точность местоопределения ИНР с точки зрения классической теории оценивания и ее применения в связи и управлении (Э. Сейдж, Дж. Мелс, Дж. Вебер, Р.Л. Стратонович, В.И. Тихонов и др.). Устранение указанной особенности составляет суть настоящей статьи.

1. Содержательная постановка задачи фильтрации

Анализ результатов [1–4; 8–10] определяет особенности для формирования обобщенной схемы фильтрации координат наземного ИНР в ССС с прямой ретрансляцией (рис. 1).

 

Рис. 1. Схема фильтрации координат наземного ИНР

Fig. 1. Filtration scheme for coordinates of ground-based INR

 

Основные объекты схемы: ИНР; комплекс радиоконтроля (КРК); $K\ge 2$ СР; $N\ge 1$ реперных станций (РС). Вектор параметров состояния k-го СР ${\overrightarrow{\lambda }}^k\left(t\right)=\left\{{\overrightarrow{r}}^k\left(t\right),{\overrightarrow{V}}^k\left(t\right)\right\}$ $\left(k=\overline{\mathrm{1,}K}\right)$ является векторным случайным процессом, который задает положение ${\overrightarrow{r}}^k\left(t\right)$ и скорость ${\overrightarrow{V}}^k\left(t\right)$ в момент времени t k-го СР. Оценка средневыборочных значений ${\overrightarrow{r}}^k\left(t\right)$ и ${\overrightarrow{V}}^k\left(t\right)$ в момент времени t осуществляется по данным телеметрии при численном решении дифференциального уравнения, описывающего возмущенное движение СР [8], при усреднении получаемых значений в момент времени t по данным измерений РС. Положения КРК и n-x $\left(n=\overline{\mathrm{1,}N}\right)$ РС являются известными и характеризуются векторами ${\overrightarrow{r}}^{\text{0}}$ и ${\overrightarrow{r}}^n$ соответственно.

Вектор параметров ${\overrightarrow{\lambda }}^{\text{н}}\left(t\right)=\left\{{\overrightarrow{r}}^{\text{н}}\left(t\right)\right\}$ является неизвестным векторным случайным процессом, который задает положение ${\overrightarrow{r}}^{\text{н}}\left(t\right)$ в момент времени t ИНР. Реперные станции n-е передают тестовый сигнал $\zeta \left(t\right)$ в k-х $\left(k=\overline{\mathrm{1,}K}\right)$ направлениях на соответствующий k-й СР с частотой $f_{nk}^{\text{вв}}\left(t\right)$. Передача сигнала выполняется синхронно с помощью меток времени спутниковых радионавигационных систем – ГЛОНАСС и/или GPS. На входе приемника k-го СР наблюдается сигнал ${\tilde{\zeta }}_{nk}\left(t\right)$, различие которого от исходного $\zeta \left(t\right)$ обуславливается воздействием помех среды распространения радиосигнала и смещения $f_{nk}^{\text{вв}}\left(t\right)$ на величину доплеровского сдвига частоты $\Delta f_{nk}^{\text{вв}}\left({\overrightarrow{\lambda }}^k,t\right)$ на участке распространения от n-й РС до k-го СР – $\left|{\overrightarrow{r}}^n-{\overrightarrow{r}}^k\left(t\right)\right|$. Ретрансляция сигнала $\zeta \left(t\right)$ n-й РС k-м спутником связи осуществляется на частоте $f_{nk}^{\text{вв}}\left(t\right)+\Delta f_{nk}^{\text{вв}}\left({\overrightarrow{\lambda }}^k,t\right)+\Delta f_k^{\text{ср}}\left(t\right)$, где $\Delta f_k^{\text{ср}}\left(t\right)$ – частота переноса k-го СР на линию вниз, которая формируется из суммарной частоты всех гетеродинов приемника и передатчика k-го СР. На входе приемника КРК сигнал $\zeta \left(t\right)$ n-й РС наблюдается на частоте $f_{nk}^{\text{вв}}\left(t\right)+\Delta f_{nk}^{\text{вв}}\left({\overrightarrow{\lambda }}^k,t\right)+\Delta f_k^{\text{ср}}\left(t\right)+\Delta f_{nk}^{\text{вн}}\left({\overrightarrow{\lambda }}^k,t\right)$, где $\Delta f_{nk}^{\text{вн}}\left({\overrightarrow{\lambda }}^k,t\right)$ – доплеровский сдвиг частоты на участке распространения от k-го СР до КРК.

Введем обозначение $\overrightarrow{\lambda }\left(t\right)=\left\{{\overrightarrow{\lambda }}^{\text{н}}\left(t\right),{\overrightarrow{\lambda }}^1\left(t\right)\mathrm{,...,}{\overrightarrow{\lambda }}^K\left(t\right)\right\}$. ИНР осуществляет передачу радиосигнала $s\left(t\right)$ через 1-й СР – направление основного излучения. Направления передачи между ИНР и остальными СР (смежные СР) являются побочными. На входе приемника k-го $\left(k=\overline{\mathrm{1,}K}\right)$ СР наблюдается сигнал ${\tilde{s}}_k\left(t\right)$, различие которого от исходного $s\left(t\right)$ обуславливается воздействием помех среды распространения радиосигнала и смещения $f_{}^{\text{нвв}}\left(t\right)$ на величину доплеровского сдвига частоты $\Delta f_k^{\text{нвв}}\left(\overrightarrow{\lambda },t\right)$ с учетом протяженности k-го интервала на линии вверх – $\left|{\overrightarrow{r}}^{\text{н}}\left(t\right)-{\overrightarrow{r}}^k\left(t\right)\right|$. Сигнал $\zeta \left(t\right)$ РС формируется таким образом, чтобы его частотный диапазон не содержал общих частот из частотного диапазона сигнала $s\left(t\right)$ ИНР.

На выходе приемника КРК формируется суммарный сигнал:

$\begin{array}{l}\eta \left(t\right)=\sum _{k=1}^K\left(\begin{array}{l}\\ \end{array}\right.A_k^s\left(t\right)s\left(t-{{\rm T}}_k^s\left(\overrightarrow{\lambda },t\right)\right)e^{2\pi i\Delta F_k^s\left(\overrightarrow{\lambda },t\right)t}+\\ +\sum _{n=1}^N{A}_{nk}^{\zeta }\left(t\right)\zeta \left(t-{{\rm T}}_{nk}^{\zeta }\left({\overrightarrow{\lambda }}^k,t\right)\right)\times \\ \left.\times e^{2\pi i\Delta F_{nk}^{\zeta }\left({\overrightarrow{\lambda }}^k,t\right)t}\right)+\upsilon \left(t\right),\end{array}$

составленный из: сдвинутых по временному ${{\rm T}}_k^s\left({\overrightarrow{\lambda }}^k,t\right)$, ${{\rm T}}_{nk}^{\zeta }\left({\overrightarrow{\lambda }}^k,t\right)$ и частотному $\Delta F_k^s\left(\overrightarrow{\lambda },t\right)$, $\Delta F_{nk}^{\zeta }\left({\overrightarrow{\lambda }}^k,t\right)$ параметрам сигналов ИНР, n-х РС соответственно; затуханий $A_k^s\left(t\right)$, $A_{nk}^{\zeta }\left(t\right)$ сигналов на участках распространения ИНР–k-й СР–КРК, n-я РС–k-й СР–КРК; шумового параметра $\upsilon \left(t\right)$ который в рассматриваемой математической модели непрерывного канала связи определяет помехи среды распространения на линиях вверх и вниз.

Величины времени прихода ${{\rm T}}_k^s\left(\overrightarrow{\lambda },t\right)$, ${{\rm T}}_{nk}^{\zeta }\left({\overrightarrow{\lambda }}^k,t\right)$ сигналов $s\left(t\right), \zeta \left(t\right)$ от ИНР на КРК, от n-й РС на КРК складываются ${{\rm T}}_k^s\left(\overrightarrow{\lambda },t\right)={\tau }_k^s\left(\overrightarrow{\lambda },t\right)+{\delta }_k^{s\tau }\left(t\right)$, ${{\rm T}}_{nk}^{\zeta }\left({\overrightarrow{\lambda }}^k,t\right)={\tau }_{nk}^{\zeta }\left({\overrightarrow{\lambda }}^k,t\right)+{\delta }_{nk}^{\zeta \tau }\left(t\right)$ из ошибок ${\delta }_k^{s\tau }\left(t\right), {\delta }_{nk}^{\zeta \tau }\left(t\right)$ измерения и истинного времени   прихода на КРК сигналов от ИНР, n-й РС соответственно. Значения ${\delta }_k^{s\tau }\left(t\right), {\delta }_{nk}^{\zeta \tau }\left(t\right)$ удовлетворяют тождествам:

$\begin{array}{l}{\tau }_k^s\left(\overrightarrow{\lambda },t\right)={с}^{-1}\left(\left|{\overrightarrow{r}}^{\text{н}}\left(t\right)-{\overrightarrow{r}}^k\left(t\right)\right|+\left|{\overrightarrow{r}}^k\left(t\right)-{\overrightarrow{r}}_0\right|\right);\\ {\tau }_{nk}^{\zeta }\left({\overrightarrow{\lambda }}^k,t\right)={с}^{-1}\left(\left|{\overrightarrow{r}}^n-{\overrightarrow{r}}^k\left(t\right)\right|+\left|{\overrightarrow{r}}^k\left(t\right)-{\overrightarrow{r}}_0\right|\right),\end{array}$ (1)

где $c$ – скорость света в вакууме.

Ошибки ${\delta }_k^{s\tau }\left(t\right), {\delta }_{nk}^{\zeta \tau }\left(t\right)$ обуславливаются ионосферной и тропосферной ошибками измерения на линиях вверх, вниз и ошибками синхронизации передачи сигналов РС по меткам времени.

Величина частотных сдвигов $\Delta F_k^s\left(\overrightarrow{\lambda },t\right), \Delta F_{nk}^{\zeta }\left({\overrightarrow{\lambda }}^k,t\right)$ принимаемых КРК сигналов $s\left(t\right), \zeta \left(t\right)$ определяется

$\Delta F_k^s\left(\overrightarrow{\lambda },t\right)=f_{}^{\text{нвв}}\left(t\right)+\Delta f_k^{\text{нвв}}\left(\overrightarrow{\lambda },t\right)+\Delta f_k^{\text{ср}}\left(t\right)+$

null,

$\Delta F_{nk}^{\zeta }\left({\overrightarrow{\lambda }}^k,t\right)=f_{nk}^{\text{вв}}\left(t\right)+\Delta f_{nk}^{\text{вв}}\left({\overrightarrow{\lambda }}^k,t\right)+$

$+\Delta f_k^{\text{ср}}\left(t\right)+\Delta f_k^{\text{вн}}\left({\overrightarrow{\lambda }}^k,t\right)+\Delta f_{}^{\text{крк}}\left(t\right)=$

$=\Delta f_{nk}^{\zeta }\left({\overrightarrow{\lambda }}^k,t\right)+{\delta }_{nk}^{\zeta f}\left(t\right),$

где $\Delta f_{}^{\text{крк}}\left(t\right)$ – суммарная частота всех гетеродинов приемника КРК; $\Delta f_k^s\left(\overrightarrow{\lambda },t\right)$, $\Delta f_{nk}^{\zeta }\left({\overrightarrow{\lambda }}^k,t\right)$ и ${\delta }_k^{sf}\left(t\right), {\delta }_{nk}^{\zeta f}\left(t\right),$ – истинный частотный сдвиг принанятых КРК сигналов $s\left(t\right), \zeta \left(t\right)$ и ошибки измерения $\Delta f_k^s\left(\overrightarrow{\lambda },t\right), \Delta f_{nk}^{\zeta }\left({\overrightarrow{\lambda }}^k,t\right)$ соответственно. Значение $\Delta f_k^s\left(\overrightarrow{\lambda },t\right), \Delta f_{nk}^{\zeta }\left({\overrightarrow{\lambda }}^k,t\right)$ удовлетворяет тождествам:

$\Delta f_k^s\left(\overrightarrow{\lambda },t\right)={с}^{-1}{\overrightarrow{V}}^k\left(t\right)\left({\tilde{f}}_{}^{\text{нвв}}\frac{{\overrightarrow{r}}^k\left(t\right)-{\overrightarrow{r}}^{\text{н}}\left(t\right)}{\left|{\overrightarrow{r}}^k\left(t\right)-{\overrightarrow{r}}^{\text{н}}\left(t\right)\right|}+{\tilde{f}}_k^{\text{нвн}}\frac{{\overrightarrow{r}}_0-{\overrightarrow{r}}^k\left(t\right)}{\left|{\overrightarrow{r}}_0-{\overrightarrow{r}}^k\left(t\right)\right|}\right);$ (2)

$\Delta f_{nk}^{\zeta }\left({\overrightarrow{\lambda }}^k,t\right)={с}^{-1}{\overrightarrow{V}}^k\left(t\right)\left({\tilde{f}}_{nk}^{\text{вв}}\frac{{\overrightarrow{r}}^k\left(t\right)-{\overrightarrow{r}}^n}{\left|{\overrightarrow{r}}^k\left(t\right)-{\overrightarrow{r}}^n\right|}+{\tilde{f}}_{nk}^{\text{вн}}\frac{{\overrightarrow{r}}_0-{\overrightarrow{r}}^k\left(t\right)}{\left|{\overrightarrow{r}}_0-{\overrightarrow{r}}^k\left(t\right)\right|}\right),$

где null, ${\tilde{f}}_{nk}^{\text{вв}}=\mathbb{E}\left[f_{nk}^{\text{вв}}\right]$

и ${\tilde{f}}_k^{\text{нвн}}=\mathbb{E}\left[f_k^{\text{нвн}}=f_{}^{\text{нвн}}+\Delta f_k^{\text{ср}}\right],$, ${\tilde{f}}_{nk}^{\text{вн}}=\mathbb{E}\left[f_{nk}^{\text{вн}}=f_{nk}^{\text{вв}}+\Delta f_{}^{\text{ср}}\right].$

Ошибки ${\delta }_k^{sf}\left(t\right), {\delta }_{nk}^{\zeta f}\left(t\right),$ обуславливаются нестабильностями гетеродинов передатчиков ИНР, n-х РС, k-х СР и приемников k-х СР, КРК.

Учитывая содержательное описание схемы фильтрации канонических параметров ИНР и введенную в [8] систему ограничений, сформулируем основные допущения, выделяемые при решении задачи фильтрации координат наземного ИНР: 1) ИНИ находится на поверхности Земли и является неподвижным объектом; 2) СР осуществляют прямую ретрансляцию сигнала ИНИ; 3) диапазон частот сигнала, передаваемого ИНИ в направлении основного излучения, известен; 4) положения объектов модели фильтрации координат наземного ИНИ, работающего через СР, не совпадают, а сами объекты модели представляются в качестве кинематически не взаимосвязанных материальных точек; 5) различие времени распространения сигналов ИНИ и РС до КРК через соответствующие СР значительно превышают обратную величину ширины спектра этих сигналов; 6) математическая модель пространственного движения k-го СР формируется в земной системе координат и задается в кинематических параметрах по правилам [8]:

$\frac{d{\overrightarrow{V}}_{ср}\left(t\right)}{dt}={\mathbf{M}}_{ct}\left(t\right)\left(\overrightarrow{\Psi }\left(t\right)-\frac{\mu }{{\left|{\overrightarrow{r}}_{ср}\left(t\right)\right|}^3}{\overrightarrow{r}}_{ср}\left(t\right)\right)-$

$-{\mathbf{M}}_{ct}\left(t\right){{\mathbf{M}}^{\mathbf{\text{'}}}}_{ct}^T\left(t\right){\overrightarrow{V}}_{ср}\left(t\right)+{\overrightarrow{\Delta }}_V\left(t\right);$

$\begin{array}{l}d{\overrightarrow{r}}_{ср}\left(t\right)/dt=\\ ={\overrightarrow{V}}_{ср}\left(t\right)-{\mathbf{M}}_{ct}\left(t\right){{\mathbf{M}}^{\mathbf{\text{'}}}}_{ct}^T\left(t\right){\overrightarrow{r}}_{ср}\left(t\right)+{\overrightarrow{\Delta }}_r\left(t\right),\end{array}$ (3)

где $\mu$ – гравитационный параметр Земли; ${\overrightarrow{\Delta }}_V\left(t\right)$ и ${\overrightarrow{\Delta }}_r\left(t\right)$ – векторные функции, которые характеризуют случайный характер изменения вектора скорости и положения СР соответственно; $\overrightarrow{\Psi }\left(t\right)$ – возмущающее ускорение, обусловленное действием сил на СР; ${{\mathbf{M}}^{\text{'}}}_{ct}^T\left(t\right)$ – транспонированная матрица ${{\mathbf{M}}^{\text{'}}}_{ct}^{}\left(t\right)=d{\mathbf{M}}_{ct}\left(t\right)/dt; {\mathbf{M}}_{ct}\left(t\right)$ – матрица преобразования небесной системы координат в земную; 7) ошибки измерений ${\delta }_k^{s\tau }\left(t\right), {\delta }_{nk}^{\zeta \tau }\left(t\right), {\delta }_k^{sf}\left(t\right), {\delta }_{nk}^{\zeta f}\left(t\right)$ и эволюций ${\overrightarrow{\Delta }}_{V_k}\left(t\right), {\overrightarrow{\Delta }}_{r_k}\left(t\right)$, характеризующие флуктуацию случайных параметров, являются случайными величинами, которые распределены по нормальному закону; 8) особенности реализации многостанционного доступа к частотно-энергетическому ресурсу СР не рассматриваются.

2. Математическая постановка задачи фильтрации

Математическая модель фильтрации координат ИНР задается при представлении $\overrightarrow{\lambda }\left(t\right)$ непрерывным векторным марковским процессом. Эволюцию элементов вектора параметров $\overrightarrow{\lambda }\left(t\right)=\left\{{\overrightarrow{\lambda }}^{\text{н}}\left(t\right),{\overrightarrow{\lambda }}^1\left(t\right)\mathrm{,...,}{\overrightarrow{\lambda }}^K\left(t\right)\right\}$ в общем виде определим уравнением состояния:

${\dot{\overrightarrow{\lambda }}}_t=\overrightarrow{\Phi }\left({\overrightarrow{\lambda }}_t,t\right)+{\overrightarrow{w}}_t,$ (4)

где ${\overrightarrow{\lambda }}_t\equiv \overrightarrow{\lambda }\left(t\right)$; ${\dot{\overrightarrow{\lambda }}}_t\equiv d{\overrightarrow{\lambda }}_t/dt$; $\overrightarrow{\Phi }\left({\overrightarrow{\lambda }}_t,t\right)$ – вектор сноса; ${\overrightarrow{w}}_t\equiv \overrightarrow{w}\left(t\right)$ – белый гауссовский вектор-шум при $\mathbb{E}\left[{\overrightarrow{w}}_t\right]=0$ и $\mathrm{cov}\left[{\overrightarrow{w}}_t\right]={\mathbf{Q}}_t,$ наблюдаемый по правилу

${\overrightarrow{\eta }}_t=\overrightarrow{\Upsilon }\left({\overrightarrow{\lambda }}_t,t\right)+{\overrightarrow{v}}_t,$ (5)

где ${\eta }_t\equiv \eta \left(t\right)$; ${\overrightarrow{v}}_t\equiv \overrightarrow{v}\left(t\right)$ – белый гауссовский вектор-шум при $\mathbb{E}\left[{\overrightarrow{W}}_t\right]=0$ и $\mathrm{cov}\left[{\overrightarrow{v}}_t\right]={\mathbf{R}}_t;$ $\overrightarrow{\Upsilon }\left({\overrightarrow{\lambda }}_t,t\right)$ – функция измерения канонических параметров СР и положения ИНР по принимаемым сигналам РС и ИНР, определяемая по методу Ньютона [12] численным решением системы уравнений, составленных тождествами (1), (2).

С учетом введенной постановки задачи фильтрации в соотношениях состояния (4) и наблюдения (5) вектор сноса $\overrightarrow{\Phi }\left({\overrightarrow{\lambda }}_t,t\right)$ и вектор-шум ${\overrightarrow{w}}_t$ для:

$\overrightarrow{\lambda }\left(t\right)=\left\{{\overrightarrow{\lambda }}^{\text{н}}\left(t\right),{\overrightarrow{\lambda }}^1\left(t\right)\mathrm{,...,}{\overrightarrow{\lambda }}^K\left(t\right)\right\}$ при ${\overrightarrow{\lambda }}^{\text{н}}\left(t\right)=\left\{{\overrightarrow{r}}^{\text{н}}\left(t\right)\right\}$;

${\overrightarrow{r}}^{\text{н}}\left(t\right)=\left\{x^{\mathrm{н}}\left(t\right), y^{\mathrm{н}}\left(t\right), z^{\mathrm{н}}\left(t\right)\right\}$;

${\overrightarrow{\lambda }}^k\left(t\right)=\left\{{\overrightarrow{r}}^k\left(t\right), {\overrightarrow{V}}^K\left(t\right)\right\}$;

${\overrightarrow{r}}^k\left(t\right)=\left\{x^k\left(t\right), y^k\left(t\right), z^k\left(t\right)\right\}$;

${\overrightarrow{V}}^k\left(t\right)=\left\{V_x^k\left(t\right), V_y^k\left(t\right), V_z^k\left(t\right)\right\}$,

определяются по правилу:

$\overrightarrow{\Phi }\left({\overrightarrow{\lambda }}_t,t\right)=\left(\begin{array}{c}\overrightarrow{0}\\ {\overrightarrow{V}}^1\left(t\right)-{\mathbf{M}}_{ct}\left(t\right){{\mathbf{M}}^{\mathbf{\text{'}}}}_{ct}^T\left(t\right){\overrightarrow{r}}^1\left(t\right)\\ \begin{array}{c}{\mathbf{M}}_{ct}\left(t\right)\left({\overrightarrow{\Psi }}^1\left(t\right)-\mu {\left|{\overrightarrow{r}}^1\left(t\right)\right|}^{-3}{\overrightarrow{r}}^1\left(t\right)\right)-\\ -{\mathbf{M}}_{ct}\left(t\right){{\mathbf{M}}^{\mathbf{\text{'}}}}_{ct}^T\left(t\right){\overrightarrow{V}}^1\left(t\right)\end{array}\\ \mathrm{...}\\ {\overrightarrow{V}}^K\left(t\right)-{\mathbf{M}}_{ct}\left(t\right){{\mathbf{M}}^{\mathbf{\text{'}}}}_{ct}^T\left(t\right){\overrightarrow{r}}^K\left(t\right)\\ \begin{array}{c}{\mathbf{M}}_{ct}\left(t\right)\left({\overrightarrow{\Psi }}^K\left(t\right)-\mu {\left|{\overrightarrow{r}}^K\left(t\right)\right|}^{-3}{\overrightarrow{r}}^K\left(t\right)\right)-\\ -{\mathbf{M}}_{ct}\left(t\right){{\mathbf{M}}^{\mathbf{\text{'}}}}_{ct}^T\left(t\right){\overrightarrow{V}}^K\left(t\right)\end{array}\end{array}\right);$ (6)

$\begin{array}{l}{\overrightarrow{w}}_t=\left(\begin{array}{cccc}{\overrightarrow{\Delta }}_{{\overrightarrow{r}}^{н}}\left(t\right)& {\overrightarrow{\Delta }}_{{\overrightarrow{r}}^1}\left(t\right)& {\overrightarrow{\Delta }}_{{\overrightarrow{V}}^1}\left(t\right)& \mathrm{...}\end{array}\right.\\ {\left.\begin{array}{cc}{\overrightarrow{\Delta }}_{{\overrightarrow{r}}^K}\left(t\right)& {\overrightarrow{\Delta }}_{{\overrightarrow{V}}^K}\left(t\right)\end{array}\right)}^T,\end{array}$

а функция измерения $\overrightarrow{\Upsilon }\left({\overrightarrow{\lambda }}_t,t\right)$ задается решением системы из 3 уравнений и $K$ систем, включающей 6 уравнений. Каждая k-я система из $K$ систем уравнений составляется из тождеств (1), (2) при обозначениях

$G_{nk}^{\tau }\left(\overrightarrow{\lambda },t\right)=\left|{\overrightarrow{r}}^n-{\overrightarrow{r}}^k\left(t\right)\right|+\left|{\overrightarrow{r}}^k\left(t\right)-{\overrightarrow{r}}_0\right|-с{\tau }_{nk}^{\zeta }\left({\overrightarrow{\lambda }}^k,t\right),$

$G_{nk}^f\left(\overrightarrow{\lambda },t\right)={\overrightarrow{V}}^k\left(t\right)\left(\frac{{\displaystyle {\tilde{f}}_{nk}^{\mathrm{вв}}}}{∆f_{nk}^{\zeta }\left({\overrightarrow{\lambda }}^k, t\right)} \frac{{\overrightarrow{r}}^k\left(t\right){\displaystyle -}{\displaystyle {{\displaystyle \overrightarrow{r}}}^n}}{\left|{\overrightarrow{r}}^k\left(t\right)-{\overrightarrow{r}}^n\right|}+\frac{{\tilde{f}}_{nk}^{\mathrm{вн}}}{∆f_{nk}^{\zeta }\left({\overrightarrow{\lambda }}^k, t\right)} \frac{{\displaystyle {\overrightarrow{r}}_0}{\displaystyle -}{\displaystyle {{\displaystyle \overrightarrow{r}}}^k}{\displaystyle \left(t\right)}}{\left|{\overrightarrow{r}}_0-{\overrightarrow{r}}^k\left(t\right)\right|}\right)-с$,

задаваемых в виде

$\begin{array}{l}\left(\begin{array}{cccc}{G}_{n_{1}k}^{\tau }\left(\overrightarrow{\lambda },t\right)& G_{n_{2}k}^{\tau }\left(\overrightarrow{\lambda },t\right)& G_{n_{3}k}^{\tau }\left(\overrightarrow{\lambda },t\right)& G_{n_{1}k}^f\left(\overrightarrow{\lambda },t\right)\end{array}\right.\\ {\left.\begin{array}{cc}{G}_{n_{2}k}^f\left(\overrightarrow{\lambda },t\right)& G_{n_{3}k}^f\left(\overrightarrow{\lambda },t\right)\end{array}\right)}^{\text{T}}={\overrightarrow{G}}_k\left(\overrightarrow{\lambda },t\right)=\mathrm{0,}\end{array}$ (7)

где $n_1,n_2,n_3\in \left[1;N\right]:n_1\ne n_2\ne n_3$.

Элементы системы из трех уравнений ${\overrightarrow{\tilde{G}}}_{k_1{k}_2}\left(\overrightarrow{\lambda },t\right)={\left(\begin{array}{ccc}{\tilde{G}}_{k_1{k}_2}^{\tau }\left(\overrightarrow{\lambda },t\right)& {\tilde{G}}_{k_1{k}_2}^f\left(\overrightarrow{\lambda },t\right)& F_{\text{З}}\left(\overrightarrow{\lambda }\right)\end{array}\right)}^{\text{T}}$

в $\overrightarrow{\Upsilon }\left({\overrightarrow{\lambda }}_t,t\right)$ составляются с учетом обобщенной постановки задачи метода FDOA&TDOA [6; 7] в следующем виде:

$\begin{array}{l}{\tilde{G}}_{k_1{k}_2}^{\tau }\left(\overrightarrow{\lambda },t\right)=\left(\left|{\overrightarrow{r}}^{\text{н}}\left(t\right)-{\overrightarrow{r}}^{k_1}\left(t\right)\right|+\left|{\overrightarrow{r}}^{k_1}\left(t\right)-{\overrightarrow{r}}_0\right|\right)-\\ -\left(\left|{\overrightarrow{r}}^{\text{н}}\left(t\right)-{\overrightarrow{r}}^{k_2}\left(t\right)\right|+\left|{\overrightarrow{r}}^{k_2}\left(t\right)-{\overrightarrow{r}}_0\right|\right)-с\Delta {\tau }_{k_1{k}_2}^s\left(\overrightarrow{\lambda },t\right);\end{array}$ (8)

$\begin{array}{l}{\tilde{G}}_{k_1{k}_2}^f\left(\overrightarrow{\lambda },t\right)={\overrightarrow{V}}^{k_1}\left(t\right)\left({\tilde{f}}_{}^{\text{нвв}}{\left[\Delta f_{k_1{k}_2}^s\left(\overrightarrow{\lambda },t\right)\right]}^{-1}\times \right.\\ \left.\times \left[{\overrightarrow{r}}^{k_1}\left(t\right)-{\overrightarrow{r}}^{\text{н}}\left(t\right)\right]{\left|{\overrightarrow{r}}^{k_1}\left(t\right)-{\overrightarrow{r}}^{\text{н}}\left(t\right)\right|}^{-1}\right)-\\ -{\overrightarrow{V}}^{k_2}\left(t\right)\left({\tilde{f}}_{}^{\text{нвв}}{\left[\Delta f_{k_1{k}_2}^s\left(\overrightarrow{\lambda },t\right)\right]}^{-1}\times \right.\\ \left.\times \left[{\overrightarrow{r}}^{k_2}\left(t\right)-{\overrightarrow{r}}^{\text{н}}\left(t\right)\right]{\left|{\overrightarrow{r}}^{k_2}\left(t\right)-{\overrightarrow{r}}^{\text{н}}\left(t\right)\right|}^{-1}\right)-c,\end{array}$ (9)

где $k_1,k_2\in \left[1;K\right]:k_1\ne k_2$; $\Delta {\tau }_{k_1{k}_2}^s\left(\overrightarrow{\lambda },t\right)$ и $\Delta f_{k_1{k}_2}^s\left(\overrightarrow{\lambda },t\right)$ – разница по времени запаздывания и величин доплеровского смещения пар сигналов, принятых КРК от СР $k_1$ и $k_2$.

Приняв дополнительное ограничение о том, что на длительности $∆t$ элементы вектора ${\overrightarrow{\lambda }}_t$ при достаточно медленном изменении постоянны, разобьем интервал наблюдения $\left[0, T\right]$ на тактовые подынтервалы $t\in \left[t_j, t_{j+1}\right]$ $\left(∆t=t_{j+1}-t_j; j=\overline{1, N_t}; N_t=⌊T/∆t⌋\right)$ и перепишем уравнения состояния (4) и наблюдения (5) в виде

${\dot{\overrightarrow{\lambda }}}_j=\overrightarrow{\Phi }\left({\overrightarrow{\lambda }}_j,t_j\right)+{\overrightarrow{w}}_j;$ (10)

${\overrightarrow{\eta }}_j=\overrightarrow{\Upsilon }\left({\overrightarrow{\lambda }}_j,t_j\right)+{\overrightarrow{v}}_j,$ (11)

где ${\overrightarrow{\lambda }}_j\equiv \overrightarrow{\lambda }\left(t_j\right)$; ${\overrightarrow{\eta }}_j\equiv \overrightarrow{\eta }\left(t_j\right)$; ${\overrightarrow{w}}_j\equiv \overrightarrow{w}\left(t_j\right)$; ${\overrightarrow{v}}_j\equiv \overrightarrow{v}\left(t_j\right)$.

Решение системы нелинейных уравнений

$\overrightarrow{G}\left(\overrightarrow{\lambda },t\right)=\left\{{\overrightarrow{\tilde{G}}}_{k_1{k}_2}\left(\overrightarrow{\lambda },t\right),{\overrightarrow{G}}_{k_1}\left(\overrightarrow{\lambda },t\right),{\overrightarrow{G}}_{k_2}\left(\overrightarrow{\lambda },t\right)\right\}=\overrightarrow{0}$

предлагается выполнять методом Ньютона [12] при составлении приближения наблюдаемого вектора состояния ${\overrightarrow{\eta }}_j$ по положению ИНР и каноническим параметром СР $k_1$ и $k_2$ итерационным вычислением

${\overrightarrow{\eta }}_j^{\left[n+1\right]}={\overrightarrow{\eta }}_j^{\left[n\right]}-{\mathbf{H}}^{-1}\left({\overrightarrow{\eta }}_j^{\left[n\right]},t_j\right)\overrightarrow{G}\left({\overrightarrow{\eta }}_j^{\left[n\right]},t_j\right),$ (12)

где $n\in \left[0; N\right]$, $N$ – максимальное число итераций; ${\overrightarrow{\eta }}_j^{\left[0\right]}={\overrightarrow{\lambda }}_j$; $\mathbf{H}\left({\overrightarrow{\eta }}_j^{\left[n\right]},t_j\right)$ – матрица Якоби для системы функций $\overrightarrow{G}\left({\overrightarrow{\eta }}_j^{\left[n\right]},t_j\right)$. Матрица $\mathbf{H}\left(\overrightarrow{\lambda },t\right)$ задается в блочной форме:

$\mathbf{H}\left(\overrightarrow{\lambda },t\right)=\left(\begin{array}{ccc}\stackrel{\mathbf{~}}{\mathbf{H}}\left(\overrightarrow{\lambda },t\right)& {\tilde{\mathbf{H}}}_{k_1}\left(\overrightarrow{\lambda },t\right)& {\stackrel{\mathbf{~}}{\mathbf{H}}}_{k_2}\left(\overrightarrow{\lambda },t\right)\\ \mathbf{O}& {\mathbf{H}}_{k_1}\left(\overrightarrow{\lambda },t\right)& \mathbf{O}\\ \mathbf{O}& \mathbf{O}& {\mathbf{H}}_{k_2}\left(\overrightarrow{\lambda },t\right)\end{array}\right),$ (13)

где ${\mathbf{H}}_{k1}\left(\overrightarrow{\lambda },t\right)$ и ${\mathbf{H}}_{k2}\left(\overrightarrow{\lambda },t\right)$ – матрицы Якоби размера $6\times 6$ для СР $k_1$ и $k_2$ составленные с использованием частных производных [8] для соотнесенных к соответствующим СР; $\mathbf{О}$ – нулевая матрица размера $6\times 6$; $\tilde{\mathbf{H}}\left(\overrightarrow{\lambda },t\right)$ – симметричная матрица размера $3\times 3$, элементы которой определяются через частные производные элементов вектора ${\overrightarrow{\tilde{G}}}_{k_1{k}_2}\left(\overrightarrow{\lambda },t\right)$ по переменным ${\overrightarrow{r}}^{\text{н}}$; ${\mathbf{H}}_{k1}\left(\overrightarrow{\lambda },t\right)$ и ${\mathbf{H}}_{k2}\left(\overrightarrow{\lambda },t\right)$ – матрицы размера $3\times 6$ элементы которых определяются через частные производные элементов вектора ${\overrightarrow{\tilde{G}}}_{k_1{k}_2}\left(\overrightarrow{\lambda },t\right)$ по каноническим параметрам ${\overrightarrow{\lambda }}^{k_1}$ и ${\overrightarrow{\lambda }}^{k_2}$ СР соответственно.

Прогнозирование вектора состояния ${\overrightarrow{\lambda }}_{j+1}$ в момент времени $t_{j+1}$ выполняется аналогично решением обыкновенного дифференциального уравнения ${\dot{\overrightarrow{\lambda }}}_j=\overrightarrow{\Phi }\left({\overrightarrow{\lambda }}_j,t_j\right)$ из (3.6) с применением явного метода Рунге – Кутты [13]:

$\begin{array}{l}{\overrightarrow{\lambda }}_{j+1}={\overrightarrow{\lambda }}_j+h\sum _{i=0}^s{b}_i{\kappa }_{ji};\\ {\kappa }_{ji}=\overrightarrow{\Phi }\left({\overrightarrow{\lambda }}_j+h\sum _{i^{\text{'}}=1}^{i-1}{a}_{i{i}^{\text{'}}}{\kappa }_{j{i}^{\text{'}}},t_j+c_{i}h\right),\\ {\kappa }_{j0}=\overrightarrow{\Phi }\left({\overrightarrow{\lambda }}_j,t_j\right),\end{array}$ (14)

где $h\le ∆t$ – величина шага по временной сетке; $s$ – число этапов численной схемы; $a_{i{i}^{\text{'}}},b_i,c_i$ – вещественные коэффициенты, которые могут быть заданы по правилам из [13].

Указанные методы приближенного вычисления функций экстраполяции $\overrightarrow{\Phi }\left({\overrightarrow{\lambda }}_j,t_j\right)$ и измерения $\overrightarrow{{\rm Y}}\left({\overrightarrow{\lambda }}_j,t_j\right)$ позволяют задать нелинейную задачу фильтрации канонических параметров ИНР в виде:

$\begin{array}{l}{\overrightarrow{\lambda }}_{j+1}={\overrightarrow{\lambda }}_j+{\overrightarrow{\Phi }}_{{\epsilon }_h}\left({\overrightarrow{\lambda }}_j,t_j\right)+{\overrightarrow{w}}_{j+1};\\ {\overrightarrow{\eta }}_{j+1}={\overrightarrow{\Upsilon }}_{{\epsilon }_N}\left({\overrightarrow{\lambda }}_{j+1},t_{j+1}\right)+{\overrightarrow{v}}_{j+1}.\end{array}$ (15)

По аналогии с [8] выполним линеаризацию сформированной модели (15) при разложении функции эволюции ${\overrightarrow{\lambda }}_j+{\overrightarrow{\Phi }}_{{\epsilon }_h}\left({\overrightarrow{\lambda }}_j,t_j\right)$ и функции измерения ${\overrightarrow{\Upsilon }}_{{\epsilon }_N}\left({\overrightarrow{\lambda }}_j,t_j\right)$ в ряд Тейлора в окрестности точки  по первому члену:

$\begin{array}{l}{\overrightarrow{\lambda }}_{j+1}{=}_j{\overrightarrow{\lambda }}_j+{\overrightarrow{B}}_j+{\overrightarrow{w}}_{j+1};\\ {\overrightarrow{\eta }}_{j+1}{=}_{j+1}{\overrightarrow{\lambda }}_{j+1}+{\overrightarrow{D}}_{j+1}+{\overrightarrow{v}}_{j+1},\end{array}$ (16)

где

${\mathit{A}}_j={\mathbf{J}}_M\left({\overrightarrow{\lambda }}_j,t_j\right)+\mathbf{E};$

${\overrightarrow{B}}_j={\overrightarrow{\Phi }}_{{\epsilon }_h}\left({\overrightarrow{\lambda }}_j,t_j\right)-{\mathbf{J}}_M\left({\overrightarrow{\lambda }}_j,t_j\right){\overrightarrow{\lambda }}_j;$

${\mathbf{J}}_M\left({\overrightarrow{\lambda }}_j,t_j\right)$ – матрица Якоби для системы функций ${\overrightarrow{\Phi }}_{{\epsilon }_h}\left({\overrightarrow{\lambda }}_j,t_j\right);$

${\overrightarrow{D}}_{j+1}=-{\mathbf{H}}^{-1}\left({\overrightarrow{\lambda }}_{j+1},t_{j+1}\right)\overrightarrow{G}\left({\overrightarrow{\lambda }}_{j+1},t_{j+1}\right)+$

$+{\mathbf{J}}_S\left({\overrightarrow{\lambda }}_{j+1},t_{j+1}\right){\overrightarrow{\lambda }}_{j+1};$

${\mathbf{C}}_{j+1}=\mathbf{E}-{\mathbf{J}}_S\left({\overrightarrow{\lambda }}_{j+1},t_{j+1}\right);$

${\mathbf{J}}_S\left({\overrightarrow{\lambda }}_{j+1},t_{j+1}\right)$ – матрица Якоби для системы функций ${\mathbf{H}}^{-1}\left({\overrightarrow{\lambda }}_{j+1},t_{j+1}\right)\overrightarrow{G}\left({\overrightarrow{\lambda }}_{j+1},t_{j+1}\right)$; $\mathbf{E}$ – единичная матрица. Элементы матриц ${\mathbf{J}}_M\left({\overrightarrow{\lambda }}_j,t_j\right)$, ${\mathbf{J}}_S\left({\overrightarrow{\lambda }}_{j+1},t_{j+1}\right)$ аналогично [8] предлагается вычислять численно с применением интерполяционного метода Лагранжа. Следует отметить, что представление матриц ${\mathbf{J}}_M\left({\overrightarrow{\lambda }}_j,t_j\right)$, ${\mathbf{J}}_M\left({\overrightarrow{\lambda }}_j,t_j\right)$ задается в блочной форме:

${\mathbf{J}}_M\left(\overrightarrow{\lambda },t\right)=\left(\begin{array}{ccc}\stackrel{\mathbf{~}}{\mathbf{E}}& \stackrel{\mathbf{~}}{\mathbf{O}}& \stackrel{\mathbf{~}}{\mathbf{O}}\\ \mathbf{O}& {\mathbf{J}}_M^{k_1}\left(\overrightarrow{\lambda },t\right)& \mathbf{O}\\ \mathbf{O}& \mathbf{O}& {\mathbf{J}}_M^{k_2}\left(\overrightarrow{\lambda },t\right)\end{array}\right);$ (17)

${\mathbf{J}}_S\left(\overrightarrow{\lambda },t\right)=\left(\begin{array}{ccc}{\stackrel{\mathbf{~}}{\mathbf{J}}}_S\left(\overrightarrow{\lambda },t\right)& {\stackrel{\mathbf{~}}{\mathbf{J}}}_S^{k_1}\left(\overrightarrow{\lambda },t\right)& {\stackrel{\mathbf{~}}{\mathbf{J}}}_S^{k_2}\left(\overrightarrow{\lambda },t\right)\\ \mathbf{O}& {\mathbf{J}}_S^{k_1}\left(\overrightarrow{\lambda },t\right)& \mathbf{O}\\ \mathbf{O}& \mathbf{O}& J_S^{k_2}\left(\overrightarrow{\lambda },t\right)\end{array}\right),$

где ${\mathbf{J}}_M^{k_1}\left(\overrightarrow{\lambda },t\right),{\mathbf{J}}_S^{k_1}\left(\overrightarrow{\lambda },t\right)$ и ${\mathbf{J}}_M^{k_2}\left(\overrightarrow{\lambda },t\right),{\mathbf{J}}_S^{k_2}\left(\overrightarrow{\lambda },t\right)$ – матрицы Якоби размера $6\times 6$ для СР $k_1$ и $k_2$ формируемые относительно уравнений состояния, наблюдения соответственно; $\tilde{\mathbf{O}}$ – нулевая матрица размера $3\times 6$; $\stackrel{\mathbf{~}}{\mathbf{E}}$ – единичная матрица $3\times 3$; $\tilde{\mathbf{H}}\left(\overrightarrow{\lambda },t\right)$ – симметричная матрица размера $3\times 3$, элементы которой определяются через частные производные элементов вектора ${\overrightarrow{\tilde{G}}}_{k_1{k}_2}\left(\overrightarrow{\lambda },t\right)$ по переменным ${\overrightarrow{r}}^{\text{н}};$ ${\mathbf{J}}_S^{k_1}\left(\overrightarrow{\lambda },t\right)$ и ${\mathbf{J}}_S^{k_2}\left(\overrightarrow{\lambda },t\right)$ – матрицы размера $3\times 6$ элементы которых определяются через частные производные соответствующих элементов вектора ${\mathbf{H}}^{-1}\left({\overrightarrow{\lambda }}_j,t_j\right)\overrightarrow{G}\left({\overrightarrow{\lambda }}_j,t_j\right)$ по каноническим параметрам ${\overrightarrow{\lambda }}^{k_1}$ и ${\overrightarrow{\lambda }}^{k_2}$ СР.

Для сформированной математической модели (15) с учетом правила (16) по ее линеаризации в калмановской модели фильтрации [14] правило оценки вектора состояния

${\widehat{\overrightarrow{\lambda }}}_{j+1\left|j+1\right.}={\widehat{\overrightarrow{\lambda }}}_{j+1\left|j\right.}+{\mathbf{K}}_{j+1}\left({\overrightarrow{\eta }}_{j+1}-{\overrightarrow{\Upsilon }}_{{\epsilon }_N}\left({\widehat{\overrightarrow{\lambda }}}_{j+1\left|j\right.},t_{j+1}\right)\right)$

в момент времени $t_{j+1}$ на интервале наблюдения $\left[t_j, t_{j+1}\right]$ где оптимальная, по Калману, матрица коэффициентов усиления ${\mathbf{K}}_{j+1}$ определится аналогичным способом, рассмотренным в [8].

3. Алгоритм фильтрации координат наземного ИНР

С учетом формализованных соотношений (15), (16) и представлений основу алгоритма фильтрации координат наземного ИНР в системе спутниковой связи с прямой ретрансляцией составляет правило рекурсивной оценки вектора состояния ${\widehat{\overrightarrow{\lambda }}}_{j+1\left|j+1\right.}$ СР, которое позволяет сформировать итоговый алгоритм для 2 СР в следующем виде.

Шаг 1. С учетом параметров принимаемых сигналов от n-х РC и ИНР k-х определить среднеквадратические ошибки оценки временных задержек ${\tau }_{nk}^{\zeta }\left({\overrightarrow{\lambda }}^k,t\right)$, $\Delta {\tau }_{k_1{k}_2}^s\left(\overrightarrow{\lambda },t\right)$ и доплеровских смещений $\Delta F_{nk}^{\zeta }\left({\overrightarrow{\lambda }}^k,t\right)$, $\Delta f_{k_1{k}_2}^s\left(\overrightarrow{\lambda },t\right)$ [3].

Шаг 2. Рекуррентным соотношением (12) вычислить приближения значений наблюдаемого по правилам (11) вектора состояния ИНР и пары спутников-ретрансляторов (основного и зеркального) ${\overrightarrow{\eta }}_j^{{\epsilon }_N}$ для $j=\overline{\mathrm{1,}{N}_t}$ и его отклоненного наблюдаемого значения ${\tilde{\overrightarrow{\eta }}}_j^{{\epsilon }_N}$ с учетом вычисленных значений среднеквадратических ошибок на первом шаге.

Шаг 3. С использованием вычислительной схемы Рунге – Кутты (14) выполнить экстраполяцию вектора состояния ИНР и пары СР ${\overrightarrow{\lambda }}_j^{{\epsilon }_h}$ для $j=\overline{\mathrm{1,}{N}_t}$ и с учетом оценок ошибок моделей прогнозирования положения Луны, Солнца и др. небесных тел, аппроксимации поверхности Земли геоидом при использовании ограниченного числа полиномов, неполноты учета возмущающих орбитальное движение СР сил [11] выполнить экстраполяцию отклоненного вектора состояния ИНР и пары СР ${\tilde{\overrightarrow{\lambda }}}_j^{{\epsilon }_h}$.

Шаг 4. Вычислить элементы $Q_{l_1,l_2}^j$, $R_{l_1,l_2}^j$ $(l_1,l_2\in \left\{\overline{\mathrm{1,15}}\right\})$ соответствующих ковариационных матриц ${\mathbf{Q}}_j, {\mathbf{R}}_j$ по правилам из [15]:

$Q_{l_1,l_2}^j=\left\{\begin{array}{l}0\text{ if }\left(l_1\le 3\right)\vee \left(l_2\le 3\right),\\ {\tilde{Q}}_{l_1,l_2}^j\text{ otherwise;}\end{array}\right.$

${\mathbf{R}}_{\mathbf{j}}\mathbf{=}\left(\begin{array}{ccc}{\tilde{R}}_j^{}& \tilde{O}& \tilde{O}\\ O& R_j^{k_1}& O\\ O& O& R_j^{k_2}\end{array}\right)\mathbf{;}$ (18)

$\begin{array}{l}{\tilde{Q}}_{l_1,l_2}^j=\frac{2}{N_t\left(N_t-1\right)}\sum _{j^{\text{'}}=1}^{N_t}\left(N_t-\left|j-j^{\text{'}}\right|\right)\times \\ \times \left({\left({\tilde{\overrightarrow{\lambda }}}_j^{{\epsilon }_h}\right)}_{l_1}-{\left({\overrightarrow{\lambda }}_j^{{\epsilon }_h}\right)}_{l_1}\right)\left({\left({\tilde{\overrightarrow{\lambda }}}_j^{{\epsilon }_h}\right)}_{l_2}-{\left({\overrightarrow{\lambda }}_j^{{\epsilon }_h}\right)}_{l_2}\right),\end{array}$

где ${\tilde{\mathbf{R}}}_j$, ${\mathbf{R}}_j^{k_1}$ и ${\mathbf{R}}_j^{k_2}$ – блочные матрицы размерности $3\times 3$, $6\times 6$ и $6\times 6$ соответственно, задаваемые из соответствующих блочных элементов матрицы:

$\begin{array}{l}{\tilde{R}}_{l_1,l_2}^j=\frac{2}{N_t\left(N_t-1\right)}\sum _{j^{\text{'}}=1}^{N_t}\left(N_t-\left|j-j^{\text{'}}\right|\right)\times \\ \times \left({\left({\tilde{\overrightarrow{\eta }}}_j^{{\epsilon }_N}\right)}_{l_1}-{\left({\overrightarrow{\eta }}_j^{{\epsilon }_N}\right)}_{l_1}\right)\left({\left({\tilde{\overrightarrow{\eta }}}_j^{{\epsilon }_N}\right)}_{l_2}-{\left({\overrightarrow{\eta }}_j^{{\epsilon }_N}\right)}_{l_2}\right).\end{array}$ (19)

Шаг 5. Вычислить ${\mathbf{A}}_j$, ${\overrightarrow{B}}_j$, ${\mathbf{C}}_j$, ${\overrightarrow{D}}_j$ в лианирезованной модели (16) фильтрации координат ИНР в ССС с прямой ретрансляцией для $j=\overline{\mathrm{1,}{N}_t}$.

Шаг 6. Положить $j=0; {\widehat{\overrightarrow{\lambda }}}_{j\left|j\right.}={\overrightarrow{\lambda }}_j^{{\epsilon }_h}; {\mathbf{P}}_{j|j}=0$.

Шаг 7. Выполнить прогнозирование оценки вектора состояния ${\widehat{\overrightarrow{\lambda }}}_{j+1\left|j\right.}={\widehat{\overrightarrow{\lambda }}}_{j\left|j\right.}+{\overrightarrow{\Phi }}_{{\epsilon }_h}\left({\widehat{\overrightarrow{\lambda }}}_{j\left|j\right.},t_j\right)$ и вычислить ковариационную матрицу $P_{j+1|j}$.

Шаг 8. Вычислить вектор отклонения ${\overrightarrow{\zeta }}_{j+1}={\overrightarrow{\eta }}_{j+1}-{\overrightarrow{\Upsilon }}_{{\epsilon }_N}\left({\widehat{\overrightarrow{\lambda }}}_{j+1\left|j\right.},t_{j+1}\right)$ и рассчитать ковариационную матрицу ${\mathbf{S}}_{j+1}={\mathbf{B}}_{j+1}{\mathbf{P}}_{j+1|j}{\mathbf{C}}_{j+1}^T+{\mathbf{R}}_{j+1}$ для ${\overrightarrow{\zeta }}_{j+1}$.

Шаг 9. Вычислить оптимальную, по Калману, матрицу коэффициентов усиления ${\mathbf{K}}_{j+1}={\mathbf{P}}_{j+1|j}{\mathbf{C}}_{j+1}^T{\mathbf{S}}_{j+1}^{-1}$.

Шаг 10. Выполнить коррекцию ${\widehat{\overrightarrow{\lambda }}}_{j+1\left|j+1\right.}={\widehat{\overrightarrow{\lambda }}}_{j+1\left|j\right.}+{\mathbf{K}}_{j+1}{\overrightarrow{\zeta }}_{j+1}$ экстраполяции вектора состояния ${\widehat{\overrightarrow{\lambda }}}_{j+1\left|j\right.}$.

Шаг 11. По правилам калмановской фильтрации рассчитать ковариационную матрицу ${\mathbf{P}}_{j+1|j+1}$ оценки вектора состояния ${\widehat{\overrightarrow{\lambda }}}_{j+1\left|j+1\right.}$.

Шаг 12. Проверить условие окончания алгоритма:

1. если $j+1<N_t$, положить $j=j+1$ и перейти к шагу 7;
2. если $j+1=N_t$ завершить итерационный процесс и вывести полученный результат фильтрации вектора состояний ${\widehat{\overrightarrow{\lambda }}}_{j\left|j\right.}$ ИНР и пары спутников ретрансляторов для $j=\overline{\mathrm{1,}{N}_t}$.

Шаг 13. Для полученных значений отфильтрованного вектора состояния ${\widehat{\overrightarrow{\lambda }}}_{j\left|j\right.}$ по правилу (11) выполнить оценку $N_t$ значений ${\widehat{\overrightarrow{r}}}_{\text{н}}^j$ возможного местоположения ИНР.

Шаг 14. Выполнить усреднение ${\widehat{\overrightarrow{r}}}^{\text{н}}=\frac{1}{N_t}\sum _{j=1}^{N_t}{\widehat{\overrightarrow{r}}}_{\text{н}}^j$ оценок $N_t$ по местоположению ИНР, а результат ${\widehat{\overrightarrow{r}}}^{\text{н}}$ считать итоговым и завершить работу алгоритма.

4. Верификация сформированных решений

Для практического обоснования результативности сформированного итогового алгоритмического решения фильтрации координат наземного ИНР в ССС с прямой ретрансляцией в Mathcad проведено математическое моделирование. Исходные данные вычислительного эксперимента выбраны следующими.

1. Начальное время и дата моделирования 2 часов 26 минут 15 секунд 7.04.2019 г. Длительность измерений 1 день.
2. Параметры СР:
   
   1. Основной СР – Ямал-402, подспутниковая точка 54,9 град в.д. (номер спутника в базе NORAD 39022); номер витка на момент эпохи 2324; эксцентриситет орбиты 0,0002895; наклон орбиты 0,00010821041362365 рад; долгота восходящего узла 4,396838687611871 рад; аргумент перигея 2,036889994198486 рад; средняя аномалия 0,9513335967015571 рад; частота обращения 1,00274382 об./день; значения первой и второй производных среднего движения по времени –0,00000082 и 0 соответственно.
   2. Зеркальный СР – KazSat 3, подспутниковая точка 58,5 град в.д. (номер спутника в базе NORAD 39728); номер витка на момент эпохи 1808; эксцентриситет орбиты 0,0000279; наклон орбиты 0,00023561944901923 рад; долгота восходящего узла 3,7041890659946533 рад; аргумент перигея 0,7237583702047645 рад; средняя аномалия 0,6355110325654273 рад; частота обращения 1,00264004 об./день; значения первой и второй производных среднего движения по времени 0,0000006 и 0 соответственно.
3. Параметры КРК: широта 52°57’54’’; долгота 36°04’42’’.
4. Параметры РС: число РС $K=3$; ${\tilde{f}}_n^{\text{вв}}=\text{14250}$ Мгц; ${\tilde{f}}_n^{\text{вн}}=\text{11000}$ МГц; $B=48$ кГц; широта РС №1 55°45’07’’; долгота РС №1 37°36’56’’; широта РС №2 54°50’37’’; долгота РС №2 38°11’10’’; широта РС №3 56°01’06’’; долгота РС №3 92^o52’01’’; ОСШ реперных станций 0 дБ; длительности тактового подынтервала оценки $∆t=90°$с.
5. Параметры модели состояния СР: поверхность Земли аппроксимируется геоидом EGM2008; возмущение параметров орбиты в модели состояния СР задано притяжением центрированного поля Земли, Луны и Солнца; возмущение параметров орбиты в реальной модели состояния СР дополнено давлением солнечного излучения и нецентрированным притяжением Земли.
6. Параметры ИНР: широта 52°57’54’’; долгота 36°00’51’’; ${\tilde{f}}_{}^{\text{нвв}}=\text{14500}$ МГц; ${\tilde{f}}_k^{\text{нвн}}=\text{10750}$ МГц; $B=256$ кГц; КУ антенны в направлении основного СР 45 дБ; КУ антенны в направлении на зеркальный СР 20 дБ.

Для выбранных параметров среднеквадратическая ошибка за время наблюдения модели состояния: 1) основного СР по положению ${\sigma }_r=70,04$ м, по скорости ${\sigma }_V=1,275·{10}^{-6}$ м/с; 2) зеркального СР по положению ${\sigma }_r=60,672$ м, по скорости ${\sigma }_V=1,044·{10}^{-6}$ м/с.

На рис. 2 приведены результаты оценок местоположения ИНР в ССС с прямой ретрансляцией при определении положения СР при фильтрации канонических параметров СР [8], предполагая что задача местоопределения ИНР решается раздельно – существующий алгоритм.

 

Рис. 2. Оценки местоположения ИНР в ССС с прямой ретрансляцией при раздельном решении

Fig. 2. Estimates of the location of the INR in the CCS with direct retransmission with a separate solution

 

Для полученных оценок местоположения ИНР средняя ошибка расчетов в случае раздельного решения – 5753,25 м, а усредненное – 307,83 м.

На рис. 3 и 4 приведены сравнительные примеры за время наблюдения ошибки фильтрации канонических параметров основного и зеркального СР соответственно по результатам работы итогового алгоритма фильтрации. Для задач фильтрации получены следующие значения среднеквадратической ошибки за время наблюдения: 1) основного СР по положению ${\sigma }_r=\mathrm{30,685}$ м, по скорости ${\sigma }_V=\mathrm{0,647}\cdot {10}^{-6}$ м/с; 2) зеркального СР по положению ${\sigma }_r=\mathrm{23,075}$ м, по скорости ${\sigma }_V=\mathrm{0,336}\cdot {10}^{-6}$ м/с.

 

Рис 3. Графики ошибки фильтрации канонических параметров основного СР: а – положение; б – скорость

Fig 3. Graphs of the filtering error of the canonical parameters of the main SR: a – position; b – speed

 

Рис. 4. Графики ошибки фильтрации канонических параметров зеркального СР: а – положение; б – скорость

Fig. 4. Graphs of the filtering error of the canonical parameters of the mirror SR: a – position; b – speed

 

На рис. 5 представлен результат фильтрации координат ИНР в ССС с прямой ретрансляцией при работе итогового алгоритма.

 

Рис. 5. Оценка местоположения ИНР в ССС с прямой ретрансляцией по разработанному алгоритму фильтрации

Fig. 5. Evaluation of the location of the INR in the CCS with direct retransmission according to the developed filtering algorithm

 

Для полученных оценок местоположения ИНР с применением итогового алгоритма фильтрации средняя ошибка расчетов от истинного составила 5037,69 м, а усредненное – 264,3 м.

Заключение

В целом результаты проведенного моделирования подтверждают работоспособность и результативность разработанного алгоритмического решения задачи фильтрации координат ИНР в ССС с прямой ретрансляцией. По данным апостериорных исследований, сформированный итоговый алгоритм фильтрации относительно известных решений в среднем позволяет снизить среднеквадратическую ошибку оценки в 1,34 раза.

Следует отметить, что полученные результаты определяют особенность построения итогового алгоритма, в котором сигнал ИНР при его неподвижности выступает в качестве дополнительной реперной станции относительно разработанной в [8] математической модели фильтрации канонических параметров СР при орбитальном движении.

About the authors

Ivan S. Polyansky

The Academy of the Federal Guard Service of the Russian Federation

Email: van341@mail.ru
ORCID iD: 0000-0002-1282-1522

Doctor of Physical and Mathematical Sciences, member of the Academy of the Federal Guard Service of the Russian Federation, Oryol, Russia. The number of scientific publications – 146.

Research interests: mathematical modeling, dynamic systems, differential equations, optimization methods, optimal control, conformal mapping, computational electrodynamics, digital signal processing.

Russian Federation, 35, Priborostroitelnaya Street, Oryol, 302015

Tuan Zap Pham

The Academy of the Federal Guard Service of the Russian Federation

Email: giapptvn1412@gmail.com

adjunct of the Academy of the Federal Guard Service of the Russian Federation, Oryol, Russia. Number of scientific publications – 8.

Research interests: radio engineering, theory of filtering coordinates of moving objects, algorithms of digital signal processing on PLD.

Russian Federation, 35, Priborostroitelnaya Street, Oryol, 302015

Alexey V. Tikhonov

The Academy of the Federal Guard Service of the Russian Federation

Email: sec@academ.msk.rsnet.ru

Candidate of Technical Sciences, member of the Academy of the Federal Guard Service of the Russian Federation, Oryol, Russia. The number of scientific publications – 78.

Research interests: radio engineering, theory of filtering coordinates of moving objects, digital signal processing, theory of electrical communication.

Russian Federation, 35, Priborostroitelnaya Street, Oryol, 302015

Boris G. Katygin

The Academy of the Federal Guard Service of the Russian Federation

Author for correspondence.
Email: katygin.b@yandex.ru

Candidate of Technical Sciences, member of the Academy of the Federal Guard Service of the Russian Federation, Oryol, Russia. The number of scientific publications – 42.

Research interests: radio engineering, theory of electrical communication, optimal control, optimization methods.

Russian Federation, 35, Priborostroitelnaya Street, Oryol, 302015

References

Supplementary files