Тема 11: АДАПТИВНАЯ ФИЛЬТРАЦИЯ ДАННЫХ

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

h_oR(M) + h₁R(M-1)+ h₂R(M-2)+ ....... + h_MR(0) = 0.

При расчете коэффициентов фильтра значение s(0) обычно принимается равным 1.

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

10.5. Фильтры прогнозирования.

Если в правой части уравнения (10.3.3) желаемым сигналом задать входной сигнал со сдвигом на величину kDt, то при этом B(m) = R(m+k) и уравнение принимает вид:

h(n) _* R(m-n) = R(m+k). (10.5.1)

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

10.6. Фильтр обнаружения сигналов.

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

Частотная характеристика. Для расчета фильтра требуется задать известную форму полезного сигнала s(k) ó S(w) и функцию автокорреляции или спектр мощности помех R_q(m) ó W_q(w). Полный входной сигнал принимается по аддитивной модели: x(t) = s(t)+q(t). На выходе проектируемого фильтра h(n) ó H(w) для составляющих выходного сигнала имеем:

y(t) = H(w) S(w) exp(jwt) dw, (10.6.1)

s² = |H(w)|² W_q(w) dw, (10.6.2)

где s - средняя квадратическая амплитуда выходной помехи. Значения (10.6.1, 10.6.2) используются для задания критерия максимума отношения сигнал/шум (10.2.3) для произвольной точки t_i:

r = [y(t_i)]²/d². (10.6.3)

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

H(w) = exp(-jwt_i) S*(w) / W_q(w), (10.6.4)

Без потери общности можно принять t_i=0:

H(w) = S*(w)/W_q(w) = |S(w)|exp(jj_s(w)) / W_q(w). (10.6.5)

При переходе во временную (координатную) область:

H(w)W_q(w) = S*(w) ó h(n) _* R_q(n-m) = s(-m). (10.6.6)

Система линейных уравнений для расчета фильтра:

h_oR_q(0)+ h₁R_q(1)+ h₂R_q(2)+ h₃R_q(3)+ ...+ h_MR_q(M) = S(-M),

h_oR_q(1)+ h₁R_q(0)+ h₂R_q(1)+ h₃R_q(2)+ ...+ h_MR_q(M-1)= S(-M+1),

h_oR_q(2)+ h₁R_q(1)+ h₂R_q(0)+ h₃R_q(1)+ ...+ h_MR_q(M-2)= S(-M+2),

h_oR_q(M)+ h₁R_q(M-1)+ h₂R_q(M-2)+ ..... + h_MR_q(0) = S(0).

Эффективность фильтра. Из выражения (10.6.5) можно видеть, что фильтр имеет максимальный коэффициент передачи на частотах доминирования сигнала и минимальный коэффициент передачи на частотах доминирования помех. Кроме того, фазовая характеристика фильтра j(w) = -j_s(w) обеспечивает синфазность всех частотных составляющих выходного сигнала и соответственно максимальную его амплитуду в заданный момент времени t_i = 0:

y(0)ó S(w) H(w) = = .

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

Согласованный фильтр. При помехах типа белого шума W_q(w) = s² и H(w) = S*(w)/s². Множитель 1/s² не влияет на отношение сигнал/помеха и может быть опущен. Частотная характеристика фильтра определяется только спектром сигнала, при этом:

h(n) = s(-n). (10.6.7)

Фильтр получил название согласованного (по частотной характеристике со спектром сигнала). Он мало эффективен при коротком импульсном или длинном гармоническом сигнале.

Обратный фильтр. Допустим, что помеха имеет такой же частотный состав, что и полезный сигнал, т.е.:

W_q = s²|S(w)|².

Выделение полезного сигнала в таких условиях весьма сомнительно. Тем не менее, определим оптимальный фильтр:

H(w) = S*(w) / [s²|S(w)|²] = 1 / [s²S(w)]. (10.6.8)

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

W_q = a²|S(w)|²+b².

Подставляя это выражение в (10.6.5), с точностью до множителя получаем:

H(w) = S*(w) / [|S(w)|²+g²], (10.6.9)

где g = b/a - отношение дисперсий шума и сигнала. Фильтр стремится к согласованному при больших g, и к обратному (фильтру сжатия) при малых.

10.7. Энергетический фильтр.

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

Критерий оптимальности. Энергия сигнала на выходе фильтра:

E_sh = S_k s_k² = S_k (S_n h_n s_k-n)² = S_k h_k S_n h_n R_s(k-n), (10.7.1)

где R_s- функция автокорреляции сигнала. В векторной форме:

E_sh = . (10.7.2)

Аналогично, выражение для энергии помех на выходе:

E_qh = S_k h_k S_n h_n R_q(k-n) = , (10.7.3)

где R_q - функция автокорреляции помех. При некоррелированной помехе E_qh = s².

Подставим (10.7.2, 10.7.3) в выражение (10.2.4):

r= / . (10.7.4)

Расчет векторов операторов фильтров. Для определения значений вектора продифференцируем r по и приравняем производную к нулю:

.

. (10.7.5)

В системе уравнений (10.7.5) неизвестны собственные значения r матрицы и значения коэффициентов h_n, при этом система имеет N+1 ненулевых решений относительно значений r и соответствующих этим значениям векторов . Для определения коэффициентов фильтра приравнивается к нулю и решается относительно r определитель матрицы , после чего максимальное значение r_max подставляется в (10.7.5) и система уравнений решается относительно коэффициентов h_i вектора . При фильтрации сигнала вектор обеспечивает выделение первой по мощности главной компоненты сигнала, т.е. составляющей сигнала, которая имеет наибольшую энергию и отношение сигнал/шум. В сложных геофизических полях такая компонента, как правило, соответствует региональному фону.

В принципе, расчет может быть продолжен и для других значений r<r_max и определены значения коэффициентов векторов , и т.д., с использованием которых могут выделяться вторая и прочие компоненты сигнала. Наиболее эффективно такой метод используется для разделения сигналов (полей) при некоррелированных помехах. В этом случае корреляционная матрица помех является единичной (единицы по диагонали, остальное - нули) и уравнение (10.7.5) имеет вид:

. (10.7.6)

В развернутой форме:

h_o(R_s(0)-r)+ h₁R_s(1)+ h₂R_s(2)+ h₃R_s(3)+ ...+ h_MR_s(M) = 0,

h_oR_s(1)+ h1(Rs(0)-r)+ h₂R_s(1)+ h₃R_s(2)+ ...+ h_MR_s(M-1) = 0,

h_oR_s(2)+ h₁R_s(1)+ h₂ (Rs(0)-r)+ h₃R_s(1)+ ...+ h_MR_s(M-2) = 0,

h_oR_s(M)+ h₁R_s(M-1)+ h₂R_s(M-2)+ ..... + h_M (Rs(0)-r) = 0.

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

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

Пусть они постараются подчинить себе обстоятельства, а не подчиняются им сами.

Гораций. Послания.

(................)

Если в этой теории Вы не увидите смысла, тем лучше. Можно пропустить объяснения и сразу приступить к ее практическому использованию.

Валентин Ровинский. Теория карточных игр.

(Киевский геофизик Уральской школы).

Содержание: 11.1. Введение. 11.2. Основы статистической группировки информации. Предпосылки метода. Задача статистической группировки. Использование априорных данных. Эффективность метода. 11.3. Статистическая регуляризация данных. Проверка теоретических положений метода. Оценка СРД. Результаты моделирования. 11.4. Статистическая группировка полезной информации. Сущность аппаратной реализации. Реализация систем СГПИ. Пример технического исполнения системы СГПИ. Литература.

11.1. Введение.

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

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

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

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

11.2. Основы статистической группировки информации.

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

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

В общем случае полезная (целевая) информация может присутствовать в нескольких энергетических интервалах спектра излучения. Рабочими интервалами измерений обычно счи­таются участки спектра, где полезная информация присутствует в "чистом" виде либо в смеси с помехами (фоном), значение которых может быть учтено при обработке результатов изме­рений. Так, например, при гамма-опробовании пород на содержание естественных радионук­лидов (ЕРН) регистрируется излучение с энергией более 250-300 кэВ, представленное в основном первичными и однократно рассеянными квантами, плотность потока которых про­порциональна массовой доле ЕРН в породах. Плотность потока излучения в низкоэнергети­ческом интервале спектра (20-250 кэВ, в основном многократно рассеянное излучение) также зависит от массовой доли ЕРН, но эта зависимость является параметрически связанной с эф­фективным атомным номером излучающе-поглощающей среды в области детектора, вариации которого по стволу скважины могут приводить к большой погрешности интерпретации ре­зультатов измерений. Между тем плотность потока полезной информации (относительно массовой доли ЕРН) в интервале 20-250 кэВ много выше, чем в интервале более 250 кэВ, осо­бенно при регистрации излучения сцинтилляционными детекторами малых объемов, которые имеют повышенную чувствительность именно к низкоэнергетической части спектра излуче­ния.

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

Обозначим потоки, а равно и частоты основного и дополнительного потоков сигналов индексами n и m (импульсов в секунду), связь потоков по частотам индексом х = m/n. Определению подлежит частота потока n. Значение х может изменяться за счет влияния дестабилизирующих факторов на поток m и в общем случае представляет собой случайную величину, распределенную по определенному закону с плотностью вероятностей Р(х), мате­матическим ожиданием , и дисперсией D_x.

На основе теоремы Байеса, плотность вероятностей распределения частоты n по измеренному за единичный интервал t числу отсчетов сигнала N определяется выражением:

P_N(n) = P(n) P_n(N) /P(N), (11.2.1)

где: P_n(N) = (nТ)^N e^-nt /N! , (11.2.2)

P(N) =P_n(N) P(n) dn, (11.2.3)

здесь P(n)- априорная плотность вероятностей частоты n, P_n(N)- апостериорное распределе­ние вероятностей числовых отсчетов N (закон Пуассона). Принимая в дальнейшем в качестве искомой величины значения отсчетов z=nt по интервалам t (экспозиция цифровых отсчетов или скользящее временное окно аналоговых данных) и подставляя (2) и (3) в (1), получаем:

P_N(z) = P(z) z^N e^-z /P(z) z^N e^-z dz. (11.2.4)

При неизвестном распределении значений z априорная плотность распределения P(z) принимается равномерной от 0 до ¥, при этом из выражения (11.2.4) следуют общеизвестные выражения:

z = D_z = N+1 @ N, (11.2.5)

d_z² = D_z /z² = 1 /(N+1) @ 1/N, (11.2.6)

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

Как следует из теории гамма-каротажа (ГК) и достаточно хорошо подтверждено практикой гамма-опробования, пространственная разрешающая способность гамма-каротажных измерений при интерпретации результатов ГК на содержание естественных радиоактивных элементов в породах по стволу скважин в среднем составляет 10 см, а в скважинах малого диаметра может даже повышаться до 5-7 см. Однако реализация такой разрешающей способности возможна только в условиях достаточно "хорошей" статистики. Коэффициент усиления дисперсии помех цифровых фильтров деконволюции, которые используются при интерпретации ГК, в среднем порядка 12 и изменяется от 4 до 25 в зависимости от плотности пород, диаметра скважин, диаметра скважинных приборов и пр. Отсюда следует, что для достижения разрешающей способности в 10 см при нормативной погрешности дифференциальной интерпретации не более 10-20 % статистическая погрешность измерений не должна превышать 3-7 %. А это, в свою очередь, определяет объем отсчета за единичную экспозицию не менее 200-1000 импульсов. При гамма-каротаже последнее возможно только для пород с относительно высоким содержанием ЕРН (более 0.001 % эквивалентного урана), при использовании детекторов больших размеров (с эффективностью регистрации более 10 имп/сек на 1 мкР/час) и при низкой скорости каротажа (не более 100-300 м/час). В той или иной мере эта проблема характерна для всех методов ядерной геофизики, и особенно остро стоить в спектрометрических модификациях измерений.

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

Простейшим способом подготовки цифровых данных для интерпретации является их низкочастотная фильтрация методом наименьших квадратов (МНК) или весовыми функциями (Лапласа-Гаусса, Кайзера-Бесселя и др.). Однако любые методы низкочастотной фильтрации данных снижают пространственную разрешающую способность интерпретации, так как кроме снижения статистических флюктуаций приводят к определенной деформации частотных составляющих полезной части сигнала, спектр которого по условиям деконволюции должен иметь вещественные значения вплоть до частоты Найквиста. В определенной мере ликвидировать этот негативный фактор позволяет метод адаптивной регуляризации данных (АРД).

Выражения (11.2.5-6) получены в предположении полной неизвестности априорного распределения P(z)для отсчетов в каждой текущей экспозиции t. Между тем, при обработке данных непрерывных измерений, и тем более каротажных данных, которые обычно являются многопараметровыми, для каждого текущего отсчета при обработке данных может проводиться определенная оценка распределения P(z). Как минимум, можно выделить два способа оценки распределения P(z).

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

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

Использование априорных данных.Допустим, что кроме основного массива данных N,подлежащего обработке (подготовке к интерпретации), мы располагаем дополнительным массивом данных M, значения которого в определенной степени коррелированы с массивом N. При отсутствии дополнительных массивов способ 2 позволяет получить массив М обработкой массива N цифровым фильтром МНК (или любым другим весовым фильтром) со скользящим временным окном T = kt при к ³ 3 (M_i = m_itсглаженного сигнала m_i = n_{i *} h_k, где h_k – оператор симметричного цифрового фильтра с k – окном). Отметим также, что 2-ой способ всегда может использоваться для регуляризации данных независимо от наличия данных для 1-го метода.

Массив М позволяют дать оценку статистических характеристик распределения P(z). Так, если для тех же интервалов времени t в массиве М имеются отсчеты М = m_it (или приведенные к ним отсчеты какого-либо другого параметра), то можно записать:

P_M(z) =, (11.2.7)

где Р(х) – априорная плотность распределения значений x_i = m_i/n_i, которые в общем случае также могут быть случайными. При равномерном распределении Р(х) от 0 до ¥ для отсчета М равновероятно любое значение z, т.е. эффект от измерений в потоке m отсутствует. Однако по исходным условиям задачи в потоке m обязательно присутствие полезной информации, а, следовательно, и суще­ствование, как минимум, определенных границ распределения Р(х) от х_min > 0 до x_max << ¥, и среднего значения по пространству измерений. При этом из выражения (11.2.7) следует, что наиболее вероятное значение z_a, "априорное" для отсчетов z=nt в потоке n по измерениям в потоке m (отсчетам М), должно быть равно:

z_a = (M+1)/@ М/. (11.2.8)

При статистической независимости величин х и М относительная средняя квадратическая погрешность определения значений z_a по отсчетам в массиве М:

d_za² = d_M² + d_x². (11.2.9)

Отсюда дисперсия распределения значений z_a:

D_za = (D_M+M²d_x²)/² = D(M) /², (11.2.10)

D(M) = D_M+M²d_x² = D_M+D_xm , (11.2.11)

D_M = М+1 @ М, D_xm = M²d_x²,

где значение дисперсии D_M определяется статистикой отсчетов в массиве М при х = const, значение D_xm представляет собой дисперсию значений М за счет флюктуаций величины х, а сумма D(M) определяет полную дисперсию отсчетов М.

Влияние Р(х) на форму распределения Р_М(z) сказывается в его "растягивании" по координате z относительно модального значения, при этом решение интеграла (11.2.7) в первом приближении может быть представлено в следующем виде:

P_M(z) @ be^-bz. (11.2.12)

Для данного распределения:

= z_a = a/b, (11.2.13)

D_za = a/b², (11.2.14)

С учетом выражений (11.2.8) и (11.2.10):

a = MD_M/(D_za²) = MD_M/D(M), (11.2.15)

b = D_M/(D_za) = D_М/D(M). (11.2.16)

Значение 'а' в выражении (11.2.15) принимается целочисленным. Выражение (11.2.12) может быть принято для распределения (11.2.4) в качестве априорного распределения вероятностей Р(z), при этом:

P_N(z) = (b+1)e^-z(b+1). (11.2.17)

Отсюда, математическое ожидание и дисперсия z:

z = (N+a)/(b+1), (11.2.18)

D_z = (N+a)/(b+1)². (11.2.19)

C использованием выражений (11.2.15-16):

z = bN+(1-b)M/, (11.2.20)

где b и (1-b) – весовые коэффициенты доверия отсчетам N и M:

b = D(M)/(D_N²+D(M)). (11.2.21)

Дисперсия и относительная средняя квадратическая погрешность отсчетов z:

D_z = D(M), (11.2.22)

d_z² =1/(N+MD_M/D(M)). (11.2.23)

Эффективность метода. Сравнение выражений (11.2.20-23) и (11.2.5-6) позволяет дать оценку эффекта использования дополнительной информации из статистически независимого от N потока М (произвольная дополнительная информация).

1. При Þ const имеет место d_х² Þ 0, D_xm Þ 0 и дисперсия отсчетов в массиве М определяется только статистикой потока:

D(M) Þ D_M = M, z = (N+M) /(+1),

d_z² Þ 1/(N+M) < d_N² = 1/N, (11.2.24)

h = d_N² /d_z² = [N+M²/D(M)] /N Þ 1+M/N,

что соответствует определению z по двум независимым измерениям и эффект использования дополнительной информации максимален. Так, при M » N, h Þ 2 и погрешность измерений уменьшается в ~1.4 раза.

2. В общем случае D_xm ¹ 0, при этом D(M) > D_М и положительный эффект снижается. В пределе: d_x Þ ¥, D_xm Þ ¥, D(M) Þ ¥, h Þ 1, z Þ N, d_z Þ d_N и положительный эффект полностью вырождается. Во всех остальных случаях h > 1 и d_z < d_N. Отсюда следует, что при наличии коррелированной информации в массиве М положительный эффект, в той или иной мере, всегда имеет место.

3. Положительный эффект тем больше, чем больше значение x = m/n, меньше флюктуации х (величина d_х), и меньше значения отсчетов N = nt. Положительный эффект увеличивается именно в тех случаях, когда особенно остро ощущается недостаток информации: при малых значениях плотности потока излучения и/или экспозиции измерений.

Аналогичный эффект будет иметь место и при формировании отсчетов M_i по окрестностям текущих точек обработки данных путем определения их среднего значения (низкочастотное сглаживание массива n). Предварительное низкочастотное сглаживание может применяться и для статистически независимого дополнительного массива m, что будет повышать достоверность прогнозных отсчетов и увеличивать глубину регуляризации, если это сглаживание при регуляризации по формулам (11.2.20 и 21) не сказывается на изменении формы основного сигнала. Последнее определяется соотношением частотных спектров основного сигнала и оператора сглаживания.

Возможны два способа реализации уравнения (11.2.20): непосредственно в процессе измерений методом статистической группировки полезной информации (СГПИ) в реальном масштабе времени, или методом статистической регуляризации данных (СРД), зарегистрированных в виде временного (пространственного) распределения в параллельных массивах отсчетов.

11.3. Статистическая регуляризация данных.

Как следует из выражения (11.2.21), для практического использования информации из дополнительных потоков данных необходимо установить значения и диспер­сию D(M), причем, исходя из задания последней по выражению (11.2.11), должно быть известно значение d_x - относительной средней квадратической флюктуации вели­чины х.

Применительно к СРД определение значений и d_x по зарегистрированным мас­сивам данных не представляет затруднений как в целом по пространству измерений, так и в виде распределений в скользящем окне усреднения данных. Последнее эквивалентно приведе­нию D_xm => 0 для текущей точки обработки данных по информации ее ближайших окрестно­стей и позволяет производить максимальное извлечение полезной информации из дополни­тельных потоков сигналов, если частотный спектр распределения величины х по пространству измерений много меньше частотного спектра полезного сигнала. Отметим, что информация о распределении х также может иметь практическое значение (в частности, при гамма-опробовании с дополнительным потоком сигналов в низкоэнергетическом диапазоне спектра излучения - для оценки эффективного атомного номера горных пород).

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

В таблице 1 приведены 4 группы результатов обработки по формулам (11.2.20-21) двух статистически независимых и постоянных по средним значениям массивов данных n и m (модели постоянных полей) при различных установках СРД по скользящему окну К_с счета текущих значений = m_i/n_i и D_i(М) по массиву m. Текущая точка обработки данных – по центру окна. Количество отсчетов в каждом массиве – 1000, распределение значений отсчетов соответствует закону Пуассона. Определение прогнозных отсчетов М_i по массиву m для использования в уравнении (11.2.20) проводилось со сглаживанием отсчетов в скользящем окне K_s низкочастотного цифрового фильтра (вариант без сглаживания при K_s = 1). В качестве низкочастотного фильтра в алгоритме СРД используется (здесь и в дальнейшем) весовое окно Лапласа-Гаусса. Теоретическое значение D_z.т. дисперсии результатов z определялось по выражению (11.2.22) с расчетом дисперсии D(M) по выражению D(M) =[1+(1/(K_c)+1/(K_c))]. При сглаживании прогнозных отсчетов значение D_M в выражении (11.2.22) принималось равным D_M. = ×H_s, где H_s – коэффициент усиления сглаживающим фильтром дисперсии шумов (сумма квадратов коэффициентов цифрового фильтра). Дополнительно в таблице приводятся зарегистрированные средние значения коэффициента снижения статистических флюктуаций h = d_n²/d_z².

Таблица 1.

Статистика результатов моделирования СРД.

(Основной массив = 9.9, D_n = 9.7, дополнительный массив = 9.9, D_m = 9.9, 1000 отсчетов.)

Как видно из данных таблицы, практические результаты фильтрации достаточно хорошо совпадают с ожидаемыми по данным теоретических расчетов. Некоторое уменьшение среднего значения z по отношению к исходному среднему значению n определяется асимметричностью пуассоновского типа модели. При малых средних значениях модельных отсчетов в массиве m это приводит к определенной статистической асимметрии в работе СРД, т.к. при (+s_m)² > (-s_m)² среднестатистическое доверие к дополнительной информации с отсчетами M_i+s меньше, чем с отсчетами M_i-s. Этим же фактором, по-видимому, вызвано и большее расхождение между теоретическими и фактическими значениями D_z при малых значениях окна К_с. Можно также заметить, что по значению коэффициента h фильтрация выходит на теоретические значения (Þ 1+M/N) только при достаточно точном определении значений и D_i(М), что требует увеличения окна К_с счета этих параметров для полного использования дополнительной информации.

Таблица 2.

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

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

При установке СРД без усреднения данных по массиву М (К_s = 1, прогноз М_i по текущим значениям массива М) при любых значениях окна К_с выходной массив Z без всяких изменений повторяет массив N, т.е. не изменяет полезный сигнал и полностью сохраняет его частотные характеристики. Естественно при условии, что массив М пропорционален массиву N.

При К_s > 1 форма выходных кривых несколько изменяется и приведена на рис. 11.3.1. В индексах выходных кривых z приведена информация по установкам окон СРД: первая цифра - окно счета дисперсии D_M и текущего значения (в количестве точек отсчетов), вторая цифра (через флеш) - окно сглаживания отсчетов М весовой функцией Лапласа-Гаусса и определения прогнозных отсчетов М_i. Для сравнения с результатами типовой низкочастотной фильтрации на рисунке приведена кривая n25 отсчетов N, сглаженных весовой функцией Лапласа-Гаусса с окном 25 точек.

Рис. 11.3.1. СРД прямоугольного импульса. Счет D_m по несглаженному массиву М.

На рис. 11.3.1а приведен результат СРД прямоугольного импульса с амплитудным значением 10 на фоне 5 при отношении m/n = 1 (равные значения отсчетов N и М). Дисперсия D_N в выражении (11.2.21) принималась равной значению отсчетов N (статистика Пуассона). Как видно на рисунке, при сохранении фронтов сигнальной функции сглаживание прогнозных значений М_i приводит к появлению искажения формы сигнала по обеим сторонам скачка, интервал которого тем больше, чем больше значение K_s. Амплитудное значение искажений, как это и следует из выражения (11.2.21), в первую очередь зависит от соотношения текущих значений D_N и D(M) и в меньшей степени от глубины сглаживания прогнозных отсчетов.

Максимальную величину искажения для точек скачка в первом приближении можно оценить из следующих соображений. Значения D(M) между точками скачка равны D(M) = А²/4, где А - амплитуда скачка, при этом значения коэффициента b для нижней и верхней точек скачка определяются выражениями b » А²/(4D_N+A²), где D_N = N точки скачка (для статистики Пуассона). Отсюда, при прогнозном значении М » N+А/2 для нижней точки скачка и M » N-A/2 для верхней точки относительная величина изменений N определится выражением d » 1/(2N/A+A), т.е. будет тем меньше, чем больше значения А и N и больше отношение N/A, что можно наглядно видеть на рис. 11.3.1в. Из этого выражения также следует, что максимальные искажения скачков, вносимые системой СРД, будут всегда в несколько раз меньше, чем статистические флюктуации непосредственных отсчетов d = 1/на краях скачков.

При увеличении глубины регуляризации введением счета дисперсии D(M) по сглаженному массиву М картина искажений несколько изменяется и приведена на рис. 11.3.2. Реакция СРД на сглаживание дисперсии D(M) проявляется в своеобразной компенсации абсолютных отклонений отсчетов непосредственно по сторонам скачка отклонениями противоположного знака в более дальней зоне от скачка. Максимальные значения искажений остаются примерно на таком же уровне, как и для работы по несглаженной дисперсии D(M), с несколько меньшей зависимостью от увеличения значений N и А.

Рис. 11.3.2. СРД прямоугольного импульса. Счет D_m по сглаженному массиву М.

В приведенных примерах значение окна счета К_с принималось равным значению окна сглаживания К_s дополнительного массива М. При К_с > K_s картина процесса практически не изменяется. При обратном соотношении размеров окон вступает в действие второй фактор - отклонение от фактических значений счета текущих значений x_i = m/n в малом окне К_с по массиву отсчетов, сглаженных с большим окном K_s. На расстояниях от скачка функции, больших К_с/2, СРД переходит в режим предпочтения сглаженных значений массива М, т.к. D(M) Þ 0, что при К_с < K_s может приводить к появлению существенной погрешности – выбросов на расстояниях ± К_с/2 от скачков. Естественно, что при практических измерениях таких условий наблюдаться не будет и эффект резко уменьшится, но для полного его исключения вариант K_c ³ K_s можно считать предпочтительным.

Рис. 11.3.3. СРД сигнала N по массиву M. Рис. 11.3.4. Коэффициент b. (Счет Dm по несглаженному массиву М). (Среднее статистическое по 50 циклам)

На рис. 11.3.3 приведен пример регистрации рандомизированного модельного сигнала в виде прямоугольного импульса амплитудой 40 на фоне 10, на котором виден принцип работы СРД. Как и следовало ожидать, СРД производит сглаживание статистических флюктуаций фона и сигнала за пределами зоны ±К_с от скачка, отдавая предпочтение сглаженным прогнозным значениям М_i, и не изменяет значения фона и сигнала в пределах этой зоны в связи с резким возрастанием текущих значений D(M) в выражении (11.2.21). Изменение коэффициента b в зоне скачка, управляющего формированием выходных отсчетов, приведено на рис. 11.3.4 (среднестатистическое по 50-ти циклам рандомизации для модельного импульса на рис. 11.3.3) и наглядно показывает принцип адаптации СРД к динамике изменения значений обрабатываемых сигналов.

Рис. 11.3.5. Статистика сигнала N Рис. 11.3.6. Статистика сигнала Z (Измерения по 50-ти циклам). (50 циклов. Счет Dm по несглаженному М)

Статистическая оценка работы СРД по прямоугольным импульсам проводилась по 50-ти циклам рандомизации исходных массивов N и M. В качестве примера на рисунках 11.3.5 и 6 приведены результаты обработки статистики массивов N и Z. Кроме статистики циклов рандомизации проводилась суммарная обработка всех циклов по общей статистике фона и вершины импульсов. Результаты обработки для тех же установок фильтров приведены в таблице 3.

Таблица 3.

Статистика значений фона и вершины импульсов (50 циклов).

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

При отсутствии дополнительных массивов М, коррелированных с регуляризируемым массивом N, формирование прогнозных значений М_i может производиться по ближайшим окрестностям текущих значений N_i в скользящем окне K_s. При строго корректном подходе текущая точка N_i не должна включаться в число счета прогнозных значений M_i, но, как показало моделирование, это практически не влияет на результаты регуляризации. При прогнозировании M_i по всем точкам окна K_s массив М формируется любым методом сглаживания из массива N, и все особенности работы СРД по сглаженным массивам М, рассмотренные выше, остаются без изменений при условии счета значений D_m в окне К_с по массиву М. Для исключения выбросов по обе стороны от скачков полезного сигнала счет D_m как дисперсии прогнозных значений M_i необходимо выполнять непосредственно по массиву N.

Рис. 11.3.7. СРД одиночного массива N (3 прохода. Счет Dm по массиву n)

Фундаментальной особенностью СРД является возможность последовательной многократной фильтрации данных, при которой может осуществляться преимущественное повышение степени регуляризации данных с минимальными искажениями формы полезного сигнала. Для выполнения последнего размер окна К_с счета x_i и D_m устанавливается минимальным (3-5 точек), а глубина регуляризации данных (степень подавления шумов) устанавливается количеством последовательных операций фильтрации (до 3-5 проходов). Пример регуляризации модельного массива N в три прохода приведен на рис. 11.3.7. Для сравнения пунктиром на рисунке приведено сглаживание массива 5-ти точечным фильтром Лапласа-Гаусса, который имеет коэффициент подавления шумов, эквивалентный 3-х проходному СРД (см. рис. 11.3.9).

Рис. 11.3.8. Статистика средних значений Рис. 11.3.9. Статистика дисперсий (25 циклов. Счет Dm по массиву n) (25 циклов. Счет Dm по массиву n)

На рисунках 11.3.8 и 11.3.9 приведены результаты статистической обработки 3-х проходной СРД для 25 циклов моделирования в сравнении с 1-м проходом и с 5-ти точечным фильтром Лапласа-Гаусса (кривая n5).

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

Рис. 11.3.10. Модули спектров модельных сигналов. Рис.11.3.11. Участок спектра. (1– входной массив N, 2– выходной массив Z, один цикл CРД, 3– выходной массив Z,три цикла CРД), 4 – массив нерандомизированного меандра).

Частотное представление работы СРД хорошо видно на рис. 11.3.10, где приведены модули спектров рандомизированного сигнала в виде меандра (средние значения в минимуме - 20, в максимуме - 100, 25 периодов по 40 отсчетов, всего 1000 отсчетов) и результатов его обработки СРД (окно К_с= 3, окно К_s= 3). Модуль спектра основного полезного сигнала (в данном случае чистого меандра) представляет собой последовательность отдельных частотных гармоник по всему диапазону спектра. В спектре рандомизированного меандра эти частотные гармоники суммируются со спектром шума, статистически равномерно распределенным по всему частотному диапазону (спектр шума на рисунке для наглядности сглажен). СРД осуществляет подавление шумовых составляющих сигнала, практически не затрагивая частотных гармоник меандра и не изменяя их по амплитуде. Последнее можно видеть на рис. 11.3.11, где представлен отрезок спектра сигналов в высокочастотной части главного диапазона в области одной гармоники меандра (частотные составляющие шума не сглажены). При 3-х цикловом СРД высокочастотные составляющие шумов подавляются практически на порядок.

Рис. 11.3.12. Диаграммы ГК.

На рис. 11.3.12 и 11.3.13 приведен пример практического использования СРД при опробовании участка скважины, пересекающей пласты каменной соли, на содержание сильвинита по гамма-излучению Калия-40. По данным геологического опробования пласты сильвинита в толще вмещающих пород (галита) имеют достаточно резкие границы и однородны по содержанию сильвинита в пределах пластов. Исходная диаграмма ГК (детектор CsJ(Tl) со свинцовым фильтром толщиной 2 мм) и результаты фильтрации исходного массива данных ГК с использованием СРД и низкочастотного фильтра с весовым окном Лапласа-Гаусса приведены на рис. 11.3.12.

Рис. 11.3.13. Результаты интерпретации диаграмм ГК.

Результаты интерпретации диаграмм ГК симметричным деконволюционным цифровым фильтром (окно 13 точек) приведены на рис. 11.3.13. Как видно на рисунке, деконволюция по несглаженной диаграмме ГК дает существенные вариации содержания сильвинита в пределах пластов. Применение низкочастотной фильтрации диаграммы ГК снимает флюктуации содержания в пределах пластов, но существенно сглаживает границы пластов. Использование СРД позволяет устранить этот недостаток.

В заключение отметим, что СРД может использоваться для регуляризации не только ядернофизических данных, но и любых других числовых массивов непрерывных измерений, если радиус их корреляции не менее 3-5 отсчетов. В качестве примера на рис. 11.3.14 приведена диаграмма акустического каротажа, зарегистрированная с шагом дискретизации данных 20 см, сглаживание которой проведено СРД без потери пространственного разрешения.

Рис. 11.3.14. Диаграмма акустического каротажа и результат ее обработки СРД

(5 циклов, K_c = K_s = 3, физическое окно 0.6 м).

11.4. Статистическая группировка полезной информации.

Что касается аппаратных способов реализации СГПИ, то он может быть выполнен в реальном масштабе времени, если информация представлена потоком импульсов и основным информативным параметром является скорость следования импульсов.

Сущность аппаратной реализации заключается в статистической (близкой к статистической) нормированной выборке импульсов из дополнительного потока m и их суммировании с основным потоком n с заданием условий выборки по отношению частоты следования импульсов в потоках. Пола­гая для непрерывного режима измерений M+1 = М, перепишем выражение (5.2.20) с подстановкой значения b в следующем виде:

z = N + (M/-N)·M/(M+D(M)). (11.4.1)

Умножим левую и правую части выражения на нормировочный коэффициент размножения выходного потока K = l+R:

Z = K·z= N + RN+(M/-N)·KM/(M+D(M). (11.4.2)

Заменим отсчеты RN выборкой сигналов из потока m:

RN = Р_вМ, (11.4.3)

где Р_в - вероятность выборки сигналов из потока m. Если вероятность выборки сигналов под­держивать равной значению

P_в = R/, (11.4.4)

то при этом будет иметь место

M/-N = Р_вM/R-N ® 0, (11.4.5)

и соответственно для выражения (11.4.2) имеем:

(M/-N)·KM/(M+D(M) ® 0, (11.4.6)

Z = N+P_вM ® N+RN. (11.4.7)

При статистической независимости величины х от частоты потоков n и m приведен­ные выражения действительны при определении значения как в целом по пространству измерений, так и для скользящих окон текущих значений по определенным интервалам предшествующих измерений. Действительно и обратное заключение: если по определенному интервалу измере­ний выражение (11.4.5) обращается в нуль, то установленная вероятность выборки соответствует условию (11.4.4). На этом принципе может проводиться аппаратная реализация СГПИ с авто­матической адаптацией к условиям измерений: управление процессом выборки импульсов из потока m и направление их на суммирование с потоком n по сигналам обратной связи с уст­ройства, следящего за обращением в нуль выражения (11.4.5).

Отметим особенности аппаратной реализации СГПИ с автоматической адаптацией под условия измерений.

Значение вероятности выборки Р_в не может быть больше 1. Отсюда из (11.4.3) следует, что для любых интервалов измерений должно выполняться условие М ≥ RN, а соответственно по всему пространству измерений должно выполняться условие ≥ R, чем и обуславливается выбор коэффициента R. Значение коэффициента R принципиально ограничивает степень положительного эффекта СГПИ (k_max ® 1+R), в отличие от СРД, где такого ограничения не имеется.

Относительная статистическая погрешность измерений выходного потока отсчетов Z соответ­ствует выражению (11.2.23) при условии постоянного значения величины Р_в, т.е. при установке значения Р_в по среднему значению величины в целом по пространству измерений. При автоматической адаптации под условия измерений значение вероятности Р_в по текущему среднему значению отношения n/m определенного предшествующего интервала измерений также является статистически флюктуирующей величиной с дисперсией распределения (без учета изменений действительного значения х):

D_p = R²(n+m)n/(m³T), (11.4.8)

где Т- интервал усреднения информации при определении текущего значения . Соответственно, дисперсия и средняя квадратическая погрешность текущих отсчетов Z:

D_z = D_N+ P_вD_M+M²D_p = N+Р_вМ+М²D_р, (11.4.9)

d_z² = (N+Р_вМ+М²D_р)/(N+Р_вМ)². (11.4.10)

При постоянной экспозиции измерений t положительный эффект возрастает с увеличением значения Т:

k = K²/(K+R²(n+m)t/mT). (11.4.11)

k_max ® 1+R, d_z² ® 1/(N+Р_вМ) при Т ® ¥. (11.4.12)

В общем случае, с учетом средней квадратической ошибки прогнозирования d_xi значе­ний x_i для текущих точек измерений по значениям в предшествующих интервалах при Т > t:

D_z = N+Р_вМ+M²(D_p+P_в² d_xi²). (11.4.13)

Формирование значения Р_в на основе информации по средним значениям интер­валов измерений, предшествующих текущим, определяет СГПИ как динамическую систему с соответствующей постоянной времени реакции на изменение условий измерений. Учитывая, что, во-первых, для любой точки пространства измерений должно выполняться условие m > nR, и, во-вторых, увеличение интервала Т приводит к возрастанию времени реакции на из­менение условий измерений, значение Т целесообразно ограничивать величиной порядка (5-10) значений текущих экспозиций. Чем меньше пространственная частота рас­пределения х по отношению к распределению n, тем большее значение Т допустимо.

Реализация систем СГПИ значительно облегчается при чисто пр