Радиотехника и электроника, 2021, T. 66, № 8, стр. 772-781

ВВЕДЕНИЕ

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

Существует множество практических приложений, требующих оценки параметров импульсного сигнала неизвестной формы и длительности, регистрируемого на фоне шума с неизвестными характеристиками. Такие задачи рассматривались, например, в работах [7, 8], посвященных исследованию источников помех в прецизионных, в том числе криогенных, измерительных системах. Работа [9] посвящена измерению времени прихода импульса артериального давления при оценке кардиоваскулярного риска. В работах [10, 11] предложены алгоритмы оценки параметров импульсных сейсмических сигналов в распределенных системах охраны протяженных периметров. При оценке параметров сигналов в рассмотренных работах авторы, как правило, используют модель аддитивного белого гауссовского шума, аппроксимацию сигнала известной функцией и стандартные методы оценки (байесовский, максимального правдоподобия и т.п.). Здесь качество оценки во многом определяется точностью аппроксимации сигнала и адекватностью применяемой модели аддитивного белого шума.

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

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

1. МОДЕЛЬ НАБЛЮДАЕМОГО ПРОЦЕССА

В качестве модели наблюдаемого на временном интервале T_x процесса x(t) примем аддитивную смесь полезного сигнала s(t) и фоновой составляющей η(t): $x(t) = s(t) + \eta (t)$, при этом сигнал s(t) имеет конечную длительность τ₀ и полностью находится в пределах интервала наблюдения T_x, т.е. ${{\tau }_{0}} \leqslant {{T}_{x}}$.

Характер фона во многом зависит от конкретного приложения. Например, в прецизионных системах низкочастотных измерений существенную роль играет фликкер-шум и электромагнитные наводки. В том случае, когда интервал наблюдения много меньше времени корреляции фликкер-шума или периода электромагнитной наводки, текущую реализацию фона можно представить в виде суммы квазидетерминированной линейной составляющей $g(t) = {{a}_{0}} + {{a}_{1}}t$ с неизвестными параметрами a₀, a₁ и белого гауссовского шума $\xi (t)$ с нулевым средним и неизвестной спектральной плотностью: $\eta (t) = {{a}_{0}} + {{a}_{1}}t + \xi (t)$.

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

Сигнал s(t), наблюдаемый на интервале T_x, характеризуется энергией, временным положением и длительностью. Введем понятие текущей энергии сигнала $E(t)$, накопленной в интервале наблюдения T_x к текущему моменту времени t:

Тогда полная энергия полезного сигнала на интервале T_x определится как

а распределение накопленной на интервале T_x энергии сигнала будет определяться функциями нормированных текущей энергии и “плотности” распределения энергии

Рассмотрим теперь понятия временного положения и длительности сигнала. В литературе имеются различные подходы к их определению. Одним из них является “энергетический” подход, при котором временное положение сигнала $\bar {t}$ определяется как координата центра тяжести текущей энергии сигнала на интервале наблюдения T_x:

Аналогично при энергетическом подходе за длительность сигнала принимают эффективную длительность – среднеквадратическое отклонение от координаты центра тяжести:

Альтернативным способом оценки временных параметров сигнала является подход, предложенный авторами в работе [12], также основанный на энергетических свойствах сигнала. За положение T₀ сигнала на оси времени принимается медиана распределения накопленной энергии на интервале T_x, т.е. решение уравнения $e({{T}_{0}}) = {{L}_{0}} = {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0em} 2}$, а длительность сигнала определяется длительностью интервала $\tau = {{T}_{2}} - {{T}_{1}}$, в течение которого энергия сигнала возрастет от ${{E}_{{{{T}_{1}}}}} = {{L}_{1}}{{E}_{s}}$ до ${{E}_{{{{T}_{2}}}}} = {{L}_{2}}{{E}_{s}}$, где L₁, L₂ – относительные пороговые уровни текущей энергии, T₁, T₂ – моменты времени, в которые происходит пересечение соответствующего порогового уровня. Эти определения временного положения и длительности сигнала используются в данной работе.

На рис. 1а приведены три импульсных сигнала: прямоугольный, двойной экспоненциальный и косинус-квадрат с гармоническим заполнением. Все импульсы имеют одинаковые энергии ${{E}_{s}} = 1$ и длительность по основанию τ₀. На рис. 1б приведены графики зависимости текущей энергии от времени для этих же импульсов (сплошные линии) и пороговые значения энергии (пунктир). Относительные пороговые значения энергии приняты равными ${{L}_{0}} = 0.5$, ${{L}_{1}} = 0.2$ и ${{L}_{2}} = 0.8$.

2. СИНТЕЗ АЛГОРИТМОВ ОЦЕНКИ ПАРАМЕТРОВ СИГНАЛA И ФОНОВОЙ СОСТАВЛЯЮЩЕЙ

В качестве исходных данных для синтеза алгоритмов оценки параметров фоновой составляющей, энергетического параметра и параметров временного положения и длительности импульсных сигналов примем M-мерный выборочный вектор $\vec {x} = {{\left\{ {{{x}_{1}},{{x}_{2}},...,{{x}_{M}}} \right\}}^{T}}$, составленный из статистически независимых отсчетов наблюдаемого процесса (1), взятых в моменты времени kΔt, $k = 1,...,M$, где Δt – шаг дискретизации, MΔt = T_x, ${{\left\{ \cdot \right\}}^{T}}$ означает транспонирование. Вектор $\vec {x}$ в общем случае представляет собой сумму векторов полезного сигнала $\vec {s} = {{\left\{ {{{s}_{1}},{{s}_{2}},...,{{s}_{M}}} \right\}}^{T}}$, линейной составляющей фона

и шума $\vec {\xi } = {{\left\{ {{{\xi }_{1}},{{\xi }_{2}},...,{{\xi }_{M}}} \right\}}^{T}}$:

Компоненты вектора $\vec {s}$, параметры a₀ и a₁ вектора $\vec {g}$ являются априорно неопределенными детерминированными величинами, а компоненты вектора $\vec {\xi }$ – статистически независимыми гауссовскими величинами с нулевым средним значением и неизвестной дисперсией σ².

Поскольку исходными данными для синтеза алгоритма является дискретная выборка из наблюдаемого процесса, то вместо выражения (2) будем пользоваться его дискретным аналогом:

Очевидно, что полная энергия сигнала на интервале наблюдения T_x, подлежащая оценке, равна E_M.

Плотность распределения вероятностей вектора $\vec {x}$ имеет следующий вид:

Если дисперсия шумовой составляющей фона σ², параметры a₀, a₁ и форма сигнала неизвестны, то число неизвестных параметров распределения (5) на 3 больше числа отсчетов M наблюдаемого процесса к моменту времени M, поэтому необходима дополнительная информация о характеристиках наблюдаемого процесса. Поскольку форма сигнала по условию является неизвестной, то единственной возможностью получения дополнительной информации является расширение интервала наблюдения за счет включения в него участка T_y, заведомо содержащего только фоновый процесс. Воспользуемся методом контраста, известным из теории обнаружения сигналов [13]. Указанный метод состоит в том, что на интервале времени ${{T}_{y}} = N\Delta t$, в котором заведомо отсутствует полезный сигнал, формируется “фоновый” выборочный вектор $\vec {y} = {{\left\{ {{{y}_{1}},{{y}_{2}},...,{{y}_{N}}} \right\}}^{T}}$, а затем на интервале T_x – выборочный вектор $\vec {x}$, содержащий отсчеты полезного сигнала и фоновой составляющей. Будем полагать, что интервалы T_y и T_x следуют непосредственно друг за другом.

Совместное распределение вероятностей выборочных векторов $\vec {x}$ и $\vec {y}$ имеет следующий вид:

Для удобства синтеза алгоритма оценки параметров сигнала и фона представим выражение (6) в следующем виде:

где – нормирующий множитель.

Полученное распределение принадлежит экспоненциальному семейству, характеризуется (M + 3) параметрами с достаточными статистиками:

для параметра ${{ - 1} \mathord{\left/ {\vphantom {{ - 1} {(2{{\sigma }^{2}})}}} \right. \kern-0em} {(2{{\sigma }^{2}})}}$ –

для параметра a₀/σ² – для параметра ${{{{a}_{1}}\Delta t} \mathord{\left/ {\vphantom {{{{a}_{1}}\Delta t} {{{\sigma }^{2}}}}} \right. \kern-0em} {{{\sigma }^{2}}}}$ – и Z_k = x_k для параметров s_k/σ², k = 1, …, M. При неизвестных σ², a₀, a₁, s₁, …, s_M вектор параметров принимает значения из области $( - \infty ,0) \times ( - \infty ,\infty ) \times $ $ \times \,\,( - \infty ,\infty ) \times ... \times ( - \infty ,\infty )$, т.е. содержит (M + 3)-мерный интервал, поэтому достаточные статистики являются полными [14].

Для оценки дисперсии шума σ², параметров линейной составляющей a₀, a₁ и текущей энергии сигнального импульса E_m воспользуемся методом, аналогичным методу моментов. Классический метод моментов предполагает выражение моментов распределения вероятностей наблюдаемых данных нескольких порядков через искомые параметры и отыскание последних решением полученной системы уравнений. Затем истинные моменты заменяются выборочными, и при больших объемах выборки полученная оценка стремится к истинному значению. В нашем случае для составления системы уравнений будем использовать выражения для начальных моментов первого и второго порядка полных достаточных статистик. В полученном решении истинные моменты заменим значениями полных достаточных статистик в соответствующей степени. Полученные таким образом оценки будут зависеть от исходной выборки только через полные достаточные статистики. Особенность оценок, зависящих от полных достаточных статистик, состоит в том, что если искомый параметр представляет собой линейную комбинацию моментов полных достаточных статистик произвольного порядка, то полученная оценка будет эффективной (несмещенной и с минимальным значением квадратичной функции потерь) [13], причем эта оценка будет существенно единственной [16]. Если же зависимость оценки от полных достаточных статистик нелинейная, то полученная оценка может иметь смещение, однако и в этом случае она будет минимизировать квадратичную функцию потерь в классе K_b оценок, имеющих смещение b, и будет существенно единственной [17].

Получим выражения для следующих начальных моментов полных достаточных статистик:

В выражениях (8)–(12) ${\mathbf{M}}\left\{ \cdot \right\}$ означает вычисление математического ожидания.

Решая систему уравнений (8)–(12) относительно σ², a₀, a₁, $\sum\nolimits_{k = 1}^M {s_{k}^{2}} $, получим

Заменяя в выражениях (13)–(16) моменты ${\mathbf{M}}\left\{ \cdot \right\}$ полными достаточными статистиками, получим выражения для оценок интересующих параметров:

Для получения оценки текущей энергии ${{\hat {E}}_{m}}$ в выражении (20) следует заменить M на m. В принципе, в выражении (20) в силу случайного характера компонентов выборочных векторов может получиться отрицательное значение энергии, особенно при малых объемах выборки, как это имеет место на начальном этапе оценивания текущей энергии ${{\hat {E}}_{m}}$. В этом случае следует принять ${{\hat {E}}_{m}} = 0$.

Окончательные выражения полученных оценок нетрудно получить, заменив ${{Y}_{0}},{{Y}_{1}},{{Z}_{0}},{{Z}_{1}},...,{{Z}_{M}}$ их значениями. Приведем окончательные выражения для оценок параметров фоновой составляющей:

Из выражений (21)–(23) видно, что параметры фона оцениваются только по фоновой выборке на интервале T_y. На интервале T_x выделить фоновую составляющую не представляется возможным, так как ни уровень, ни форма сигнала, ни его длительность и положение на временной оси неизвестны.

Оценки ${{\hat {a}}_{0}}$ и ${{\hat {a}}_{1}}$ – несмещенные и эффективные, так как они являются линейными комбинациями полных достаточных статистик. Оценки $\widehat {{{\sigma }^{2}}}$ и ${{\hat {E}}_{m}}$ представляют собой нелинейные функции полных достаточных статистик, поэтому требуют отдельного исследования на несмещенность.

Для оценки временного положения T₀ и длительности τ сигнала выберем моменты времени, когда оценка текущей энергии ${{\hat {E}}_{m}}$ достигает значений ${{E}_{{{{T}_{0}}}}} = {{L}_{0}}{{\hat {E}}_{M}}$, ${{E}_{{{{T}_{1}}}}} = {{L}_{1}}{{\hat {E}}_{M}}$ и ${{E}_{{{{T}_{2}}}}} = {{L}_{2}}{{\hat {E}}_{M}}$. Соответствующие моменты времени ${{\hat {T}}_{0}}$, ${{\hat {T}}_{1}}$ и ${{\hat {T}}_{2}}$ будут оценками временного положения, начала и окончания сигнала, а величина $\hat {\tau } = {{\hat {T}}_{2}} - {{\hat {T}}_{1}}$ – оценкой его длительности.

3. ОЦЕНКА ЭФФЕКТИВНОСТИ АЛГОРИТМОВ

Как уже отмечалось, оценки ${{\hat {a}}_{0}}$ и ${{\hat {a}}_{1}}$ являются несмещенными, их дисперсии $\sigma _{{{{{\hat {a}}}_{0}}}}^{2}$ и $\sigma _{{{{{\hat {a}}}_{1}}}}^{2}$ легко вычисляются аналитически:

где ${\mathbf{D}}( \cdot )$ означает вычисление дисперсии, σ² – дисперсия гауссовского шума.

Анализ эффективности оценок оставшихся параметров выполнялся методом статистических испытаний на ЭВМ по 10 000 экспериментов. Для этого на интервалах T_y и T_x в соответствии с формулой (3) моделировались выборочные векторы $\vec {y}$ и $\vec {x}$ размерностью соответственно N и M при следующих значениях параметров: ${{\sigma }^{2}} = 1$, ${{a}_{0}} = - 2\sigma $, ${{a}_{1}} = {{{{{10}}^{{ - 3}}}\sigma } \mathord{\left/ {\vphantom {{{{{10}}^{{ - 3}}}\sigma } {\Delta t}}} \right. \kern-0em} {\Delta t}}$, $\Delta t = 1$. В качестве сигналов использовались рассмотренные ранее прямоугольный, двойной экспоненциальный и косинус-квадрат с гармоническим заполнением импульсы (см. рис. 1) одинаковой длительности по основанию ${{\tau }_{0}} = 120$ отсчетов. Энергия сигналов принималась одинаковой, и ее значение зависело от задаваемого отношения сигнал/шум (ОСШ) ${{q}^{2}} = {{{{E}_{M}}} \mathord{\left/ {\vphantom {{{{E}_{M}}} {({{\tau }_{0}}{{\sigma }^{2}})}}} \right. \kern-0em} {({{\tau }_{0}}{{\sigma }^{2}})}}$.

На рис. 2 показаны полученные в результате моделирования зависимости относительного среднеквадратического отклонения (СКО) оценки $\widehat {{{\sigma }^{2}}}$ дисперсии гауссовского шума (${{\delta }_{{\widehat {{{\sigma }^{2}}}}}}$, кривая 1), оценок коэффициентов линейной составляющей фона ${{\hat {a}}_{0}}$ (${{\delta }_{{{{{\hat {a}}}_{0}}}}}$, кривая 2) и ${{\hat {a}}_{1}}$ (${{\delta }_{{{{{\hat {a}}}_{1}}}}}$, кривая 3) от объема шумовой выборки. Нормировка производилась для оценки $\widehat {{{\sigma }^{2}}}$ – к дисперсии гауссовского шума σ², для оценки ${{\hat {a}}_{0}}$ – к среднеквадратическому значению σ, для ${{\hat {a}}_{1}}$ – к параметру σ/Δt.

Как уже отмечалось, оценки ${{\hat {a}}_{0}}$ и ${{\hat {a}}_{1}}$ являются несмещенными. Их характеристики, полученные в результате моделирования, хорошо совпали с результатами аналитических расчетов. Оценка $\widehat {{{\sigma }^{2}}}$ также является практически несмещенной, величина смещения по результатам моделирования оказалась более чем на порядок меньше СКО оценки.

На рис. 3 приведены зависимости относительного СКО ${{\delta }_{{{{{\hat {E}}}_{M}}}}}$ оценки энергии сигналов трех типов от объема шумовой выборки N при фиксированном ОСШ, равном 10 дБ, объеме сигнальной выборки M = 600 и длительности сигнала τ₀ = = 120 отсчетов. Нормировка производилась к истинной энергии импульса E_M. При объеме шумовой выборки N < 200 относительное СКО оценки энергии становится больше 0.3, практически не зависит от формы сигнала и определяется в основном только объемом выборки. При больших объемах шумовой выборки различие в относительных СКО оценки энергии для рассматриваемых сигналов не превышало 0.01.

На рис. 4 представлены зависимости относительного СКО ${{\delta }_{{{{{\hat {E}}}_{M}}}}}$ оценки энергии сигналов трех типов от ОСШ при фиксированных длительности сигнала τ₀ = 120 отсчетов и объемах шумовой выборки N = 600 и 5000. Объем сигнальной выборки M = 1000, 600 и 300. Видно, что при малых ОСШ СКО оценки зависит от формы импульса, соотношения длительности импульса и длины временного интервала, на котором формировалась сигнальная выборка, а также от объема шумовой выборки. Чем меньше доля импульса в сигнальном интервале и чем меньше шумовая выборка, тем менее точными получаются оценки. При больших отношениях сигнал/шум и рассматриваемых объемах выборки эта разница для рассматриваемых трех сигналов становится незначительной.

На рис. 5а–5в приведены зависимости относительных СКО ${{\delta }_{{\hat {T}}}}$ оценок временного положения точек пересечения текущей энергией сигнала соответствующих пороговых уровней от величины этих порогов для трех типов сигналов: прямоугольный импульс, двойной экспоненциальный импульс и косинус-квадрат импульс с гармоническим заполнением. Значения СКО получены при ОСШ = 5, 10 и 15 дБ. Длительность сигнала τ₀ = 120 отсчетов и объем шумовой выборки N = 5000 оставались неизменными, объемы сигнальных выборок принимались равными M = 300, 600 и 1000. Нормировка осуществлялась к длительности импульса.

Из рисунка видно, что оценка временного положения существенно зависит от формы импульсов, при этом наименьшая ошибка оценки времени пересечения порога наблюдается на участках сигнала, имеющих наибольшую крутизну. Поэтому, если это возможно, значения порогов L₁, L₂, используемых при оценке длительности сигнала, следует выбирать таким образом, чтобы крутизна кривой текущей энергии в окрестности уровней ${{E}_{{{{T}_{1}}}}} = {{L}_{1}}{{E}_{M}}$ и ${{E}_{{{{T}_{2}}}}} = {{L}_{2}}{{E}_{M}}$ была максимальной. При ОСШ более 10 дБ для прямоугольного импульса эти значения находятся в диапазоне 0.1…0.3 и 0.7…0.9, а для косинус-квадрат импульса с гармоническим заполнением оптимальными являются пороги в диапазоне 0.2…0.4 и 0.6…0.8. При малых ОСШ (5 дБ) при пороговых уровнях больше 0.7E_M алгоритм становится неработоспособным, ошибки превышают 0.2 и только в случае, когда сигнал занимает не менее 40% от длительности сигнального интервала, относительное СКО не превышает 0.1.

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

Видно, что размер шумовой выборки существенно влияет на точность оценки при малых ОСШ, а также при больших пороговых уровнях, что имеет особое значение при оценке длительности импульса. При ОСШ = 5 дБ относительные СКО различаются при пороге 0.9 на 70…100%, при пороге 0.8 – на 20…40%, при пороге меньше 0.7 – меньше чем на 10%. Но уже при ОСШ 10 дБ объем шумовой выборки существенен только для порога 0.9. Таким образом, если при построении алгоритма оценки временного положения не представляется возможным использовать достаточно объемную шумовую выборку для оценки параметров шума, то следует выбрать пороги, не превышающие 0.7.

Для сравнения приведем характеристики алгоритмов, предложенных в работе [18] для оценки временного положения прямоугольного импульса с известной длительностью и полностью известным статистическим описанием наблюдаемых данных и в работе [19] – для импульса известной формы, наблюдаемого на фоне гауссовского шума с неизвестными параметрами. По данным работы [18], при ОСШ = 10 дБ относительное СКО оценки временного положения прямоугольного импульса составляет величину 0.425%. Это значительно меньше значения, полученного предложенным методом (1.4%), однако требует полной информации о форме сигнала и шуме. При неизвестных параметрах шума, но известной форме сигнала [19] СКО оценок временного положения имеют тот же порядок (2…8%), что и приведенные в табл. 1. Существенным преимуществом предлагаемого метода является то, что он допускает присутствие линейно изменяющейся помехи в фоновой составляющей наблюдаемого процесса и не требует априорной информации ни о параметрах фона, ни о параметрах и форме импульсного сигнала.

ЗАКЛЮЧЕНИЕ

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

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