Coordinates estimation of the unmanned aerial vehicles by using distributed system of base stations

Full Text

Введение

В последнее время беспилотные летательные аппараты (БПЛА) нашли применение в гражданских задачах, таких как слежение, управление движением колонны автомобилей, фотосъемка, доставка посылок и т. д. [1]. На данный момент перспективы применения систем достаточно широки.

Для управления БПЛА требуется высокая точность определения местоположения (позиционирования), для достижения которой используют системы, основанные на комбинации глобальной навигационной спутниковой системы (GNSS) и дополнительных датчиков, таких как радары, камеры, инерциальные датчики [1]. Однако, в случаях кратковременного пропадания сигнала GNSS, БПЛА требуют дополнительной информации о собственном местоположении. Для повышения точности и робастности позиционирования предлагается применять сигналы наземных базовых станций (БС) [1].

Известно четыре метода для определения координат БПЛА на основе сигналов, принимаемых от БС [2]:

• метод, использующий время прибытия сигналов (англ. time of arrival, TOA);
• метод, использующий разницу времени прибытия сигналов (англ. time difference of arrival, TDoA);
• метод, использующий углы прибытия сигналов (англ. angle of arrival, AoA);
• метод, использующий мощности принятых сигналов (англ. receive strength signal, RSS).

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

Методы, основанные на измерении времени прибытия сигнала, более простые и универсальные, вследствие чего они часто используются для позиционирования в современных БПЛА [2].

1. Алгоритм оценки координат методом TOA

Метод определения координат TOA основан на измерении расстояния от нескольких передатчиков до приемника. В качестве передатчиков выступает БС, в качестве приемника – БПЛА. Расстояние вычисляется путем оценки разности между временем начала передачи сигнала БС $T_0$ и временем приема сигнала $T_i$ в БПЛА:

$d_i=с(T_i-T_0),$ (1)

где $d_i$ – расстояние от БПЛА до i – базовой станции; с – скорость света.

Вектор оценок дальностей $d_i$ определяет точки равного расстояния от БПЛА до $Б{С}_i.$ В двумерном случае каждая такая оценка дальности эквивалентна окружности с радиусом $d_i$ и центром в месте расположения $Б{С}_i.$ Координаты БПЛА вычисляются как точка пересечения нескольких окружностей.

На рис. 1 изображена геометрия измерения координат методом TOA.

 

Рис. 1. Геометрия измерения координат методом TоA

Fig. 1. Geometry of coordinate measurement by ToA method

 

Используя вектор измерений дальностей от БС до БПЛА, можно записать систему из N нелинейных уравнений:

$\sqrt{{(х-x_1)}^2 + {(y-y_1)}^2}=d_1+{\eta }_{\mathrm{1,}\text{prop}},$ (2)

$\begin{array}{l}\sqrt{{(х-x_2)}^2 + {(y-y_2)}^2}=d_2+{\eta }_{\mathrm{2,}\text{prop}},\\ \dots ,\\ \sqrt{{(х-x_i)}^2 + {(y-y_i)}^2}=d_i+{\eta }_{i,\text{prop}},\end{array}$

где $(x_i, y_i)$ – координаты i-й БС; ${\eta }_{i, \text{prop}}$ – независимая гауссова случайная величина с нулевым математическим ожиданием и дисперсией ${\left\{{\sigma }_{\text{prop}}\right\}}^2.$

Одним из способов решения системы уравнений (2) является преобразование ее в линейную форму с помощью алгоритма TOA [4].

Для линеаризации системы (2) возведем каждое ее уравнение в квадрат:

${(x-x_i)}^2+{(y-y_i)}^2=$ (3)

$=d_i^2+2{\eta }_{i, \text{prop}}{d}_{i }+{\eta }_{i, \text{prop}}^2,i=\mathrm{1,} \dots ,N.$

Раскроем скобки в выражении (3). Тогда уравнение (3) преобразуется к виду

$x^2-2x{x}_i+x_i^2+y^2–2y{y}_i+y_i^2=$ (4)

$=d_i^2+2{\eta }_{i, \text{prop}}{d}_{i }+{\eta }_{i, \text{prop}}^2,i=\mathrm{1,} \dots ,N.$

Введем замену переменных $R^2=x^2+y^2.$ Тогда выражение (4) можно представить как (5).

$-2x{x}_i–2y{y}_i+R^2=$ (5)

$=d_i^2-x_i^2-y_i^2+2{\eta }_{i, \text{prop}}{d}_{i }+{\eta }_{i, \text{prop}}^2,i=\mathrm{1,} \dots ,N.$

Используя результат преобразований (5), система уравнений (2) будет:

$\left[\begin{array}{ccc}{x}_1& y_1& -\mathrm{0,5}\\ x_2& y_2& -\mathrm{0,5}\\ ⋮& ⋮& ⋮\\ x_N& y_N& -\mathrm{0,5}\end{array}\right] \left[\begin{array}{c}x\\ y\\ R\end{array}\right] =$ (6)

$=-\frac{1}{2} \left[\begin{array}{c}{d}_1^2-x_1^2-y_1^2\\ d_2^2-x_2^2-y_2^2\\ ⋮\\ d_N^2-x_N^2-y_N^2\end{array}\right]-\frac{1}{2}\left[\begin{array}{c}2{\eta }_{\mathrm{1,} \text{prop}}{d}_{1 }+{\eta }_{\mathrm{1,} \text{prop}}^2\\ 2{\eta }_{\mathrm{2,} \text{prop}}{d}_{2 }+{\eta }_{\mathrm{2,} \text{prop}}^2\\ ⋮\\ 2{\eta }_{N, \text{prop}}{d}_{N }+{\eta }_{N, \text{prop}}^2\end{array}\right].$

Другим методом приведения системы уравнений (2) к линейному виду является модифицированный вариант алгоритма TOA [2]. В статье предложено после раскрытия скобок в выражении (3) вычесть первое уравнение из остальных уравнений и после этого перегруппировать переменные:

$-2x(x_i-x_1)-2y( y_i-y_1)+x_i^2-x_1^2+y_i^2-y_1^2=$ (7)

$={d_i}^2-{d_1}^2+2{\eta }_{i, \text{prop}}{d}_{i }+{\eta }_{i, \text{prop}}^2-$

$-2{\eta }_{\mathrm{1,} \text{prop}}{d}_1-{\eta }_{\mathrm{1,} \text{prop}}^2,i=\mathrm{2,} \dots , N.$

Введем замену переменных ${r_i}^2=x_i^2+y_i^2.$ Тогда система уравнений (2) преобразуется в:

$2\left[\begin{array}{cc}{x}_2-x_1& y_2-y_1\\ x_3-x_1& y_3-y_1\\ ⋮& ⋮\\ x_N-x_1& y_N-y_1\end{array}\right]\left[\begin{array}{c}x\\ y\end{array}\right]=-\left[\begin{array}{c}{r}_2^{2 }-d_2^2-r_1^2+{d_1}^2 \\ r_3^{2 }-d_3^2-r_1^2+{d_1}^2\\ ⋮\\ r_N^{2 }-d_N^{2 }-r_1^2+d_1^{2 }\end{array}\right]+$ (8)

$+\left[\begin{array}{c}2{\eta }_{\mathrm{1,}\text{prop}}{d}_{1 }+{\eta }_{\mathrm{1,}\text{prop}}^2-2{\eta }_{\mathrm{2,}\text{prop}}{d}_2-{\eta }_{\mathrm{2,}\text{prop}}^2\\ 2{\eta }_{\mathrm{1,}\text{prop}}{d}_{1 }+{\eta }_{\mathrm{1,}\text{prop}}^2-2{\eta }_{\mathrm{3,} \text{prop}}{d}_3-{\eta }_{\mathrm{3,} \text{prop}}^2\\ ⋮\\ 2{\eta }_{\mathrm{1,}\text{prop}}{d}_{1 }+{\eta }_{\mathrm{1,}\text{prop}}^2-2{\eta }_{N, \text{prop}}{d}_N-{\eta }_{N,\text{prop}}^2\end{array}\right].$

Система уравнений (8) имеет на одну строку и один столбец меньше, чем в способе (6). В матричном виде выражения (6) и (8) запишем:

$Ax=b,$ (9)

где матрица $A$ содержит координаты опорных точек, относительно которых происходит позиционирование, вектор $b$ – измеренные значения дальности до базовых станций, вектор $x$ – содержит координаты, которые необходимо найти.

Для решения матричного выражения (9) можно использовать метод наименьших квадратов (МНК).

$x_{ls}={\left(A^{T}A\right)}^{-1}{A}^{T}b.$ (10)

2. Чувствительность алгоритмов к ошибкам округления

Для оценки чувствительности алгоритмов (11) и (13) к ошибкам округления проанализируем число обусловленности матрицы $A,$ которое можно определить как произведение нормы матрицы на норму обратной к ней матрицы [5]:

$cond\left(A\right)=||A||||A^{-1}||.$ (11)

Анализ числа обусловленности осуществлялся с применением метода математического моделирования (см. рис. 2). В процессе моделирования был использован набор из 19 базовых станций (БС), размещенных в соответствии со стандартом 38.901 [6], который определяет расположение БС в плотной городской застройке. Расположение БС организовано в гексагональной структуре, при этом расстояние между БС задается параметром ISD (англ. Inter-Site Distance, ISD), который в данном исследовании составляет 500 м.

 

Рис. 2. Сравнение чисел обусловленности для двух методов решения системы уравнений (2)

Fig. 2. Conditionality numbers comparison for two methods of solving a system of equations (2)

 

Результаты моделирования указывают на более низкую устойчивость алгоритма TOA к ошибкам округления по сравнению с его модифицированным вариантом. Это объясняется присутствием дополнительного столбца, состоящего из постоянных значений, в матрице  что делает ее более близкой к вырожденной по сравнению с модифицированным вариантом.

3. Точность измерения координат

Для оценки точности измерения координат рассмотрим выражения (5) и (8). Из уравнения (5) дисперсия ошибки измерения координат при использовании алгоритма TOA имеет вид

$D\{2{\eta }_{i, prop}{d}_i+{\eta }_{i, prop}^2\}=4d_i^{2}D\left\{{\eta }_{i, prop}\right\}+D\left\{{\eta }_{i, prop}^2\right\}.$ (12)

При использовании модифицированного варианта алгоритма TOA (8) дисперсия ошибки измерения координат может быть вычислена следу ющим образом:

$D(2{\eta }_{\mathrm{1,} prop}{d}_1+{\eta }_{\mathrm{1,} prop}^2-2{\eta }_{i, prop}{d}_i-{\eta }_{i, prop}^2=$ (13)

$=4d_1^{2}D\left\{{\eta }_{\mathrm{1,}prop}\right\}+D\left\{{\eta }_{\mathrm{1,}prop}^2\right\}+$

$+4d_i^{2}D\left\{{\eta }_{i, \text{prop}}\right\}+D\left\{{\eta }_{i, \text{prop}}^2\right\}.$

В случае равноточных измерений дальностей от БС до БПЛА дисперсия ошибки оценки координат вторым методом в 2 раза выше, чем первым. На рис. 3 показан график зависимости точности оценки координат двумя методами в зависимости от числа БС. Для моделирования использовался сценарий Uma, описанный в стандарте 38.901 [6]. Считалось, что мощность излучаемого опорного сигнала составляет 49 дБм, ширина спектра частот опорного сигнала равна 16 МГц, коэффициент шума приемника БПЛА составляет 12 дБ. Стандартное отклонение ошибки (СКО) измерения координат оценивалось по 5000 ансамблей реализаций.

Из рис. 3 видно, что СКО измерения координат алгоритмом TOA ниже модифицированного варианта в 1,41 раза.

 

Рис. 3. Точность определения координат БПЛА алгоритмом TOA и его модифицированной версией в зависимости от количества БС

Fig. 3. Accuracy of the coordinate estimation of the UAV by the TOA algorithm and its modified version, depending on the number of BS

 

4. Вычислительная сложность алгоритмов оценки координат

Основную вычислительную сложность метода МНК составляют операции умножения матриц $A^{T}A$ и обращения матрицы ${\left(A^{T}A\right)}^{-1}.$ Существуют различные способы решения данной проблемы, выбор которых зависит от множества факторов, таких как размер матрицы, стабильность решения к ошибкам округления и т. д. [5].

В качестве критерия оценки сложности обычно используется количество флопов (англ. FLoating-point OPerations, FLOP) [5], требуемых на выполнение операций. Наиболее популярными методами решения уравнения (10) являются метод на основе разложения Холецкого $(2m{n}^2+n^3/3+$$mn+2n^2$ флопов) [5], метод симметричного QR-разложения $(2m{n}^2+4n^3/3+$$mn+2n^2$ флопов) [5] и метод асимметричного QR-разложения $\left(2n^2\right(m-n/3)+$$m^2+n^2$ флопов) [5].

Вычислительная сложность алгоритмов оценки координат БПЛА в зависимости от числа БС показана на рис. 4.

 

Рис. 4. Сравнение вычислительной сложности различных методов позиционирования

Fig. 4. Comparison of computational complexity of different positioning methods

 

Метод симметричного QR-разложения требует большего числа вычислений по сравнению с методом на основе разложения Холецкого. Однако метод QR более устойчив к малым погрешностям в данных, таким как шум в исходных данных. Это обусловлено наличием ортогональной матрицы Q, которая обеспечивает стабильность вычислений с плавающей точкой [5]. Основным преимуществом метода на основе разложения Холецкого является то, что операция факторизации осуществляется над матрицей меньшего размера $\left(A^{T}A\right)$ [5]. Выигрыш от данного подхода тем выше, чем больше БС используется для оценки координат БПЛА.

Метод асимметричного QR-разложения применяется к полной матрице А, размерность которой растет с увеличением числа БС. Однако, как показано на рис. 4, при малом числе БС данный метод может быть более эффективным.

Заключение

В данной исследовательской работе был про анализирован наиболее популярный метод позиционирования БПЛА, основанный на оценке дальностей до базовых станций. Для решения задачи вычисления координат были применены методы, основанные на алгоритме TOA. В ходе исследования были рассмотрены полная и модифицированная версии данного алгоритма. Результаты статьи показывают, что в сценарии Uma стандартное отклонение ошибки оценки координат составляет приблизительно 1 м. При этом точность оценки координат полной версии алгоритма превышает точность модифицированной версии в 1,41 раза, в то время как вычислительная сложность полной версии вдвое выше по сравнению с модифицированной. Выбор между использованием полной или модифицированной версии алгоритма зависит от характеристик вычислительной платформы, ресурсов, доступных вычислителю, требований к точности и скорости выполнения вычислений, а также от конкретных целей и ограничений задачи.

About the authors

Dmitry V. Tyurin

Moscow Technical University of Communications and Informatics

Email: turin88@bk.ru
ORCID iD: 0009-0001-4048-1959

lecturer of the Department of Infocommunication and Professional Disciplines, Volga-Vyatsky branch

Russian Federation, 15, Mendeleev Street, Nizhny Novgorod, 603011

Sergei V. Shishanov

Moscow Technical University of Communications and Informatics

Email: tribott@mail.ru
ORCID iD: 0009-0005-2955-8790
SPIN-code: 7344-5310
ResearcherId: J-3656-2017

Candidate of Technical Sciences, lecturer of the Department of Infocommunication and Professional Disciplines, Volga-Vyatsky branch

Russian Federation, 15, Mendeleev Street, Nizhny Novgorod, 603011

Vyacheslav V. Kazakov

Moscow Technical University of Communications and Informatics

Author for correspondence.
Email: vvfmtuci@mts-nn.ru
SPIN-code: 9550-7102

Candidate of Technical Sciences, associate professor, director of Volga-Vyatsky branch

Russian Federation, 15, Mendeleev Street, Nizhny Novgorod, 603011

References

Supplementary files