Введение

Интерполяция случайных процессов – один из важнейших аспектов обеспечения надежной и точной передачи информации. Развитие науки и техники, компьютерных технологий стимулирует современных ученых разрабатывать новые, более точные подходы к интерполированию случайных процессов [1–3]. Данная технология цифровой обработки сигналов широко применяется при разработке современных радиотехнических устройств [4; 5]. В области обработки цифровых изображений и графической информации интерполяция играет ключевую роль при восстановлении утерянных пикселей, устранении размытия и повышении разрешения видеосигнала. Применение этих методов позволяет формировать более реалистичные и точные визуальные модели объектов [6; 7]. Интерполяционные методы находят обширное применение в задачах радиофизики, например, для обработки сигналов радиолокаторов с синтезированной апертурой (SAR) [8]; для восстановления электромагнитных полей [9].

Впервые задачу интерполяции случайных процессов сформулировал А. Н. Колмогоров [10], который решил ее для произвольных стационарных последовательностей в случае, когда интервал интерполирования 0 < s < T состоит из единственной точки. Выделение класса случайных процессов с рациональной спектральной плотностью, для которого возможно построение явных формул интерполяции, связано с работами Дж. Дуба. Соответствующие аналитические выражения были предложены Н. Винером, а их строгое математическое обоснование представлено А. М. Ягломом [11]. В отличие от классических результатов для рациональных спектральных плотностей, в работе [12] авторами получены в явном виде интерполяционные формулы для процессов с квазирациональными спектральными плотностями и представлен пример численного расчета для L-марковского процесса с квазирациональной спектральной плотностью в интервале 1< τ <2.

Научная новизна данной работы заключается в разработке комплексного аналитико-численного подхода к синтезу вероятностной модели и расчету оператора интерполирования для L-марковского процесса с другим вариантом квазирациональной спектральной плотности, отличающимся от предыдущей работы авторов, на расширенном интервале, составляющем порядка 30 % выборки, с околонулевыми временными отсчетами. Новый результат, полученный в данной работе, является развитием предыдущего. Практическая реализация метода обеспечивается созданием виртуального прибора в среде LabVIEW для имитационного моделирования процесса и специализированного программного обеспечения на языке Python, реализующего сквозную векторизацию каскадных систем уравнений и автоматический анализ метрик качества модели. Устойчивость и высокая скорость вычислений достигаются за счет разложения подзнаковых квазирациональных слагаемых стохастического интеграла на простейшие дроби, проверки условий аналитичности, а также оптимизации сборки комплексного отклика через численные свертки тригонометрического полинома с ядрами Коши в точках спектральных полюсов.

Цель исследования – синтезирование вероятностной модели и расчет числовых значений оператора интерполирования для L-марковского процесса с квазирациональной спектральной плотностью для расширенного интервала интерполирования с околонулевыми временными отсчетами.

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

Постановка задачи: Найти наилучшее линейное приближение ξ(t;τ) к случайной величине ξ(t;τ), 0 < τ < T , зная значения процесса ξ(t) на двух полупрямых s ≤ t, s ≥ t + T.

В основу построения спектральной характеристики интерполирования в настоящей работе положен метод А. М. Яглома [11, с. 36], определяющий достаточные условия того, что некоторая функция комплексного переменного Фτ,Т(ω) является спектральной характеристикой интерполирования. Ключевыми из этих условий являются аналитичность функций Ф(1)τ,Т(ω) и Ф(2)τ,Т(ω), составляющих функции

(1)

в нижней (H–) и, соответственно, в верхней (H+) полуплоскостях (условие А); их возрастание в этих полуплоскостях не быстрее, чем некоторая степень |ω| (условие Б) и аналитичность функции

(2)

во всей плоскости (условие В).

В работе [11] на основе сформулированных выше достаточных условий были получены явные формулы интерполяции для процессов с рациональной спектральной плотностью

. (3)

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

, (4)

где Q0(ω) = ω + i, Q1(ω) = ω + 2i – многочлены первой степени, g > 0 – действительное число. Исследуемый стационарный гауссовский процесс обладает спектральной плотностью, заданной выражением (4). Согласно теореме Левинсона – Маккина – Питта, случайный процесс с такой спектральной плотностью относится к классу L-марковских процессов. В работе [13] обосновано, что если 1 / f(ω) – спектр однородного обобщенного поля {ξ*, Н*}, то условие L-марковости поля {ξ, Н} равносильно представлению спектральной плотности f(ω) в виде:

,

где μk – конечные меры с носителями в (L–L), . В частном случае классического марковского процесса L = {0}, а функция f(ω) является полиномом. Поскольку знаменатель G(ω) спектральной плотности (4) представляет собой целую функцию экспоненциального типа, исследуемый процесс гарантированно является L-марковским.

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

Соотношения между Т, g и τ существенно влияют на окончательную формулу спектральной характеристики интерполирования. В связи с этим целесообразно отдельно рассмотреть следующие случаи 1) Т = 2; 0 < τ < 1 и 2) Т = 2; 1 < τ < 2; для обоих вариантов g = 1.

Формулы спектральных характеристик экстраполирования и фильтрации, полученные авторами в работах [14; 15] для L-марковских процессов с квазирациональной спектральной плотностью, позволяют предположить, что в случае 1 (Т = 2; 0< τ <1) спектральную характеристику интерполирования для процесса со спектральной плотностью (4) следует искать в виде

(5)

где М(ω), М0(ω) – неизвестные на данном этапе многочлены, а W(ω), Г(ω) и Z(ω) – целые функции, рассчитываемые таким образом, чтобы для функций Фτ,Т(ω) и Ψτ,Т(ω)

(6)

выполнялись требования А), Б) и В).

Ключевым из этих требований будет условие В, согласно которому функция Ψτ,Т(ω) является целой. Этого можно добиться тогда, когда нули функции М(ω) – Г(ω) в формуле (٦) будут совпадать с нулями квазиполинома |G(ω)|2, которых, как известно, бесконечное множество. Очевидно, что это возможно только в том случае, когда

. (7)

Учитывая требование Б) о возрастании Ф(1)τ,Т(ω) и Ф(2)τ,Т(ω) в H– и, соответственно, H+ не быстрее некоторой степени |ω|, приходим к заключению, что числа αk при всех целых k должны удовлетворять неравенствам

. (8)

Для выполнения условия Б) на целые функции W(ω), Г(ω) и Z(ω) из представления (5) спектральной характеристики Фτ,Т(ω) достаточно принять следующие требования:

а) функция содержит только слагаемые с множителями exp(iωα) при α ≤ – τ или α ≥ Т – τ;

б) для целых функций W(ω) и Z(ω) имеют место представления: ; , в которых многочлены W0(ω) , W1(ω), Z0(ω), Z1(ω), D0(ω), D1(ω) и D2(ω) должны удовлетворять системам линейных разностных уравнений

, (9)

, (10)

. (11)

В системах (9)–(11) приняты следующие обозначения:

, , . (12)

Решением системы (9) в классе многочленов являются полиномы:

; ; ; . (13)

Решением системы (10) являются рациональные дроби

; . (14)

Решением системы (11) являются рациональные функции

; , (15)

где С1(ω) , С2(ω) –неизвестные на данном этапе многочлены.

Выбрав в качестве многочлена М0(ω) в знаменателях выражений (5) и (6) функцию вида

, (16)

можно обеспечить совпадение общих знаменателей дробей (١٤) и (١٥) с М0(ω). Тогда в качестве полиномиальных решений для функций W0(ω) , W1(ω), Z0(ω), Z1(ω) можно принять следующие выражения:

, , (17)

, . (18)

Подстановка и многочленов (13), (17) и (18) в выражения (5) и (6) с учетом равенств (9)–(11) позволяет получить окончательные формулы для функций Фτ,Т(ω) и Ψτ,Т(ω):

, (19)

. (20)

Чтобы выделить функции Ф(1)τ,Т(ω) и Ф(2)τ,Т(ω) из функции Фτ,Т(ω) для дальнейшего их исследования на аналитичность в нижней и, соответственно, верхней полуплоскостях, перепишем спектральную характеристику в виде

. (21)

Осталось удовлетворить требованию А) аналитичности функций Ф(1)τ,Т(ω) в H– и Ф(2)τ,Т(ω) в H+ и требованию В) целости функции Ψτ,Т(ω).

Из формул (20) и (21) следует, что знаменатели функций Фτ,Т(ω) и Ψτ,Т(ω) имеют четыре корня, совпадающих с корнями функции M0 = V2 – LU, а также двукратные корни функции V2:

; ;

; ;

– двукратные корни V2 из ∈ H+;

– двукратные корни V2.

Полученные восемь уравнений относительно восьми неизвестных коэффициентов многочленов третьей степени и представляют собой условия, обеспечивающие аналитичность функции Ф(1)τ,Т(ω) в H– (с ١-го по ٤-е уравнение системы) и функции Ф(2)τ,Т(ω) в H+ (с ٥-го по ٨-е уравнение).

Полученная система линейных алгебраических уравнений порядка ٨×٨ относительно неизвестных коэффициентов b0, b1, b2, b3, d0, d1, d2, d3) имеет следующий вид:

; ;

;

;

; ;

;

. (22)

В системе приняты обозначения: , , , , k = 1, 2, 3, 4, 5, 6; , .

Структура системы (22) гарантирует, что ее определитель отличен от нуля, следовательно, система имеет единственное решение.

В результате решения системы линейных уравнений (22) определены коэффициенты bk и dk (k = 0, 1, 2, 3) многочленов С1(ω) и С2(ω).

Для завершения доказательства необходимо показать выполнение требования В) о целости функции Ψτ,Т(ω). Целость данной функции непосредственно вытекает из системы уравнений (22), обеспечивающей аналитичность функций Ф(1)τ,Т(ω) в H– и Ф(2)τ,Т(ω) в H+ рассматриваемой области. Действительно, если v – корень многочлена М0(ω), то числитель функции Ψτ,Т(ω) обращается в нуль в точке ω = v в силу одновременного выполнения 1-го и 5-го (либо 2-го и 6-го) уравнений системы (22).

Таким образом, все требования (А, Б, В) к спектральной характеристике выполнены, следовательно, функция (21) является искомой спектральной характеристикой интерполирования.

Следующим шагом является построение оператора интерполирования . Как известно [11], данный оператор задается спектральным представлением вида

, (23)

где – случайная функция с некоррелированными приращениями, участвующая в спектральном представлении изучаемого случайного процесса

. (24)

Подстановка формулы (21) в выражение (23) приводит к следующему соотношению:

. (25)

Подынтегральные выражения в ١, ٢ и ٦-м интегралах являются неправильными рациональными дробями, в остальных шести – правильными. Выделив целую часть в интегралах J1, J2 и J6, разложим все получившиеся правильные дроби в выражениях (25) на простейшие. Это необходимо для того, чтобы свести каждое слагаемое к интегралу от исследуемого процесса ξ(s).

Проиллюстрируем этот шаг на примере последнего интеграла из системы (٢٥):

. (26)

Приведя дроби в подынтегральном выражении последнего интеграла к общему знаменателю и приравнивая коэффициенты при одинаковых степенях ω в числителях предпоследнего и последнего интегралов в формулы (26), получаем систему шести линейных уравнений относительно шести неизвестных x9k (неопределенных коэффициентов), k = 1, 2, 3, 4, 5, 6:

ω5 : ;

ω4 : ;

ω3 : ;

ω2 :;

ω : ;

Свободный член : . (27)

В системе (27) приняты следующие обозначения:

; ; ; ;

; ; ;

; ; ; ; ;

; ; ; ; ; ;

; ; ; . (28)

Таким образом, интеграл J9 принимает следующий вид:

(29)

Воспользуемся далее представлением

, ν = 1,2,3… (30)

и выразим интеграл J9 через случайный процесс ξ(s), приняв в формуле (30) ν = 1:

. (31)

Согласно спектральному представлению (٢٤) случайного процесса ξ(u), выражение в квадратных скобках формулы (٣١) представляет собой значение процесса в точке , потому интеграл J9 можно записать в виде

. (32)

Аналогичные рассуждения позволяют получить соответствующие представления для остальных восьми интегралов из формулы (25). Объединение всех найденных выражений приводит к окончательному виду оператора интерполирования:

(33)

В формуле (33) приняты следующие обозначения:

; ; ; ; ; ; ;

; ; ; ; ; r = 1;2;3;4;5;6.

; ; ; , r = ٥;٦, j = 4;5;6.

; ; j = 4;5;6, r = ٥;٦. (34)

Из формулы (33) следует, что оператор интерполирования исследуемого L-марковского процесса представляет собой сумму линейной комбинации значений процесса в точках ; ; и интегралов от значений процесса в точках zr (r = 1, 2, …, 8).

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

Моделирование процесса осуществлялось с помощью виртуального прибора, разработанного в среде LabVIEW. Формирование цифрового массива данных исследуемого процесса производилось на основе заданной спектральной плотности (4) посредством обратного преобразования Фурье. Данная процедура концептуально эквивалентна получению экспериментальных данных с помощью измерительного комплекса, состоящего из анализатора спектра и осциллографа (например, на базе приборов R&S®FSU и R&S®RTO).

Для обеспечения строгой методологической корректности вычислительного эксперимента и полного исключения эффекта утечки данных процедура построения аппроксимационной модели была жестко фрагментирована. Все отсчеты процесса, попадающие во внутренний интервал моделируемого пропуска (0; 2), были предварительно удалены из исходного массива, переведены в разряд тестовых данных и использовались авторами строго на финальном этапе для независимой оценки метрик качества алгоритма. Настройка структуры модели, включая расчет коэффициентов тригонометрического полинома Ak, Bk, а также минимизацию информационных критериев AIC/BIC для определения рационального числа гармоник, осуществлялась исключительно на основе доступной обучающей выборки, то есть вне интервала пропуска. На рис. 1 приведен график реализации процесса с искусственно удаленным фрагментом данных на интервале 0 < τ < 2.

Рис. 1. Исследуемый процесс ξ(u) и аппроксимирующий его тригонометрический полином Примечание: составлен авторами по результатам данного исследования

Из графика процесса ξ(u) (рис. 1) видно, что он близок к линейной комбинации тригонометрических функций sinx и cosx. В связи с этим целесообразно аппроксимировать ξ(u) тригонометрической регрессией вида

ξтеоретич(u) , (35)

где ее параметры А0, Аk и Вk (k = 1, 2, 3,…, N) определяются формулами

, , . (36)

Пошаговое увеличение числа гармоник N в формуле (35) и последовательный расчет RMSE проводился по формуле

. (37)

Для определения оптимального числа гармоник N исследован характер изменения RMSE, AIC и BIC в диапазоне от 1 до 40. Все три метрики одновременно достигают минимума при N = 15 (рис. 2 и 3) (RMSE = 1,2966; AIC = 1935,8067; BIC = 2432,8173), что подтверждает оптимальность модели с точки зрения баланса точности и сложности. При N > 15 RMSE возрастает и стабилизируется в районе 1,58, а AIC и BIC монотонно увеличиваются.

Принципиальное отличие предложенного подхода от классического МНК заключается в способе расчета коэффициентов. Коэффициенты ряда вычисляются по прямым аналитическим формулам (36) как интегральные оценки на доступном множестве точек, а не через минимизацию суммы квадратов остатков. Коэффициенты (N−1)-го порядка вычисляются один раз и фиксируются; переход к порядку N происходит путем удлинения ряда без пересчета предшествующих членов. Поскольку коэффициенты определяются независимо, их усечение не эквивалентно приравниванию к нулю в рамках МНК, а представляет собой исключение слагаемых с уже фиксированными весами.

Вследствие этого кривая RMSE имеет принципиально иной характер. После исчерпания полезного сигнала добавление избыточных гармоник вносит шум, который фиксированная модель не может сгладить пересчетом параметров. Ошибка начинает возрастать, образуя четко выраженный минимум при N = 15, который служит объективным критерием останова. В отличие от МНК, где RMSE на обучающей выборке монотонно убывает и не имеет внутреннего минимума, предложенный подход позволяет использовать саму кривую RMSE как прямой индикатор оптимальности модели без привлечения штрафных критериев.

Рис. 2. Зависимость RMSE от количества гармоник Примечание: составлен авторами по результатам данного исследования

Рис. 3. График зависимости BIC_AIC от числа гармоник Примечание: составлен авторами по результатам данного исследования

Следует подчеркнуть, что тригонометрическая аппроксимация по формулам (36) не является окончательным результатом восстановления. Она служит промежуточным приближением – аргументом в операторе интерполяции , который удовлетворяет граничным условиям на краях пропуска и корректирует искажения, вызванные неполнотой данных. Таким образом, погрешности прямых коэффициентов Фурье, связанные с применением ортогональных формул к неполной выборке, не переходят в финальный результат, что позволяет использовать простые и устойчивые формулы (36) в качестве входного приближения.

На основании вышеизложенного для всех последующих расчетов принято N = 15.

Аппроксимирующий тригонометрический полином при N = 15 гармоник использовался как базовый уровень (baseline) и промежуточный элемент численной схемы для вычисления значений в узлах оператора восстановления.

Подставим выражение (35) вместо ξ(u) в формулу (33), предварительно сделав в интегралах очевидную замену: u = zj – s.

Интегрирование произведений экспоненциальных и тригонометрических функций на интервале (0; ∞) позволяет получить окончательное выражение для оператора интерполирования:

, (38)

где

, (39)

, (40)

, (41)

где – выражение в квадратных скобках Y1.

‧

‧

. (42)

В формулах (39)–(42) параметры zj (j = 1, 2…, 9), αjr, Θjr и εjr определяются выражениями (34).

В работе выполнено табулирование значений оператора интерполирования при T = 2 на интервале 0 < τ < 1 с шагом hτ = 0,01. На основе полученных данных построены графики зависимости от параметра τ а) действительной и б) мнимой частей оператора интерполирования (рис. 4).

а) б)

Рис. 4. График зависимости а) действительных и б) мнимых частей интерполятора от параметра интерполирования τ Примечание: составлен авторами по результатам данного исследования

Мнимая часть полученного результата определяет фазовые изменения процесса и используется при анализе качества его фильтрации [15].

Разработанная авторами в среде Python программа позволяет вычислять значения оператора интерполирования для любого заданного момента времени τ (в рассматриваемом случае из диапазона (0, 1)).

Например, при τ = 0,06 .

Аналогично ;

.

Рассмотрим теперь случай Т = 2, 1 < τ < 2.

Аналогичные рассуждения, что приводились при выводе формулы (21) для спектральной характеристики исследуемого случайного процесса при Т = 2, 0 < τ < 1, привели к следующей формуле для спектральной характеристики Фτ,Т(ω) при 1 < τ < 2:

. (43)

В формуле (43) и – многочлены второй степени с неизвестными коэффициентами bk, dk, k = 0, 1, 2.

Следует отметить, что знаменатель спектральной характеристики в этом случае отличается от знаменателя первого случая, когда Т = 2, 0 < τ < 1.

Из представлений для Ф(1)τ,Т(ω) и Ф(2)τ,Т(ω) видим, что условие Б) выполнено, так как в первой из этих функций все слагаемые имеют множителями χ(ω)еiωα при α < 0, а во второй Ф(2)τ,Т(ω) – множителями χ(ω)еiωα при α > 0, так как 1 < τ < 2.

Требования А) аналитичности функций Ф(1)τ,Т(ω) в H– и Ф(2)τ,Т(ω) в H+ будут соблюдены, если будут выполнены следующие условия:

;

;

;

;

;

(44)

В формуле (44) , , , ωk, k = 1, 2, 3, 4, 5, 6 – корни знаменателя спектральной характеристики V·(V2 – LU).

Из системы (44) определяются неизвестные коэффициенты bk и dk (k = 0, 1, 2) многочленов С1(ω) и С2(ω).

Выполнение требования В) о целости функции

следует из условий (44) аналитичности функций Ф(1)τ,Т(ω) и Ф(2)τ,Т(ω). Докажем этот факт на примере корня . В данной точке выполняется последнее равенство системы (44), которое совпадает с числителем функции Ψτ,Т(ω) после его умножения на е-iωg. Первое и третье слагаемые числителя Ψτ,Т(ω) обращаются в нуль в этой точке, так как V(-2,5i) = 0. Таким образом, все требования А), Б) и В) к спектральной характеристике выполнены, следовательно, функция Фτ,Т(ω), заданная формулой (43), является спектральной характеристикой L-марковского процесса со спектральной плотностью (4) при Т = 2 и 1 < τ < 2.

Подставляя формулу (43) в спектральное представление интерполятора (23), получим

. (45)

Разложим подынтегральные дроби во всех интегралах в (45) на простейшие дроби подобно тому, как это было сделано в случае Т = 2 и 0 < τ < 1, и далее рассчитаем неопределенные коэффициенты xjk (j – номер интеграла, k – номер неизвестного корня знаменателя) , j = 1, 2…8; k = 1, 2…6 – числители простейших дробей в случае первых четырех интегралов J1, J2 , J3, и J4 из четырех систем четырех уравнений с четырьмя неизвестными xjk , где j = 1, 2, 3, 4; k = 1, 2, 3, 4.

; ; ; . (46)

В случае интегралов J5 и J6 xjk рассчитываются из двух следующих систем с шестью неизвестными xjk , где j = 5, 6; k = 1, 2, 3, 4, 5, 6.

1) ;

2) ;

3) ;

4) ;

5) ;

6)

и, наконец, для интегралов J7 и J8 – из двух систем линейных уравнений с шестью неизвестными xjk, где j = 7, 8 определяет номер интеграла, а k = 1, 2, 3, 4, 5, 6 – номер неизвестного корня знаменателя.

; .

Правые части σjk данных систем уравнений выражаются через коэффициенты bk и dk, k = 0, 1, 2 многочленов С1(ω) и С2(ω), определенных ранее из системы (44). Соответствующие расчетные соотношения, полученные методом неопределенных коэффициентов, имеют вид

σ11 = 2i; σ12 = –2; σ13 = 4i; σ14 = –3; σ21 = b1+ib2; σ22 = b0+ib1–3b2; σ23 = ib0+2b1; σ24 = 2b0–7b2;

σ31 = –2i; σ32 = – 1; σ33 = –5i; σ34 = –4; σ41 = d1–id2; σ42 = d0–id1–3d2; σ43 = –id0+2d1; σ44 = 2d0–7d2;

σ51 = d1+2id2; σ52 = d0+2id1–9d2/2; σ53 = 2id0+3d1+4id2; σ54 = 3d0+ 4id1– 31d2/2; σ55 = 4id0+4d1;

σ56 = –4d0–35d2/2; σ61 = b1–2ib2; σ62 = b0–2ib1–9b2/2; σ63 = –2ib0+3b1–4ib2; σ64 = 3b0–4ib1–31b2/2;

σ65 = –4ib0+4b1 ; σ66 = 4b0–35b2/2; σ75 = b1 ; σ76 = b0–5b2/2; σ85 = d1 ; σ86 = d0–5d2/2.

Параметры AA, AB, AC, AD, AF, AG, φ1, φ2, φ3, φ4, а также μk, Nk и æk (k = 1, 2, 3, 4) определены выражениями (28), полученными при рассмотрении случая Т = 2 ; 0 < τ < 1.

Приняв далее во всех восьми интегралах и аппроксимируя стохастические интегралы от ξ(u) тригонометрическим полиномом (35) с гармониками, после их вычисления придем к окончательной формуле для оператора интерполирования при Т = 2 и 1 < τ < 2: , где ;

,

где – выражение в квадратных скобках формулы для Y1;

; .

В этой формуле ; ; ; ; ;

; ; ; ; ; ; ; j = 1,2,3,4;

; ; j = 1,2,3,4,5,6; ; ; j = 5,6.

; ; ; ; ; ; ; .

В работе показано, что если для марковского процесса с рациональной спектральной плотностью вида Sрац(ω) = (где К(ω) – многочлен) «оператор интерполирования зависит от значений процесса в концах интервала интерполирования, то есть в точках s = t и s = t + T» [11], то для исследуемого L-марковского процесса с квазирациональной спектральной плотностью (4) ситуация меняется. При T = 2 оператор зависит от значений процесса в следующих точках:

в шести точках ; ; ; ;

; , общих для обоих диапазонов параметра τ;

в двух дополнительных точках и при 0 < τ < 1;

в двух дополнительных точках и при 1 < τ < 2.

С помощью разработанной авторами программы на языке Python выполнено табулирование значений оператора интерполирования при T = 2 на интервале 1 < τ < 2 с шагом hτ = 0,01. Это позволяет вычислять численные значения оператора для любого заданного момента времени из указанного диапазона. Например, при τ = 1,15 значение оператора составляет ; аналогично при τ = 1,28 оно будет равно .

На заключительном графике (рис. 5) продемонстрирована эффективность полученных формул спектральных характеристик при 0 < τ < 1 и 1 < τ < 2, а также соответствующих им операторов интерполирования.

Рис. 5. Исходные, в том числе «скрытые», данные L-марковского процесса и решение задачи интерполирования Примечание: составлен авторами по результатам данного исследования

Метрики оценки качества аппроксимационной модели

Примечание: составлена авторами на основе полученных данных в ходе исследования.

Для независимой проверки точности восстановления данных на тестовом интервале 𝜏 ∈ (0;2) был рассчитан комплекс взаимодополняющих критериев. Выбранный пакет метрик позволяет всесторонне оценить величину погрешности, структуру ошибок и геометрическое соответствие моделируемой кривой исходному процессу. Результаты расчетов приведены в таблице.

Анализ результатов вычислений позволяет сделать следующие выводы:

– Разработанная модель обладает высокой объясняющей способностью: коэффициент детерминации R2 = 0,908 подтверждает, что построенная модель успешно описывает 90,8 % дисперсии амплитуды исследуемого процесса на участке лакуны.

– Взвешенная относительная ошибка wMAPE составила 29,1%. Использование данной метрики взамен классической MAPE методологически обосновано знакопеременной природой стационарного гауссовского процесса. Для крупного моделируемого пропуска, составляющего 31,83 % от общей длины реализации (1591 отсчет из 4998), полученный уровень погрешности свидетельствует об устойчивости алгоритма к сингулярностям.

– Средняя абсолютная ошибка составляет MAE = 0,802. Соотношение RMSE/MAE ≈ 1,19. Близость этого коэффициента к единице указывает на отсутствие на тестовом интервале локальных аномальных всплесков, паразитных осцилляций и краевых эффектов (явлений Рунге или Гиббса). Отклонения модели от истинных значений распределены гомогенно по всему скрытому фрагменту, что подтверждает гладкость и вычислительную устойчивость предложенной интерполяции. Поскольку исследуемый стационарный гауссовский процесс обладает свойством эргодичности, полученные на данной реализации оценки являются статистически устойчивыми и инвариантными относительно конкретной траектории.

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

Эксперимент 1: устойчивость к выбору числа гармоник. Была выполнена серия численных экспериментов с варьированием числа гармоник N в диапазоне от 10 до 20 при фиксированном количестве наблюдаемых данных. Оптимальные аналитические показатели достигнуты при N = 15 (R2 = 0,908, RMSE = 0,954, MAE = 0,802). В диапазоне N ≥ Nopt качество остается стабильным: R2 ∈ [0,773;0,908], а относительное приращение ошибки RMSE не превышает 14,07 %. При N < Nopt наблюдается резкое снижение качества вплоть до отрицательных значений R2, что объясняется недостаточностью количества гармоник для аппроксимации сигнала и автоматически исключается критерием минимума RMSE на наблюдаемых данных. Таким образом, устойчивость метода к переопределению модели количественно подтверждена тремя независимыми метриками и охарактеризована через CV(RMSE), а не через единичный показатель.

Эксперимент 2: инвариантность к локальным флуктуациям. Для исследования робастности алгоритма к 2 % наблюдаемых точек добавлялся импульсный шум амплитудой до 2σ, после чего коэффициенты тригонометрического полинома и оператор ξ(τ) пересчитывались полностью. Численный эксперимент повторялся 5 раз с генерацией различных реализаций шума с последующим усреднением результатов. Отклонение ξnoisy от ξclean оценивалось по четырем показателям: MAE = 0,037, RMSE = 0,046, max|Δξ| = 0,106, средняя относительная погрешность равна 1,11 % при максимальной 3,18 %. Выявленная инвариантность объясняется тем, что коэффициенты тригонометрического полинома являются интегральными характеристиками сигнала: локальные выбросы усредняются при вычислении коэффициентов и не оказывают существенного влияния на результат интерполяции.

Использование нескольких независимых метрик (MAE, RMSE, max|Δξ|, относительная погрешность) исключает зависимость выводов от какой-либо одной характеристики. Полученные количественные результаты подтверждают как вычислительную устойчивость метода в окрестности N = 15, так и его инвариантность к локальным флуктуациям наблюдаемых данных.

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

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

Выводы

В работе предложена и исследована вероятностная модель линейного интерполятора для L-марковского процесса с квазирациональной спектральной плотностью, обобщающей класс классических рациональных функций. В основе разработанного математического аппарата лежит метод спектральных представлений А. М. Яглома.

По результатам теоретических исследований и вычислительного эксперимента сформулированы следующие ключевые выводы:

Разработана аналитическая модель оператора интерполяции L-марковского процесса, базирующаяся на спектральной теории случайных процессов, аппарате разностных уравнений и теории функций комплексного переменного. На основе теоремы Левинсона – Маккина – Питта доказана L-марковость исследуемого стационарного гауссовского процесса. Установлено его фундаментальное отличие от классического марковского процесса: если в традиционном случае число граничных точек строго ограничено двумя, то для исследуемого L-марковского процесса с квазирациональной спектральной плотностью их количество возрастает до восьми. Данный эффект обеспечивает более детализированное и точное восстановление пропущенных фрагментов.

Имитационное моделирование L-мар-ковского процесса по заданной спектральной плотности реализовано в виде виртуального прибора в среде LabVIEW. Вычислительное ядро алгоритма интерполяции и аналитические процедуры реализованы в виде специализированного программного обеспечения на языке Python.

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

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