CC BY
50
Обработка результатов натурных испытаний с учетом траекторных наблюдений маневров
Ю.И. Юдин, Р.Г. Степахно
Судоводительский факультет МГТУ, кафедра УС и ПР
Аннотация. Предлагается новый подход к обработке результатов натурных экспериментов по маневрированию судов. Он состоит в более полном использовании траекторной информации маневра, что позволяет получать все кинематические характеристики и, тем самым, упростить последующую идентификацию параметров уравнения вращения судна вокруг вертикальной оси. Обработка использует сплайновую аппроксимацию траектории, параметризованной временем, и позволяет в некоторой мере сгладить влияние погрешностей результатов наблюдений. Однако метод предъявляет довольно жесткие требования к уровню этих погрешностей, достижение которого связано с дальнейшим совершенствованием средств обсервации места судна (например, использование DGPS или P-кода GPS).
Abstract. The new approach to processing of full-scale experiments' results on vessel manoeuvring has been proposed. It includes more complete using of the manoeuvre trajectory information, that allows to obtain all kinematic characteristics. The processing uses a spline-approximation of the trajectory parametrized by time, and allows to smooth the influence of the observation result errors to some extend. However, the method claimes a rather rigid requirements to the level of these errors and it leads to the further perfecting of the vessel observation tools (for example, using of DGPS or GPS P-code).
1. Введение
Одной из основных задач при постановке натурных экспериментов по маневрированию судна и последующей обработке их результатов является оценка параметров выбранной модели судна. Это весьма сложная задача, если учесть, что сама модель существенно нелинейна по параметрам состояния. Погрешности измерений наблюдаемых параметров при проведении экспериментов достигают значительных величин, а часть параметров в принципе ненаблюдаема и находится численным дифференцированием. Поэтому любое изменение методики обработки результатов таких экспериментов, ведущее к повышению точности конечных результатов, есть явление положительное.
Обычно (Соболев, 1976) при идентификации коэффициентов в уравнении управляемости, т.е. в уравнении для угловой скорости поворота судна вокруг вертикальной оси а>, используют только данные по измерениям курса судна в определенные моменты времени t. По этим измерениям находят угловую скорость как производную от курса по времени, используя конечные разности первого порядка. Известно, что точность такого численного дифференцирования крайне низка, и ее можно повысить несколькими способами (Ланцош, 1961). Это интерполяция опытных данных по нескольким точкам и последующее дифференцирование, или аппроксимация опытных точек, скажем, полиномами и их дифференцирование, или, наконец, применение сплайнов (Носач, 1994), которые представляют собой также набор полиномов, но "склеенных" особым образом в промежуточных точках. К сожалению, все эти способы чувствительны к погрешностям исходных данных, и нам приходится выбирать из них оптимальный в некотором смысле, например, в смысле минимума средней квадратической погрешности. Именно этим мы и займемся ниже.
2. Основные направления решения задачи
Проведем сравнительный анализ методов интерполяции, полиномиальной аппроксимации и сплайн-аппроксимации данных. При этом вначале в качестве данных будем рассматривать не данные реального опыта, а сгенерированные, "чистые" данные, погрешности которых можно устанавливать на любом уровне и, следовательно, контролировать. Это даст возможность оценивать погрешности полученных в обработке величин, что невозможно сделать с данными натурных экспериментов, когда отсутствуют эталонные результаты. И лишь после выбора лучшего способа обработки применим его для натурных результатов.
Первой рассмотрим интерполяцию в конечных разностях Ньютона, используя для нее пять точек. Тем самым строится интерполяционный полином четвертого порядка. Этот полином является скользящим, его коэффициенты изменяются при сдвиге пятерки точек по массиву экспериментальных точек на одну позицию. Это дает возможность последовательно определять первую и вторую
производные в первой точке из текущей выбранной пятерки. Пять точек - это традиционное количество обрабатываемых за один раз точек, дающее хорошую точность при сравнительно низких затратах вычислительных ресурсов - времени и объема памяти.
Интерполяционная формула Ньютона для 5-ти точек функции у = у(х) может быть записана в форме полинома 4-го порядка:
у(х) = У1 + (х-х1)4(х1,хг) + (х-х1)(х-х2)^2(х1,х2,хз) + (х-х1)(х-х2)(х-хз)/1з(хьх2,хз,х4) +
+ (х-х1)(х-х2) (х-х3)(х-х4)Л4(х1,х2,х3,х4,х5). (1)
Здесь хк, Ук - значения аргумента и функции в выбранных 5-ти точках, а Лк - конечные разности различных порядков. Они строятся рекуррентно, начиная с разности первого порядка:
4(хЬх2) = (У2 - У1) /(х2 - х1).
Лк(х1,...х+1) = (Лк-1 (хъ...,хк) - 4к-1(хь-• -,хк-1))/(хк-хО. (2)
Поскольку при обработке данных нас обычно интересует не только сама функция, но и ее производные, то выпишем выражения для двух первых производных в головной точке х1:
у'(х) = 4(хьх2) + (х1-х2)4>(хьх2,хз) + (х^Хх^хз^х^хзЛ) +
+ (х1-х2)(х1-хз)(х1-х4)^4(х1,х2,хз,х4,х5), (з)
у"(х) = 24>(х1,х2,хз) + 2[(х1-х2) + (х1-хз)]Лз(х1,х2,хз,х4) + + 2[(х1-х2)(х1-хз) + (х1-х2)(х1-х4) + (х^зХ^^/Цх^хз^хз). (4)
Этих формул достаточно, чтобы организовать вычисление самого интерполяционного полинома и его двух производных для произвольного массива точек, что и будет сделано ниже.
Вторым рассмотрим использование ортогональных полиномов на множестве из тех же пяти точек. Такой набор из четырех полиномов на отрезке [-2,+2] хорошо известен:
Ро(0 = 1, = 1/2, Р2(Ь = (г2 - 1)/2, Рз(0 = (513 - 17)/6. (5)
Их парные произведения обладают свойством ортогональности на совокупности равномерно распределенных на указанном отрезке из пяти точек -2, -1, 0, +1, +2:
2кР,(4)Р,<4) = 0 1Ф] , к = 1,... ,5, (6)
но они не нормированы к единице, т.е. сумма (6) не равна 1, когда / = ]. Однако их нормировка известна и нормирующие делители Ыр] найдены по общей формуле и равны, соответственно: 5, 5/2, 7/2, 10 и 70. Чтобы воспользоваться приведенными полиномами для обработки произвольных данных хк, используется преобразование переменных г = 2(2х - а - Ь)/(Ь - а), где а, Ь начальная и конечная точки интервала, т.е. в нашем случае а = х\(г = -2), Ь = х5. Правда, требуется равномерная расстановка значений аргумента хк на отрезке [а, Ь].
Теперь мы можем применить эти полиномы для аппроксимации функции, заданной наблюдаемыми значениями ук. Если представить функцию в виде линейной комбинации полиномов
у(г(х)) = 2,а,Р,«х)), (7)
то придется найти неопределенные коэффициенты такого представления. Их находят методом наименьших квадратов, что приводит к следующему выражению:
а} = ЪуРА 4) / ^Р] 4) = 2уР,( 4) / Ыр]. (8)
Найденное таким образом представление функции существенно отличается от представления с помощью интерполяционного полинома Ньютона (1). Если полином (1) дает кривую, которая проходит через все пять заданных точек, то представление (7) определяет кривую, которая может не пройти ни через одну из них. Однако она наилучшим образом представляет заданную совокупность опытных точек, приводя к минимуму сумму квадратов отклонений опытных точек от этой кривой, и потому обладает сглаживающим эффектом, чем не обладает полином Ньютона (1).
Две первые производные в головной точке пятерки в данном случае выглядят следующим образом:
у'(х}) = 0.5(х5 - х1)[а1/2 - 2а2 - 4заз/6],
у"(х0 = 0.25(х5 - х02[а2 - 10аз/6]. (9)
Наконец, сплайновое представление с помощью кубического сплайна. Это также представление функции по набору ее значений, но в таком представлении участвуют сразу все точки набора, а не выбранный пул точек (скажем, пять точек). Сплайн "склеивается" из полиномов третьей степени таким образом, что в точках склейки остаются непрерывными первая и вторая производные. Другая важная для нас особенность в применении сплайна состоит в том, что определенная методика позволяет сделать сплайн как интерполяционным, так и аппроксимирующим. Эта методика состоит в том, что для сплайна (f(x) минимизируется функционал Лагранжа следующего вида:
i [^"(x)]2dx + ZkA(^Xk)-y)2. (10)
Первый интеграл минимизирует среднюю кривизну сплайна (fix), а вторая сумма минимизирует взвешенную сумму квадратов отклонений сплайна от опытных точек. Изменяя множители Лагранжа -веса рк, мы можем делать сплайн ближе к интерполяционному (pk ^ 0), либо в той или иной мере аппроксимирующим (pk ^ да).
Если коэффициентами сплайна считать вторые производные в опытных точках Щ, то условие минимума (10) приводит к системе уравнений для коэффициентов
С\К\ + Ь\К2 + a1K3 = d\ b\K\ + c2K2 + b2K3 + a2K4 = d2
а].2К].2 + Ъ-ХКЯ + с3К3 + Ъ / + а/ = (11)
+ ЪЫ-2Км-2 + сЫ-\Км-\ + ЪЫ-\Км = dN_l + ЪЫ-\К;М-\ + сыКм =
В правых частях первого и последнего уравнений системы (11) участвует первая производная, соответственно, в первой и последней точке обрабатываемого набора данных. Поскольку эти производные неизвестны, то их следует найти каким-либо способом, например, с помощью обычной первой разности. Это может сказаться на результатах в виде краевого эффекта, т.е. только на точках, близких к крайним.
Это система линейных уравнений с пятидиагональной матрицей, метод ее решения хорошо известен как метод прогонки (прямой и затем обратной) и достаточно легко программируется. После определения коэффициентов сплайна Щ можно оперировать сплайном в виде кусочного полинома
<р(х) = Щл(хГх)ъЩ + К/х-хГ-1)3/6кГ + (хгх)/к/угл - К/12/6) +
+ (х-хя)%(у - Г/6), (12)
где х е (хг-1, хг), = хг-хг-1,/ = 2,..
Используя представление (12), можно записать явно выражения для первой и второй производных в произвольной точке х соответствующего интервала (х^ь х;):
<р\х) = -К]Л(х]-х)112Н] + К/х-х]Л)212Н / - (уя - К^2,/6)% + (у - Г/б)/^
^"(х) = К]Л(хгх)/к, + Г-хя)//г (13)
3. Сравнительный анализ численных результатов решения
Проведем сравнительный анализ трех методов, описанных выше. Для этого сгенерируем идеальные данные, т.е. данные не содержащие погрешностей наблюдений. Выберем характерный режим установившейся циркуляции судна с угловой скоростью 1 град/с на диаметре 200 м, т.е. с линейной скоростью 1.745 м/с или примерно 3.49 узла. Сгенерировано множество из 60 точек с фиксацией курса через каждые 6 град.
Применение всех трех методов к этому идеальному множеству данных дает одинаковые и практически точные результаты. А именно, сам курс и его первая производная - угловая скорость -находятся совершенно точно, а вторая производная отличается от точного нулевого значения только в методе аппроксимации не более чем на 10-6 град.
Однако если сгенерировать данные с некоторым уровнем погрешности, картина существенно изменяется. В этом случае ни один из методов не может восстановить точные значения ни самого курса, ни одной из двух его производных. Но теперь преимущества метода сплайн-аппроксимации очевидны.
Данные расчетов с помощью всех трех моделей - интерполяционной, аппроксимирующей и сплайновой, приведены ниже в табл. 1. Обработан один и тот же набор сгенерированных данных, причем при генерации заложена среднеквадратическая погрешность (СКП) курса в 2 град., а СКП
координат (они понадобятся далее для решения еще одной задачи обработки) в 20 м. В табл. 1 приведены результаты по восстановлению значений самого курса (столбец з) и его первой (столбец 4) и второй (столбец 5) производных. Напомним при этом, что точные значения двух последних параметров в нашем случае равны 1 и 0 соответственно. Данные приведены с четырехкратным разрежением - 15 вместо сгенерированных 60 - с целью сокращения материала. Они не подвергались какому-либо специальному отбору и остались полностью репрезентативными.
Таблица 1
Время, с 1 Курс, град (факт) 2 Курс, град (модель) 3 1-я производная (модель) 4 2-я производная (модель) 5
интер аппрок сплайн интер аппрок сплайн интер аппрок сплайн
6 8.з55 8. з55 8.з47 7.504 0.501 0.557 1.з01 0.1з1 0.120 0.00з
з0 з5.066 з5.066 з4.794 з4.610 -0.826 1.066 1.091 0.з20 -0.058 -0.002
54 57.597 57.597 57.604 57.200 1.58з 1.540 1.08з -0.182 -0.174 -0.001
78 82.606 82.606 82.605 81.575 0.8з9 0.846 1.186 0.01з 0.118 -0.002
102 104.748 104.748 104.68з 104.920 0.691 1.145 0.970 0.07з -0.017 0.001
126 128.679 128.679 128.921 129.164 з.252 1.572 1.006 -0.502 -0.166 0.000з
150 155.056 155.056 154.777 15з.904 -1.з96 0.545 1.174 0.510 0.121 -0.002
174 178.5з6 178.5з6 178.649 178.49з 1.024 0.2з8 0.911 0.052 0.210 0.001
198 201.77з 201.77з 201.651 202.452 -0.з47 0.50з 0.815 0.з49 0.179 0.002
222 225.425 225.425 225.467 227.166 2.57з 2.282 0.9з2 -0.з24 -0.266 0.00з
246 252.148 252.148 252.554 251.629 з.006 0.186 1.090 -0.з54 0.210 -0.001
270 272.598 272.598 272.911 275.062 4.175 1.999 0.801 -0.6з0 -0.195 0.004
294 297.4з1 297.4з1 297.509 297.з54 1.191 0.652 0.906 -0.04з 0.065 0.0005
з 18 з25.з04 з25.з04 з25.146 з2з.669 -1.з91 -0.29з 1.090 0.466 0.246 -0.005
з42 з42.87з з42.87з з4з.142 0.888 0.767 0.062 0.005
Хорошо видно, что методы интерполяции и аппроксимации полиномами не выдерживают конкуренции со стороны сплайн-аппроксимации. При этом следует отметить, что сплайн-аппроксимация проведена с весами рк = 10000. Как говорилось выше, это означает, что сплайн намного ближе к аппроксимирующему, он очень близко проходит к заданным при генерации точкам. Дальнейшее увеличение весов уже не дает практически никакого эффекта. Однако их уменьшение существенно ухудшает результаты. Так, даже при весах, равных 1000, сплайновая аппроксимация не намного лучше, чем аппроксимация полиномами или интерполяция. Особенно это верно для 1-ой и 2-ой производных. Итак, будем считать доказанным экспериментально, что сплайн-аппроксимация при весах порядка 10000 дает вполне удовлетворительные результаты по определению двух первых производных курса, заданного набором значений. Найденный результат вполне естественен с общих позиций. Так как в образовании сплайна участвуют одновременно все точки, то такая полнота информации ведет и к большей точности результатов.
4. Обработка траекторных наблюдений
Предыдущий раздел показал путь к достаточно точному дифференцированию экспериментальных данных по маневрированию. При этом использовалась лишь часть информации, которая наблюдается в натурных экспериментах (курсовой угол). В настоящее время при использовании таких средств обсервации как спутниковые навигационные системы (СНС) можно получить достаточно точные координаты траектории движения маневрирующего судна. Особенно верно это, если натурные эксперименты происходят в районах уверенного функционирования дифференциальных опорных пунктов СНС.
Эта дополнительная траекторная информация с хорошей точностью 5 м) при
соответствующей обработке дает возможность получить те кинематические параметры состояния объекта-судна, которые или непосредственно не наблюдаются (угол дрейфа и его производная) или наблюдаются недостаточно точно (скорость хода, угловая скорость поворота).
Скорость судна определяется как первая производная от пройденного пути по времени, т.е. это проблема численного дифференцирования, которая описана выше, и лишь источник информации другой: координаты, а не курс. Но траекторные наблюдения позволяют получать и угол дрейфа судна, его
производную и текущую кривизну траектории. Действительно, угол дрейфа есть разность между курсом судна и путевым углом. Но путевой угол есть угол наклона касательной к траектории и, следовательно, может быть найден как производная координаты У по координате X. Однако лучше обе проблемы - и со скоростью, и с углом дрейфа - решать в параметрическом представлении траектории. Параметром траектории естественным образом является время /. При этом скорость перемещения судна будет выражаться с помощью первых производных по времени от координат X и У:
V = ^[(Х/)2 + (У/)2]. (14)
Угол дрейфа определяется как разность между курсом и путевым углом, выраженным через те же производные:
р= К - ПУ = К - агС£[У//Х/]. (15)
Для целей параметрической идентификации модели судна необходима также и первая производная от угла дрейфа по времени. Ее также можно найти, используя параметрическое представление траектории. Но при этом появятся вторые производные от координат по времени:
ёр/Л = со-¥/Як = со-[УЙ"Х/ - УХа]/У2. (16)
Наличие в формуле вторых производных не усложняет нашей задачи, т.к. используя сплайн-аппроксимацию мы находим в качестве коэффициентов сплайна именно вторые производные. В формуле (16) использовано выражение для радиуса кривизны траектории
Як = V3 / [УЙ"Х/ - У/Хй"]. (17)
Для проверки предлагаемого подхода к обработке траекторией информации маневра снова сгенерируем файл данных установившейся циркуляции, но теперь, кроме курса, он будет содержать и координаты траектории. Генерация осуществлялась с помощью следующих формул для координат:
Х = хо + Як 8т(К) У = >-о - Як С08(К), (18)
где х0, >>0 - координаты центра циркуляции, К - путевой угол. В файл при этом выводится не сам путевой угол, а измененный на некоторую постоянную величину, которая и является определяемым далее углом дрейфа (в нашем случае 4 град.). Координаты (17) "портятся" случайными погрешностями с нормальным распределением и заданной СКП и затем заносятся в файл данных. Были выбраны два уровня погрешностей: 1) СКП курса 1 град., СКП координат 5 м; 2) СКП курса 0.5 град., СКП координат 2 м.
Таблица 2
, С 1 К факт 2 К модель 3 ЛгС^У/Х") 4 V 5 Як 6 Р = (3)-(4) 7
6 9.59 11.19 338.59 1.83 48.82
30 34.27 33.52 29.56 1.68 64.73 3.96
54 57.90 57.56 55.47 1.81 116.65 2.09
78 82.15 81.75 77.76 1.76 98.34 3.99
102 105.69 105.80 103.35 1.77 98.65 2.45
126 129.67 129.81 126.11 1.73 100.71 3.70
150 154.26 153.95 149.69 1.70 95.08 4.26
174 178.13 178.14 174.86 1.76 105.74 3.28
198 201.94 202.24 198.05 1.74 99.00 4.19
222 225.86 226.35 221.99 1.73 99.09 4.36
246 250.54 250.41 244.79 1.73 99.96 5.60
270 273.65 274.26 269.95 1.76 94.61 4.31
294 297.86 297.90 294.34 1.79 101.90 3.56
318 322.83 321.62 319.36 1.72 77.44 2.26
342 345.22 346.21 7.48 1.47 29.05
Для обработки этих сгенерированных данных использовалась сглаживающая сплайн-аппроксимация с весовыми коэффициентами рк = 10000. Она, как и в случае обработки курса, показала лучшие по точности результаты. Разумеется, идеально точные результаты дала обработка файла, не
содержащего погрешностей. Обработка же данных с погрешностями ведет и к погрешностям в конечных результатах. В табл. 2 приведены с четырехкратным прореживанием результаты обработки для второго варианта погрешностей исходных данных (0.5 град. по курсу и 2 м по координатам). Это высокие требования к точности, пожалуй, недостижимые в настоящее время на судах без специального оборудования для обсерваций. Заметим также значительное действие краевого эффекта - влияния задания краевых условий для определения коэффициентов сглаживающего сплайна. Здесь они заданы в виде первых производных, находимых по двум краевым точкам исходных данных. Хорошо видно, что все кинематические характеристики (за исключением нескольких краевых значений) найдены вполне удовлетворительно. Для этого достаточно сравнить полученные значения линейной скорости 5-ой колонки табл. 2 с точным ее значением - 1.745 м/с и значения радиуса кривизны 6-ой колонки с точным значением - 100 м. Угол дрейфа находится простым вычитанием из значения курса наблюдаемого (колонка 2) или смоделированного сплайном курса (колонка 3) величины аг^^/Х1) (колонка 4). Именно последняя разность и занесена в колонку 7. Точное значение, заложенное в данные и равное 4 град., определяется также достаточно хорошо.
Как уже было сказано выше, определение указанных кинематических характеристик движения судна во время маневрирования позволит проще идентифицировать параметры модели криволинейного движения судна. Обычно (Справочник по теории корабля, 1985) рассматривают систему уравнений криволинейного движения судна в следующем дифференциальном виде (Н(а) = со + у2 о}):
рекомендованным XIV Международной конференцией опытовых бассейнов (МКОБ) для задач по управляемости судов. Не имея данных по углу дрейфа и его производной, нет возможности напрямую использовать уравнение первого порядка (19). Поэтому для идентификации параметров этой модели приводят систему из двух уравнений к одному уравнению второго порядка и идентифицируют в этой эквивалентной модели необходимые 8 параметров (у1, у2, q\2, q3l, Г12, г31, ¿21, ¿31) Естественно, что такая идентификация для уравнения второго порядка значительно сложнее (см. описание этой непростой процедуры на с.240-248 указанного источника, причем точность подобной идентификации совершенно непредсказуема).
Возможность определения угла дрейфа и его первой производной позволяет идентифицировать с помощью непосредственно уравнения (19) сразу пять параметров, а остальные 3 параметра определить с помощью уравнения (20). При этом следует определить часть параметров из установившихся движений, например, трех циркуляций с разными кладками руля. Для остальных параметров придется использовать любое неустановившееся движение при маневрировании, в котором с необходимой точностью фиксируются курс судна и его координаты. Ими могут быть и классические маневры типа "Зигзаг" и "Спираль".
5. Заключение
Показано, что при обработке траекторных данных маневрирования судна наилучшие по точности результаты дает сплайн-аппроксимация при больших (~10000) коэффициентах веса. Однако даже такой аппарат требует высокой точности получения траекторных данных (~2 м). Найденные предлагаемым способом кинематические параметры движения судна могут быть использованы при параметрической идентификации модели судна. Именно это будет предметом дальнейших исследований авторов по данной тематике.
Корн Г., Корн Т. Справочник по высшей математике для научных работников. М., Наука, 831 е., 1978
Ланцош К. Практические методы прикладного анализа. М., Физматгиз, 423 е., 1961.
Носач В.Б. Решение задач аппроксимации с помощью ПК. СПб, Бином, 317 е., 1994.
Соболев Г.В. Управляемость корабля и автоматизация судовождения. Л., Судостроение, 477 е., 1976.
Справочник по теории корабля. Под редакцией Войткунского Я.И. Л., Судостроение, т.3, 544 е., 1985.
ёр/А + q12p + г12а-qH(aí)/q31 + = 0, ёа/А + qъlfi + г31®+ = 0,
(19)
(20)
Литература
