Секция 7 - MATLAB в образовании и Интернете (1250002), страница 19
Текст из файла (страница 19)
Их влияниеможно исключить с помощью того же аппарата передаточных функций,примененного к иной группе рядов наблюденных данных, связанных сэтим сторонним процессом. Связи вида (1) и (2) не универсальны. Ониимеют место, когда наблюдаемые данные линейно связаны с основнымпроцессом и связи между их рядами стационарны (хотя сами ряды могут ине обладать свойством стационарности). Передаточные функции зависят,как правило, от физико-геологического строения среды в окрестностипункта наблюдения. Если это строение изменяется в силу протекающих всреде процессов, то изменяются и сами передаточные функции.
Их изменение несет дополнительную информацию об изучаемых процессах.Адаптивные методы решения задач геофизикиПередаточные функции gi(t) могут быть найдены путем решения интегральных уравнений свертки (1) или (2) на отрезке времени, на которомзаданы значения рядов xi(t) и yi(t). После того, как эти функции определены, соотношения (1) и (2) могут использоваться для предсказания поведения ряда y(t) по известному поведению рядов xi(t) (но только, если передаточные функции неизменны во времени).
Если значения остаточного ряда∆y(t) превышают измерительные шумы и погрешности вычислений, то онможет нести независимую от xi(t) информацию. Безотносительно величины∆y(t) полезную независимую информацию могут давать также измененияво времени передаточных функций gi(t).Решение интегральных уравнений свертки (1) или (2) может бытьнайдено хорошо известными методами [4], однако для целей мониторинганаиболее подходящими представляются итерационные адаптивные методы[5]. В таких методах последовательно получаемые решения уравненийдостигаются при постоянно обновляющихся с течением времени значени1920Секция 7.
MATLAB в образовании и Интернетеях временных рядов y(t) и xi(t). Если связи между рядами не изменяются, топередаточные функции gi(t) принимают после ряда итераций значенияблизкие к истинным, и в дальнейшем колеблются в их малой окрестности.Если стационарность связей между рядами нарушается, то адаптивные методы позволяют отслеживать изменения передаточных функций во времени. Это, конечно, возможно лишь при более медленных их изменениях посравнению с вариациями самих рядов.Наиболее простым и в то же время эффективным методом решениясистем алгебраических уравнений, к которым после дискретизации сводятся уравнения свертки, зарекомендовал себя метод наименьших квадратовУидроу-Хоффа [5].
Итерационный алгоритм этого метода чрезвычайнопрост:(3)g i (k + 1, l ) = g i (k , l ) + µ ⋅ ∆y (k ) ⋅ xi (k − l ) ,где индекс k обозначает дискретизированное текущее время наблюдений,l — время задержки импульсной переходной характеристики, а µ — величину параметра сходимости. Значения ∆y(k) находятся путем вычитания изнаблюденного значения y(k) его синтезированного значения yS(k):NLy (k ) = ∑∑ g i (k , l ) ⋅ xi (k − l ) ⋅ ∆tS(4)i =1 l =1Последнее выражение представляет собой дискретизированнуюформу соотношения свертки (2) с определенными на k-ом отрезке временизначениями передаточных функций, ∆t — интервал дискретизации.Характерным примером является обработка данных наблюдений естественных электромагнитных полей Земли с применением описанноговыше аппарата передаточных функций.
Как известно, эти вариации, в основном, порождаются вариациями токов в ионосфере Земли. В силу удаленности этих токов от точек наблюдения поля на поверхности Земли, ихможно рассматривать в первом приближении как токи на плоскости, случайным образом изменяющие свою интенсивность и направление. В соответствии с выбранным представлением систему таких токов можно охарактеризовать двумя независимыми субпроцессами — изменениями составляющих этих токов по двум координатным направлениям на плоскости.
Эти токи порождают электромагнитное поле в форме вертикально падающей на Землю плоской волны. Если целью мониторинга является слежение за ионосферными токами, то в наблюдаемом на поверхности Землиполе следует выбрать также две независимо изменяющиеся во временикомпоненты электромагнитного поля [6]. В качестве таких компонент естественно выбрать две горизонтальные компоненты геомагнитного поля(Hx, Hy), на которые электрическое строение Земли оказывает наименьшеевоздействие — ввиду большой контрастности по сопротивлению атмосферы и литосферы эти компоненты просто удваиваются по сравнению с падающим полем.
Остальные компоненты электромагнитного поля тесно1921Труды II научной конференции «Проектирование инженерных и научных приложений в среде MATLAB»связаны с геоэлектрическим разрезом (электрическое поле Ex, Ey) или появляются в силу его горизонтальной неоднородности (вертикальная компонента магнитного поля Hz). В рамках сформулированной модели вариации этих компонент в соответствии с выражением (2) могут быть описанысоотношением:H Z (t ) = H X (t ) ∗ I ZX (t ) + H Y (t ) ∗ I ZY (t ) + ∆H Z (t ) ,(5)где Izx(t) и Izy(t) — компоненты индукционного вектора или, в более общейтрактовке, — импульсные переходные характеристики соответствующихпередаточных функций, отражающие электрическое строение Земли в окрестности точки наблюдения.
В вариациях остаточного поля (∆Hz) в значительной степени ослаблена зависимость от ионосферных токов. Существенный вклад в них вносят поля внутреннего, геодинамического происхождения, возникающие в недрах Земли в результате разнообразных механоэлектрических преобразований. Если в результате каких-либо геодинамических процессов не только генерируются электромагнитные поля, но иизменяется геоэлектрическое строение среды, то соответственно изменяются и передаточные функции Iij.
Их изменения несут независимую посравнению с остаточными полями информацию о протекающих внутриЗемли процессах, которая может быть эффективно использована для целейгеодинамического мониторинга.Примеры реализации вычислительного алгоритмаБлок-схема вычислительного алгоритма метода наименьших квадратов Уидроу-Хоффа представлена на рис. 1.Рис. 1. Блок-схема вычислительного алгоритма.1922Секция 7. MATLAB в образовании и ИнтернетеНа этой схеме блоки «Адаптивный фильтр» 1 и 2 реализуют итерационный алгоритм (3). В соответствии с выражением (5) в качестве опорных сигналов на входах блоков используются две горизонтальные компоненты геомагнитного поля (Hx, Hy), а в качестве основного — вертикальнаякомпонента Hz. Разностный сигнал на выходе находится путем вычитанияиз наблюденного значения Hz его синтезированного значения в соответствии с выражением (4).Рассмотренный алгоритм адаптивного фильтра был применен намидля анализа данных магнитотеллурического мониторинга, проводимого наБишкекском геофизическом полигоне в начале 90-х годов [7].
С помощьюпакета программ, написанных на Паскале, были обработаны ряды данныхнепрерывной в течение нескольких лет регистрации электромагнитногополя Земли, и графики, характеризующие динамику передаточных функций во времени, были сопоставлены с данными о сейсмической активностирегиона.
На основании полученных результатов были отработаны методики мониторинга сейсмической активности.В пакете расширения Filter Design, начиная с версии 2.1, входящей впоставку MATLAB 6.1, содержатся функции, реализующие алгоритмыадаптивной фильтрации [8]. Однако, по целому ряду соображений, намибыли использованы функции собственной разработки. Код основнойMATLAB-программы, реализующей алгоритм метода наименьших квадратов Уидроу-Хоффа, приведен ниже:k=1:L;Ys = sum(sum(W(:,k).*X(:,i-k+1)));E(i) = S(i)-Ys;for j=1:NrefW(j,k)=W(j,k)+2*U(j)*E(i)*X(j,i-k+1);end;Здесь: L — длина фильтра; Nref — число опорных сигналов; E — разностный сигнал; YS — синтезированный сигнал; W — передаточные функции.Разработанная нами программа позволяет использовать любое число опорных сигналов и произвольное количество циклов адаптации.Новые возможности, предоставляемые пакетом MATLAB, позволилипровести дополнительную обработку данных.
На Рис. 1. показан примервычисления передаточных функций между компонентами магнитного поляЗемли в соответствии с выражением (5) продолжительностью в 24 часа. Вкачестве основного сигнала принята вертикальная компонента вариацийгеомагнитного поля, в качестве опорных — горизонтальные компоненты.На графиках представлены исходный основной временной ряд (Raw data),разностный временной ряд (Residuals), а также передаточные функции между вертикальной и горизонтальными компонентами.1923Труды II научной конференции «Проектирование инженерных и научных приложений в среде MATLAB»Рис.
1. Пример вычисления передаточных функций между компонентамиэлектромагнитного поля Земли: 1 — основная компонента HZ; 2 — остаточное поле;3–4 — передаточные функции по компонентам HX, HY.Из графика видно, как в течение суток меняется характер вариацийэлектромагнитного поля, а также характер связи между компонентами, чтоотражается в изменениях передаточных функций. Детальный анализ физического смысла этих изменений выходит за рамки данной работы.Спектрально-временной анализ узкополосным адаптивным фильтромОпределенный интерес представляет анализ изменений спектра колебаний в разностном (остаточном) поле. Учитывая специфику геофизических исследований, а именно, интерес к низкочастотной части спектра вдиапазоне периодов от единиц минут до часа, пришлось отказаться от традиционных методов спектрального анализа.
Для спектрально-временногоанализа рядов данных нами был использован метод выделения гармонических составляющих с помощью узкополосных адаптивных режекторныхфильтров [9]. Их преимущества заключаются в простоте перестройки полосы пропускания, практически неограниченном подавлении соседнихгармоник и точном слежении за частотой. Принцип их действия можетбыть описан той же блок-схемой рис.1. Отличие заключается лишь в том,что на входы адаптивных фильтров 1 и 2 в качестве опорных сигналов подаются две квадратурные составляющие моногармонического сигнала с1924Секция 7. MATLAB в образовании и Интернетечастотой ω 0 , равной частоте выделяемой гармоники.
Анализируемый временной ряд S i поступает на вход основного сигнала. Формируемые в процессе вычислений квадратурные составляющие x1i и x2i в каждой точкеанализируемого временного ряда принимают значения:x1i = C ⋅ cos(ω 0 i∆t + ϕ )(6)x2i = C ⋅ sin(ω 0 i∆t + ϕ )где C — амплитуда колебания, ∆t — шаг дискретизации. Из этих значений и весовых коэффициентов адаптивного фильтра W1i и W2i формируется синтезированный сигнал S si = x1iW1i + x2iW2i . Разность между исходными синтезированным сигналами ε i = S i − S si используется в качестве параметров, управляющего обновлением весовых коэффициентов фильтра накаждом шаге вычисления в соответствии с выражениями:W1,i +1 = W1,i + 2 µε i x1i(7)W2,i +1 = W2,i + 2 µε i x2iгде µ — параметр, характеризующий скорость сходимости алгоритмаадаптации коэффициентов фильтра.
Добротность такого фильтра определяется выражением Q = ω 0 ∆t /( 2 µC 2 ) . Синтезированный сигнал S si является выделяемой гармоникой с частотой ω 0 , огибающая которой находитсяпо формуле:Ai = C W1i2 + W22i(8)Код MATLAB-программы узкополосного адаптивного режекторногофильтра в соответствии с выражениями (6) и (7), приведен ниже:Xs1 = Amax.*sin(pii.*step);Xs2 = Amax.*cos(pii.*step);step = step+1;Ys = Ws1.*Xs1+Ws2.*Xs2;Eps = AvErr(i)-Ys;Ws1 = Ws1 + (2*Eps).*Xs1;Ws2 = Ws2 + (2*Eps).*Xs2;При помощи такого фильтра были обработаны данные временногоряда остаточного поля из предыдущего примера. Результат их обработки ввиде динамического спектра показан на рис.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
