Содержание к диссертации
Введение
1. Становление электромагнитного поля в горизонтально-слоистой среде 19
1.1. Математическая модель поля вертикального магнитного диполя 19
1.2. Решение задачи методом конечных элементов 21
1.3. Вычисление значений dBJdt 23
1.4. Применение методики выделения основного поля 30
1.5. Вычисление ЭДС 38
1.6. Проблема вычисления ЭДС на ранних временах и больших разносах и ее решение 40
2. Восстановление параметров горизонтально-слоистой среды 57
2.1. Построение трансформаций результатов геофизических наблюдений 58
2.2. Автоматическое восстановление параметров горизонтально-слоистой среды 62
2.2.1. Построение начального приближения 63
2.2.2. Уточнение параметров 68
3. Трехмерная обратная задача для технологии с закрепленным источником 74
3.1. Основные трудности применения традиционных методов интерпретации для технологии ПЗС-ЗИ 75
3.2. Анализ возможности построения системы интерпретации данных для технологии ПЗС-ЗИ 77
3.3. Основные идеи предлагаемого способа интерпретации практических данных 91
3.4. Итерационная процедура восстановления свойств среды с использованием регуляризации 93
3.5. Пример применения процесса восстановления свойств среды для модели с кимберлитовым объектом 96
3.6. Анализ работоспособности процедуры интерпретации при несовпадении размеров «пробных» объектов с размерами поисковых 98
3.7. Анализ работоспособности процедуры интерпретации при неточных параметрах вмещающей модели 100
3.8. Примеры решения трехмерных обратных задач на синтетических данных 105
3.8.1. Проводящий объект в слоистой среде 105
3.8.2. Высококонтрастные объекты в слабопроводящей среде 106
3.8.3. Непроводящий объект в слоистой среде 109
3.8.4. Два проводящих объекта в слоистой среде 110
3.9. Апробация системы интерпретации на экспериментальных данных 111
3.9.1. Поиск слабоконтрастных объектов 111
3.9.2. Поиск высококонтрастных объектов 113
4. Трехмерная обратная задача для технологии с коаксиальной установкой 118
4.1. Одномерная интерпретация 118
4.1.1. Интерпретация с использованием одномерных трансформаций 118
4.1.2. Интерпретация на основе слоистой модели среды 122
4.2. Трехмерная интерпретация 124
4.3. Примеры решения трехмерных обратных задач на синтетических данных 130
4.3.1. Профилирование над неоднородностью в пласте 130
4.3.2. Профилирование над кимберлитовой трубкой 131
4.3.3. Три проводящих объекта в слоистой среде 134
5. Программный комплекс 139
5.1. Программный комплекс HORIZON 139
5.2. Сравнение возможностей программного комплекса HORIZON с существующими разработками 144
5.3. Программный комплекс IFSTEM 147
5.3.1. Трехмерная интерпретация 150
5.3.2. «Ручное» восстановление параметров объектов 152
5.3.3. Трехмерная визуализация 153
Заключение 156
Список использованных источников 158
- Применение методики выделения основного поля
- Автоматическое восстановление параметров горизонтально-слоистой среды
- Анализ возможности построения системы интерпретации данных для технологии ПЗС-ЗИ
- Примеры решения трехмерных обратных задач на синтетических данных
Введение к работе
Одним из ведущих электроразведочных методов, применяемых при исследовании земных недр, является индукционное электромагнитное зондирование методом переходных процессов (Ml 111). Главным параметром, изучаемым в этом методе, является электропроводность среды. Задача состоит в том, чтобы выяснить распределение этого параметра внутри земли по результатам измерений характеристик электромагнитного поля на поверхности земли (так называемой дневной поверхности) или в воздухе. Ml 111 основан на изучении неустановившегося электромагнитного поля источника с импульсно изменяющимся током. Метод применяется для решения широкого круга геологических задач, связанных с исследованиями как больших, так и малых глубин, а также с поисками слабоконтрастных объектов. Основы применения неустановившихся полей в электроразведке содержатся в классических теоретических работах А.Н. Тихонова, СМ. Шейнмана, Дж.Р. Уэйта и других российских и зарубежных ученых [2-4, 7, 8, 12-16, 21, 26,50,77,79,89,92,100-102].
Для решения многих задач исследования структуры среды с использованием Ml ill актуальной является задача моделирования электромагнитного поля в среде, состоящей из нескольких однородных по проводимости а горизонтальных слоев. Как правило, первым этапом интерпретации результатов является получение оценки геоэлектрического разреза на основе горизонтально-слоистой модели.
Эффективность различных процедур восстановления свойств среды по данным экспериментальных наблюдений и процедур выделения аномалий в горизонтально-слоистой среде во многом определяется возможностью быстро вычислять значения характеристик электромагнитного поля в заданных точках расчетной области.
В настоящее время для интерпретации данных электромагнитных зондирований методом переходных процессов существует ряд программных ком- плексов, довольно характерными примерами которых являются комплексы ЭРА (Л.А. Табаровский, М.И. Эпов, И.Н. Ельцов - Россия), ПОДБОР (B.C. Могилатов, А.В. Злобинский - Россия), TEMIX (Ch.H. Stoyer, W.L. Anderson - США), EMVision (Австралия). Расчеты электромагнитных полей в такого рода комплексах базируются в основном на аналитических методах решения краевых задач [1-5, 7, 8, 14, 21, 26-34, 48, 50, 51, 75, 77, 79, 89, 92, 100-108].
Трудности, связанные с применением аналитических методов решения многомерных краевых задач, привели к необходимости поиска новых подходов к решению задачи становления поля, основанных на применении одного из наиболее универсальных методов решения дифференциальных уравнений в частных производных - метода конечных элементов (МКЭ). Однако, применение стандартных конечноэлементных схем [10, 18, 24, 42, 45, 57-60, 62-64, 74, 91, 99] не приводит к быстрым алгоритмам вычисления характеристик электромагнитного поля. Использование же методики выделения основного поля, предложенной в работах [37—41, 52-54, 56, 57, 64-71, 111-113], позволяет получать решение с высокой точностью и относительно малыми вычислительными затратами. В предлагаемой диссертационной работе будет построена основанная на таком подходе и адаптированная для одномерного случая схема моделирования электромагнитного поля в горизонтально-слоистой среде на базе постановки, предложенной А.Н.Тихоновым [79].
Большинство существующих методов решения обратной задачи в рамках горизонтально-слоистой среды сводятся к уточнению параметров некоторым методом многомерной оптимизации и, следовательно, требуют задания начального приближения, из которого в дальнейшем осуществляется спуск в точку минимума некоторого функционала. Это обстоятельство осложняет использование таких методов, так как при недостаточно удачном выборе начального приближения методы многомерной оптимизации могут останавливаться в локальных минимумах. Поэтому большой интерес вызывает задача автоматического построения начального приближения с заданной точностью.
На практике при проведении электромагнитных исследований становлением поля наиболее часто используется установка с соосными генераторной и приемной петлями, перемещаемыми по некоторой площади или отдельным профилям. Однако технология с закрепленной генераторной петлей и перемещаемыми приемными петлями гораздо более экономична и обладает не меньшей разрешающей способностью [70, 83-87]. Более широкому внедрению такой технологии препятствуют в основном трудности, связанные с интерпретацией измеренных характеристик электромагнитных полей. Для закрепленного источника ЭМ поля характерны знакопеременные процессы становления поля в пространственно-временной области. Сигнал переходного процесса за пределами генераторного контура может иметь несколько переходов через ноль даже в случае горизонтально-слоистой среды [26, 30]. Поэтому все попытки получить устойчивую процедуру трансформации в кажущиеся параметры одномерной среды не дали положительного результата [4, 14, 22, 23, 26, 30]. Более того, попытки использовать методы интерпретации, основанные на одномерных моделях строения среды, для установки с закрепленным источником чаще всего не позволяют локализовать поисковые объекты из-за несоответствия аномалий проводимости в исследуемой среде аномалиям измеренных на дневной поверхности значений ЭДС как по месту их расположения в плане, так и по знаку аномалии. Так, для объекта с повышенным сопротивлением в каждый момент времени на дневной поверхности есть точки в окрестности этого объекта как с положительными аномальными значениями ЭДС в них, так и с отрицательными, причем максимумы абсолютных значений аномальных ЭДС могут быть довольно сильно смещены относительно границ поискового объекта в плане [17, 70, 72, 73, 86, 87]. Поэтому подходы к интерпретации данных для установки с закрепленным источником, основанные на одномерных моделях среды, вряд ли приведут к положительному результату- поле влияния трехмерного объекта, смещенного относительно центра генераторной петли, не может быть описано с использованием одномерной модели среды. Таким обра- зом, необходимость разработки трехмерного подхода к обработке данных индуктивных импульсных зондирований очевидна.
Способ обработки данных от закрепленного источника, основанный на упрощенной модели физического процесса, а именно - на приближенном (линеаризованном) трехмерном подходе применялся, например, в работах [25, 26, 30, 35, 36]. Разработанный в них метод основан на построении линеаризованного решения трехмерной прямой задачи и дальнейшей линейной инверсии (обращении линейной системы, связывающей экспериментальные данные и возмущения проводимости относительно вмещающей среды). Для этого при расчете поля влияния трехмерной неоднородности (которая рассматривается в виде сектора кольца) полное поле в соответствующем интегральном уравнении заменяется нормальным. В результате такого подхода поле влияния неоднородности оказывается линейно связанным с ее аномальной проводимостью.
Такой упрощенный подход к решению прямой трехмерной задачи наталкивается на серьезные трудности, в частности, при значительных отклонениях проводимости локальной неоднородности от проводимости вмещающей среды. Это обусловлено самим подходом к вычислению аномального поля (то есть поля влияния изучаемой неоднородности)- линеаризованное решение будет близким к истинному только в достаточно малой окрестности базовой модели. Поэтому этот метод может давать приемлемые результаты только для не слишком контрастных по проводимости поисковых объектов.
Другой проблемой метода из работ [25, 26, 30, 35, 36] являются трудности, связанные с оценкой точности получаемого решения обратной задачи-после вычисления значений проводимости путем обращения линейной системы не ясно, каким образом можно оценить правильность полученного распределения проводимости в изучаемой среде. Эта проблема также связана со способом получения решения прямой трехмерной задачи в рассматриваемом подходе: используется только линеаризованное решение, хотя реальный электромагнитный процесс зависит от проводимости не линейно.
Поэтому актуальной является разработка новых подходов к интерпретации снятых на дневной поверхности данных для технологии с закрепленным источником, позволяющих восстанавливать распределение проводимости в исследуемом объеме, то есть разработка новых подходов к решению трехмерной обратной задачи. При этом способ решения прямой трехмерной задачи должен быть не только достаточно быстрым, но и достаточно точно описывать поле влияния реального трехмерного объекта даже при больших контрастах прово-димостей объекта по отношению к проводимости вмещающей среды, то есть должен учитывать нелинейную зависимость получаемого решения от проводимости объекта. Всем этим требованиям вполне удовлетворяет способ аппроксимации поля влияния трехмерного аномального объекта системой токовых колец и токового диска, рассматриваемый в данной диссертационной работе.
Наиболее распространенный способ исследований становлением поля предполагает использование симметричных приемно-генераторных конструкций. В нем измерительный датчик (петля) находится в центре генераторного контура. Эта технология получила широкое применение вследствие ограниченных возможностей математического аппарата интерпретации. До сих пор интерпретация данных этих методик основана на одномерных математических моделях, что не только существенно снижает точность интерпретации, но и иногда делает ее результаты совершенно неприемлемыми для практики. Поэтому актуальной является проблема создания новых методик интерпретации данных соосных профильно-площадных зондирований, основанных на трехмерных моделях строения среды.
Таким образом, актуальность данной диссертационной работы определяется: - необходимостью разработки максимально автоматизированного способа решения обратной задачи в рамках одномерной модели среды по экспериментальным данным, зарегистрированным в широком диапазоне времен и раз- носов, использующего новые быстрые методы математического моделирования электромагнитного поля в горизонтально-слоистой среде; потребностью изучения и создания новых подходов к интерпретации данных площадных электромагнитных зондирований, основанных на быстрых, но достаточно точных аппроксимациях полей влияния трехмерных объектов в горизонтально-слоистой среде; необходимостью создания удобных для пользователя программных средств трехмерной интерпретации профильно-площадных данных электроразведки становлением поля, позволяющих восстанавливать распределение проводимости в исследуемом объеме и локализовывать поисковые объекты.
В предлагаемой работе будут рассмотрены новые подходы к разработке быстрых процедур численного решения задачи становления электромагнитного поля в горизонтально-слоистой среде, с использованием которых будет разработана процедура автоматического построения начального приближения и последующего уточнения параметров при решении обратной задачи в рамках горизонтально-слоистой среды.
Много внимания в данной работе будет уделено исследованию эффективности аппроксимации полей влияния трехмерных объектов, получаемых с использованием трехмерного математического моделирования, системой токовых колец и токовых дисков, которые позволяют довольно точно вычислять отклики трехмерных объектов с помощью численного анализа двумерных (осе-симметричных) полей.
Основное внимание в диссертационной работе будет уделено разработке новых методик интерпретации данных, основанных на трехмерных моделях строения среды, как для технологии с закрепленным источником, так и для со-осных профильно-площадных зондирований. Эффективность разрабатываемых методик интерпретации, а также достоверность получаемых результатов будет проверена как на модельных, так и на практических задачах.
Таким образом, основной научной проблемой, решению которой посвящена данная диссертационная работа, является проблема построения эф- фективных процедур интерпретации данных импульсной индуктивной электроразведки с контролируемым источником электромагнитного поля, основанных на трехмерных моделях строения среды и использовании быстрых численных схем расчетов полей влияния трехмерных объектов на базе их двумерных приближений.
Цель исследования состоит в разработке методик и реализации быстрых и достаточно точных алгоритмов конечноэлементного моделирования электромагнитного поля в горизонтально-слоистой среде; в создании процедур автоматического построения начального приближения и последующего уточнения параметров при решении обратной задачи в рамках горизонтально-слоистой модели среды; в подтверждении эффективности аппроксимации полей влияния трехмерных объектов полями токовых колец и токовых дисков, получаемых как решения осесимметричных задач; в разработке методик и реализации эффективных алгоритмов интерпретации данных импульсной индуктивной электроразведки с контролируемым источником электромагнитного поля.
Научная новизна работы состоит в следующем:
Предложен новый подход к решению задачи о становлении электромагнитного поля магнитного дипольного источника в горизонтально-слоистой среде, базирующийся на использовании конечноэлементных аппроксимаций и методики выделения основного поля.
Для решения проблемы получения значений ЭДС на ранних временах и больших разносах предложен подход, основанный на переходе к решению двумерной задачи в осесимметричной постановке на тех временах, для которых не удается достаточно точно вычислять ЭДС на основе решения одномерных задач.
Разработаны процедуры автоматического построения начального приближения для решения задачи восстановления параметров горизонтально-слоистой среды.
Исследована возможность аппроксимации поля влияния трехмерного объекта полями модельных токовых колец и дисков, которые могут быть получены как решения соответствующих осесимметричных задач.
Разработана процедура восстановления проводимости среды по данным импульсной индуктивной электроразведки как для технологии с закрепленным источником, так и для соосных профильно-площадных зондирований, основанная на аппроксимации полей влияния трехмерных объектов полями модельных токовых колец и дисков.
Практическая ценность работы и реализация результатов
Разработанные методы и алгоритмы конечноэлементного моделирования электромагнитного поля в горизонтально-слоистой среде, а также созданные процедуры автоматического построения начального приближения и последующего уточнения параметров при решении обратной задачи в рамках горизонтально-слоистой модели среды реализованы в программном комплексе HORIZON. Разработанные процедуры восстановления проводимости среды по данным импульсной индуктивной электроразведки как для технологии с закрепленным источником, так и для соосных профильно-площадных зондирований, основанные на аппроксимации полей влияния трехмерных объектов полями модельных токовых колец и дисков, реализованы в программном комплексе IFSTEM. Программные комплексы HORIZON и IFSTEM использовались при решении исследовательских и практических задач индукционной электроразведки.
Достоверность результатов подтверждена решением ряда модельных задач, сравнением полученных результатов решения некоторых задач с результатами их решения с помощью программных комплексов, успешно зарекомендовавших себя при решении такого рода задач, а также с результатами натурных экспериментов.
Защищаемые положения
1. Предлагаемая процедура решения задачи о становлении электромагнитного поля магнитного дипольного источника в горизонтально-слоистой ере- среде, базирующаяся на методике выделения основного поля, в комбинации с численным методом решения этой задачи в двумерной постановке позволяет получать решение в любом диапазоне времен и для любых разносов с высокой точностью и относительно небольшими вычислительными затратами.
Разработанная процедура автоматического построения начального приближения и последующее уточнение параметров горизонтально-слоистой модели с использованием модификации метода Ньютона позволяет достаточно быстро и с хорошей точностью восстанавливать параметры горизонтально-слоистой среды.
Схема аппроксимации полей влияния трехмерных объектов полями токовых колец и дисков, которые могут быть вычислены как решения сформулированных специальным образом двумерных задач, позволяет добиться хорошего совпадения реального поля влияния трехмерного объекта с модельным полем токовых колец и дисков в широком диапазоне времен, в котором поле влияния поискового объекта проявляется достаточно четко на фоне поля вмещающей среды.
Разработанная итерационная процедура восстановления проводимости среды по данным импульсной индуктивной электроразведки, основанная на аппроксимации полей влияния трехмерных объектов полями модельных токовых колец и дисков, позволяет достаточно быстро и с высокой точностью определять местоположение поисковых объектов во вмещающей горизонтально-слоистой среде как для технологии с закрепленным источником, так и для со-осных профильно-площадных зондирований.
Личный вклад
Разработаны и программно реализованы схемы конечноэлементного моделирования электромагнитного поля магнитного диполя в горизонтально-слоистой среде, основанные на решении задачи в постановке А.Н. Тихонова [79]. Для решения проблемы получения значений ЭДС на ранних временах и больших разносах разработан и реализован подход, основанный на переходе к решению двумерной задачи в осесимметричной постановке на тех временах, где не удается достаточно точно вычислять ЭДС полуаналитическим методом. Разработана процедура автоматического построения начального приближения при решении задачи восстановления параметров слоистой среды. Для уточнения полученных на начальном этапе параметров горизонтально-слоистой модели разработана модификация метода Ньютона. Исследована эффективность схемы аппроксимации полей влияния трехмерных объектов полями модельных токовых колец и дисков, которые получаются двумерным математическим моделированием. Для решения трехмерной обратной задачи разработана итерационная процедура восстановления свойств среды с использованием регуляризации, основанная на аппроксимации полей влияния трехмерных объектов полями модельных токовых колец и дисков.
По результатам исследований, представленных в работе, разработаны программные комплексы HORIZON и IF STEM.
Апробация работы
Основные результаты работы были представлены на IV международной конференции «Актуальные проблемы электронного приборостроения» АПЭП-98 (Новосибирск, 1998г.); The Third Russian-Korean International Symposium on Science and Technology KORUS'99 (Новосибирск, 1999г.); Четвертом сибирском конгрессе по прикладной и индустриальной математике, посвященном памяти С.Л.Соболева (Новосибирск, 2000г.); на международной геофизической конференции и выставке EAGE (Нидерланды, 2001г.); региональной научной конференции «Наука, Техника, Инновации» (Новосибирск, 2001,2002гг.); на научно-практической конференции «Инновационные технологии в области поисков, разведки и детального изучения месторождений нефти и газа» (Москва, 2002г.); на семинарах НГТУ и СНИИГГиМС. Результаты проведенных исследований включались в отчеты по НИР СНИИГГиМС.
Публикации
По результатам выполненных исследований опубликовано 14 печатных работ.
Структура работы
Диссертация состоит из введения, пяти глав, заключения и списка использованных источников (114 наименований). Работа изложена на 172 страницах, включая 65 иллюстраций.
Краткое содержание работы
В первой главе рассматривается задача о становлении электромагнитного поля в горизонтально-слоистой среде. В п. 1.1 приводится математическая модель поля вертикального магнитного диполя, предложенная Тихоновым А.Н. [79], а в п. 1.2 рассматривается решение соответствующей одномерной начально-краевой задачи с использованием стандартного МКЭ. В п. 1.3 рассматривается алгоритм численного вычисления значений дВ2 / dt. Для этого анализируется поведение поля на примере двухслойной среды и формируются требования для алгоритма вычисления значений dBzldt с целью уменьшения суммарных вычислительных затрат. В п. 1.4 рассматривается подход к решению одномерной краевой задачи, к необходимости решения которой сводится постановка Тихонова А.Н., основанный на выделении основного поля [52-54, 57, 69, 70], который позволяет существенно снизить вычислительные затраты на получение численного решения задачи. В качестве основного предлагается использовать поле от полупространства с фиксированной проводимостью. При этом, для обеспечения максимальной близости основного поля к искомому, разработан алгоритм выбора значения проводимости полупространства, через которую в свою очередь вычисляются значения основного поля простой заменой переменных. Рассмотрены часто встречающиеся на практике ситуации, когда приповерхностные слои являются слабо проводящими по сравнению с более глубоко залегающими, а также когда значение проводимости подстилающего основания значительно меньше проводимости слоев, лежащих выше, то есть когда задан непроводящий фундамент. В этих случаях использование в качестве основного поля значений поля, рассчитанного для полупространства, становится малоэффективным. Показано, что решение первой проблемы естественным образом описывается краевым условием, заданным на поверхности среды. В результате изменяется начальное приближение задачи на добавочное поле. Для решения второй проблемы предложено использование в качестве основного поля значений поля, рассчитанного для среды, состоящей из однородного слоя некоторой толщины и бесконечного фундамента. Показано, что предлагаемый подход в несколько десятков раз эффективнее стандартной ко-нечноэлементной схемы решения этой задачи. В п. 1.5 приводится достаточно эффективный способ вычисления ЭДС для приемно-генераторной установки с петлями произвольной формы путем интегрирования значений dBz/dt по площади генератора и приемника. В п. 1.6 для решения проблемы получения значений ЭДС на ранних временах и больших разносах предлагается подход, основанный на решении двумерной задачи в осесимметричной постановке на тех временах, где не удается достаточно точно вычислять ЭДС полуаналитическим методом. Использование комбинации полуаналитического и численного методов позволяет сочетать преимущества этих методов на разных временах становления поля. Для решения двумерной задачи с некруглым генераторным контуром предлагается два способа аппроксимации генератора: при помощи суперпозиции полей от магнитных диполей и при помощи набора круглых петель достаточно малого радиуса. В результате проведенного исследования показано, что более эффективным является способ аппроксимации некруглого генератора с использованием набора круглых петель. Для этого способа (в случае, когда генераторный контур является квадратом) выведены значения оптимальной расчетной площади приемника (с точки зрения минимизации ошибки аппроксимации квадратного генератора набором круглых).
Во второй главе рассматривается задача восстановления параметров горизонтально-слоистой среды. В п.2.1 проводится анализ методик трансформации результатов наблюдения в кажущиеся параметры среды на примере одного из наиболее распространенных способов трансформации измеряемого сигнала в кажущееся удельное сопротивление. В п.2.2 приводится способ автоматического построения начального приближения, реализованный в виде последовательности обратных задач в разных классах моделей, сходящуюся к исходной задаче. Проводится исследование сходимости метода при построении начального приближения с различной точностью. Для уточнения параметров модели используется модификация метода Ньютона, позволяющая с высокой скоростью улучшать параметры разреза (в смысле среднеквадратичного отклонения от экспериментальной ЭДС), построенные на первом этапе решения обратной задачи. Основное ускорение достигается при вычислении производных по параметрам с использованием методики выделения основного поля. Проводится исследование сходимости метода Ньютона из разных начальных приближений. В третьей главе разрабатывается новый подход к решению задачи восстановления проводимости среды по данным площадных зондирований, полученным с использованием технологии с закрепленным источником. В п.3.1 рассматриваются некоторые трудности применения методик интерпретации, основанных на одномерной модели строения среды, по данным, полученным трехмерным математическим моделированием поля кимберлитового объекта, залегающего в горизонтально-слоистой среде. В п.3.2 предлагается новый способ аппроксимации поля влияния трехмерного объекта, основанный на решении ряда осесимметричных задач с токами в дисках и усеченных полых цилиндрах. На примерах численного моделирования полей влияния трехмерных объектов показывается эффективность предлагаемого способа аппроксимации. В п.3.3 рассматриваются основные идеи предлагаемого способа интерпретации данных для технологии с закрепленным источником, а в п.3.4 приводится итерационная процедура восстановления проводимости среды с использованием регуляризации. В п.3.5 приводится пример восстановления проводимости среды на данных, полученных трехмерным математическим моделированием поля кимберлитового объекта. В пп.3.6-3.7 проводится анализ работоспособности предлагаемой процедуры интерпретации при несовпадении размеров и положения «пробных» объектов с размерами и положением поисковых, а также при неточных параметрах вмещающей горизонтально-слоистой среды. Показано, что размеры и положение «пробных» объектов, а также точность параметров вмещающей горизонтально-слоистой среды достаточно существенно влияют как на получаемое распределение проводимости, так и на значение функционала ошибки. Однако минимальное значение функционал ошибки принимает при наилучшем совпадении размеров и положения «пробных» объектов с размерами и положением искомого объекта, а также при максимально точных параметрах вмещающей среды. В п.3.8 приводится серия примеров применения предлагаемой процедуры интерпретации на данных, полученных трехмерным математическим моделированием. Демонстрируется высокая достоверность получаемого трехмерного распределения проводимости. В п.3.9 приводятся примеры решения практических задач. Достоверность распределений проводимости, полученных по экспериментальным данным, подтверждается данными бурения.
В четвертой главе разрабатывается подход к решению задачи восстановления проводимости среды по данным профильно-площадных зондирований, полученным с использованием технологии с коаксиальными установками, когда измерительная петля перемещается вместе с генераторной петлей. В п.4.1 проводится исследование общепринятых методик интерпретации, основанных на одномерных моделях строения среды, на данных трехмерного математического моделирования. Показано, что использование этих методик не позволяет различать эквивалентные по отклику среды. В п.4.2 для решения этой проблемы предлагается использовать методику интерпретации, основанную на трехмерной модели строения среды, с использованием разработанной в гл.З итерационной процедуры с регуляризацией. Показывается возможность аппроксимации полей влияния трехмерных объектов полями модельных токовых колец и дисков, которые могут быть получены как решения соответствующих осесимметричных задач. В п.4.3 на серии примеров демонстрируется достоверность получаемого трехмерного распределения проводимости в результате применения разработанной процедуры интерпретации на данных, полученных трехмерным математическим моделированием.
В пятой главе описываются разработанные по результатам предложенных в диссертационной работе подходов программные комплексы, предназна- ченные для интерпретации профильно-площадных данных импульсной индуктивной электроразведки с контролируемым источником электромагнитного поля HORIZON и IFSTEM. В п.5.1 приводится описание возможностей комплекса HORIZON для обработки, исследования и интерпретации данных в классе горизонтально-слоистых моделей. В п.5.2 проводится сравнение возможностей программного комплекса HORIZON с существующими разработками. В п.5.3 приводится описание возможностей программного комплекса IFSTEM для решения трехмерной обратной задачи электроразведки становлением поля.
В заключении приводятся основные результаты диссертационной работы.
Применение методики выделения основного поля
При решении задачи (1.3)-(1.6) методом конечных элементов напрямую для достижения необходимой точности требуется довольно подробная сетка по пространству и времени. Для ускорения расчета применим подход, основанный на выделении основного поля [52-54, 57, 69, 70]. Представим функцию Z(m,z,t), являющуюся решением начально-краевой задачи (1.3)-(1.6) при фиксированном значении параметра т, в следующем виде: где функция Z является решением задачи (1.3)-(1.6) при O = J . Подставим (1.32) в (1.3)-(1.6), в результате получим краевую задачу для Z: Отметим, что при решении задачи (1.33)-(1.36) значение Z считается известным. Так как a(z) - кусочно-постоянная функция, будем считать а также кусочно-постоянной. Более того, примем a = const, z 0, то есть будем считать, что функция Z задана для полупространства с проводимостью a . Выбор СУ в виде a = const, z 0 обусловлен простотой вычисления Z и ее производной по времени, а также универсальностью и простотой использования значения а в качестве приближения функции c(z). Наилучший выбор значения а будет зависеть как от значений параметров среды, так и от значения параметра т, что делает невыгодным априорное задание a . Фактически, требуется наилучшим образом приблизить кусочно-постоянную функцию a(z) константой а в расчетной области для каждого т, причем нужно выбрать такое значение a , которое обеспечило бы максимальную близость Z к Z. При этом следует учитывать то обстоятельство, что при больших т на функцию Z оказывают большее влияние приповерхностные слои среды, а с уменьшением т - более глубокие и более мощные. Для этого введем понятие «эффективной глубины» Н при фиксированном значении параметра т, которая должна показывать до какой глубины следует учитывать функцию a(z) в процессе выбора а (с необходимой точностью).
После проведения ряда расчетов было установлено, что для полупро странства, то есть для O[Z) = JS= const, z О, значение Н « 3/ т и практиче р ски не зависит от значения os. Однако, в многослойных средах значение Н невозможно задать априорно. Поэтому поступим следующим образом для оп ределения Н , а следовательно и ст . На первом этапе примем в качестве Н значение Н =Ъ1т (для полупространства) и т как среднее о по слоям. Затем преобразуем мощности слоев hj как а и определим какой преобразованной глубине соответствует Н , то есть сколь-ко преобразованных слоев попадает в Н , что равносильно выполнению где Ш -І) - число слоев, полностью лежащих в 0,Я , Ah Е - часть мощности последнего слоя, частично лежащего в О,Я (рис. 1.4). р Окончательное значение Н (через обратное преобразование hj): Таким образом, мы приблизительно определили Н для заданной слоистой среды. Значение о вычислим как среднюю продольную проводимость, соответствующую глубине И :
Теперь обратимся к вычислению функции Z . Покажем, что Z может быть получена из функции Z, являющейся решением задачи (1.3)-(1.6) при а = 1, заменой переменной / = от. Действительно, в этом случае справедливы соотношения: причем для функции Z справедливо уравнение (получаемое из (1.3) заменой а на о0 и Z на Z ): Учитывая (1.39) и (1.40), преобразуем (1.41) к виду: то есть функция Z удовлетворяет уравнению (1.3) при a = 1. Легко видеть, что Z удовлетворяет также (1.4)-(1.6). Преимущество постановки (1.33)-(1.36) заключается в том, что, зная вы численную заранее с высокой точностью функцию Z, при вычислении Z можно использовать гораздо менее подробные пространственную и временную сетки, чем это необходимо для решения с той же точностью краевой задачи (1.3)-(1.6) напрямую, и это позволит существенно сократить время решения основной задачи. В то же время значения функции Z и ее производной dZ Idt могут быть получены по значениям функции Z и ее производной dZldt на основа нии соотношений (1.39) и (1.40), а значения Z и dZ/dt могут быть заранее вычислены, на подробной сетке и сохранены в файле. Заметим однако, что достаточно хранить значения Z и dZldt только для одного значения т (например, для т = 1). Это связано с тем, что функция Z при тф\ может быть получена из функции Z, являющейся решением (1.42) при т -1, заменой переменных м = т t и x = mz. Действительно, в этом случае справедливы соотношения то есть Z удовлетворяет (1.46) при m = 1. Таким образом, для вычисления dBzldt прежде всего необходимо создать файл, в котором будут сохранены полученные заранее значения функции Z в точке z = 0 (так как нас интересует значение dBzldt на дневной поверхности) и производной dZldt во всей расчетной области. Затем для каждого значения т, которое используется при численном интегрировании в (1.16) нуж но получить значения функции Z в точке z = 0 и производной dZldt по формулам (1.43) и (1.44); выбрать такое значение ст , которое обеспечило бы максимальную близость Z к Z, то есть воспользоваться (1.37)-(1.38).
После этого по формулам (1.39) и (1.40) могут быть вычислены Z и а затем для каждого т решена задача (1.33)-(1.36) и получены при z = 0 значения функции Z и ее производной dZldt, которые используются в интеграле (1.16). В тех случаях, когда приповерхностные слои являются слабо проводящими по сравнению с более глубоко залегающими, например, а =0.01 См/м, а2 =1 См/м, решение задачи в постановке (1.33)-(1.36) осложняется тем обстоятельством, что при использовании основного поля, рассчитанного для полупространства, для достижения необходимой точности решения нужно значительно увеличивать пространственную и временную сетки. Это связано с тем, что функция T(Z) плохо приближается константой а (слишком большой контраст проводимостей у слоев), а следовательно при решении задачи в постановке (1.33)-(1.36) в правой части уравнения (1.33) возникает мощный источник, для учета которого потребуется довольно подробная пространственно-временная сетка. В этом случае поступим следующим образом.
Автоматическое восстановление параметров горизонтально-слоистой среды
В основе разработанного автоматического построения начального приближения лежит метод послойного восстановления параметров слоистой модели среды [93, 95]. Выбор именно этого метода обусловлен, прежде всего, характерной особенностью зависимости ЭДС от времени (отклика) для горизонтально слоистой среды: влияние более глубоких слоев на суммарный отклик очень сильно зависит от свойств приповерхностных слоев. Поэтому если параметры приповерхностных слоев определены с высокой погрешностью, восстанавливать параметры подстилающих слоев не имеет смысла. Итак, пусть имеется экспериментальная ЭДС Є (7), заданная в моменты времени tfc, к = l,...,K, полученная при помощи некоторой (известной) прием-но-генераторной установки. Задача построения начального приближения состоит в нахождении параметров среды, включающих число слоев N, мощности / и проводимости а,-(или сопротивления р,) каждого из слоев, i = \,...,N, с некоторой заданной точностью с„. Точность определения начального приближения гр может быть выбрана гораздо хуже, чем требуемая точность einv совпадения экспериментальной кривой с восстановленной, которая должна быть получена в результате решения обратной задачи. Однако начальное приближение должно быть достаточно близким к истинному решению, чтобы в дальнейшем метод многомерной оптимизации смог сойтись в глобальный минимум с точностью гіт. На первом этапе построения начального приближения получим зависимость кажущегося удельного сопротивления px(t) от времени для заданной экспериментальной кривой Є (ґ) способом, описанном выше. Далее примем искомую модель среды в виде полупространства (N = l, pj = рх (/j), /zj =со) и уточним значение pj одномерной минимизацией по параметру pj. После этого определим, нужно ли добавлять в модель следующий (второй) слой. Если нужно - добавим его и начнем оптимизацию по параметрам новой модели (двухслойной) по следующей схеме. Вначале подберем fy, затем р2, затем вернемся к P] и так далее, пока не получим удовлетворительный результат.
В качестве минимизируемого функционала при оптимизации по одному из параметров модели будем использовать где pf - текущий оптимизируемый параметр, n\zp\- функция, определяемая соотношением Фактически, функция п(е„] определяет число точек по времени, восстановленных на текущем этапе с точностью є „. Для определения необходимости добавления следующего слоя в модель достаточно провести следующие вычисления. В текущую восстановленную модель из N слоев добавим следующий слой с сопротивлением PJY и вычислим аномальную ЭДС от этого нового слоя. Если переход через ноль полученной зависимости аномальной ЭДС от времени лежит в текущем восстановлен ном временном диапазоне, то есть в диапазоне , то будем считать, что новый слой действительно необходим. Значение рдг приблизительно определим по кривой отличия восстанов то pjy должно быть выставлено больше сопротивления последнего восстановленного слоя pN, иначе - меньше. Процесс добавления новых слоев, то есть увеличения значения N, и оптимизации по тройкам параметров {hN_\, pN, р#_і} повторяется до тех пор, пока не будет достигнута требуемая точность є „. В случае, когда определено, что добавлять в модель следующий слой нет необходимости (его влияние незначимо), происходит возврат к оптимизации ранее полученных параметров, начиная с первого слоя. Если и после этой процедуры не достигается требуемая точность с„, значение ер увеличивается до некоторого sp, и процесс повторяется (с использованием полученной модели)
Этим достигается предотвращение «зацикливания» алгоритма вследствие увеличения значения функции n\Zp\, а следовательно усиления влияния на значение минимизируемой функции более глубоких слоев. После построения модели с увеличенным значением є„, происходит возврат к старому значению гр и оптимизация повторяется на основе модели, восстановленной с точностью є„. Опишем процесс выбора значений zp. Если пользователь указал желае мую точность построения начального приближения гр, то для ее достижения поступим следующим образом. На первом этапе (после вычисления зависимо и так далее, пока не дойдем до є „. Фактически, такой процесс означает регуляризацию решения обратной задачи, состоящий в ограничении на каждом этапе, соответствующем значению zJp, класса моделей, в котором происходит восстановление параметров среды [19, 43, 50, 76, 77, 90, 92, 102]. Это ограничение происходит вследствие изменения значения zJp от большего к меньшему, так как большему значению zJp соответствует более узкий класс моделей (с меньшим числом слоев). В результате мы получаем последовательность обратных задач в разных классах моделей, соответствующих значениям Бд, сходящуюся к исходной, то есть к задаче с требуемой точностью гр. Скорость сходимости этой последовательности определяется значением коэффициента регуляризации а [19, 43, 50, 76]. Важное место в процедуре построения начального приближения занимает выбор параметра регуляризации а.
Очевидно, что если а слишком мал, эффект стабилизации может оказаться недостаточным и практическая неустойчивость решения будет еще значительна. Если а слишком велик, то получим слишком длительный процесс решения всей обратной задачи, то есть будет низкая скорость сходимости последовательности задач [19, 43, 50, 76]. На рис. 2.4 представлены графики отличия (в процентах) восстановленных с разной точностью ЭДС от экспериментальной. В качестве экспериментальной использована теоретическая ЭДС для среды с параметрами \ =100м, Р! = 100 Ом л/, h2 = 20м, р2 = 10 Ом -м, / = , р3 = 500 Ом м . Генератор -квадрат со стороной 200 м, приемник - точечный. Различные графики соответствуют разным значениям zJp. В таблице 2.1 представлены параметры среды, восстановленные описанным методом с разной погрешностью. Заметим, что на практике экспериментальная ЭДС задана с некоторой погрешностью 5. Поэтому точность решения обратной задачи (и построения начального приближения в частности) не может быть выше точности начальных данных. Так как процесс построения начального приближения основан на методе послойного восстановления параметров слоистой среды, это довольно длительный процесс. Более того, в случае необходимости восстановления параметров модели по зондированиям в близких по геологическим параметрам средах (с условием, что количество слоев одинаково и параметры слоев достаточно близки), использование описанного выше алгоритма становится неэффективным.
Примером такого рода ситуации может служить обработка профильно-площадных данных, полученных с использованием соосной установки, когда необходимо получить параметры одномерной модели в близких (по расположению на поверхности Земли) точках [14, 26, 30, 81-83].
Анализ возможности построения системы интерпретации данных для технологии ПЗС-ЗИ
Будем считать, что электромагнитное поле возбуждается импульсным изменением тока в круговой петле. Такое поле может быть описано системой уравнений [64, 69, 70, 113] где A - вектор-потенциал электромагнитного поля, V — скалярный потенциал электрического поля, J— сторонние токи, а- проводимость среды, JJ.Q- магнитная проницаемость вакуума. Индукция магнитного поля В и напряженность электрического поля Ё связаны с потенциалами AYLV соотношениями Если обозначить через а = а (z) проводимость нормальной (референтной) среды (то есть под нормальной (референтной) средой понимается вмещающая горизонтально-слоистая среда, описываемая кусочно-постоянной функцией a (z)), то поле А, удовлетворяющее системе уравнений (3.1)-(3.2), может быть найдено в виде суммы нормального поля А в референтной среде с проводимостью а и добавочного поля А+ [64, 69, 70, 113]. У нормального поля А в цилиндрической системе координат отлична от нуля только р я компонента. Эта компонента является функцией пространственных координат г и z и времени t (т.е. A, = Ay(r,z,t)) и удовлетворяет уравнению [64, 69, 70, 113] Добавочное же поле А+ удовлетворяет системе уравнений т.е. источником добавочного поля является величина (а-а )Ё , где Ё = , и вектор-потенциал А определяется функцией А (г,z,t), являю щейся решением уравнения (3.3). Таким образом, аномальное поле Ы+ М фактически возбуждается токами Ja = -а agradF + (а-а \Ё , которые можно считать аномальными токами, поскольку они являются разностью сум марных токов Js =-ст agradF, текущих в исследуемой среде с проводи мостью а, и нормальных токов J =а Ё , текущих в референтной среде с проводимостью а . Рассмотрим структуру вектора плотности аномальных токов Ja. Аномальные токи порождаются членом (а-а \Ё , являющимся произведением напряженности нормального поля Е на разность проводимостеи исследуемой и референтной сред.
Этот член также является составной частью аномальных токов внутри объектов с измененной относительно референтной среды проводимостью, но вне этих объектов аномальные токи формируются только за счет во времени, и компоненты -agradF, порождаемой электрическими зарядами на поверхностях разрыва проводимости а. Эти заряды делают поле аномальных токов Ja удовлетворяющим уравнению неразрывности (3.5). Таким образом, аномальные токи Ja фактически формируются пропорциональными нормальному полю Ё внутри объектов с аномальной проводимостью.
Направлено ние этих токов определяется в основном направлением поля Е внутри аномального объекта и направлением поля -gradF вне его. При этом, как уже было сказано выше, компонента -agradF поля аномальных токов формируется в основном таким образом, чтобы вектор плотности полных аномальных токов удовлетворял уравнению неразрывности, то есть на любой границе его нормальная (перпендикулярная этой границе) составляющая была непрерывна. А поскольку член (СУ-а \Ё , являющийся составной частью аномального тока, имеет разрывную нормальную составляющую на границах объекта с аномальной проводимостью, то основная часть поля -gradK формируется так, чтобы компенсировать этот разрыв в нормальной составляющей поля (а-СУ \Ё на границах аномального объекта (некоторая часть поля -agradF компенсирует разрыв нормальной составляющей поля -о , и принцип ее формирования и воздействия точно такой же).
В результате для объектов с повышенной аномальной проводимостью поле -gradF создает такие токи, которые в окрестности границы объекта с перпендикулярными к ней токами о Ё имеют совпадающее с Ё направление вне объекта и противоположное Ё направление внутри объекта. То есть основная часть поля -gradF формируется таким образом, чтобы обеспечить непрерывный выход суммарных токов за границу объекта. Таким образом, основная часть поля -gradF внутри объекта с повышенной аномальной проводимостью фактически ослабляет вклад компоненты в поле суммарных токов. Для объектов же с пониженной аномальной проводимостью в окрестности границы объекта с перпендикулярными ей нормальными токами а Ё основная часть поля -gradF, наоборот, совпадает с направлением нормального поля Ё внутри объекта и имеет противоположное Ё направление вне объекта, и фактически эта часть электрического поля компенсирует разрыв в компоненте (а-а )Ё аномальных токов на рассматриваемой границе усилением суммарных токов внутри объекта и их ослаблением вне объекта.
Примеры решения трехмерных обратных задач на синтетических данных
В качестве тестовой модели рассмотрим проводящее и непроводящее замещение проводящего слоя в слоистой среде - имитация изменения типа флюидонасыщения, одного из важнейших признаков при поиске углеводородов, решения промысловых и экологических задач. На рис. 4.9а,б приведены две модельные среды с противоположным типом замещения. На рис. 4.9а -контрастное непроводящее замещение, а на рис. 4.96 - проводящее.
Анализ глубинного разреза с использованием одномерной трансформации сигнала в кажущуюся удельную проводимость [14] дает весьма приблизительную ориентировку о присутствии локальных неоднородностей вблизи инспектируемого слоя - рис. 4.9в,г. 3D реконструкция с весьма высокой точностью позволяет повторить геометрию этих неоднородностей - рис. 4.9д,е. Рассмотрим ситуацию, когда искомый кимберлитовый объект имеет размеры 100 х 100 м2 в плане и залегает на глубине 20 м в слоистой среде с пара метрами: h} = 20 м, GI = 0.001 См/м, h2 = да, а2 = 0.002 См/м. На рис. 4.10а приведена структура описанной модели. Объект является слабо проводящим относительно вмещающей среды (стобъекта = 0.005 См/м). Профилирование происходит с соосной установкой (генератор - радиусом 141м, приемник - радиусом 0.564 м) с шагом 50 м вдоль профиля через эпицентр объекта. На рис. 4.10б,в представлены результаты интерпретации данных профилирования над проводящим кимберлитовым объектом: на рис. 4.106- результат одномерной трансформации измеряемого сигнала в кажущуюся проводимость, а на рис. 4.1 Ов- результат восстановления проводимостей «пробных» метрами: h} = 20 м, GI = 0.001 См/м, h2 = да, а2 = 0.002 См/м.
На рис. 4.10а приведена структура описанной модели. Объект является слабо проводящим относительно вмещающей среды (стобъекта = 0.005 См/м). Профилирование происходит с соосной установкой (генератор - радиусом 141м, приемник - радиусом 0.564 м) с шагом 50 м вдоль профиля через эпицентр объекта. На рис. 4.10б,в представлены результаты интерпретации данных профилирования над проводящим кимберлитовым объектом: на рис. 4.106- результат одномерной трансформации измеряемого сигнала в кажущуюся проводимость, а на рис. 4.1 Ов- результат восстановления проводимостей «пробных» Теперь рассмотрим ситуацию, когда искомый кимберлитовый объект имеет те же размеры, но является непроводящим относительно вмещающей среды (субъекта= 0.001 См/м) и залегает на глубине 20 м в слоистой среде с параметрами: П! = 20м, Gi = 0.001 См/м, п2 = оо, 02 = 0.05 См/м. На рис. 4.11а приведена структура описанной модели. Как и ранее профилирование происходит с соосной установкой (генератор - радиусом 141 м, приемник - радиусом 0.564 м) с шагом 50 м вдоль профиля через эпицентр объекта. На рис. 4.11б,в представлены результаты интерпретации данных профилирования над непроводящим кимберлитовым объектом: на рис. 4.116- результат применения одномерной трансформации измеряемого сигнала в кажущуюся проводимость, а на рис. 4.11в - результат восстановления проводимости при помощи «пробных» объектов. Как и ранее, объект слабо выделяется на фоне вмещающей среды при использовании одномерной трансформации.
В результате применения трехмерной интерпретации расположение объекта отслеживается с достаточной точностью. При этом, как и в случае с проводящим объектом, основная аномальная проводимость восстановлена в центре искомого объекта по глубине. Значение функционала ошибки при этом уменьшилось более чем в 3 раза. Рассмотрим ситуацию, когда три объекта с проводимостями СТобьектаІ = 0.1 См/м, СГ0бъекта2 = 0.1 См/м И ОобъеетаЗ = Ю См/м Залегают В ГОрИЗОН тально слоистой среде с параметрами: hi = 50 м, Сті = 0.02 См/м, 1 = 850 м, а2 = 0.0005 См/м, h3 = 300 м, а3 = 0.01 См/м, 1ц = да, а4 = 0.002 См/м. Первый и второй объекты (помехи) с размерами 1000x1000x300 м залегают на глубине 900 м (замещают третий слой). Эпицентр первого объекта расположен в точке (500, 1000), эпицентр второго - в точке (1500,-2000). Третий объект с размера-ми 3000x1000x20 м залегает на глубине 1630 м и моделирует поисковую рудную залежь. Профилирование проводится с соосной установкой: генераторный контур радиусом 564 м (с моментом 10 А-м ), приемник - радиусом 56.4 м. На рис. 4.12 приведена структура описанной модели. На рис. 4.13 приведены разрезы кажущейся удельной проводимости, полученные путем одномерной трансформации суммарного сигнала. На рис. 4.14 приведены результаты восстановления проводимости с помощью «пробных» объектов. Как видно из представленных результатов, объект, моделирующий рудную залежь, восстановлен несколько искаженным в плане: размер (в плоскости XY) центральной части объекта немного (примерно на 20 %) больше, чем размеры боковых. Это объясняется тем, что на поле в центре объекта оказывают достаточно сильное влияние его боковые части, а также тем, что мы не учитываем взаимоиндукцию электромагнитных полей, порождаемых различными объектами. Однако полученной точности восстановления этого объекта вполне достаточно для практики. Таким образом, представленные на рис. 4.14 результаты подтверждают эффективность предлагаемой процедуры восстановления проводимости исследуемой среды при наличии в ней нескольких поисковых объектов. Значение функционала ошибки при этом уменьшилось более чем в 5 раз.
