Nonlinear stochastic estimation of the navigation parameters of the antenna of a mobile radio engineering complex by inertial-satellite measurements

Full Text

1. Введение. Для современных подвижных радиотехнических комплексов (РТК) точность пространственной ориентации их антенн, размещенных на мачтах различной конструкции, расположенных на объекте, является одним из основных факторов, определяющих эффективность функционирования подвижных РТК. Это, в свою очередь, выдвигает в качестве одной из центральных задач повышение точности определения ориентации антенн как в условиях возмущенного движения мачты, так и при произвольном характере движения объекта, причем, при неизбежных помехах измерения параметров движения как антенны, так и объекта [1$–$3]. Актуальность этой задачи привела к интенсивным исследованиям в данном направлении и разработке ряда методов оценки пространственной ориентации антенны в условиях, как внешних возмущений, так и внутренних помех измерительного комплекса [4$–$9]. Одним из широко распространенных является метод, использующий спутниковые измерения, применение которого позволяет решать одновременно как задачу пространственной ориентации, так и задачу определения текущих координат подвижного объекта [7$–$13]. Главным преимуществом данного метода является отсутствие операции непрерывного интегрирования измерений чувствительных элементов (ЧЭ), приводящей к накоплению ошибок с течением времени (характерной для непосредственной обработки инерциальных измерений). В то же время низкая частота спутниковых сообщений наряду с высокой интенсивностью помех при их приеме существенно затрудняют использование подобных методов для решения задачи определения ориентации антенны, находящейся на высокодинамичном подвижном основании [14$–$18].

Альтернативным направлением определения ориентации антенн является использование показаний бесплатформенных инерциальных систем ориентации (БИСО) [18$–$23]. Но в большинстве современных методов, использующих БИСО, не учитывается динамика изменения параметров ориентации, измеряемых в условиях интенсивных помех, что не позволяет достичь требуемой точности и устойчивости процесса определения ориентации антенны подвижного РТК [1$–$3]. В работе [22] был построен алгоритм оценки ориентации антенны стационарного РТК с учетом помех измерения БИСО в предположении, что антенна может изменять (в общем случае случайным образом) свою ориентацию относительно мачты. Реализация данного алгоритма, решающего задачу определения ориентации антенны на качающейся мачте, требует привлечения в составе БИСО трех датчиков угловой скорости (ДУС) и трех акселерометров. Но при расположении мачты с антенной на подвижном объекте определить ее ориентацию с данным составом ЧЭ БИСО и с использованием алгоритма [22] уже не удается.

В связи с этим рассмотрим далее для мачтовой антенны подвижного РТК синтез алгоритма оценки ее параметров ориентации, инвариантного к характеру движения, как подвижного объекта, так и мачты относительно объекта, и обеспечивающего устойчивость и требуемую точность оценивания ориентации при самых общих предположениях о характере помех ЧЭ, используемых в данном случае в составе БИСО.

2. Постановка задачи. Для последующего решения задачи автономной ориентации антенны на подвижном основании по измерениям БИСО в общей постановке полагаем далее, что центр масс (ЦМ) антенны соединен жестким стержнем длиной R (моделирующим мачту) с ЦМ подвижного объекта и может вращаться вокруг него с произвольной угловой скоростью во всех направлениях под действием внешних возмущений (рис. 1). В свою очередь, объект перемещается по сфере Земли с переменной высотой и совершает произвольное вращение относительно его ЦМ.

 

Рис. 1. Системы координат

 

Также введем следующие системы координат (СК) (рис. 1):

• первую приборную СК J₁ (ПСК1) Ox₁y₁z₁ с началом в ЦМ объекта, оси которой направлены по взаимно ортогональным осям чувствительности ЧЭ первой группы, входящих в общий состав БИСО;
• вторую приборную СК J₂ (ПСК2) Ox₂y₂z₂ с началом в ЦМ антенны, оси которой направлены по взаимно ортогональным осям чувствительности ЧЭ второй группы, входящих в состав БИСО;
• сопровождающую СК S (ССК) OX_SY_SZ_S с началом в ЦМ объекта (одновременно $–$ в точке крепления стержня длиной R), ось OY_S которой лежит в плоскости местного меридиана и направлена на Север, ось OZ_S направлена от центра Земли, а ось OX_S дополняет СК до правой;
• опорную СК Q (ОСК) OXYZ с началом в ЦМ антенны, ось OZ которой направлена вдоль стержня длиной R, направления осей OX, OY, OZ в начальный момент времени совпадают с направлениями соответствующих осей ССК OX_S, OY_S, OZ_S;
• инерциальную СК I (ИСК) Oξηξ с началом в центре Земли, ось Oη которой направлена по оси вращения Земли от ее центра, ось Oξ в начальный момент времени лежит в плоскости нулевого меридиана, а ось Oξ дополняет СК до правой;
• геоцентрическую СК G (ГСК) OX_GY_GZ_G с началом в центре Земли, ось OY_G которой направлена по оси вращения Земли от ее центра, ось OZ_G лежит в плоскости нулевого меридиана, а ось OX_G дополняет систему координат до правой.

В соответствии с введенными СК, под задачей ориентации антенны, расположенной на качающейся мачте, размещенной, в свою очередь, на подвижном объекте, далее понимается текущая оценка параметров разворота (в качестве которых далее рассматриваются параметры Родрига$–$Гамильтона [23]) ПСК2 J₂ относительно ССК S при одновременном определении текущих координат (долготы λ и широты φ) ЦМ подвижного объекта в ГСК.

Полагаем также, что в состав БИСО входят две группы ЧЭ, состоящих из трех акселерометра и трех ДУС, расположенных ортогонально: первая группа $–$ в ЦМ объекта, вторая группа $–$ в ЦМ антенны. С целью сохранения общности решения в качестве моделей помех ЧЭ БИСО выберем аддитивные белые гауссовские шумы (БГШ) с нулевыми математическими ожиданиями и известными интенсивностями как наиболее адекватные практике использования БИСО. Учет корреляции помех или наличия в них регулярных составляющих (в том числе с неизвестными параметрами) легко обеспечивается соответствующим расширением вектора оцениваемых параметров и не влияет на существо предлагаемого далее подхода [24]. В этом случае модели выходных сигналов ЧЭ БИСО имеют вид [19$–$21]:

$–$ для акселерометров:

$Z_{ai}=a_J{}_i+W_{ai},$

где i = 1,2 $–$ номер группы ЧЭ, $Z_{ai}={\left|\begin{array}{ccc}{Z}_{axi}& Z_{ayi}& Z_{azi}\end{array}\right|}^T$ $–$ вектор выходных сигналов трех ортогональных акселерометров i-й группы; $a_J{}_i={\left|\begin{array}{ccc}{a}_{xi}& a_{yi}& a_{zi}\end{array}\right|}^T$ $–$ вектор ускорений ЦМ (объекта или антенны) в i-й ПСК; $W_a{}_i={\left|\begin{array}{ccc}{W}_{xi}& W_{yi}& W_{zi}\end{array}\right|}^T$ $–$ вектор помех измерения акселерометров i-й группы (центрированный БГШ с матрицей интенсивностей ai);

$–$ для ДУС:

$Z_{di}={\omega }_{Ji}+W_{di},$

где $Z_{di}={\left|\begin{array}{ccc}{Z}_{xi}& Z_{yi}& Z_{zi}\end{array}\right|}^T$ $–$ вектор выходных сигналов трех ортогональных ДУС i-й группы; ${\omega }_J{}_i={\left|\begin{array}{ccc}{\omega }_{xi}& {\omega }_{yi}& {\omega }_{zi}\end{array}\right|}^T$ $–$ вектор абсолютной угловой скорости вращения i-й ПСК; $W_d{}_i={\left|\begin{array}{ccc}{W}_{xi}& W_{yi}& W_{zi}\end{array}\right|}^T$ $–$ вектор помех измерения ДУС i-й группы (центрированный БГШ с матрицей интенсивностей di).

Таким образом, окончательно поставленную задачу можно сформулировать как задачу стохастического оценивания текущей ориентации ПСК2 J₂ относительно ССК S при одновременном определении текущих координат (долготы и широты) ЦМ подвижного объекта по зашумленным измерениям ЧЭ БИСО при apriori неопределенном характере изменения векторов угловой скорости ЦМ антенны относительно начала ЦМ объекта (точки крепления мачты длиной R), антенны относительно ее ЦМ и самого подвижного объекта относительно его ЦМ при неизвестных углах начального рассогласования ПСК2 J₂ и ССК S (т.е. неточно решенной задаче начальной выставки БИСО).

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

Текущая ориентация трехгранника S ССК относительно триэдра G ГСК описывается известными уравнениями [18, 20, 21, 25]:

 $\left|\begin{array}{l}\dot{\lambda }\\ \dot{\phi }\end{array}\right|=\left|\begin{array}{cc}0& {\left(\cos \phi \right)}^{-1}\\ 1& 0\end{array}\right|=\left|\begin{array}{l}{V}_{Y_S}\\ V_{X_S}\end{array}\right|{\left(r+H\right)}^{-1};$ $\lambda \left(0\right)={\lambda }_0,\phi \left(0\right),V_{Z_S},$ (3.1)

где $\lambda$$–$ долгота, $\phi$$–$ широта объекта, $V_{X_S}$, $V_{Y_S}$ $–$ проекции линейной скорости объекта на оси ССК, r $–$ радиус Земли, $H=H_0+{\displaystyle \underset{0}{\overset{t}{\int }}{V}_{Z_S}}\left(\tau \right)d\tau$ $–$ текущая высота объекта, $H_0$$–$ начальная высота, $V_{Z_S}$ − вертикальная скорость объекта в ССК.

Для синтеза уравнений неизвестных проекций $V_{X_S}$, $V_{Y_S}$ и $V_{Z_S}$ обратимся к основному уравнению инерциальной навигации [23, 25], записанному в ССК S для вектора ускорений $A_S={\left|\begin{array}{ccc}{A}_{X_S}& A_{Y_S}& A_{Z_S}\end{array}\right|}^T$, измеряемых акселерометрами при движении объекта по поверхности сферы Земли:

$A_S={\dot{V}}_S+\left(2{\Omega }_S+{\omega }_S\right)\times V_S-g_S,$ (3.2)

где $V_S={\left|\begin{array}{ccc}{V}_{X_S}& V_{Y_S}& V_{Z_S}\end{array}\right|}^T$, ${\Omega }_S={\left|\begin{array}{ccc}{\Omega }_{X_S}& {\Omega }_{Y_S}& {\Omega }_{Z_S}\end{array}\right|}^T$ − вектор угловой скорости вращения Земли, проекции которого на оси выбранной ССК имеют вид:

${\Omega }_{X_S}=\mathrm{0,}{\Omega }_{Y_S}=\Omega \text{cos}\phi ,{\Omega }_{Z_S}=\Omega \sin \phi ,$

$\Omega$ $–$ угловая скорость вращения Земли, ${\omega }_S={\left|\begin{array}{cccc}{\omega }_{X_S}& {\omega }_{Y_S}& {\omega }_{Z_S}& \text{tg}\phi \end{array}\right|}^T$ $–$ вектор угловой скорости ССК, обусловленной движением объекта относительно Земли, ${\omega }_{Y_S}=V_{X_S}{\left(r+H\right)}^{-1}$${\omega }_{X_S}=-V_{Y_S}{\left(r+H\right)}^{-1}$, $g_S$ $–$ вектор ускорения силы тяжести.

Для рассматриваемой ориентации осей ССК проекции вектора $g_S={\left|g_{X_S}{g}_{Y_S}{g}_{Z_S}\right|}^T$ на оси ССК определяются как:

$g_{X_S}=0$, $g_{Y_S}=-{\Omega }^2\left(r+H\right)\cos \phi \sin \phi ,$ $g_{Z_S}={\Omega }^2\left(r+H\right){\cos }^2\phi -g,$

$g$ $–$ гравитационное ускорение.

Система уравнений (3.2) в проекциях на оси ССК S, с учетом приведенных проекций векторов $g_S$ и ${\Omega }_S$, трансформируется к виду:

$\begin{array}{c}{\dot{V}}_{X_S}=A_{X_S}+\left(2\Omega \sin \phi +V_{X_S}\text{tg}\phi {\left(r+H\right)}^{-1}\right)V_{Y_S}-\left(2\Omega \cos \phi +V_{X_S}{\left(r+H\right)}^{-1}\right)V_{Z_S}\\ {\dot{V}}_{Y_S}=A_{Y_S}-\left(2\Omega \sin \phi +V_{X_S}\text{tg}\phi {\left(r+H\right)}^{-1}\right)V_{X_S}-V_{Y_S}{\left(r+H\right)}^{-1}{V}_{Z_S}-{\Omega }^2\left(r+H\right)\cos \phi \sin \phi ,\\ {\dot{V}}_{Z_S}=A_{Z_S}+\left(2\Omega \cos \phi +V_{X_S}{\left(r+H\right)}^{-1}\right)V_{X_S}+V_{Y_S}^2{\left(r+H\right)}^{-1}+{\Omega }^2\left(r+H\right){\cos }^2\phi -g\end{array}$ (3.3)

В свою очередь, вектор ускорений $A_S$, измеряемых акселерометрами, может быть представлен в ССК следующим образом:

$A_S=C^T\left(Z_{a1}-W_{a1}\right)$, (3.4)

где $Z_{a1}$ $–$ вектор выходных сигналов акселерометров первой группы, $W_{a1}$ $–$ вектор помех акселерометров первой группы, $C=D{B}^T$ $–$ матрица поворота (текущей ориентации) ПСК1 относительно ССК, D $–$ матрица поворота [21, 25] ПСК1 относительно ИСК, B $–$ матрица текущей ориентации ССК относительно ИСК (Приложение 1).

Соотношения (3.3), (3.4) позволяют записать дифференциальные уравнения, описывающие изменение проекций скорости $V_{X_S},V_{Y_S}$ и $V_{Z_S}$, в следующей векторной форме Ланжевена:

$\begin{array}{c}\left|\begin{array}{c}{\dot{V}}_{X_S}\\ {\dot{V}}_{Y_S}\\ {\dot{V}}_{Z_S}\end{array}\right|=C^T{Z}_{a1}+\left(\left(2\left|\begin{array}{c}0\\ \Omega \cos \phi \\ \Omega \sin \phi \end{array}\right|+{\left(r+H\right)}^{-1}\left|\begin{array}{c}-V_{Y_S}\\ V_{X_S}\\ V_{X_S}\text{tg}\phi \end{array}\right|\right)\times \left|\begin{array}{c}{V}_{X_S}\\ V_{Y_S}\\ V_{Z_S}\end{array}\right|\right)+\\ +\left|\begin{array}{c}0\\ -{\Omega }^2\left(r+H\right)\cos \phi \sin \phi \\ {\Omega }^2\left(r+H\right){\cos }^2\phi -g\end{array}\right|-C^T{W}_{a1},\end{array}$ (3.5)

где высота объекта, входящая в уравнения (3.5), описывается уравнением

$\dot{H}=V_{Z_S}$

В свою очередь, текущую ориентацию трехгранника ПСК1 J₁ относительно трехгранника I ИСК зададим, используя параметры Родрига$–$Гамильтона ${\mu }_1,{\mu }_2,{\mu }_3,{\mu }_4$, определяющие матрицу $D=D(\mu )$ (Приложение 2):

$\dot{\mu }=\frac{1}{2}\Phi \left(\mu \right){\omega }_{J1},{\mu }_1\left(0\right)={\mu }_{1_0},{\mu }_2\left(0\right)={\mu }_{2_0},{\mu }_3\left(0\right)={\mu }_{3_0},{\mu }_4\left(0\right)={\mu }_{4_0},$ (3.6)

где

$\mu ={\left|\begin{array}{cccc}{\mu }_1& {\mu }_2& {\mu }_3& {\mu }_4\end{array}\right|}^T,\Phi \left(\mu \right)=\left|\begin{array}{ccc}-{\mu }_2& -{\mu }_3& -{\mu }_4\\ {\mu }_1& -{\mu }_4& {\mu }_3\\ {\mu }_4& {\mu }_1& -{\mu }_2\\ -{\mu }_3& {\mu }_2& {\mu }_1\end{array}\right|,$

${\omega }_{J1}$ $–$ вектор абсолютной угловой скорости вращения первого приборного трехгранника, который может быть выражен через вектор показаний $Z_{d1}$ трех ДУС первой группы:

${\omega }_{J1}=Z_{d1}-W_{d1}$, (3.7)

где $W_{d1}$$–$ вектор помех измерения ДУС первой группы.

С учетом (3.7) угловое движение объекта (3.6) относительно ИСК может быть представлено следующим образом:

$\dot{\mu }=\frac{1}{2}\Phi \left(\mu \right)\left(Z_{d1}-W_{d1}\right)$

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

$\dot{\mu }=\frac{1}{2}\Phi \left(\mu \right)\left(Z_{d1}-W_{d1}\right)$ (3.8)

$\left|\begin{array}{c}\dot{\lambda }\\ \dot{\phi }\end{array}\right|=\left|\begin{array}{cc}0& {\left(\cos \phi \right)}^{-1}\\ 1& 0\end{array}\right|\left|\begin{array}{c}{V}_{Y_S}\\ V_{X_S}\end{array}\right|{\left(r+H\right)}^{-1}$

$\begin{array}{c}\left|\begin{array}{c}{\dot{V}}_{X_S}\\ {\dot{V}}_{Y_S}\\ {\dot{V}}_{Z_S}\end{array}\right|=C^T{Z}_{a1}+\left(\left(2\left|\begin{array}{c}0\\ \Omega \cos \phi \\ \Omega \sin \phi \end{array}\right|+{\left(r+H\right)}^{-1}\left|\begin{array}{c}-V_{Y_S}\\ V_{X_S}\\ V_{X_S}\text{tg}\phi \end{array}\right|\right)\times \left|\begin{array}{c}{V}_{X_S}\\ V_{Y_S}\\ V_{Z_S}\end{array}\right|\right)+\\ +\left|\begin{array}{c}0\\ -{\Omega }^2\left(r+H\right)\cos \phi \sin \phi \\ {\Omega }^2\left(r+H\right){\cos }^2\phi -g\end{array}\right|-C^T{W}_{a1},\end{array}$

$\dot{H}=V_{Z_S},$

где $C=C\left(\mu ,\lambda ,\phi ,t\right)=D\left(\mu \right)B^T\left(\lambda ,\phi ,t\right)$.

Для дальнейшего построения уравнений текущей ориентации антенны рассмотрим предварительно уравнения ее углового движения относительно мачты $–$ разворота ПСК2 J2 относительно ОСК Q, описываемого вектором параметров Родрига$–$Гамильтона $\eta ={\left|\begin{array}{cccc}{\eta }_1& {\eta }_2& {\eta }_3& {\eta }_4\end{array}\right|}^T$ [23, 25]:

$\dot{\eta }=\frac{1}{2}\Phi \left(\eta \right){\omega }_{JQ};\eta \left(0\right)={\eta }_0,$ (3.9)

где $\eta \left(0\right)={\eta }_0$ $–$ неизвестные параметры начального рассогласования ПСК2 J₂ и ОСК Q; $w_{JQ}$$–$ вектор случайной угловой скорости ПСК2 J₂ (второго триэдра БИСО) относительно ОСК Q , аппроксимируемый центрированным БГШ с матрицей интенсивностей D_J.

Аналогично рассмотрим динамику изменения ориентации мачты относительно объекта $–$ разворота ОСК Q относительно ПСК1 J₁, описываемую вектором параметров Родрига$–$Гамильтона $\lambda ={\left|\begin{array}{cccc}{\lambda }_1& {\lambda }_2& {\lambda }_3& {\lambda }_4\end{array}\right|}^T$:

$\dot{\lambda }=\frac{1}{2}\Phi \left(\lambda \right){\omega }_Q;\lambda \left(0\right)={\lambda }_0,$ (3.10)

где ${\omega }_Q={\left|\begin{array}{ccc}{\omega }_X& {\omega }_Y& 0\end{array}\right|}^T$ $–$ вектор угловой скорости вращения ОСК Q относительно ПСК1 J₁.

Для описания вектора ω_Q угловой скорости движения ОСК Q относительно ПСК1 J1 воспользуемся выражением для вектора ускорения, возникающего при движении материальной точки (МТ) по сфере радиуса R, записанным в ОСК Q:

$A_Q={\dot{V}}_Q+\left(2{\Omega }_Q+{\omega }_Q\right)\times V_Q-g_Q,$ (3.11)

где $A_Q={\left|\begin{array}{ccc}{A}_X& A_Y& A_Z\end{array}\right|}^T$ $–$ вектор ускорений МТ в ОСК Q, $V_Q={\left|\begin{array}{ccc}{V}_X& V_Y& V_Z\end{array}\right|}^T$ $–$ вектор скорости МТ в ОСК; ${\Omega }_Q={\left|\begin{array}{ccc}{\Omega }_X& {\Omega }_Y& {\Omega }_Z\end{array}\right|}^T$ − вектор угловой скорости вращения Земли в ОСК; $g_Q={\left|\begin{array}{ccc}{g}_X& g_Y& g_Z\end{array}\right|}^T$ $–$ вектор ускорения силы тяжести в ОСК.

В соответствии с приведенными выше выражениями проекций вектора угловой скорости вращения Земли и вектора ускорения силы тяжести на оси ССК S проекции векторов g_Q, Ω_Q на оси ОСК Q определяются как

$\begin{array}{c}{\Omega }_Q\left(\mu ,\lambda ,\lambda ,\phi ,t\right)={\left|\begin{array}{ccc}{\Omega }_X& {\Omega }_Y& {\Omega }_Z\end{array}\right|}^T=D\left(\lambda \right)C\left(\mu ,\lambda ,\phi ,t\right)\left|\begin{array}{c}0\\ \Omega \cos \phi \\ \Omega \sin \phi \end{array}\right|\\ g_Q\left(\mu ,\lambda ,\lambda ,\phi ,t\right)={\left|\begin{array}{ccc}{g}_X& g_Y& g_Z\end{array}\right|}^T=D\left(\lambda \right)C\left(\mu ,\lambda ,\phi ,t\right)\left|\begin{array}{c}0\\ -\Omega \left(r+H\right)\cos \phi \sin \phi \\ {\Omega }^2\left(r+H\right)\cos \phi -g\end{array}\right|,\end{array}$

где $D\left(\lambda \right)$ $–$ матрица поворота (направляющих косинусов) ОСК Q относительно ПСК1 J₁.

Система уравнений (3.11) в проекциях на оси ОСК с учетом очевидных равенств

 ${\omega }_X=-V_Y{R}^{-1},{\omega }_Y=V_X{R}^{-1},$

а также приведенных выше проекций векторов g_Q, Ω_Q, трансформируется к виду:

$\begin{array}{c}{A}_X={\dot{\omega }}_{Y}R-2{\Omega }_Z\left(\mu ,\lambda ,\lambda ,\phi ,t\right)R{\omega }_X-g_X\left(\mu ,\lambda ,\lambda ,\phi ,t\right)\\ A_Y={\dot{\omega }}_{X}R-2{\Omega }_Z\left(\mu ,\lambda ,\lambda ,\phi ,t\right){\omega }_{X}R-g_Y\left(\mu ,\lambda ,\lambda ,\phi ,t\right)\\ A_Z\left(2{\Omega }_Y\left(\mu ,\lambda ,\lambda ,\phi ,t\right)+{\omega }_Y\right)R{\omega }_Y+\left(2{\Omega }_X\left(\mu ,\lambda ,\lambda ,\phi ,t\right)+{\omega }_X\right)R{\omega }_X-g_Z\left(\mu ,\lambda ,\lambda ,\phi ,t\right)\end{array}$ (3.12)

Вторая группа акселерометров измеряет как проекции вектора ускорений A_Q на оси ПСК2: $D\left(\eta \right)A_Q$, так и проекции вектора ускорений ЦМ объекта, измеряемые первой группой акселерометров в ПСК1 (т.е. $Z_{a1}-W_{a1}$ ), на оси ПСК2: $D\left(\eta \right)D\left(\lambda \right)Z_{d1}-W_{d1}$. Таким образом, с учетом принятой модели измерений акселерометров, справедливо равенство:

$D\left(\eta \right)A_Q+D\left(\eta \right)D\left(\lambda \right)\left(Z_{a1}-W_{a1}\right)+W_{a2}=Z_{a2}$,

где $Z_{a2}$ $–$ вектор выходных сигналов акселерометров второй группы, $W_{a2}$ $–$ вектор помех акселерометров второй группы, откуда имеем выражение вектора ускорений A_Q:

$A_Q=D^T\left(\eta \right)\left(Z_{a2}-W_{a2}\right)-D\left(\lambda \right)\left(Z_{a1}-W_{a1}\right)$

Данное соотношение совместно с уравнениями (3.12) позволяет сформировать стохастические уравнения, описывающие вектор ω_Q угловой скорости движения ОСК относительно ПСК1:

$\begin{array}{c}{\dot{\omega }}_Q=-R^{-1}\left|\begin{array}{ccc}0& 1& 0\\ 1& 0& 0\\ 0& 0& 0\end{array}\right|\left(D^T\left(\eta \right)Z_{a2}-D\left(\lambda \right)Z_{a1}\right)+2{\Omega }_Z\left(\mu ,\lambda ,\lambda ,\phi ,t\right)\left|\begin{array}{ccc}0& -1& 0\\ 1& 0& 0\\ 0& 0& 0\end{array}\right|{\omega }_Q+\\ +R^{-1}\left|\begin{array}{c}-g_Y\left(\mu ,\lambda ,\lambda ,\phi ,t\right)\\ -g_X\left(\mu ,\lambda ,\lambda ,\phi ,t\right)\\ 0\end{array}\right|-R^{-1}\left|\begin{array}{ccc}0& 1& 0\\ 1& 0& 0\\ 0& 0& 0\end{array}\right|\left|D\left(\lambda \right)⋮-D^T\left(\eta \right)\right|\left|\begin{array}{c}{W}_{a1}\\ W_{a2}\end{array}\right|\end{array}$ (3.13)

Здесь важно отметить, что полученная система уравнений, описывающая текущую ориентацию ОСК в ПСК1, в соответствии с поставленной задачей оказывается полностью инвариантна к характеру динамики движения основания (мачты), качающегося относительно начала ПСК1.

Объединяя системы уравнений (3.8)$–$(3.10), (3.13), стохастические уравнения полного вектора параметров текущей ориентации БИСО на подвижном основании запишем следующим образом:

$\dot{Y}=F\left(Y,t\right)+F_1\left(Y,t\right)W,$ (3.14)

где

$Y={\left|\begin{array}{cccccccc}{\mu }^T& \lambda & \phi & V_S& H& {\eta }^T& {\lambda }^T& {\omega }_Q{}^T\end{array}\right|}^T,W={\left|\begin{array}{cccc}{W}_{d1}^T& W_{a1}^T& W_{a1}^T& W_{J\text{Q}}^T\end{array}\right|}^T$

$F\left(Y,t\right)=\left|\begin{array}{c}\frac{1}{2}\Phi \left(\mu \right)Z_{d1}\\ \left|\begin{array}{cc}0& \left(\cos {\phi }^{-1}\right)\\ 1& 0\end{array}\right|\left|\begin{array}{c}{V}_{Y_S}\\ V_{X_S}\end{array}\right|{\left(r+H\right)}^{-1}\\ C^T{Z}_{a1}+\left(\left(2\left|\begin{array}{c}0\\ \Omega \cos \phi \\ \Omega \sin \phi \end{array}\right|+{\left(r+H\right)}^{-1}\left|\begin{array}{c}-V_{Y_S}\\ V_{X_S}\\ V_{X_S}\text{tg}\phi \end{array}\right|\right)\times V_S\right)+\left|\begin{array}{c}0\\ -{\Omega }^2\left(r+H\right)\cos \phi \sin \phi \\ {\Omega }^2\left(r+H\right)\cos \phi -g\end{array}\right|\\ V_{Z_S}\\ 0\\ \dot{\lambda }=\frac{1}{2}\Phi \left(\lambda \right){\omega }_Q\\ -R\left|\begin{array}{ccc}0& 1& 0\\ 1& 0& 0\\ 0& 0& 0\end{array}\right|\left(D^T\left(\eta \right)Z_{a2}-D\left(\lambda \right)Z_{a1}\right)+2\Omega \left(\mu ,\lambda ,\lambda ,\phi ,t\right)\left|\begin{array}{ccc}0& 1& 0\\ 1& 0& 0\\ 0& 0& 0\end{array}\right|{\omega }_Q+\\ +R^{-1}\left|\begin{array}{c}-g_Y\left(\mu ,\lambda ,\lambda ,\phi ,t\right)\\ g_X\left(\mu ,\lambda ,\lambda ,\phi ,t\right)\\ 0\end{array}\right|\end{array}\right|$

где 0 $–$ нулевая матрица соответствующей размерности.

Для стохастической оценки состояния нелинейных динамических систем вида (3.9) наиболее эффективным подходом в настоящее время является использование методов теории стохастической фильтрации [24, 26], из которых самым широко известным и общеупотребительным является нелинейный (расширенный или обобщенный) фильтр Калмана. Но для его применения необходим предварительный синтез уравнения наблюдателя компонентов вектора Y (т.е. информационной модели сигнала измерения, явно зависящей от составляющих вектора Y).

В рассматриваемом случае в качестве сигналов наблюдения вектора Y можно выбрать выходные сигналы трех ортогональных ДУС второй группы. Действительно, вектор ${\omega }_{J2}$ абсолютной угловой скорости вращения ПСК2 J₂, измеряемый ДУС второй группы, определяется суммой вектора ω_JQ, проекций векторов ${\omega }_Q$, ${\Omega }_Q$ в ПСК2 $J_2\left({\omega }_{\text{Q}J2},{\Omega }_{\text{Q}J2}\right)$ и проекции вектора абсолютной угловой скорости объекта ${\omega }_{J1}$, измеряемой первой группой ДУС в ПСК1, на оси ПСК2 $\left({\omega }_{J12}\right)$:

${\omega }_{J2}={\omega }_{JQ}+{\omega }_{\text{Q}J2}+{\Omega }_{\text{Q}J2}+{\omega }_{J12}$,

где ${\omega }_{QJ2}=D\left(\eta \right){\omega }_Q,{\Omega }_{QJ2}=D\left(\eta \right)D\left(\lambda \right)C\left(\mu ,\lambda ,\phi ,t\right)\left|\begin{array}{c}0\\ \Omega \cos \phi \\ \Omega \sin \phi \end{array}\right|$

${\omega }_{J12}=D\left(\eta \right)D\left(\lambda \right)\left(Z_{d1}+W_{d1}\right)$

Это позволяет, исходя из приведенного выражения для ${\omega }_{J2}$ и уравнения вектора выходных сигналов ДУС второй группы, представить стохастическую модель вектора наблюдения следующим образом:

$\begin{array}{c}{Z}_{d2}={\omega }_{JQ}+D\left(\eta \right)D\left(\lambda \right)C\left(\mu ,\lambda ,\phi ,t\right)\left|\begin{array}{c}0\\ \Omega \cos \phi \\ \Omega \sin \phi \end{array}\right|+\\ +D\left(\eta \right)D\left(\lambda \right)\left(Z_{d1}-W_{d1}\right)+W_{d2}=H_d\left(Y,t\right)+H_{d1}\left(Y\right)\left|\begin{array}{c}{\omega }_{JQ}\\ W_{d1}\\ W_{d2}\end{array}\right|,\end{array}$

где

$\begin{array}{c}{H}_d\left(Y,t\right)=D\left(\eta \right){\omega }_Q+D\left(\eta \right)D\left(\lambda \right)C\left(\mu ,\lambda ,\phi ,t\right)\left|\begin{array}{c}0\\ \Omega \cos \phi \\ \Omega \sin \phi \end{array}\right|+D\left(\eta \right)D\left(\lambda \right)Z_{d1}\\ H_{d1}\left(Y\right)=\left|E_3⋮-D\left(\eta \right)D\left(\lambda \right)⋮E_3\right|,\end{array}$ (3.15)

Е₃ $–$ единичная матрица размерности 3.

Несмотря на то, что данный наблюдатель позволяет в явном виде наблюдать подавляющее большинство компонентов вектора состояния Y, его особенностью является невозможность наблюдения вектора линейной скорости объекта V_S, что существенно влияет на сходимость и устойчивость процесса оценивания всего вектора состояния. Для формирования сигнала наблюдения вектора скорости объекта V_S можно привлечь или измерения доплеровского датчика скорости, или доплеровские измерения спутника. Рассмотрим второй вариант как более технологичный и дешевый, полагая при этом частоту поступления спутниковых измерений высокой (в настоящее время $–$ до 100 Гц), позволяющей считать характер спутниковых измерений по отношению к динамике изменения навигационных параметров рассматриваемого объекта непрерывным [1, 4]. В этом случае информационный сигнал доплеровских измерений (псевдоскорости) Z_V в ССК S может быть представлен, как показано в [20], следующим образом:

$\begin{array}{c}{Z}_V=\left[\left({\xi }_c-\left(r+H\right)\cos \phi \sin \lambda \right)\right.\left(V_{\xi c}-B_{\left(1\right)}^T\left(\phi ,\lambda \right)V_S\right)+\left({\eta }_c-\left(r+H\right)\sin \phi \right)\left(V_{\eta c}\right.-\\ \left.\left.-B_{\left(2\right)}^T\left(\phi ,\lambda \right)V_S\right)+\left({\varsigma }_c-\left(r+H\right)\cos \phi \cos \lambda \right)\left(V_{\varsigma c}-B_{\left(3\right)}^T\left(\phi ,\lambda \right)V_S\right)\right]\times \\ \times {\sqrt{{\left({\xi }_c-\left(r+H\right)\cos \phi \sin \lambda \right)}^2+{\left({\eta }_c-\left(r+H\right)\sin \phi \right)}^2+{\left({\zeta }_c-\left(r+H\right)\cos \phi \cos \lambda \right)}^2}}^{-1}+\\ +W_{Z_V}=H_V\left(Y,t\right)+W_V,\end{array}$ (3.16)

где ${\xi }_c$, ${\eta }_c$, ${\varsigma }_c$ $–$ известные декартовы координаты спутника в ГСК, $V_{\xi c},V_{\eta c},V_{\varsigma c}$$–$ известные проекции вектора скорости спутника на оси ГСК, WV $–$ БГШ с нулевым средним и известной интенсивностью DV, обусловленный погрешностями измерения, $B_{\left(i\right)}^T\left(\phi ,\lambda \right)$ $–$ i-я строка матрицы $B^T\left(\phi ,\lambda \right)\left(B\left(\phi ,\lambda \right)=B{\left.\left(\lambda ,\phi ,t\right)\right|}_{t=0}\right).$

Функция наблюдения $H_V\left(Y,t\right)$ в (2.16) явно зависит от всех параметров линейного движения объекта, в т.ч. и от вектора скорости объекта V_S, что позволяет в совокупности с наблюдателем (3.15) сформировать наблюдатель всех компонентов вектора Y. Т.к. спутниковый навигационный приемник кроме доплеровских измерений принимает еще и кодовые измерения, то для повышения информативности наблюдения их также целесообразно включить в состав комплексного наблюдателя. Учитывая, что модель информационного сигнала кодовых измерений имеет вид [13, 20]:

$\begin{array}{c}{Z}_R=\sqrt{{\left({\xi }_c-\left(r+H\right)\cos \phi \sin \lambda \right)}^2+\left({\eta }_c-\left(r+H\right)\sin \phi \right)+{\left({\varsigma }_c-\left(r+H\right)\cos \phi \cos \lambda \right)}^2}+\\ +W_{Z_R}=H_R\left(Y,t\right)+W_{Z_R},\end{array}$

где W_R $–$ центрированный БГШ с известной интенсивностью D_R, обусловленный алгоритмически нескомпенсированными ошибками часов спутников и приемника, задержками сигнала при прохождении ионосферы и тропосферы, ошибками многолучевости и другими погрешностями, уравнения комплексного наблюдателя можно записать следующим образом:

$\begin{array}{c}Z=\left|\begin{array}{c}{Z}_{d2}\\ Z_V\\ Z_R\end{array}\right|=\left|\begin{array}{c}{H}_d\left(Y,t\right)\\ H_V\left(Y,t\right)\\ H_R\left(Y,t\right)\end{array}\right|+\left|\begin{array}{ccc}{H}_{d1}\left(Y\right)& ⋮& 0\\ 0& ⋮& E_2\end{array}\right|\left|\begin{array}{c}{\omega }_{JQ}\\ W_{d1}\\ W_{d2}\\ W_V\\ W_R\end{array}\right|=H\left(Y,t\right)+H_1\left(Y\right)\left|\begin{array}{c}{\omega }_{JQ}\\ W_{d1}\\ W_{d2}\\ W_V\\ W_R\end{array}\right|\\ H\left(Y,t\right)=\left|\begin{array}{c}{H}_d\left(Y,t\right)\\ H_V\left(Y,t\right)\\ H_R\left(Y,t\right)\end{array}\right|,H_1\left(Y\right)=\left|\begin{array}{ccc}{H}_{d1}\left(Y\right)& ⋮& 0\\ 0& ⋮& E_2\end{array}\right|,\end{array}$ (3.17)

где Е₂ $–$ единичная матрица размерности 2.

Уравнения (3.14), (3.17) «объект$–$наблюдатель» позволяют построить оценку вектора состояния Y в виде расширенного фильтра Калмана $–$ наиболее эффективного на сегодняшний день алгоритма оценивания для динамических нелинейных стохастических систем. Особенностью здесь является наличие корреляции шумов объекта (3.14) и наблюдателя (3.17), которую необходимо далее учитывать при последующем построении фильтра.

Расширенный фильтр Калмана, построенный по уравнениям «объект$–$наблюдатель» (3.14), (3.17) и обеспечивающий принципиальное решение поставленной задачи, имеет следующий вид [21, 24]:

$\begin{array}{c}\widehat{\dot{Y}}=F\left(\widehat{Y},t\right)+K\left(\widehat{Y},t\right)\left[Z-H\left(\widehat{Y},t\right)\right]\\ K\left(\widehat{Y},t\right)=\left(R\frac{\partial H^T\left(\widehat{Y}\right)}{\partial \widehat{Y}}+\Theta \right){\left(H_1\left(\widehat{Y}\right)D_H{H}_1^T\left(\widehat{Y}\right)\right)}^{-1}\\ \dot{R}\left(\widehat{Y},t\right)=\left\{\frac{\partial F^T\left(\widehat{Y},t\right)}{\partial \widehat{Y}}-\Theta {\left(H_1\left(\widehat{Y}\right)D_H{H}_1^T\left(\widehat{Y}\right)\right)}^{-1}\frac{\partial H^T\left(\widehat{Y}\right)}{\partial \widehat{Y}}\right\}R\left(\widehat{Y},t\right)+\\ +R\left(\widehat{Y},t\right)\left\{\frac{\partial F^T\left(\widehat{Y},t\right)}{\partial \widehat{Y}}-\Theta {\left(H_1\left(\widehat{Y}\right)D_H{H}_1^T\left(\widehat{Y}\right)\right)}^{-1}\frac{\partial H^T\left(\widehat{Y}\right)}{\partial \widehat{Y}}\right\}+F_1\left(\widehat{Y},t\right)D_0{F}_1^T\left(\widehat{Y},t\right)-\\ -R\frac{\partial H^T\left(\widehat{Y}\right)}{\partial \widehat{Y}}{\left(H_1\left(\widehat{Y}\right)D_H{H}_1^T\left(\widehat{Y}\right)\right)}^{-1}\frac{\partial H^T\left(\widehat{Y}\right)}{\partial \widehat{Y}}R-\Theta {\left(H_1\left(\widehat{Y}\right)D_H{H}_1^T\left(\widehat{Y}\right)\right)}^{-1}{\Theta }^T,\end{array}$ (3.18)

где $\widehat{Y}$ $–$ вектор текущей оценки вектора состояния Y(t); $R\left(\widehat{Y},t\right)$ $–$ апостериорная ковариационная матрица;

${\widehat{Y}}_0=M\left(Y_0\right);R_0=M\left\{\left(Y_0-{\widehat{Y}}_0\right){\left(Y_0-{\widehat{Y}}_0\right)}^T\right\}$

$D_0=\left|\begin{array}{l}{D}_{d1}0\\ 0D_{a1}0\\ 0D_{a2}0\\ 0D_J\end{array}\right|,D_H=\left|\begin{array}{l}{D}_{J}0\\ 0D_{d1}0\\ 0D_{d2}0\\ 0D_{V}0\\ 0D_R\end{array}\right|,\Theta =F_1\left(\widehat{Y},t\right)\left|\begin{array}{l}0D_{d1}0\\ 0\\ D_{J}0\end{array}\right|H_1{}^T\left(\widehat{Y}\right)$

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

4. Результаты имитационного моделирования. Для анализа устойчивости и сходимости процесса оценивания вектора Y(t) с использованием фильтра (3.18) было выполнено численное моделирование процесса оценки ориентации антенны на высокодинамичном подвижном основании. Моделирование осуществлялось на временно́м интервале [0, 1000] секунд с использованием при интегрировании уравнений оценки метода Рунге$–$Кутты 4-го порядка с шагом 0.01 с.

Движение объекта задавалось вдоль местного меридиана из точки с долготой 30°, широтой 45° и высотой 3 м. Проекции скорости объекта на оси CСК задавались в функции времени как:

$V_X=\mathrm{0,}{V}_Y=100+\sin 0.02t,V_Z=1.5\cos t,$

проекции вектора угловой скорости ${\omega }_{J1}$, соответственно:

${\omega }_{x1}={10}^{-1}\cos 3t,{\omega }_{y1}={10}^{-3}\cos 4t,{\omega }_{z1}=1.5\cdot {10}^{-2}\cos 7t$

Истинные текущие координаты объекта формировались путем интегрирования данных проекций скорости в соответствии с уравнениями (3.1), (3.6).

Угловые скорости, определяющие динамику движения антенны относительно мачты, были заданы центрированными случайными гауссовскими последовательностями с с.к.о. 10^-2 рад/с, а динамику движения мачты относительно объекта, соответственно, как:

${\omega }_X=0.43\cos \left(2.45t\right)\mathrm{рад}\text{/}с\text{,}{\omega }_Y=0.3\cos \left(2.7t\right)\mathrm{рад}\text{/}с$

Определение истинного углового положения антенны осуществлялось путем интегрирования уравнений (3.9), (3.10) с учетом выбранных моделей углового движения.

Компоненты векторов помех измерения акселерометров W_ai и ДУС W_di моделировались центрированными случайными гауссовскими последовательностями с соответствующими с.к.о., приведенными ниже. Формирование реальных показаний ДУС осуществлялось путем наложения данных случайных гауссовских последовательностей на соответствующие векторы угловых скоростей.

На рис. 2$–$5 представлены графики изменения погрешностей оценивания параметров текущей ориентации антенны и самого объекта относительно рассмотренных выше СК. Обобщенный анализ результатов моделирования показал:

1. Предложенный алгоритм оценки ориентации обеспечивает быструю сходимость и высокую устойчивость процесса оценивания. Погрешности оценивания углового положения антенны по окончании переходного процесса по всем углам ориентации не превысили 10^-5 рад, объекта, соответственно, $–$ 1.5·10^-6 рад, что соответствует не только современным, но и перспективным, требованиям к системам ориентации.
2. Увеличение погрешностей начальной оценки параметров η, λ, µ приводит к увеличению времени определения пространственной ориентации. Вариации частоты и амплитуды колебаний мачты в пределах 50$–$80% приводят к незначительному росту погрешностей оценивания (не более 6%), т.е. их влияние на точность ориентации оказывается существенно меньше, нежели погрешностей начальной оценки и с.к.о. шумов измерений.

На рис. 2 приведены погрешности оценивания ориентации ПСК1 относительно ИСК, выраженные вектором параметров Родрига$–$Гамильтона µ. Погрешности их начальной оценки $–$ при пересчете в углы Эйлера, были заданы в пределах 0.01 рад, с.к.о. помех измерения ДУС $–$ 5·10^-8 рад/с, с.к.о. помех измерения акселерометров $–$ 10^-4 м/с2.

 

Рис. 2. Погрешности оценивания ориентации ПСК1 относительно ИСК, выраженные вектором параметров Родрига$–$Гамильтона µ.

 

При данном уровне помех погрешности оценки компонентов вектора µ с 50-й секунды не превышали 10^-4, что при расчете ориентации в углах Эйлера определяет ошибку оценивания в интервалах: для курсового угла: $–$2.7·10^-4 $–$ 2.8·10^-4 рад, для угла тангажа: $–$2.5·10^-4 $–$ 2.5·10^-4 рад, для угла крена: $–$2.5·10^-5 $–$ 2.5·10^-5 рад.

На рис. 3 приведены погрешности оценивания ориентации ПСК2 относительно ОСК, выраженные вектором параметров Родрига$–$Гамильтона η. Как видно из рис. 3, погрешности оценки компонентов вектора η после 200-й секунды не превышают 2.5·10^-5 , что в углах Эйлера определяет ошибку оценивания в интервалах: для курсового угла: $–$2.1·10^-5 $–$ 2.6·10^-5 рад, для угла тангажа: $–$2.2·10^-5 $–$ 2.4·10^-5 рад, для угла крена: $–$2.3·10^-5 $–$ 2.4·10^-5 рад.

 

Рис. 3. Погрешности оценивания ориентации ПСК2 относительно ОСК, выраженные вектором параметров Родрига$–$Гамильтона η.

 

На рис. 4 приведены погрешности оценивания ориентации ОСК относительно ПСК1, выраженные вектором параметров Родрига$–$Гамильтона λ. В данном случае погрешности оценки вектора λ после 50-й секунды не превышают величины 2·10^-5, что в углах Эйлера определяет ошибку оценивания в интервалах: для курсового угла: $–$1.1·10^-5 $–$ 1.6·10^-5 рад, для угла тангажа: $–$1.3·10^-5 $–$ 1.7·10^-5 рад, для угла крена: $–$1.4·10^-5 $–$ 1.8·10^-5 рад.

 

Рис. 4. Погрешности оценивания ориентации ОСК относительно ПСК1, выраженные вектором параметров Родрига$–$Гамильтона λ.

 

Очевидно, что в целом полученные результаты удовлетворяют требованиям по точности определения ориентации антенн, предъявляемым не только к современным, но и перспективным БИСО. На рис. 5 показаны ошибки оценки координат ЦМ подвижного РТК. Как видно из полученных результатов, ошибки оценки долготы λ не превысили 4·10^-7 рад, а ошибки определения широты φ $–$ 1.5·10^-6 рад, что вполне соответствует требованиям, предъявляемым к современным системам навигации подвижных РТК.

 

Рис. 5. Ошибки оценки координат подвижного РТК

 

Заключение. В целом, результаты имитационного моделирования позволяют сделать вывод о том, что устойчивость предложенного алгоритма и его высокая точность обеспечивают возможность его использования для решения задачи оперативной ориентации мачтовых антенн подвижных РТК с использованием средне- и высокоточных БИСО без коррекции в течение достаточно длительного времени.

Приложение 1

$B=B(\lambda ,\phi ,t)=\left|\begin{array}{ccc}\text{cos(}\lambda \text{+}\Omega \text{t)}& 0& \text{-sin(}\lambda \text{+}\Omega \text{t)}\\ \text{-sin(}\lambda \text{+}\Omega \text{t)sin}\phi & \cos \phi & \text{-cos(}\lambda \text{+}\Omega \text{t)sin}\phi \\ \sin \text{(}\lambda \text{+}\Omega \text{t)}\cos \phi & \sin \phi & \text{cos(}\lambda \text{+}\Omega \text{t)}\cos \phi \end{array}\right|$

Приложение 2

$D\left(\mu \right)=\left|\begin{array}{ccc}2{\mu }_1^2+2{\mu }_2^2-1& 2\left({\mu }_2{\mu }_3+{\mu }_1{\mu }_4\right)& 2\left({\mu }_2{\mu }_4-{\mu }_1{\mu }_3\right)\\ 2\left({\mu }_2{\mu }_3-{\mu }_1{\mu }_4\right)& 2{\mu }_1^2+2{\mu }_3^2-1& 2\left({\mu }_3{\mu }_4+{\mu }_1{\mu }_2\right)\\ 2\left({\mu }_2{\mu }_4+{\mu }_1{\mu }_3\right)& 2\left({\mu }_3{\mu }_4-{\mu }_1{\mu }_2\right)& 2{\mu }_1^2+2{\mu }_4^2-1\end{array}\right|$

About the authors

S. V. Sokolov

Moscow Technical University of Communications and Informatics

Author for correspondence.
Email: s.v.s.888@yandex.ru
Russian Federation, Moscow

V. A. Pogorelov

Moscow Technical University of Communications and Informatics

Email: vadim-pva@narod.ru
Russian Federation, Moscow

I. V. Reshetnikova

Moscow Technical University of Communications and Informatics

Email: irina_reshetnikova@mail.ru
Russian Federation, Moscow

References

Supplementary files