О пеленгации источников излучений

ПРИКЛАДНАЯ МАТЕМАТИКА И МЕТОДЫ МАТЕМАТИЧЕСКОГО МОДЕЛИРОВАНИЯ

УДК 621.371

А. А. Грешило в, Б. П. Назаренко, П. А. Плохута

О ПЕЛЕНГАЦИИ ИСТОЧНИКОВ ИЗЛУЧЕНИЙ

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

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

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

Математическая модель задачи соответствует модели из работ [1,2]. Пусть набор детекторов содержит М ненаправленных детекторов, расположенных на одинаковом расстоянии d друг от друга вдоль оси х (рис. 1). Фазовый центр набора помещен в начало координат, тогда координата т-го детектора рт = (т — (М + 1) /2) d, т Е {1, ...,М}.

y1

d

Рис. 1. Иллюстрация геометрии пеленгации излучателей:

излучатели uk

Предположим, что неизвестное число пеленгуемых излучателей находится на большом расстоянии от набора детекторов. Для этих излучателей кривизна фронта волны незначительна по сравнению с расстоянием между детекторами и допустима аппроксимация волны плоским фронтом. Решение волнового уравнения для одного источника, излучающего сигнал f (t), имеет вид f (t — pTa), где p — вектор положения и а — так называемый вектор задержки, ориентированный в соответствии с направлением распространения волны; его модуль равен 1/c, т.е. величине, обратной скорости распространения сигнала. Ослабление сигнала с расстоянием не учитывается, так как оно постоянно для всего набора детекторов.

Принимаемый узкополосный сигнал может быть представлен выражением u (t) exp (ju0t), где u (t) — модулирующий сигнал; ш0 — несущая частота, которая должна быть намного выше, чем полоса пропускания u (t) для регистрируемого узкополосного сигнала. Для того чтобы избежать пространственного искажения, расстояние между детекторами не должно превышать половины длины волны, т.е. d ^ Л/2 = 2пе/(2ш0). Если не оговорено иное, выбирается d = Л/2. Выход m-го детектора равен ym(t) = u(t—Tcenter) exp(j(u0(t — —Tcenter) — kTp)), где Tcenter — задержка прохождения сигнала между излучателем и фазовым центром, волновое число задано вели-

чиной k = w0a. Предположение о том, что сигнал узкополосный, позволяет не учитывать задержку kTp в показателе экспоненты распространения сигнала u (t) между детекторами; она имеет место только при модуляции. Комплексная огибающая выхода детектора m (т.е. выход после демодуляции) может быть записана так: ym(t) = u(t - Teenier) exp(—j(woTcenter + kTp)). Если измерять время относительно момента поступления сигнала в фазовый центр, зависимость от Tcenter можно не учитывать. Таким образом, для единичного излучателя, комплексная огибающая выходов детекторов может быть записана как y (t) = a (0) u (t). Вектор a (0) = exp (—jkTp) включает в себя информацию о фазовой задержке для излучателя, пеленг которого равен 0. Параметризация вектора a (0) по 0 может быть выполнена, так как kTpm = — (w0/c) (m — (M + 1)/2) d cos 0.

Поскольку система линейна, применим принцип суперпозиции, тогда модель для K узкополосных сигналов с одинаковой несущей частотой может быть записана как y (t) = A (0) u (t). Матрица A (0) размерности M x K состоит из векторов a (0) для различных излучателей, входящих в нее в качестве столбцов: A (0) = [a(0^, a(02), ..., a(0K)]. Вектор сигнала детектора y (t) — это вектор-столбец, m-м элементом которого является ym (t), и аналогично u (t) — это вектор-столбец, состоящий из сигналов uk (t), получаемых от всех K излучателей. Вектор в состоит из пеленгов всех K излучателей: в = [ 0!, ..., 0K ]T. Принимая во внимание неизбежное присутствие помех и дискретизуя сигналы, получаем окончательный вариант модели:

y

Для простоты предполагаем, что помехи являются постоянным в пространстве и времени белым шумом, некоррелирующим с излучателями и являющимся кругообразно симметричным. Ковариационная матрица принимает следующий вид: математическое ожидание E [n (ti) nT (t2)] = a2I£ (ti — t2), где £ () — кронекеровская дельта-функция, I — единичная матрица. Кругообразная симметрия помех приводит к тому, что E [n (t1) nT (t2)] = 0.

В случае, когда излучатели находятся на небольшом расстоянии от набора детекторов, решение сферического волнового уравнения на расстоянии r от единичного источника f (t) имеет вид f (r, t) = (1/r)f (t — r/c). Кроме того, поскольку узкополосный сигнал может быть записан как f (t) = u (t) exp (jw0t), комплексная огибающая выхода набора детекторов принимает вид (1/rm)u(t — rm/c) x x exp(—j(w0rm/c)). Здесь rm — расстояние от излучателя до m-го детектора. Пусть rc — расстояние от излучателя до фазового центра массива, а pm — местоположение m-го детектора. Принимая

во внимание предположение об узкополосности сигналов и выбор начала отсчета времени, соответствующего прибытию сигнала в фазовый центр, выход т-го детектора можно записать как ут (г) = (1/гт)и(Ь)ехр(—(ш0(гт - г с)/с = а(рт)п(г). Здесь т-й компонент вектора а (рт) содержит фазу и коэффициенты затухания для излучателя, сигнал от которого приходит на детектор т. С учетом принципа суперпозиции модель для К излучателей эквивалентна (1), за исключением того, что столбцы матрицы А (9) содержат векторы а (9) для излучателей, находящихся на небольшом расстоянии от набора детекторов, вместо векторов а (9) для сильно удаленных излучателей. В случае, когда излучатели находятся на небольшом расстоянии, отклик набора детекторов зависит не только от направлений на излучатели, но и от расстояний до них. Таким образом, в пеленгации участвуют оба эти параметра.

В случае широкополосных сигналов модель для узкополосных сигналов с помощью преобразования Фурье может быть применена для каждой частоты:

у (ш) = А (9, ш) и (ш) + п (ш), ш € {ш1,..., } ,

где у (ш) и и (ш) — преобразования Фурье от у (г) и и (г) соответственно. Следует обратить внимание на то, что для узкополосных сигналов имеет место только одна матрица А (9), тогда как в случае широкополосных сигналов каждой частотной составляющей ш соответствует своя матрица А (9, ш).

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

Аналогично, если набор детекторов установлен под углом 9, сигнал на т-м детекторе форсирован или задержан на время —тт (9) (отрицательная задержка относительно фазового центра). Максимум мощности достигается ориентированием набора детекторов в направлении, вдоль которого происходит распространение сигнала, при усло-

вии отсутствия отражений. Задержки равны фазовым сдвигам, которые могут быть представлены комплексными весами w детекторов. Таким образом, выход набора детекторов можно записать как г (£) = wTy (£) = wTa (0) и (£). Для того чтобы "повернуть" набор детекторов на угол 0, веса должны быть установлены так: w = а (0). Поскольку система линейна, аналогичный подход может быть использован для нахождения суперпозиции плоских волн, приходящих с разных направлений и имеющих одинаковые несущие частоты. Классический спектр может быть представлен в виде

L

Pf

t=i

Классический алгоритм — это очень простой и робастный подход, который широко применяется на практике. Однако он имеет недостаток, связанный с ограничением релеевского разрешения [1], который может быть смягчен только увеличением числа детекторов; улучшение соотношения сигнал-шум или увеличение времени наблюдения не изменяют разрешения.

Кроме классического алгоритма, применяются Capon's (MVDR) [3], MUSIC [4], алгоритмы, основанные на методе максимума правдоподобия (ММП) [2, 5] и другие. Наиболее приемлемы алгоритмы, основанные на методе максимума правдоподобия. Они принадлежат к классу параметрических алгоритмов, где оцениваются параметры падающей волны. Помеха моделируется как стационарный гауссовский сигнал. Для ММП целью является поиск в, u (t) и -2, доставляющих максимум функции правдоподобия:

T 2

ЬПм^в, u

t=i

где в — вектор пеленгов излучателей. Логарифмическая функция правдоподобия имеет вид

1 T / 2

lDML

t=i ^ ' =i

Вычислительная сложность данного алгоритма значительно превышает сложность классического алгоритма. Преимущество алгоритмов ММП-семейства состоит в возможности различать когерентные сигналы и в обеспечении более хороших статистических свойств [6]. Кроме того, с помощью теоремы Крамера-Рао ММП позволяет получить ковариационную матрицу оценок, что в ряде случаев бывает принципиально важным [12, 14].

Задача оценки пеленга — типичная некорректная обратная задача. Теоретически сумма гармонических сигналов на одной частоте дает гармонический сигнал. Цель в задачах пеленгации — найти "слагаемые" по зарегистрированной "сумме".

Рассмотрим линейные некорректные задачи. Позже задача пеленгации излучателя будет сведена к задаче данного типа. В обратных задачах [7-9] функция от неизвестной величины, которую необходимо определить, известна. Цель — найти имеющую физический смысл обратную функцию. Имеем у = Т (х), где вектор х € х известен, а у € Т — вектор измерений. Обычно Т — это удобный для анализа непрерывный оператор, и решение прямой задачи (нахождение у при заданном х) не сталкивается с серьезными трудностями. Обратное отображение от у к х в интересующих нас задачах гораздо сложнее, что связано с недоопределенностью задачи, неоднозначностью решений или разрывной зависимостью решения от результатов измерений. Наличие любой из перечисленных проблем делает задачу некорректной.

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

у = Тх, у € См, х € , Т € См.

Положим, что задача (5) — некорректная. Классический способ ее решения — вычисление псевдообратного оператора Т^. Применяя псевдообратный оператор, находим решение с минимальным квадратом нормы. Поскольку Т — линейная функция в конечномерном пространстве, то оператор обязательно непрерывен. Однако в некоторых практических задачах число обусловленности N оператора Т может быть очень большим, что делает псевдообратный оператор неприменимым на практике, так как в этом случае добавление даже небольшой помехи к измерениям может привести к физически неприемлемому решению.

Эффективное решение некорректных обратных задач основано на методах регуляризации. Регуляризация используется для решения некорректных задач посредством объединения априорных знаний об х для стабилизации задачи и предоставления корректных и практически пригодных решений. Метод состоит в минимизации некоторой меры 31 (х) близости у и множества значений Тх и, в то же время, удовлетворении, насколько это возможно, условия близости х и априорной информации о нем (посредством минимизации некоторой меры 32 (х)). Таким образом, задача сводится к многокритериальной задаче математического программирования. Две целевые функции не

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

7

Скаляр А — это параметр регуляризации, регулирующий соотношение между 71 (х) и 72 (х). Если А = 0, то получается нерегуляризованное решение. Если А велико, решение определяется априорной информацией об х. В общем случае выбор подходящего А связан с особенностями решаемой задачи и является нетривиальным.

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

Здесь норма разности £2 есть условие близости к данным измерений 71 (х), а ||х|2 играет роль 72 (х) из соотношения (6). Значение оценки х, доставляющее минимум целевой функции (7), может быть вычислено по следующей формуле [10]:

х = ±

^ V-2 + V -

Другие алгоритмы квадратичной регуляризации основаны на соот-

K

к u/

ношении из работы [11] x = } v^. Они могут рассматриваться

л

г=1

как взвешенные псевдообращения с весами г^. Регуляризация Тихонова — это частный случай при г = ( х ). Идея всех указанных

V — + А/

алгоритмов заключается в том, что большинство больших собственных чисел оставляют неизменными и ограничивают влияние маленьких собственных чисел матрицы Т, поскольку помеха усиливается в раз в направлении собственного вектора, соответствующего наименьшему собственному числу (здесь N — отношение наибольшего собственного числа матрицы Т к наименьшему). Алгоритмы квадратичной регуляризации очень удобны во многих практических задачах, и их преимущество состоит в том, что решение имеет короткую аналитическую запись. Однако, так как обратный оператор всегда является линейной функцией измерений, имеют место некоторые ограничения. В частности, в связи с линейностью невозможно восстановить

часть x, которая принадлежит нуль-пространству T (острые всплески сигнала x). Данный факт показывает необходимость разработки более общего нелинейного вида регуляризации, который сможет обеспечить серьезное улучшение характеристик в данном отношении, сохраняя острые пики и другие важные элементы, а также обеспечить сверхразрешение при восстановлении решения [11]. Вычислительная сложность в последнем случае значительно возрастает.

Две известные не квадратичные регуляризирующие функции — это полная вариация и энтропия [9]. Полная вариация накладывает штраф на сумму вариаций сигнала J2 = ||Ax|1 = ^^ |[Ax]J, где А — дис-

i

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

Регуляризация посредством максимальной энтропии использует следующую функцию, подобную энтропии: J2 (x) = ^^ |xi| log |xi|.

i

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

Другая регуляризующая функция, которая обладает таким же эффектом некорректности, — это £1-штраф: J2 = ||x|1. Данная функция подобна полной вариации, за исключением того, что берется £1 -норма значений x вместо их производных. Полная вариация допускает некорректные скачки градиента x, тогда как £1 -штраф выделяет некорректные значения x. Более того, для снижения штрафа, накладываемого на большие по значению составляющие, несколько невыпуклых функций также нашли применение, например, £р-квазинорма с p < 1

J2 = ||x||p = у |xi|p обеспечивает более сильную концентрацию p z—^i

энергии в восстановленном решении. Свойство сильной концентрации энергии делает i1- и £р-регуляризации очень удобными для восстановления некорректных сигналов. Существуют и другие методы регуляризации [12].

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

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

Априорная некорректность удобна при поиске некорректных сигналов х. Некорректность вектора х определяется как наличие в нем небольшого числа больших по значению элементов и равенство нулю остальных элементов. Важной линейной обратной задачей, для решения которой хорошо подходит априорная некорректность, является задача представления сигнала при помощи переопределенных базисов. Их основой являются работы Маллата [15], Донохо [16], Рао [17] и других, касающиеся аппроксимации функций и оптимального выбора базиса.

В течение последнего десятилетия задаче выбора подходящего базиса для семейства сигналов было уделено много внимания и введено множество новых базисов, таких как wavelets-базис, ridgelets-базис, еигуе1е18-базис и многие другие. Несмотря на то, что посредством любого базиса, определенного на интервале в конечномерном пространстве, можно хорошо представить любой сигнал из данного пространства, тем не менее для конкретных задач некоторые базисы обладают лучшими аппроксимационными свойствами, чем какие-либо другие.

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

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

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

Для построения переопределенного базиса следует оценить число источников в зарегистрированном сигнале по алгоритмам, описанным в работах [2, 12, 13].

При использовании переопределенных базисов теряется очень важное свойство — уникальность представления. Для восстановления уникальности ищется наиболее некорректное решение среди множества возможных решений. Математически при отсутствии шума задача формулируется так: дан сигнал y £ CM и переопределенный базис T £ CM xN; необходимо найти x £ CN, такой, что y = Tx и вектор x — некорректен. Определим ||x||^ как число ненулевых элементов x. Необходимо найти минимум ||x||0 при условии y = Tx, что является сложной комбинаторной задачей. Можно показать, что при некоторых условиях, накладываемых на T и x, оптимальное решение данной проблемы может быть найдено точно при помощи решения связанной задачи: найти min ||x||p при условии y = Tx, где 0 < p ^ 1 (отдельно рассматриваются случай p =1 и общий случай 0 < p ^ 1).

Естественное расширение для практического применения — добавление белого гауссовского шума n, т.е.

y = Tx + n.

Уравнение (9) может быть решено следующим образом [18, 19]:

Априорная часть J2 (x) оказывает эффект принудительной регуляризации. На рис. 2 показано, почему £р-регуляризация при p ^ 1 лучше обеспечивает некорректное представление x. На графике (б) — значения £р-нормы для p-й степени (||x||p) для p = 0,5, p =1 и p = 2 двумерного вектора. По сравнению с фиксированной £2-нормой, т.е. всеми векторами, которые лежат на окружности с фиксированным радиусом, £р-нормы при p ^ 1 являются минимизированными на координатных осях, т.е. соответствуют тому, что некоторые коэффициенты равны нулю, в то время как другие имеют большие значения. Данное утверждение может быть обобщено на векторы больших размерностей. График (а) показывает £р-нормы для тех же значений p

5 TU 15 20 25 35 40

Рис.2. Одномерный график ||х||р прир = 50,5; 1; 2

в одномерном случае, т.е. то, что штраф на большие составляющие (большие х) уменьшается с уменьшением р. Мощные составляющие штрафуются намного меньше при .р-штрафовании с р ^ 1, чем при .^-штрафовании (регуляризация Тихонова). Этим объясняется сглаживающий эффект .2-штрафования.

Другой вывод, который можно сделать из рис.2 — это то, что .1-норма — выпуклая, тогда как при р ^ 1 .р-норма уже не выпукла. При р ^ 1 неравенство треугольника не выполняется, и более справедливо применить термин "квазинорма" вместо термина "норма". Однако термин .р-норма будет использован для любого значения р. Вычислительная процедура минимизации нескольких невыпуклых регуляризирующих функций (в частности .р) может быть упрощена при использовании полуквадратичного алгоритма регуляризации [17]. Ключевая идея состоит в введении добавочного вектора б и расширенной функции стоимости Ц (х, б), которая является квадратичной по х при фиксированном б и минимум которой Ш1пв Ц (х, б) = 7 (х) для любого х. Если Ц (х, б) легко минимизируется по б (или если существует точное решение), то расширенная функция Ц (х, б) может быть оптимизирована с приемлемой эффективностью при помощи итеративного координатного спуска.

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

возможных решений. Регуляризация — очень гибкая процедура, позволяющая использовать соответствующую априорную информацию и £р-норму при p ^ 1, которая обеспечивает некорректность.

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

Кроме представления сигнала при помощи переопределенного базиса некорректная регуляризация имеет множество приложений во многих областях, таких как статистика, анализ данных и обучение машин. Важной задачей во всех трех областях является выбор подмножества. Предположим, наблюдаемая величина y зависит от множества параметров x = [x1, ..., xN], но влияние малого подмножества параметров намного больше, чем влияние всех остальных. Для того чтобы надлежащим образом построить модель для y посредством элементов x, необходимо найти малое подмножество {x^}, которое хорошо прогнозирует y.

Сначала рассмотрим вариант задачи представления сигнала посредством переопределенного базиса без помех и опишем ее решение с использованием ^-штрафования. Хотя вариант без помех имеет малое практическое значение, он тесно связан с формулировкой задачи с учетом помех, и большинство теоретических результатов получено именно для такого варианта.

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

K

y = Е Xi&i, где x = x+ — x-, x+ > 0, x- > 0, A = [a1 ... aN].

i=1

Найдем весовые коэффициенты xi, используя доступную информацию: y и A. Так как A — переопределенный базис, то нахождение xi — это некорректная обратная задача (существует огромное множество решений). Логично найти решение, решая следующую задачу:

Поскольку | x| 00 представляет число ненулевых элементов x, эта сложная комбинаторная задача заменяется задачей

min ||x|1 при ограничении y = Ax.

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

Сначала рассмотрим численное решение для случая, когда данные вещественные. Процедура переформулирования задачи приводит к задаче минимизации ^-нормы, которая является широко известной задачей линейного программирования. Введем две переменные x+ и x-, определенные как x+ = max (x^ 0}, ж- = max {—x^ 0}, при помощи которых x может быть легко восстановлен следующим образом: x = x+ — x-. Переменные x+ и x- ограничены положительным октантом и должны удовлетворять условию ж+ж- = 0, Vi. Оно автоматически выполняется, когда рассматривается задача минимизации. Переформулировав задачу (12), получим

• V

min II x- I при ограничениях

Численное решение задачи (13) легко найти при помощи симплекс-метода или метода внутренних точек при линейном программировании.

Задачу (12) без учета помех можно представить в следующем виде:

2

С учетом помех уравнение (14) принимает вид:

22

где в — параметр регуляризации, который характеризует желаемую интенсивность помехи. Также может быть предложено несколько альтернативных формулировок задачи. Наиболее широко распространенная формулировка (она, в частности, использована в работах [15, 19]) — это штрафная форма

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

Следует обратить внимание на то, что можно заменить ограничение ||y — Ax|2 ^ в2 ограничением ||y — Ax||2 ^ в, что приводит к эквивалентной штрафной форме

Отличие выражения (18) от (16) заключается в том, что первое слагаемое имеет первую степень, а не вторую. Эта форма записи может быть представлена более наглядно при помощи математического программирования второго порядка.

Численный метод для ^-оптимизации с учетом помех

был разработан в работе [20]. Метод особенно эффективен при p < 1. Использование £р-алгоритмов с p < 1 ограничено в связи с возникающими трудностями, так как £р-норма не является выпуклой для p < 1 (и даже ее множества уровня не являются выпуклыми). Таким образом, не гарантируется достижение глобального минимума при помощи методов локальной оптимизации. Вместо использования алгоритмов глобальной оптимизации, чаще используют алгоритмы локальной оптимизации, имея в виду, что есть хорошее начальное приближение.

Использование £р-норм с p < 1 для обеспечения некорректности может быть обосновано точно так же, как и использование ^-норм. Фактически, когда p < 1, необходимые условия глобального экстремума ip-задачи без помех (min ||x||p при ограничении y = Ax) становятся эквивалентными условиям глобального экстремума £0-задачи (min ||x||0 при ограничении y = Ax), которые являются менее жесткими. Однако нахождение глобального решения £р-задачи без помех для p < 1 является более сложным, чем для случая, когда p =1 (tp--задача без помех недостаточно выпукла).

Перейдем к формулировке задачи с помехами: J (x) = = min ||y — Ax||2 + Л ||x||p. Кроме невыпуклости возникает еще одна трудность при использовании £р-целевой функции, которая заключается в том, что она не дифференцируема в 0. Фактически, целевые функции не дифференцируемы для обоих случаев: p = 1 и p < 1. Но для p =1 существует способ обойти эту трудность: в случае вещественных данных представлением x = x+ — x-, где x+ ^ 0, x- ^ 0, осуществляется переход к эквивалентной задаче линейного программирования. Для p < 1 это невозможно, и вместо этого рассматривают дифференцируемые аппроксимации £р-целевой функции; эти аппроксимации можно также использовать для £1 -целевой функции. Дифференцируемые аппроксимации обычно содержат параметр, который отвечает за компромисс между гладкостью аппроксимации и близостью к не дифференцируемой функции, для которой производится аппроксимация. Другими словами, дифференцируемая аппроксимация — это семейство функций.

Для p = 1 известна следующая аппроксимация:

N

xiii ~ Y1 р (x^

где

Р

7 — параметр. При 7 ^ 0, аппроксимация равномерно стремится к ^-норме. Однако для любого 7 > 0 она дифференцируема. Данная аппроксимация может быть переформулирована в задачу квадратичного программирования [21]. Другая аппроксимация [22], которая работает при р < 1, также как и при р = 1, имеет вид

Здесь е > 0 — параметр гладкости. Данная дифференцируемая аппроксимация хорошо подходит для поставленных целей и дает модифицированную регуляризирующую функцию

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

Для решения данной задачи оптимизации в [20] используется алгоритм полуквадратичной регуляризации [23]. Полуквадратичная регуляризация переводит задачу неквадратичной оптимизации в совокупность квадратичных задач. Для наших целей достаточно рассмотреть квазиньютоновский алгоритм [14]. Однако доказательство локальной сходимости из любой начальной точки основано на полуквадратичных решениях алгоритма.

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

N

N

J

где n — номер итерации;

H (x) 4 2AtA + ЛЛ (x);

diag {•} — диагональная матрица, 1-й диагональный элемент которой вычисляется по формуле, заключенной в скобки. Здесь важным отличием от таких методов, как гребневая регрессия [27], является зависимость Н от Л (х) — варьируемого штрафа.

Алгоритм может также интерпретироваться как квазиньютоновский алгоритм с аппроксимацией гессиана. Полный гессиан определяется выражением

2

VJ = H (x) + Л diag

p

2 . \2-p/2 + £)

Для р < 2 второе слагаемое всегда отрицательно и может привести к бесконечности гессиана. Оставляя только первое слагаемое, мы получаем положительно определенную аппроксимацию гессиана. Кроме того, для р = 2 аппроксимация становится точной.

Итеративный квазиньютоновский алгоритм имеет следующий вид:

¿-

— ßH x

>

1

vi jf

>

При выборе шага в = 1 сокращение слагаемых в правой части приводит данное выражение к виду (24). Градиент ,]£ (х) выражается формулой

VxJ£ = 2AT (Ax — y) + Л vec

pxi

1—p/2

= 2 AT Ax + Л diag

p

1—p/2

x — 2ATy = H (x) x — 2ATy.

Отсюда получаем

x(n+1) = x(n) - H (x

r

1

^ XJF.

= 2H

что приводит к выражению (24).

Продолжаем итеративный процесс в (24) до тех пор, пока

x

< s,

где 8 > 0 — малая константа. Сходимость алгоритма из любой начальной точки к локальному минимуму гарантирована [28, 29].

В рассматриваемой задаче пеленгации сигналов необходимо дополнительно учесть влияние погрешностей всех величин (кроме по-

2

2

2

2

мехи зарегистрированного сигнала), входящих в уравнение падающей волны, которые не фигурировали в предыдущем анализе:

— погрешность определения расстояния между детекторами;

— точность установки направлений;

— точность определения частоты ш0.

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

ММП позволяет получить и ковариационную матрицу (матрицу

рассеяния) Ю найденных оценок 0 по матрице М, составленной из вторых производных логарифмической функции правдоподобия 1п Ь по 0, вычисленных при найденных оценках 0, т.е. Ю = М-1, д21п Ь

где М = • ^ = •

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

п

п = / & ©) = Е

г=1

В процессе наблюдений из-за влияния погрешностей измерений 8г и е^ имеем х^ = ^ + 8^; уг = п + ег, я = 1, 2,..., п. Таким образом, оцениваются параметры линейной модели с погрешностями в матрице системы и в правой части уравнений. Будем считать, что ошибки измерений суть независимые нормально распределенные случайные величины с нулевыми средними значениями и известными дисперсиями а2(х^) и а2(у^).

В этом случае логарифмическая функция правдоподобия 1п Ь, для нашего описания это функция 71 (х), с точностью до констант будет иметь вид

\2

1 ^^J

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

Xij Cij

< 3a

Здесь для упрощения вида функционала (30) ошибки измерений считаются статистически независимыми.

Задача минимизации по в функционала (30) при фиксированных значениях £ есть обычная задача регрессионного анализа. Исследования показали, что для обеспечения минимального времени счета и большей точности лучше всего на первом шаге решать систему линейных уравнений при £ц = хц методом Гаусса с выбором максимального элемента, в дальнейшем минимизируя по в функционал (30) методом сопряженных градиентов. Учитывая вид функционала (30), задачу минимизации по в легко свести к задаче поиска минимума квадратичной формы

С(в) = ^вт Ав + атв, 2

которая решается методом сопряженных градиентов. Здесь в — вектор искомых оценок параметров; А — матрица квадратичной формы, элементы которой вычисляются по формуле

т 1

Агр = / , о~? Г хгз хрЗ; г = ~1, <2, . . . , П; р = 1, 2, . . . ,n,

°2

т

хгц — значение г-й переменной, полученное в ]-м измерении; ат — вектор, компоненты которого равны

т

T X Л У3 л о

&r — / о"/ Xrj; Г — 1, 2, . . . , П.

а

3=

дГ

Пересчет точных значений £ на основании условий —— = 0 сводится к решению т несвязанных между собой систем из п линейных уравнений вида

^ Qr Qp ^ + Cpj = xpj + Qp • yj =1 а2

E

p = 1, 2,... ,n; j = 1, 2,... ,m.

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

слагаемых функционала

«Sj )

по сравнению с предыдущей итерацией, заменяются значениями ^ с предыдущего шага.

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

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

Продемонстрируем результаты модельных расчетов оценки пеленгов для двух излучателей гармонических сигналов. Полагаем, что на детекторы падают плоские волны с известной одинаковой несущей частотой ш0. Расстояние между вибраторами равно половине длины падающей волны, помеха подчиняется нормальному закону распределения с нулевым математическим ожиданием и ковариационной матрицей Е = а21. Элементы приемной антенной решетки имеют широкую диаграмму направленности.

Следует отметить, что модель (1) является нелинейной относительно пеленгов 0, что вносит определенные трудности в процедуру пеленгации. Для линеаризации (1) применен метод расширенного базиса. Пусть |0Ь 01, ..., 0мв | — дискретная сетка пеленгов, соответствующая необходимой предельной точности пеленгации. Представим матрицу А следующим образом:

Введем вектор s (t), dim [s (t)] = N x 1, причем n-й его элемент sn (t) = uk (t), если пеленг k-го излучателя равен 0n; в противном случае sn (t) = 0. Можно записать следующую модель:

Математическая модель (33) — линейная.

В примерах выбран диапазон пеленгов от 0 до п, введены координаты расширенного базиса через 0,1 рад, т.е. 32 координаты. Другими

y (t) = As (t) + n (t), t G {1, ..., L} .

словами, определяется интенсивность сигналов, поступивших на антенну под углами от 0 до п рад с шагом 0,1 рад. Таким образом, будет зарегистрировано столько сигналов, сколько их поступит на антенну (в реальных условиях шаг в выборе координат должен определяться релеевским разрешением). Соотношение числа пеленгуемых сигналов и числа элементов антенной решетки, зависимость разрешения пеленгов и отношения амплитуд регистрируемых сигналов от числа элементов антенной решетки определяются только уровнем помехи. При гауссовских помехах со средним квадратическим отклонением до 5... 10% от амплитуды сигнала достаточно 3-5 элементов для надежной пеленгации сигналов.

Исследованы различные методы регуляризации. Они отличаются разрешающей способностью сигналов, видом решения и его дисперсией. В соответствии со своими свойствами каждый метод имеет свою область применения. Для решения задачи (33) наиболее подходящим с точки зрения разрешения сигналов и минимальных погрешностей определения амплитуд сигналов является метод £1-регуляризации, с использованием которого получены приведенные ниже результаты. Они сравнивались с возможностями широко применяемого метода пеленгации MUSIC [4].

Пример 1. Пеленгация двух излучателей с близкими пеленгами (рис. 3, 4). Пеленгация производится в момент времени t = 0. Параметры набора детекторов: уравнение выходного сигнала для m-го детектора

к

Ут = ^2 U (t)exp {j (ш0/о) (m — (M + 1)/2) d cos вк},

i=1

m = 1, 2,...,M; M = 50; d = A/2.

Параметры излучателей: ш0 = 2п/0, /0 = 400 МГц, u1 (0) = 10 мВ, u2 (0) = 7 мВ, #2 — ~ 6°. Параметры помехи: аддитивная гауссовская помеха n (t) с математическим ожиданием E [n (t)] =0 и средним квадратическим отклонением а = 3 мВ.

Доверительная вероятность приведенного на рис. 3 доверительного интервала равна 0,95.

На рис. 4 круглые маркеры — точные значения сигналов (согласно модели), квадратные маркеры — значения сигналов с помехой.

При соотношении амплитуд 10:7 и СКО шума 3 мВ MUSIC устойчиво различает пики лишь при разнице пеленгов не менее 0,5 рад.

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

Рис. 3. Результат пеленгации

Номер детектора

Рис. 4. Сигналы, поступающие с выходов детекторов антенной системы

момент времени t = 0. Параметры набора детекторов: уравнение выходного сигнала для m-го детектора

к

ym = ^ иг (t)exp {j (Vе) (m - (M+1)/2) dcos 9k}, m = 1, 2,..., M;

i=1

M = 12; d = Л/2.

Параметры излучателей: = 2nf0, f0 = 400 МГц, ui (0) = 10 мВ, и2 (0) = 1 мВ, в2 — Oi ~ 11°. Параметры помехи: аддитивная гауссов-ская помеха n (t) с математическим ожиданием E [n (t)] =0 и средним квадратическим отклонением а = 2 мВ.

Рис. 5. Результат пеленгации

RetyCO)], мВ

с < р ) с ]

о Ï i < )

<> « [] j) V / 5

[ < ] [ ) ] <  о

[ <

< с ) h

0 2 4 6 8 10 12

Номер детектора

Рис. 6. Сигналы, поступающие с выходов детекторов антенной системы

Доверительная вероятность приведенного на рис. 5 доверительного интервала равна 0,95.

На рис. 6 круглые маркеры — точные значения сигналов (согласно модели), квадратные маркеры — значения сигналов с помехой.

При соотношении амплитуд 10:1 и СКО шума 2мВ MUSIC устойчиво различает пики лишь при разнице пеленгов не менее 0,8 рад.

Выводы. Исследованы различные алгоритмы методов регуляризации, применяемые в задачах пеленгации источников излучений. Особое внимание уделено методам и ^-регуляризации, которые дают возможность надежно оценивать параметры сигналов, содержащих редкие, резкие импульсы. Предлагаемые алгоритмы позволяют найти не только точечные оценки параметров сигналов, но и ковариационную матрицу оценок (матрицу рассеяния). В примерах демонстрируются пеленги двух гармонических сигналов с постоянной амплитудой. В первом случае обеспечивается разрешение соизмеримых сигналов лучше 6°, во втором сигналы различаются в 10 раз по амплитуде. В обоих случаях пеленги получены методом ^-регуляризации с применением переопределенного базиса при высоких уровнях помехи.

СПИСОК ЛИТЕРАТУРЫ

1. Johnson D. H. and Dudgeon D. E. Array signal processing - Concepts and Techniques. Prentice Hall, 1993.

2. KrimH. and VibergM. Two decades of array signal processing research. The parametric approach // IEEE Signal Proc. Mag., July 1996. - Vol. 13, №.4. -P. 67-94.

3. C a p o n J. High resolution frequency-wavenumber spectrum analysis // Proc. IEEE. - 1969. - Vol. 57, №. 8. - P. 1408-1418.

4. S c h m i d t R. O. A signal subspace approach to multiple emitter location and spectral estimation // Ph.D. thesis, Stanford Univ., 1981.

5. Stoica P. and Sharman K. C. Maximum likelihood methods for direction-of-arrival estimation // IEEE Trans. Signal Processing. July, 1990. - Vol. 38, -№.7. - P. 1132-1143.

6. Stoica P. and NehoraiA. MUSIC, maximum likelihood, and Cramer-Rao bound // IEEE Trans. May, 1989. ASSP. - Vol. 37. - №. 5. - P. 720-741

7. E n g l H. W., Hanke M. and Neubauer A. Regularization of Inverse Problems. Kluwer, 1996.

8. N e u m a i e r A. Solving ill-conditioned and singular linear systems: a tutorial on regularization. // SIAM Review. - 1998. - Vol. 40, - №. 3. - P. 636-666.

9. H a n s e n P. C. Regularization tools: A Matlab package for analysis and solution of discrete ill-posed problems // Numer. Algorithms. - 1994. - Vol. 6. - P. 1-35.

10. Тихонов А. Н., А р с е н и н В. Я. Методы решения некорректных задач. -М.: Наука, 1979.

11. Karl W. C. Regularization in image restoration and reconstruction // Handbook of Image and Video Processing. A. Bovik. Ed. Academic Press, 2000.

12. Грешилов А. А. Некорректные задачи цифровой обработки информации и сигналов. - М.: Радио и связь, 1984.

13.

14,

15.

16.

17.

18.

19.

20.

21

22

23

24

25

26

27

28

29

30

26

StoicaP. and NehoraiA. MUSIC, maximum likelihood, and Cramer-Rao bound // IEEE Trans. ASSP. May, 1989. - Vol. 37. - №. 5. - P. 720-741.

Грешилов А. А. Математические методы принятия решений: Учеб. пособие для вузов. - М.: Изд-во МГТУ им. Н.Э. Баумана, 2006.

M a 11 a t S. A wavelet tour of signal processing. - Academic Press, 1998.

Donoho D. L., Johnstone I. M., K o c h J. C., S t e r n A. S. Maximum entropy and the nearly black object // J. R. Statist. Soc. B. - 1992. - Vol. 54, №. 1. -P. 41-81.

Gorodnitsky I. F. and Rao B. D. Sparse signal reconstruction from limited data using focus: a re-weighted minimum norm algorithm // IEEE Trans. Signal Processing. Mar. 1997. - Vol. 45, №. 3. - P. 600-616.

Chen S. S., Donoho D. L. and Saunders M. A. Atomic decomposition by basis pursuit // SIAM J. Scientific Computing. - 1998. - Vol.20, №. 1. - P. 33-61.

F u W. J. Penalized regressions: the bridge versus the LASSO // Journal of Computational and Graphical Statistics. - Sept. 1998. - Vol. 7, №. 3.

GemanD. and Yang C. Nonlinear image recovery with half-quadratic regularization // IEEE Trans. Image Processing. - July, 1995. - Vol. 4, №. 7. P. 932946.

Miller A. J. Subset Selection in Regression. - Chapman and Hall, 2002.

Tibshirani R. Regression shrinkage and selection via the LASSO // Journal of Royal Statistical Society, Series B. - Nov., 1996. - Vol. 58. - P. 267-288.

CetinM. and KarlW. C. Feature-enhanced synthetic aperture radar image formation based on nonquadratic regularization // IEEE Trans. Image Processing. -Apr. 2001. - Vol. 10. - №.4. - P. 623-631.

Mangasarian O. L. and Musicant D. R. Robust linear and support vector regression // IEEE Trans. Pattern Anal. Machine Intell. - 2000. - Vol. 22, №. 9. -P. 950-955.

Vogel C. R. and Oman M. E. Fast, robust total variation-based reconstruction of noisy, blurred images // IEEE Trans. Image Processing. Jun, 1998. -Vol. 7, №.6. - P. 813-824.

Geman D. and Reynolds G. Constrained restoration and the recovery of discontinuities // IEEE Trans. Pattern Anal. Machine Intell. Mar., 1992. - Vol. 14, №.3. P. 367-383.

Hoerl A. E. and Kennard R. W. Ridge regression: Biased estimation for nonorthogonal problems // Technometrics. - 1970. - Vol. 12, №. 3. - P. 55-67.

C h a r b o n n i e r P., B l a n c - F e r a u d L., A u b e r t G. and BarlaudM. Deterministic edge-preserving regularization in computed imaging // IEEE Trans. Image Processing. - Feb. 1997. - Vol. 6, №. 2. - P. 298-310.

Delaney A. H. and Bresler Y. Globally convergent edge-preserving regularized reconstruction: an application to limited-angle tomography // IEEE Trans. Image Processing. - Feb., 1998. - Vol. 7, №.2. - P. 204-221.

Malioutov D. M. A sparse signal reconstruction perspective for source localization with sensor arrays // Master of Science thesis, Massachusetts Institute of Technology, 2003.

Статья поступила в редакцию 19.01.2007

Анатолий Антонович Грешилов родился в 1939г., окончил в 1964 г. Московский инженерно-физический институт. Д-р техн. наук, профессор кафедры "Вычислительная математика и математическая физика" МГТУ им. Н.Э. Баумана. Автор более 150 научных работ, в том числе 19 монографий, 20 авторских свидетельств и патентов в области разработки математических методов строгого учета неопределенности исходной информации в задачах математической физики, распознавания образов, прогнозирования и других технических приложений.

АА. Greshilov (b. 1939) graduated from the Moscow Institute for Engineering and Physics in 1964. D. Sc. (Eng.) professor of "Computational Mathematics and Mathematical Physics" department of the Bauman Moscow State Technical University. Author of more than 150 publications including 19 monographs, 20 author's certificates and patents in the field of development of mathematical methods for strict account of source data uncertainty in problems of mathematical physics, image identification, forecast and other technical applications.

Борис Петрович Назаренко родился в 1943 г., окончил в 1966 г. МВТУ им. Н.Э. Баумана. Канд. техн. наук, руководитель научно-учебного комплекса, декан факультета "Фундаментальные науки" МГТУ им. Н.Э. Баумана. Автор более 70 научных работ, в том числе трех монографий и трех авторских свидетельств в области математического моделирования транспортных систем, теоретической механики и технических приложений информатики.

B.P. Nazarenko (b. 1943) graduated from the Bauman Moscow Higher Technical School in 1966. Ph. D. (Eng.), chief of scientific and educational complex, dean of "Fundamental Sciences" faculty of the Bauman Moscow State Technical University. Author of more than 70 publications including 3 monographs and 3 author's certificates in the field of mathematical simulation of transport systems, theoretical mechanics and technical applications of information technology.

Павел Анатольевич Плохута родился в 1983 г., окончил МГТУ им. Н.Э. Баумана в 2006 г. Аспирант кафедры "Вычислительная математика и математическая физика" МГТУ им. Н.Э. Баумана.

P.A. Plokhuta (b. 1983) graduated from the Bauman Moscow State Technical University in 2006. Post-graduate of "Computational Mathematics and Mathematical Physics" department of the Bauman Moscow State Technical University.