Постановка задачи в классическом подходе

 

  В классическом подходе полагается, что производятся наблюдения некоторой скалярной величины d, которая связана с оцениваемыми величинами х₁, х₂,... х _N линейным уравнением следующего вида, названным уравнением наблюдения (измерения):

 

         а _I * х₁ + b _I * х₂  + c _I * х ₃ +.... = d _I ,                   (4.11)

 

где I - порядковый номер измерения,

    а _I,b _I, c _I  - коэффициенты уравнения наблюдения.

 

     Для оценки вектора Х = [ х₁, х₂, х ₃ ] проведем  m   измерений, причем m > n    Для каждого из N измерений составим уравнение наблюдения вида (4.11). Набор этих уравнений в совокупности даст нам избыточную систему уравнений следующего вида:

 

           а ₁ * х₁ + b ₁ * х₂  + c ₁ * х ₃ +.... = d ₁   

           а ₂ * х₁ + b ₂ * х₂  + c ₂ * х ₃ +.... = d ₂ 

                                                                                                  (4.12)

           ..............................................................

           а _m * х₁ + b _m * х₂  + c _m * х ₃ +.... = d _m

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

      Оценка вектора Х определяется методом наименьших взвешенных квадратов (МНВК). В этом методе для каждого из  m уравнений наблюдений рассчитываются невязки между измерением и прогнозом оценки:

 

         d _I = а _I * х₁ + b _I * х₂  + c _I * х ₃ +.... - d _I         (4.13)

 

и формируется критерий оптимальности оценки:

 

      J = S _I=1 ^m k _I * d _I ²,                                                (4.14)

 

Искомая оценка вектора Х получается из условия минимума критерия оптимальности:

 

         min J (X):     ¶ J (X) / ¶ X = 0,                     (4.15)

           X

 

где производная от скаляра J по вектору Х даст нам следующую систему N уравнений (N - размерность вектора Х):

 

¶ J / ¶ x₁ = 0 = 2* S _I=1 ^m k _I * d _I *  a _I ,

¶ J / ¶ x ₂ = 0 = 2* S _I=1 ^m k _I * d _I * b _I,

                                                                                            (4.16)

.............................................................

¶ J / ¶x_N = 0 = 2* S _I=1 ^m k _I * d _I * n _I,

 

В уравнениях (4.16) a _I - коэффициент компоненты x₁ в I-ом уравнении наблюдения,  n_I - коэффициент компоненты x_N в I-ом уравнении наблюдения.

Решение системы уравнений (4.16) и даст нам искомые оценки вектора Х. Данная система уравнений является линейной относительно составляющих вектора Х и может быть представлена в следующем виде:

 

           A * X = B.                                                     (4.17)

 

Решением (4.17) является выражение (4.18):

 

          X = A ^-1 *B,                                                         (4.18)

 

где А и В - матрицы коэффициентов.

Решение имеет место, когда детерминант матрицы А не равен нулю.

 

ТОЧНОСТЬ ОЦЕНИВАНИЯ

 

  Точность оценивания вектор Х определяется ковариационной матрицей его оценок:

 

           Сov { X } = E{ [ X- E [X] ] *[ X- E [X] ] ^T },     (4.19)

 

где математическим ожиданием оценки Х в силу несмещенности оценки является истинное значение вектора Х, обозначенное как Х^ист.

 

Из уравнения наблюдения (4.1) и уравнения фильтра МНВК (4.10) получим следующее выражение:

 

X = (H^T * R ^-1 * H) ^-1  * H ^T * R ^-1 * Z =

       (H^T * R ^-1 * H) ^-1  * H ^T * R ^-1 * H * X ^ист -

       (H^T * R ^-1 * H) ^-1 * H ^T * R ^-1 * V.                           (4.20)

 

Используя последнюю формулу, получим уравнение для разности оценки и истинного значения вектора Х:

 

Х - Х ^ист = (H^T * R ^-1 * H) ^-1 * H ^T * R ^-1 * V.               (4.21)

 

Подставляя (4.21) в (4.19) получим выражение для расчета ковариационной матрицы оценок вектора Х:

 

Сov { X } = E{ (H^T * R ^-1 * H) ^-1 * H ^T * R ^-1 * V *

                    V ^T * R ^-1* H* (H^T * R ^-1 * H) ^-1 }             (4.22)   

 

В силу того, что матрицы Н, R являются детерминированными, они выносятся за операцию математического ожидания. Случайным является вектор V, математическое ожидание произведения

 Е{V * V ^T} которого, вследствие условия несмещенности измерений: Е{V} =0, равно ковариационной матрице вектора V -

 

                                    Е{ V * V ^T } = R.  

 

В результате сокращений произведения R *R ^-1  и

                                   (H^T * R ^-1 * H)^-1*H^T * R ^-1 * H получим 

 

     Сov { X } = (H^T * R ^-1 * H) ^-1.                                  (4.23)

 

Отметим, что при составлении алгоритма в соответствии с классическим подходом выражение H^T * R ^-1 * H фактически является матрицей А из формулы (4.17), поэтому в этом случае

 

      Сov { X } = А ^-1.                                                    (4.24)

 

РАСЧЕТ ПАРАМЕТРОВ ДВИЖЕНИЯ МЕТОДОМ НВК

ПО ДАННЫМ ЛОКАЦИИ

 

 

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

В процессе определения параметров движения объектов управления: их курсов и скоростей, необходимо вести счисление пути своего судна. Проекции этого пути на оси координат < X,Y> вычисляются следующим образом:

 

S_{x, i} = S_{x, i-1} + V_c * DT * sin K c,

                                                                                               (4. 25)

S_{y, i} = S_{y, i-1} + V_c * DT * cos K c,

 

где S_{x, i}, S_{y, i}, S_{x, i-1}, S_{y, i-1},- проекции пути судна на оси X,Y в моменты i и i-1 соответственно, V_c, K c - курс и скорость судна,

DT - интервал времени между i и i-1 замерами пеленга и дистанции до объекта наблюдения.

     

Дисперсии Var { П } и Var {D } измерений полагаются известными.

     Для произвольного i -ого измерения координаты объекта x _i, y _i описываются следующей формулой:

 

 

x _i = D _i * sin П _i + S _{x, i},

                                                                                           (4.26)

y _i = D _i * cos П _i + S _{y, i},          

 

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

 

 

x _i = x ₁ + V _x * (T - T ₁),

                                                                                           (4.27)

y _i = y ₁ + V _y * (T -  T ₁ )     _,

 

где V _x,V _y - проекции скорости объекта на оси координат.

 

   Объединим уравнения (4.26) и (4.27). При этом, умножим уравнение проекций на ось Х на величину sin П _i , а уравнение проекций на ось Y на величину cos П _i  . Таким образом мы получим следующее итоговое уравнение:

 

 

y ₁ * cos П _i + x ₁ * sin П _i + V _x * (T - T ₁₎ * sin П _i +

 

+V _y * (T - T ₁) * cos П _i = D _i + S _{x, i} * sin П _i + S _{y, i}* cos П _i  .   

         (4.28)

 

  Это уравнение является линейным уравнением относительно переменных y ₁, x ₁, V _x, V _y   и может быть представлено в следующем виде:

 

 

        a _i * y ₁ + b _i * x ₁ + c _i * V _x + d _i * V _y = z _i,                     (4.29)

 

где a _i, b _i, c _i, d _i, z _i - коэффициенты уравнения, определяемые следующими выражениями:

 

 

a _i= cos П _i,b _i = sin П _i , c _i =(T - T ₁₎ * sin П _i , 

d _i = (T - T ₁₎ * cos П _i ,                                                                    (4.30)

z _i = D _i + S _{x, i} * sin П _i + S _{y, i}* cos П _i

  

Составляя избыточную систему уравнений для m измерений, получим следующую запись:

 

a ₂ * y ₁ + b ₂ * x ₁ + c ₂ * V _x + d ₂ * V _y = z ₂

a ₃ * y ₁ + b ₃ * x ₁ + c ₃ * V _x + d ₃ * V _y = z ₃

                                                                                                                                                 (4.31)

  ...................................................................

a _m * y ₁ + b _m * x ₁ + c _m * V _x + d _m * V _y = z _m

 

 Систему уравнений (4.31) можно представить в векторном виде:

 

H * X = Z,                                                                              (4.32)

 

где X = [ y ₁, x ₁, V _x, V _y ] ^T,

 

               a ₂, b ₂, c ₂, d ₂

               a ₃, b ₃, c ₃, d ₃

   H =

               a _m, b _m, c _m, d _m  ,

 

   Z = [ z ₂, z₃,...............z _m ] ^T.

 

Тогда в соответствии с рассмотренным выше математическим аппаратом МНВК, получим следующее выражение для выработки оценок вектора Х:

 

  X _k = (H_k ^T * R _k ^-1 * H_k) ^-1  * H_k ^T * R _k^-1 * Z_k = P_k * B _k,            

                                                                                          (4.33)

 

где P_k = Сov { X_k } = (H_k ^T * R _k ^-1 * H_k) ^-1  = А _k ^-1, k - номер замера и номер соответствующего ему шага итерации процедуры оценивания. 

Матрица А _k имеет следующий вид:

 

A_k =

S _i=2 ^m a _i ²*k _i; S _i=2 ^m a _i*b _i*k _i; S _i=2 ^m a _i*c _i* k _i; S _i=2 ^m a _i*d _i* k _i

 

S _i=2 ^m a _i*b _i*k _i; S _i=2 ^m b _i²*k _i; S _i=2 ^mb_i*c _i* k _i; S _i=2 ^m b _i*d _i*k _i

    

S _i=2 ^m a_i* c_i* k_i; S _i=2 ^mb_i*c_i*k_i; S _i=2 ^m c_i ²* k_i;   S _i=2 ^m c_i* d _i* k _i

       

S _i=2 ^m a_i *d_i* k_i; S _i=2 ^m b_i*d_i*k _i; S _i=2 ^m c_i*d_i*k _i; S _i=2 ^m d _i ² * k _i;

(4.34)

 

Матрица В _к вычисляется следующим образом:

               S _i=2 ^m a_i*h_i*k_i

 

               S _i=2 ^m b_i*h_i*k _i

B _k =                                                                                  (4.35)

               S _i=2 ^m c_i*h_i*k _i

 

                        S _i=2 ^m d_i*h_i*k _i

 

 

После расчета вектора Х по формуле (4.33), используя формулы (2.15),(2.16), (2.17)-(2.20), вычисляем значения курса, скорости объекта и их дисперсий. Списки исходных, используемых алгоритмом параметров, и параметров, определяемых в результате решения задачи, которые необходимы для составления алгоритма по данному математическому описанию, приведены в таблицах 4.1 и 4.2.

 

Таблица 4.1                               Используемая информация

	Величина	Обозначение
1.	Замеры пеленга на объект, град.	ПI
2.	Замеры дистанций до объекта, м.	DI
3.	СКВО измерений пеленга, град.	sПI
4.	СКВО измерений дистанции, м.	sDI
5.	Собственный курс управляющего судна, град	КS
6.	Собственная скорость управляющего судна, м/с.	VS
7.	Время I- ого замера, сек.	t I

 

 

Таблица 4.2                                      Результаты решения задачи

№	Величина	Обозначение
1.	Начальное значение координаты Х объекта, м	x1
2.	Начальное значение координаты Y объекта, м	y1
3.	Проекция скорости на ось Х, м/c	V X
4.	Проекция скорости на ось Y, м/c	V Y
5.	Скорость, м/c	V
6.	Курс, рад	K
7.	Ковариационная матрица вектора X = [ y 1, x 1, V x, V y ] ^T	Сov { Xk }
8.	Дисперсия скорости,, м^2/c^2	Var { V^*}
9.	Дисперсия курса, рад^2	Var { K^*}

 

 

ЗАДАНИЕ К КУРСОВОЙ РАБОТЕ

 

1. По приведенным выше соотношениям составить блок-схему модуля расчета оценок параметров движения объекта наблюдения и дисперсий ошибок оценивания.

 

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

 Блок-схема должна включать в себя модули:

- движения наблюдателя (см. приложение 2),

- движения объекта наблюдения (см. приложение 2),

- расчета истинных значений дистанций и пеленгов,

- упрощенных моделей функционирования систем гидролокации или радиолокации, обеспечивающих имитацию случайных гауссовых ошибок измерения (см. приложение 3) и наложение ошибок на истинные значения измеряемых параметров для моделирования измеренных значений П и D,  

- модуль расчета оценок ПДО и их расчетных дисперсий,

- модуль вычисления статистических ошибок  оценивания при повторении N раз моделируемого процесса (см. приложение 3).

Структурная схема модели приведена в приложении 1.1

 

3. Разработать программу статистического исследования алгоритма на языке программирования, заданном преподавателем.

 

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

 

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