Solving rician data analysis problems: theory and numerical modeling using computer algebra metods in Wolfram Mathematica

Full Text

1. Введение

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

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

Задача совместного оценивания сигнала и шума при двухпараметрическом анализе райсовских данных состоит в совместном расчете параметров распределения Райса. Традиционно используемый однопараметрический подход к решению задачи восстановления полезного сигнала на фоне шума предполагает априорную известность одного из райсовских параметров, а именно – дисперсии шума, и, следовательно, расчет лишь одного неизвестного параметра [2–4]. Однако условие априорной известности дисперсии шума никогда не выполняется на практике. В противоположность упрощенному однопараметрическому приближению, двухпараметрический подход не ограничен никакими априорными предположениями и тем самым обеспечивает значительно более корректное оценивание искомых величин сигнала и шума.

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

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

2. Теоретические и численные аспекты решения задачи совместного расчета райсовских параметров

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

В задачах анализа райсовского сигнала измеряемой величиной является амплитуда

$x=\sqrt{x_{\mathrm{Re}}^2+x_{\mathrm{Im}}^2}$

комплексной переменной величины с действительной x_Re и мнимой x_Im составляющими, характеризуемыми математическим ожиданием ν и искажаемыми нормально распределенным гауссовским шумом с дисперсией σ² [1–5]. Амплитуда

$x=\sqrt{x_{\mathrm{Re}}^2+x_{\mathrm{Im}}^2}$

подчиняется статистическому распределению Райса с функцией плотности вероятности, определяемой следующим выражением:

$P\left(x\left|\nu ,{\sigma }^2\right.\right)=\frac{x}{{\sigma }^2}\exp \left(-\frac{x^2+{\nu }^2}{2{\sigma }^2}\right)\cdot I_0\left(\frac{x\nu }{{\sigma }^2}\right),$ (1)

где I_a(z) — модифицированная функция Бесселя первого рода порядка a. Поставленная задача двухпараметрического анализа райсовских данных состоит в совместном расчете параметров ν и σ² распределения Райса на основе данных, полученных в выборках измерений. Именно эти параметры определяют величину исходного, неискаженного сигнала ν и дисперсию искажающего его шума σ².

На рис.1 представлена иллюстрация формирования райсовского сигнала из исходно детерминированной величины A = ν амплитуды вектора, характеризующего некоторый процесс и искажаемого под неизбежным воздействием гауссовского шума, описываемого неким вектором $\overrightarrow{r}$.

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

Рис. 1. Иллюстрация поведения функции плотности вероятности, характеризующей распределение райсовской величины x, формируемой из исходной детерминированной величины A под воздействием гауссовского шума.

Величины математического ожидания x и квадрата стандартного отклонения σ_x² даются следующими формулами:

$\begin{array}{c}\overline{x}=\sigma \sqrt{\pi /2}\cdot L_{1/2}(-{\nu }^2\text{}/2{\sigma }^2),\\ {\sigma }_x^2=2{\sigma }^2+{\nu }^2-\frac{\pi \sigma {}^2}{2}\cdot L_{1/2}(-{\nu }^2\text{}/2{\sigma }^2),\end{array}$  (2)

где L_1/2(z) — полином Лаггера. Как легко видеть, математическое ожидание райсовской величины и ее стандартное отклонение не совпадают с параметрами ν и σ², что характеризует особые, нелинейные свойства этого статистического распределения, в силу которых анализ райсовской случайной величины требует развития особых методов и соответствующего математического аппарата.

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

Ниже приведены основные выкладки, лежащие в основе данного метода. Будем использовать следующие обозначения: x_i — величина сигнала, полученная как результат i-го измерения в выборке; ν — количество элементов в выборке, называемое также длиной выборки. Рассчитанное на основе выборочных измерений $x_i\left(i=\mathrm{1,...,}n\right)$ значение ⟨x ^k⟩ стремится к значению соответствующего k-го момента $\overline{x^k}$ случайного райсовского сигнала при бесконечно большой длине ν выборки измерений:

$\overline{x^k}=\lim \frac{1}{n}\sum _{i=1}^n{x}_i^k.$

Рис. 2. Иллюстрация нелинейных свойств распределения Райса: (a) – нелинейная зависимость математического ожидания райсовской величины x от райсовского параметра ν; (б) – нелинейная зависимость квадрата стандартного отклонения райсовского сигнала от дисперсии гауссовского шума σ², формирующего райсовскую случайную величину.

Рассмотрим выборку из ν измерений величины райсовской случайной величины x. Тогда функция правдоподобия F(ν,σ²), т.е. функция совместной плотности вероятности событий, состоящих в том, что результатом i-го измерения является величина $x_i\left(i=\mathrm{1,...,}n\right),$ выражается как произведение функций плотности вероятности для каждого измерения из выборки:

$F(\nu ,{\sigma }^2)=\prod _{i=1}^{n}P\left(x_i\right|\nu ,{\sigma }^2),$  (3)

где функция P(x | ν,σ²) определяется формулой (1). При полученных в результате измерений конкретных данных выборок функция правдоподобия является функцией неизвестных статистических параметров ν и σ². Как известно, метод максимума правдоподобия состоит в определении тех значений параметров ν и σ², которые максимизируют функцию правдоподобия (3), или, что эквивалентно, ее логарифм, называемый логарифмической функцией правдоподобия распределения Райса:

$\begin{array}{c}\ln F\left(\nu ,{\sigma }^2\right)=\\ \sum _{i=1}^n\left\{\ln x_i-\ln {\sigma }^2-\frac{x_i^2+{\nu }^2}{{\sigma }^2}+\ln I_0\left(\frac{x_i\nu }{{\sigma }^2}\right)\right\}\end{array}$ . (4)

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

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

Представляемый в настоящей работе метод определения искомых параметров ν и σ² распределения Райса основан на использовании двух известных подходов к решению задач математической статистики: принципа максимума правдоподобия и метода моментов [6–8]. В качестве одного из уравнений нового комбинированного метода будем использовать уравнение метода максимума правдоподобия, полученное приравниванием нулю частной производной логарифмической функции правдоподобия по параметру сигнала ν:

$\frac{\partial }{\partial \nu }\ln F(\nu ,{\sigma }^2)=0.$

Тогда, принимая во внимание выражение (4), получаем:

$\sum _{i=1}^n\frac{\partial }{\partial \nu }\ln I_0\left(\frac{x_i\nu }{{\sigma }^2}\right)-\frac{n\cdot \nu }{{\sigma }^2}=0.$ (5)

В качестве второго уравнения комбинированного метода будем использовать выражение для первого момента райсовской случайной величины [9]:

$\overline{x}=\sigma \sqrt{\pi /2}\cdot L_{1/2}(-{\nu }^2\text{}/2{\sigma }^2).$ (6)

Формулы (5) и (6) фактически представляют собой систему двух уравнений для двух неизвестных ν и σ².

Используя известное соотношение

$\frac{d}{dz}{I}_0\left(z\right)=I_1\left(z\right)$ ,

[10] и проводя несложные преобразования, получаем следующую систему уравнений рассматриваемого метода:

$\left\{\begin{array}{l}\nu =\frac{1}{n}\sum _{i=1}^n{x}_i\cdot \tilde{I}\left(\frac{x_i\nu }{{\sigma }^2}\right),\\ \frac{1}{n}\sum _{i=1}^n{x}_i=\sigma \sqrt{\frac{\pi }{2}}\cdot L_{1/2}\left(-\frac{{\nu }^2}{2{\sigma }^2}\right).\end{array}\right.$ (7)

В первом уравнении системы (7) обозначение $\tilde{I}$ используется для функции

$\tilde{I}\left(z\right)=\frac{I_0\left(z\right)}{I_1\left(z\right)}$,

равной отношению модифицированных функций Бесселя первого и нулевого порядков. Уравнения (7) представляют собой систему двух существенно нелинейных уравнений для двух неизвестных: ν и σ².

Введем переменную r = ν²/(2σ²), которая характеризует величину отношения сигнала к шуму. Переходя от пары переменных (ν,σ) к переменным (r,σ), получим из (7) следующую систему уравнений:

$\left\{\begin{array}{l}\sigma \sqrt{2r}=\frac{1}{n}\sum _{i=1}^n{x}_i\cdot \tilde{I}\left(\frac{x_i}{\sigma }\sqrt{2r}\right),\\ 〈x〉=\sigma \sqrt{\pi /2}\cdot L_{1/2}(-{\nu }^2\text{}/2{\sigma }^2),\end{array}\right.$  (8)

где  $〈x〉=\frac{1}{n}\sum _{i=1}^n{x}_i$ — измеренное в выборках значение первого момента, или среднего значения райсовской величины x. Из второго уравнения системы (8) получаем для переменной σ следующее выражение:

$\sigma =\frac{〈x〉}{\sqrt{\frac{\pi }{2}}\cdot L_{1/2}\left(-r\right)}.$  (9)

Подставляя (9) в первое уравнение системы (8), получаем уравнение для переменной r:

$2\sqrt{\frac{r}{\pi }}\frac{〈x〉}{L{}_{1/2}\left(-r\right)}=\frac{1}{n}\sum _{i=1}^n{x}_i\tilde{I}\left(\frac{x_i\sqrt{\pi r}}{〈x〉}{L}_{1/2}\left(-r\right)\right)$  . (10)

Таким образом, задачу решения системы (8) двух уравнений для двух неизвестных (r,σ) удалось свести к задаче решения одного уравнения (10) для одной неизвестной величины r, что позволяет существенно упростить решения данной задачи посредством символьных вычислений в системе Wolfram Mathematica.

Поставленная задача расчета величины полезного сигнала и дисперсии шума состоит в том, чтобы, подставляя выборочные данные для величины сигнала $x_i\left(i=\mathrm{1,...,}n\right)$ в уравнение (10), найти решение этого уравнения для неизвестной величины r и затем рассчитать параметр шума σ, используя формулу (9). Используя полученные значения для σ и r, а также формулу r = ν²/(2σ²), нетрудно получить значение искомого параметра полезного сигнала $\nu =\sqrt{2{\sigma }^{2}r}$.

На рис. 4 представлены графики расчета райсовских параметров, полученные представленным выше методом в ходе символьных вычислений в системе Wolfram Mathematica.

Рис. 4. Результаты расчета райсовских параметров сигнала ν (а) и шума σ (б) в системе Wolfram Mathematica посредством представленного комбинированного метода.

На рис. 4 а представлена зависимость рассчитанного параметра сигнала от величины, пропорциональной отношению сигнала к шуму, при заданном значении второго параметра, причем сплошная, пунктирная и штрихпунктирная кривые соответствуют следующим значениям параметра: =0.5 (сплошная линия); (пунктирная линия); =1,5 (штрихпунктирная линия). На рис. 4 б представлены зависимости рассчитанного параметра шума при различных заданных значениях параметра сигнала, а именно: =3.0, =2.0 и =1.0, соответственно. Прямые линии на обоих графиках соответствуют исходно заданным значениям параметров. Длина выборки для проводимых вычислений составляла с усреднением по 25 выборкам (в реальных системах обработки райсовских данных число усредняемых выборок составляет, как правило, величины порядка 10³÷10⁴, что обеспечивает гораздо более высокую точность расчетов. Представленные графики демонстрируют также ожидаемое повышение точности расчетов с увеличением величины отношения сигнала к шуму.

3. Заключение

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

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

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

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

Благодарности

Автор выражает глубокую признательность Абрамову Сергею Александровичу и Рябенко Анне Андреевне за полезные обсуждения и помощь в подготовке данной работы.

About the authors

T. V. Yakovleva

Author for correspondence.
Email: tan-ya@bk.ru
ORCID iD: 0000-0003-2401-9825
Russian Federation, ul. Vavilova 44/2, Moscow, 119333

References

Supplementary files

Supplementary Files

Action

1. JATS XML

Download

2. Fig. 1. An illustration of the behavior of the probability density function characterizing the distribution of the Rice value x, formed from the initial deterministic value A under the influence of Gaussian noise.

Download (75KB)

Indexing metadata

3. Fig. 2. Illustration of the nonlinear properties of the Rice distribution: (a) – the nonlinear dependence of the mathematical expectation of the Rice value x on the Rice parameter r; (b) – the nonlinear dependence of the square of the standard deviation of the Rice signal on the dispersion of the Gaussian noise a2, forming the Rice random variable.

Download (87KB)

Indexing metadata

4. Fig. 3. Three-dimensional graphs of the likelihood function of the statistical Rice distribution, constructed in the Wolfram Mathematica system for various ratios of the values of the Rice parameters.

Download (176KB)

Indexing metadata

5. Fig. 4. The results of calculating the Rice parameters of the signal v (a) and noise a (b) in the Wolfram Mathematica system using the presented combined method.

Download (131KB)

Indexing metadata