Анализ разброса погрешностей магнитотеллурических данных в зависимости от времени суток на Балтийском щите по результатам статистической обработки данных bear Текст научной статьи по специальности «Науки о Земле и смежные экологические науки»

ВЕСТНИК САНКТ-ПЕТЕРБУРГСКОГО УНИВЕРСИТЕТА

Сер. 4. 2009. Вып. 2

УДК 53.088-550.370 А. К. Зароченцев

АНАЛИЗ РАЗБРОСА ПОГРЕШНОСТЕЙ МАГНИТОТЕЛЛУРИЧЕСКИХ ДАННЫХ В ЗАВИСИМОСТИ ОТ ВРЕМЕНИ СУТОК НА БАЛТИЙСКОМ ЩИТЕ ПО РЕЗУЛЬТАТАМ СТАТИСТИЧЕСКОЙ ОБРАБОТКИ ДАННЫХ BEAR

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

Метод магнитотеллурического зондирования (МТЗ) основан на использовании модели Тихонова-Каньяра [6], подразумевающей зависимость компонентов тензора импеданса (КТИ) только от характера проводимости Земли и от периода вариации, в предположении о горизонтальной однородности первичного поля, обусловленной удалённостью источника, а также малостью длины волны в среде по сравнению с расстоянием до источника.

К недостаткам метода относиться его малая помехоустойчивость. С ростом периода и, соответственно, длинны волны минимальное расстояние до источника, удовлетворяющее условиям модели Тихонова-Каньяра, увеличивается. Кроме того, на длину волны в среде и, соответственно, на минимальное расстояние до источника может влиять и геологические характеристики местности, такие как проводимость пород. Это приводит к тому, что на низких частотах далеко не каждый выделенный сигнал соответствующей длины волны удовлетворяет модели Тихонова-Каньяра. Кроме расстояния до источника на качество МТЗ-данных может влиять так же и тип источника, будь то продольный линейный ток, магнитный или электрический диполь [5]. Эти обстоятельства приводят к тому, что на больших периодах значение КТИ перестает быть независимым от источника. На малых периодах похожая картина может наблюдаться в случае наличия искусственных источников, но такие источники легко локализуются и отфильтровываются. Для периодов, начиная с сотен секунд, затруднена как сама локализация источника, так и определение, в случае локализации, качества этого источника. Поэтому на низких частотах разумнее определять не тип источника, а устойчивость самих элементов тензора импеданса и соответствие их модели Тихонова-Каньяра. Но такая обработка возможна только при наличии большой статистики (большого набора данных). Такой подход уже оправдал себя на примере обработки данных эксперимента BEAR [3], представлявшего собой синхронную запись 5-ти компонентов электромагнитного поля Ex, Ey, Hx, Hy, Hz с шагом в 2 с в 50 пунктах Балтийского щита в течение 2-х месяцев

© А. К. Зароченцев, 2009

Baltic Electromagnetic Array Research

10 o

40 o

200 км I________I

Рис. 1. Распределение исследовательских станций в эксперименте BEAR

(рис. 1). На основе этих данных, путём сглаживания и использования робастных оценок построены импедансные кривые для периодов до 10 000 с [10].

Но BEAR достаточно дорогой эксперимент, а как уже сказано, одно из достоинств МТЗ - его низкая стоимость. Если заранее выделить какие либо классы вариаций электромагнитного поля, на которых разброс значений вычисленных КТИ будет минимальным, но охватывающие весь изучаемый спектр частот, то впоследствии можно ограничить измерения и вычисления только выделенным классом вариаций и тем самым сократить и удешевить эксперимент. Для определения таких классов можно воспользоваться данными эксперимента BEAR, на основе которых оценить поведение величин, характеризующих разброс значений.

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

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

Данные эксперимента BEAR были уже неоднократно обработаны и интерпретированы. Для примера возьмём обработку данных, проведённую Смирновым и Варенцовым [10], результаты которой были любезно предоставлены автору М. Ю. Смирновым. В результате этой обработки были получены значения КТИ для периодов от 101 до 105 с, с относительной погрешностью от 0,5 % до 10 %. На рис. 2 в двойном логарифмическом масштабе представлена зависимость относительной погрешности о = &.Z/Z элемента Zxy тензора импеданса точки В23 от корня периода %/Т. В области периодов VT « « Зс1/2, а/Т» 1000 с погрешность имеет минимум и далее начинает возрастать, не имея более чётко выраженных, повторяемых в различных точках, экстремумов. Похожую картину можно наблюдать и для других точек и элементов. В данной работе мы не будем заострять внимание на возможных причинах такого поведения погрешностей в зависимости от периода, но отметим сам факт возрастания погрешностей начиная с периода в 1000 с.

В предположении суточной периодичности колебаний искомого класса и, следовательно, суточной периодичности увеличения отклонений значений КТИ от истинного значения, построим временной ряд значений КТИ в зависимости от времени суток. Для создания этого временного ряда будем использовать программу Астапенко САМИ [1]. Эта программа, имея в качестве исходных данных синхронизированный набор значений Ex(t), Ey(t), Hx(t), Hy(t) электромагнитного поля и искомую частоту Ю, с помощью преобразования Фурье и узкополосной цифровой фильтрации строит временной ряд Ex(t, ю), Ey(t, ю), Hx(t, ю), Hy(t, ю). На основе этого временного ряда составляется избыточная система уравнений (1) с постоянной Ю и с помощью метода рекуррентных наименьших квадратов находит элементы тензора импеданса Z:

(Ex(t, ю)\ = (Zxx(^) Zyx(rn)\ (Hx(t, ю)\ (1)

\Ey (t, ю)/ \Zxy(ю) Zyy(®>) ) \Hy (t, ю) )

щ 1

0,1 0,01 0,001

1 10 100 Г112, с1/2

Рис. 2. Поведение относительной погрешности вычисления элемента тензора импеданса р в точке Б23 в зависимости от корня периода

Рис. 3. Распределение значений КТИ Z,

ху

за разные часы в точке В23

N

В этом случае появляется возможность выделить в течении суток такие временные интервалы, где значения КТИ имеют систематическую ошибку, что может говорить о наличии близко расположенного локального источника ЕЭМП.

Процедура составления суточного временного ряда заключается в измерении КТИ за короткие промежутки ДЬ через равный интервал времени т, который составляет 1/М от суток, где М целое число:

т = 1/М сут = 24/М ч.

После этого рассматривается распределение значений элементов тензора импеданса привязанных к определённому времени суток = кт, где к < М, по всему времени измерения. Для примера, на рис. 3 показано распределение значений Zxy для точки В42 за весь период наблюдения. Значения подсчитаны за временные промежутки Ь = 10 240 с + 3 чи привязаны к чётным часам мирового времени. Пунктирной линией отмечено значение импеданса, рассчитанное М. Ю. Смирновым и И. М. Варенцовым за весь период наблюдения. Значения, подсчитанные за разные сутки, отмечены разной формой. Для каждого часа наблюдается большой разброс значений. Задача заключается в том, чтобы найти характеристику этому разбросу, оценить его параметры, и выделить закономерность изменений этих параметров в течение суток для разных точек. Найденные характеристики должны быть устойчивы к малым изменениям данных, таких как добавление или изъятие нескольких дней.

Рассмотрим различные характеристики распределения значений на примере точки В42. Для этого построим график выборочной функции распределения ) значений модуля КТИ Zxy, рассчитанные за промежуток времени ДЬ = 3 ч, начиная с 22-х часов мирового времени (иТ) на периоде Т = 1024 с для различных дней. На рис. 4а для каждого значения Zxy, отложенных по оси абсцисс, по оси ординат поставлена в соответствие величина, равная количеству значений в выборке меньших, чем Zxy, нормированная на общее количество значений в выборке N. Эта величина и определяет выборочную функцию распределения Г. На рис. 4б представлена выборочная плотность функции распределения ]^) значений модуля Zxy в точке В42, рассчитанных начиная с 22 ч на том же периоде. По оси У отмечено количество отсчётов п^ху) модуля КТИ, зафиксированных в промежутке ^ху — ДZxy/2, Zxy + ДZxy/2), где ДZxy = 3,78 • 10~4 Ом, нормированное на N и делённое на величину промежутка ДZxy. Эта величина и определяет выборочную плотность функции распределения ]^).

Я2у)

2 , Ом

ХУ

Рис. 4- Плотность распределения значений 2ху в точке В40 в 18, 20 и 22 часа по Гринвичу:

Т = 1024 с, AZxy = 0, 3 мВ/км-у, Zxy (18) = 0, 75 мВ/км-у, Zxy(20) = 1, 05 мВ/км-у, АZxy (22) = 0, 75 мВ/км-у

а

Простейшие характеристики такого распределения - это среднее арифметическое Z, которое является оценкой математического ожидания, и выборочное квадратичное отклонение 5, которое является оценкой корня дисперсии В:

Б =

N _

Е(г-г>)2

г=1

N - 1 :

где N - число значений в выборке, Z^ - г-е значение. Но такие характеристики оказались неустойчивыми для выборки состоящей из 20-50 наблюдений. Появление одного значения, отличающегося от остальных на порядок (достаточно частое явление), приводит к разбросу найденных величин как минимум в 1,2 раза, что сравнимо с теми вариациями, которые мы хотели бы выделить [4].

Можно отметить, что величина, называемая «вероятным отклонением» [3], как мера разброса более устойчива. Вероятное отклонение определяется, как половина разности квартилей Хх/4 и Х3/4 - значений X, в которых функция распределения переходит границы:

Р (Х1/4) = 1/4 Р (Х3/4) = 3/4,

где X - значения непрерывного числового ряда, описывающего распределение искомой величины.

Выборочная функция распределения определяется:

Р (^) = N

где Nz<z. - количество значений Z меньших Z^ (0 < г ^ N), а N - всё количество значений в выборке. Величина вероятного отклонения оценивается выборочным вероятным отклонением по формуле:

щг) = ^3/4 ~ ^1/4,

где Z1/2 и Zз/4 такие элементы выборки, в которых выборочная функция распределения переходит соответствующие границы:

Р(Z1/2) ^ 1/4 Р(Z3/4) ^ 3/4-

На рис. 4а горизонтальные прямые отмечают значение выборочной функции вероятности Р^ху) 1/4 и 3/4, а вертикальные прямые отмечают соответствующие квартили. На рис. 4б эти же квартили отмечены на графике выборочной плотности функции вероятности f ^ху).

На рис. 5 представлены квадратичное отклонение Б, выборочное вероятное отклонение Ш, среднее арифметическое и значение модуля КТИ, рассчитанное в работе Смирнова и Варенцова ZN для точки В42, элементы с индексом ху в зависимости от времени суток по Гринвичу. Такое относительное поведение дисперсии и вероятного отклонения

2 , Ом

ху

ИТ, ч

Рис. 5. Поведение основных характеристик распределения КТИ 2ху для точки В23:

— значение КТИ, вычисленное М. Ю. Смирновым, 2Ху — мода КТИ, вычисленная по гистограмме с шагом 0,1 мВ/км-у, 0(2Ху) — дисперсия КТИ, 2\/4 и 2з/4 — квартили

2-1/4 и 23/4

нехарактерно для нормального распределения даже с учётом точности оценки. Оценка математического ожидания может превышать величину ZN, потому что в предложенном варианте обработки не отбрасывались заведомо ложные - завышенные значения импеданса, в отличие от работы Смирнова-Варенцова, но подобное поведение так же может объяснено не «нормальным» характером распределения величины Z, когда мода - наиболее вероятное значение - отличается от математического ожидания.

Для оценки характера распределения величины Z мы имеем следующие предварительные данные: величина Z неотрицательна; распределение величины Z имеет асимметрию; значение Z может превышать среднее значение более чем в два раза. Указанным характеристикам удовлетворяет приближение логнормальным распределением |ц,о) с параметрами (о,ц) [8, 9]:

1 1п2(х/ц)

, О) = ---у= е 2а2 .

хО\/ 2п

Примем параметр распределения о за характеристику разброса значений КТИ в различное время суток. При таком подходе приближение логнормальным распределением является достаточно точным [8]. В логнормальном распределении Ь^|ц, о) величина о соответствует дисперсии нормально распределённой величины 1п Z/ц. Так как нас не интересует значение Ц, характеризующее натуральный логарифм медианы распределения - величины, делящей область определения случайной величины

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

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

Применение метода. Как уже было сказано, нас интересуют периоды, начиная с 600 с. Для первых исследований был взят период T = 1024 с. Интервал измерения At брался равным 10T = 10 240 с, т. е. примерно 3 ч, а временной промежуток между измерениями т равным 2 ч.

После предварительной обработки из 50 точек BEAR было выбрано 22 точки. На рис. 6 показано поведение параметра D ln Z для выбранных точек в зависимости от времени суток и жирной линией выведено среднее арифметическое всех точек. На рисунке наблюдается увеличение отклонений КТИ с отметки 10 ч и до отметки 18 ч. Для перевода отметок на шкале измерений в местное время напомним, что интервал измерения у нас приблизительно равен 3 ч, разумно брать середину этого интервала, т. е. прибавлять к отметке его начала 1,5 ч. Рассматриваемая область лежит между 15 и 30 параллелями восточной долготы, т. е. для приведения значений к местному времени нам надо прибавить ещё 1,5 ч. В итоге, найденный, наиболее возмущённый, интервал соответствует интервалу местного времени с 13 до 21 ч.

В приведённых примерах упоминается только компонента Zxy. Это объясняется тем, что на большинстве точек эксперимента BEAR импеданс имеет северную поляризацию и, соответственно меньшую восточную компоненту, оценка которой происходит с большей относительной погрешностью. В итоге поведение погрешностей восточной компоненты гораздо более хаотично, но и в этом случае среднее значений дисперсии логарифма импеданса в различных точках даёт тот же возмущённый интервал: 10-18 ч по используемой шкале, т. е. 13-21 ч местного времени.

Были проведены измерения тех же величин для периода 609 с. В этом случае At брался равным 10T = 6 080 с, т. е. примерно 1,5 часам, а временной промежуток т между измерениями равный 1 ч.

На рис. 7 показано поведение дисперсии логарифма КТИ для T = 609 с в 9 точках и их среднего. Как и в случае с периодом 1024 с, наблюдаются увеличение отклонений КТИ с 10 до 20 ч используемой шкалы.

Заключение. По результатам данной работы можно говорить об интервале суток, для которого характерно увеличение погрешности измерения КТИ. На данном этапе рано говорить о характере источников этих погрешностей - автор не обнаружил никаких закономерностей в поведении значений КТИ в разные интервалы суток, ни постоянного смещения в большую или меньшую сторону, ни постоянных значительных отклонений от значения, принимаемого за истинное. Следовательно нельзя пока говорить о том, что в найденном временном интервале значения КТИ отличаются от истинного для идеальной модели. Можно только говорить, что погрешность определения КТИ в данном интервале больше, и определение «истинного значения» требует большего количества данных. Такой интервал найден на конкретной территории и только для изучаемого времени года (эксперимент BEAR проводился в июне-июле) для определённых периодов (600-1000 с), но сам метод определения такого периода показал себя

М — V

Ь01 ----------

Ь02-----------

Ы1------------

ЫЗ------------

Ы5------------

Ы8 ..........

Ь20---------

Ы2 -----------

Ь23 - ■ -

Ь24 ----------

Ь27 ---------

Ь28 ----------

ЬЗО ----------

Ь41 -------

Ь05 --------

ЫО ----------

Ы4 -----------

Ы6 ----------

Ы9 - - -Ь36 ■ ■ ■

Ь42-----------

Ь45 --------

иТ, ч

Рис. 6. Поведение разброса значений элемента тензора импеданса 2ху, нормированного на «истинное» значение ,

оценённых в «спокойных» точках для периода Г = 1 024 с

Рис. 7. Поведение разброса значений элемента тензора импеданса 2ху, нормированного на «истинное» значение

оценённых для периода Г = 609 с

достаточно устойчивым, и автор надеется использовать его в дальнейших исследованиях.

В завершение подчеркнём следующее:

а) величина выборочного квадратичного отклонения логарифма модуля импеданса достаточно устойчиво характеризует разброс значений КТИ;

б) летом, днём с 13 до 21 ч местного времени на Балтийском щите снятие МТЗ-данных будет происходить с большими погрешностями и экономически менее выгодно ночного зондирования.

Автор благодарит профессора Ковтун Аиду Андреевну и Смирнова Максима Юрьевича за предоставленные материалы и помощь в работе

Литература

1. Астапенко В. Н., Емельянов А. П. Система анализа магнитотеллурической информации «САМИ» // Алгоритмы и программы анализа данных геофизических исследований: сб. статей / под ред. Ж. П. Хотько, А. П. Емельянова. Минск, 1979. C 58-86.

2. Ван дер Варден Б. Л. Математическая статистика / пер. с англ. М., 1960. C. 108-109.

3. Варенцов И. М., Соколова Е. Ю., Мартанус Е. Р. и др. Система передаточных операторов электромагнитного поля для массива синхронных зондирований BEAR // Физика Земли. 2003. № 2. C 30-61.

4. Зароченцев А. К. Анализ погрешностей магнитотеллурических данных во времени и пространстве по результатам международного эксперимента BEAR // Тр. 5-й междунар. на-учн.-техн. конф. молодых учёных и специалистов «Геофизика 2005». СПб., 2005. C. 103-105.

5. Ковтун А. А. Строение коры и верхней мантии на северо-западе Восточно-Европейской платформы по данным магнитотеллурических зондирований. Л. 1989. C. 30-44.

6. Ковтун А. А. Использование естественного электромагнитного поля при изучении электропроводности Земли. Л. 1980. C. 11-26.

7. Ковтун А. А., Вагин С. А., Варданянц И. Л. и др. Анализ магнитотеллурических и магнитовариационных результатов в интервале периодов суточных вариаций по данным BEAR и определение «нормального» разреза Балтийского щита // Физика Земли. 2002. № 11. C. 34-53.

8. Мальцев В. А. К вопросу о получении несмещённых оценок при геостатистическом моделировании с использованием нелинейных нормализующих преобразований // Горн. ин-форм.-аналит. бюл. 1999. Вып. 5.

9. http://algolist.manual.ru/maths/matstat/logNormal/index.php.

10. Смирнов М. Ю. Обработка магнитотеллурических данных с использованием робастных статистических процедур // Вопр. геофизики. 1998. Вып. 35. C. 198-205.

Принято к публикации 7 октября 2008 г.