Scientific journal
Modern high technologies
ISSN 1812-7320
"Перечень" ВАК
ИФ РИНЦ = 0,940

THE CHOICE OF BICUBIC SPLINE SMOOTHING PARAMETERS IN PROBLEMS OF NONPARAMETRIC IDENTIFICATION

Voskoboynikov Yu.E. 1, 2 Boeva V.A. 3
1 Novosibirsk State University of Architecture and Civil Engineering (SIBSTRIN)
2 Novosibirsk State Technical University
3 Novosibirsk State University of Architecture and Civil Engineering (Sibstrin)
In the last two decades, the so-called Voltaire series have been used to describe the dynamics of nonlinear systems in terms of “input-output”. The most commonly used is the quadratic term of the series, which generates a two-dimensional kernel and the corresponding two-dimensional Volterra equation. To identify the two-dimensional impulse transition function, an approach based on the inversion of the Volterra equation and using second-order partial and mixed derivatives of the system output signal are proposed. The fundamental problem of implementing this approach is the stable calculation of second-order partial and mixed derivatives from noisy data. To solve this problem, a two-dimensional smoothing bicubic spline (abbreviated SBS) can be used. However, the construction of SBS in a practice faces a fundamental difficulty – the choice of smoothing parameters. In the case of SBS, the problem of choosing a parameter becomes more complicated due to the construction of a whole ensemble (set) of cubic (one-dimensional) splines when smoothing is performed over one of the two variables of the initial data. Therefore, two types of SBS smoothing parameters are introduced in the work: scalar and vector. The scalar parameter is the same for all cubic splines in the given variable. The vector parameter is a vector whose projections are defined as the smoothing parameters for each cubic spline. The paper proposes algorithms for estimating the optimal values of both scalar and vector parameters. The computational experiment performed showed a high accuracy in estimating the optimal smoothing parameters by the proposed algorithms.
the problem of identifying the quadratic kernel of a nonlinear system
smoothing bicubic splines
scalar and vector smoothing parameters
estimating optimal smoothing parameters

Два последних десятилетия ведутся интенсивные научные исследования методов идентификации нелинейных динамических систем, представленных различными математическими моделями. Весьма перспективной в этом отношении является интегральная модель «вход-выход», состоящая из нескольких уравнений Вольтерра, ядра которых образованны из соответствующих слагаемых ряда Вольтерра [1, 2]. Особенно часто используется ряд, ограниченный квадратичным членом ряда Вольтерра [3]. В этом случае связь между входным сигналом φ(τ) моделируемой стационарной системы и ее выходным сигналом f(t) можно представить следующим интегральным соотношением [3]:

missing image file

missing image file,(1)

где k1(τ), k2(τ,s) – линейное и квадратичное ядра Вольтерра соответственно. Нужно идентифицировать эти ядра по измеренным значениям входного и выходного сигналов идентифицируемой системы.

Очевидно, что сигнал f(t) есть сумма двух неизвестных сигналов f1(t) и f2(t) – первое и второе слагаемые в (1), которые можно интерпретировать как выходные сигналы линейной и квадратичной «подсистем». Следовательно, на первом этапе идентификации необходимо выделить из f(t) сигналы f1(t), f2(t) по отдельности. Для произвольного входного сигнала такая задача неразрешима. В работах [4, 5] предложена эффективная методика такого разделения, которая подразумевает проведение активной идентификации, когда на вход системы подаются прямоугольные сигналы заданной формы и заданной амплитуды. В результате получаем следующие интегральные соотношения:

missing image file (2)

missing image file (3)

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

missing image file (4)

а из (3) – формула обращения [4, 6]:

missing image file (5)

missing image file

где используются следующие обозначения:

missing image file

missing image file (6)

К сожалению, реализация на практике полученных формул обращения сталкивается с принципиальной трудностью – некорректностью операции дифференцирования [6]. Одно из проявлений некорректности заключается в больших ошибках вычисления производной даже при очень малых погрешностях задания значений дифференцируемой функции. Таким образом, устойчивое дифференцирование зашумленных данных становится актуальной задачей для реализации формулы (5) на практике.

В работе [7] был построен алгоритм идентификации на основе формулы (4), где для устойчивого вычисления первой производной использовался сглаживающий кубический сплайн (СКС) дефекта единица с выбором параметра сглаживания из условия минимума среднеквадратической ошибки сглаживания. Алгоритм показал приемлемую точность идентификации при решении практических задач [8]. В случае идентификации квадратичного ядра k2(τ,s) использование сглаживающих сплайнов существенно усложняется. Во-первых, для вычисления смешанной производной второго порядка missing image file нужно строить уже сглаживающий бикубический сплайн (СБС), являющийся функцией двух переменных t, v. Во-вторых, из-за разной «гладкости» функции f2(t, v) по разным переменным необходимо выбрать уже несколько параметров сглаживания из условия минимума ошибки сглаживания. В-третьих, в большинстве случаев дисперсия шума (погрешностей) измерения функции f2(t, v) неизвестна.

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

− введение нескольких типов параметров сглаживания – скалярного и векторного параметров для СБС;

− раздельный выбор параметров сглаживания (по каждой переменной сплайна) из условия минимума ошибки сглаживания при неизвестной дисперсии шума измерения на основе метода L-кривой;

− исследование эффективности предложенных алгоритмов выбора с точки зрения минимизации ошибки сглаживания.

Материалы и методы исследования

Предположим, что значения двумерной функции f2(t, v) измеряются при значениях аргументов ti , i = 1…Nt , vj , j = 1…Nv, т.е. в узлах прямоугольной сетки {ti, vj}, i = 1…Nt, j = 1…Nv. Заметим, что узлы ti, и τj могут иметь неодинаковый и неравный шаг, что встречается в реальных экспериментах.

Для учета возможных погрешностей (шумов) измерений принимается следующая модель измеренных значений missing image file:

missing image file

missing image file (7)

где ηij – случайный шум измерения с нулевым средним и дисперсией сетки {ti, vj} (равноточные измерения). Требуется по исходным данным missing image file вычислить значения производных missing image file в заданных узлах сетки {ti , vj}, i = 1…Nt, j = 1…Nv.

Замечание. Из вида интеграла (3) следует, что функция f2(t, v) принимает ненулевые значения для аргументов, удовлетворяющих условию v ≤ t. Для отрицательных значений v, t функция равна нулю в силу условия технической реализуемости системы, т.е.: k2(t, v) ≡ 0, если v < 0, t < 0. Так как значения производных в (5) используются только при v ≤ t, то для устранения разрыва первого рода при значениях v = t предлагается дополнить значения функции f2(t, v) для v > t по правилу:

missing image file

missing image file (8)

Дополненную таким образом функцию будем обозначать как missing image file.

Для устойчивого вычисления производных missing image file обратимся к сглаживающим кубическим сплайнам (СКС) [9], широко используемым при обработке экспериментальных данных [10, 11], включая задачи идентификации [8]. Алгоритм построения двумерного сглаживающего бикубического сплайна как построение ансамбля (набора) одномерных СКС был изложен в работе [12]. Поэтому здесь кратко опишем основные шаги этого алгоритма, необходимые для понимания места и способов выбора параметра сглаживания. Сделаем это для нахождения смешанной производной missing image file, так как вычисление частной производной второго порядка missing image file сводится к двукратному дифференцированию одномерного сплайна с аргументом v.

Шаг 1. Для каждого j = 1…Nv формируется набор данных (фиксируется значение vj):

missing image file (9)

выбирается параметр сглаживания α1(j), и по исходным данным (9) строится СКС missing image file, допускающий на отрезках [ti, ti+1), i = 1…Nt – 1 представление:

missing image file

missing image file (10)

и имеющий непрерывную вторую производную (по переменной t) на всем интервале missing image file. По сплайну (10) вычисляются значения первой производной missing image file (оценка для производной missing image file). Шаг 1 повторяется для vj , j = 1…Nv и таким образом строится ансамбль (набор) из Nv одномерных СКС, каждому из которых необходимо задать параметр сглаживания α1(j), j = 1…Nv.

Шаг 2. Для каждого i = 1…Nt формируется набор данных (фиксируется значение ti):

missing image file (11)

выбирается параметр сглаживания α2(j) и строится СКС missing image file, первая производная которого является оценкой missing image file для смешанной производной missing image file, где missing image file – коэффициент сплайна missing image file в представлении:

missing image file

missing image file (12)

Шаг 2 повторяется для ti, i = 1…Nt, и таким образом строится ансамбль (набор) из Nt одномерных СКС, каждому из которых необходимо задать параметр сглаживания α2(j), i = 1…Nt.

Введем для СБС два новых понятия: скалярный и векторный параметр сглаживания. Под скалярным параметром будем понимать параметр, величина которого одинакова для всех сплайнов данного ансамбля. Так, на шаге 1 скалярный параметр определяется как

α1(j) = α1, j = 1…Nv , (13)

Векторный параметр сглаживания – это вектор, каждая проекция которого задает свой параметр сглаживания соответствующему СКС из ансамбля сплайнов. Так, на шаге 2 векторный параметр missing image file имеет i-ю проекцию, равную

missing image file. (14)

Возникает главный вопрос при построении сплайнов на практике: как выбрать скалярный или векторный параметры сглаживания? При ответе на этот вопрос существенное значение играет наличие достоверной информации о дисперсии шума измерений missing image file (7).

Первоначально предположим, что дисперсия missing image file шума измерения известна с точностью 5–10 %. В этом случае для выбора скалярного параметра сглаживания можно использовать алгоритм, основанный на проверке статистической гипотезы об оптимальности параметра сглаживания и предложенный в работе [11]. Для проверки этой гипотезы (предположим, что для шага 1) вводится критерий

missing image file (15)

где missing image file – невязка i-го измерения в j-м сплайне missing image file(9). В качестве скалярного параметра сглаживания α1 принимается значение α1W, для которого выполняется неравенство

missing image file (16)

где величины missing image file – квантили χ2-рас-пределения с N = Nt · Nv степенями свободы уровней β/2 и 1 – (β/2) соответственно. Величина β определяет вероятность ошибки первого рода при проверке гипотезы об оптимальности (величина среднеквадратической ошибки сглаживания) оценки α1W и, как правило, β = 0.05. Если N > 30, то для вычисления квантилей χ2-распределения при β = 0.05 можно использовать простые формулы:

missing image file

missing image file (17)

Вычисление α1W сводится к решению нелинейного уравнения

ρWS(α1) = N, (18)

итерационными алгоритмами. В качестве α1W принимается приближённое решение уравнения (18), которое удовлетворяет неравенству (16). Вычислительный эксперимент показал [12], что рассмотренный алгоритм выбора скалярного параметра позволяет с приемлемой точностью (увеличение ошибки сглаживания по сравнению с минимальной 6–10 %) оценить оптимальный скалярный параметр сглаживания α1opt, минимизирующий среднеквадратическую ошибку (СКО) фильтрации.

Такой же алгоритм выбора можно построить для вычисления проекций векторного параметра сглаживания. В качестве j-й проекции вектора missing image file принимается значение missing image file, удовлетворяющее неравенству

missing image file (19)

где

missing image file (20)

Вычисление missing image file сводится к решению нелинейного уравнения

ρWV(α1(j)) = N, (21)

итерационными алгоритмами. В качестве missing image file принимается приближённое решение уравнения (21), удовлетворяющее неравенству (19). Результаты сравнения ошибок сглаживания при использовании этих двух типов параметров сглаживания приводятся ниже.

Заметим, что вычисление как скалярного, так и векторного параметров сглаживания на основе проверки статистической гипотезы об оптимальности параметра требует задания достоверной величины дисперсии шума измерения. При заниженном значении дисперсии вычисляется заниженное значение (по сравнению с оптимальным) параметра сглаживания, что приводит к «недостаточной» фильтрации шума измерения – высокий уровень «остаточного» шума. При завышенной дисперсии наблюдается обратный эффект – сглаживаются «тонкие» структуры обрабатываемого сигнала. К сожалению, ситуация достоверного знания дисперсии шума измерения в реальном физическом эксперименте скорее исключение, чем правило. Кроме того, при вычислении смешанной производной на шаге 2 сплайны строятся по результатам дифференцирования на шаге 1, которые содержат «шум дифференцирования» неизвестной дисперсии. Все это говорит об актуальности разработки алгоритмов выбора параметров сглаживания для случая неизвестной дисперсии шума измерения, которые и составляют предмет дальнейших исследований.

В работе [13] был построен алгоритм выбора параметра сглаживания одномерного СКС на основе метода L-кривой, позволяющий с приемлемой точностью оценить оптимальный параметр сглаживания как в случае некоррелированного, так и для коррелированного шума измерений. Сам метод L-кривой иногда используется для выбора параметра регуляризации в алгоритмах решения некорректных задач (например, [14–15]), когда неизвестны характеристики погрешностей исходных данных. Поэтому попытаемся модифицировать алгоритм выбора параметра сглаживания работы [13] для вычисления скалярного и векторного параметров сглаживания.

В случае одномерного СКС в качестве параметра сглаживания принималась величина αL, являющаяся решением вариационной задачи:

missing image file (22)

Кривизна kL(α) L-кривой определяется по формуле

missing image file(23)

где

missing image file,

missing image file;

missing image file.

Одномерный СКС Sn,α(ti) строился по набору данных missing image file, а весовые множители pi в случае равноточных измерений задаются одинаковыми, например pi = 1. Для эффективного вычисления функционала γ(α) (куда входит вторая производная сплайна) в работе [13] предложена формула

missing image file (24)

где hti = ti+1 – ti, i = 1…n–1, ci , di – коэффициенты СКС, вычисленные при заданном параметре α.

Такой же алгоритм выбора можно построить для вычисления проекций векторного параметра сглаживания (для определенности на шаге 2). В качестве i-й проекции вектора missing image file принимается значение missing image file, доставляющее максимум кривизны kL(α) L-кривой (23), вычисленной для сплайна missing image file. В этом случае функционалы ρ(α), γ(α) для этого сплайна определяются выражениями

missing image file;

missing image file,

где hvj = vj+1 – vj, j = 1…Nv – 1, missing image file – коэффициенты сплайна missing image file.

Для вычисления скалярного параметра сглаживания (для определенности на шаге 2) определим функционалы ρ(α), γ(α) следующим образом:

missing image file,

missing image file,

где pi,j – весовые множители для i, j-го измерения (для равноточных измерений задаются одинаковыми), missing image file – коэффициенты сплайна missing image file, построенного при параметре α2, hvj = vj+1 – vj, j = 1…Nv – 1.

Далее в соответствии с (23) находится значение кривизны и, решая вариационную задачу (22), находим значение α2L, при котором кривизна L-кривой максимальна, и это значение является скалярным параметром сглаживания бикубического сплайна.

Результаты исследования и их обсуждение

В качестве «тестового» квадратичного ядра Вольтерра k2(τ,s) (уравнение (1)) была принята функция, изображенная на рис. 1 и часто используемая для моделирования энергетических объектов [5]. Был выполнен обширный вычислительный эксперимент по идентификации этого ядра. В силу ограниченности объема статьи кратко рассмотрим только результаты эксперимента, связанные с вычислением смешанной производной missing image file функции f2(t, v) (уравнение (3)) по бикубическому сплайну со скалярным и векторным параметрами сглаживания при неизвестной дисперсии шума измерения. Эта ситуация является более сложной и наиболее типичной в реальном эксперименте. При этом основное внимание уделяется исследованию эффективности скалярного и векторного параметров сглаживания, которая характеризует, насколько увеличивается ошибка сглаживания при этих параметрах по сравнению с оптимальными скалярным и векторным параметрами сглаживания. На рис. 2 показана смешанная производная missing image file, имеющая достаточно сложную форму. Параметры используемой сетки {ti, vj}: Nt = 80, Nv = 80, ti = 0.03·(i – 1), i = 1…Nt , vj = 0.03·(j – 1), j = 1…Nv . В качестве краевых условий на всех четырех границах прямоугольной области задавались нулевые вторые производные.

missing image file

Рис. 1. Квадратичное ядро k2(τ,s)

missing image file

Рис. 2. Смешанная производная missing image file

Точные значения «дополненной» функции missing image file (замечание 1) в узлах сетки {ti, vj} (матрица F* с элементами missing image file) искажались нормально распределенным шумом с нулевым средним и дисперсией missing image file (7). Относительный уровень δη шума измерения определялся выражением

missing image file,

где || || означает евклидову норму матрицы. Сформированная таким образом матрица missing image file с элементами missing image file missing image file является исходной для построения СБС. Так как параметры кубических сглаживающих сплайнов задаются на шаге 1 (параметр сглаживания α1(j)) и шаге 2 (α2(i)), то здесь же определялись и относительные ошибки сглаживания:

• на шаге 1:

скалярная

missing image file;

векторная

missing image file;

• на шаге 2:

скалярная

missing image file;

векторная

missing image file.

Очевидно, что на шаге 2 определяется уже относительная ошибка сглаживания первой производной. Обозначим через missing image file, missing image file, missing image file, missing image file – наименьшие относительные ошибки сглаживания при использовании оптимальных параметров сглаживания. Аналогично величины missing image file, missing image file, missing image file, missing image file определяют относительные ошибки при выборе параметров на основе метода L-кривой описанными выше алгоритмами. Для количественного сравнения относительных ошибок сглаживания введем следующие коэффициенты эффективности:

missing image file; missing image file;

missing image file; missing image file.

Очевидно, что эти коэффициенты меняются в интервале от 0 до 1 и их значения являются случайными величинами. Поэтому по выборке объемом 30 вычислялись выборочные средние missing image file. Чем больше эти средние отклоняются от 1 в меньшую сторону, тем больше проигрыш по точности сплайна, построенного с параметрами missing image file. Средние значения коэффициентов эффективности приведены в таблице для трех уровней шума: 0,02; 0,04; 0,06. Анализ этой таблицы показывает, что:

− векторный параметр сглаживания имеет меньший коэффициент эффективности по сравнению со скалярным параметром при прочих равных условиях построения сглаживающих сплайнов;

− с увеличением уровня шума измерения наблюдается уменьшение значений коэффициентов эффективности;

− коэффициент эффективности при сглаживании функции больше соответствующего коэффициента эффективности при сглаживании производных.

Средние значения коэффициентов эффективности

Средние значения

коэффициентов

эффективности

Относительный уровень шума измерений

0,02

0,04

0,06

missing image file

0,989

0,972

0,949

missing image file

0,902

0,873

0,848

missing image file

0,908

0,863

0,791

missing image file

0,816

0,774

0,694

Заключение

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

Исследование выполнено при финансовой поддержке РНФ в рамках научного проекта № 22-21-00409.