Повышение точности определения поправок за влияние рельефа при гравиметрической съемке

Размещено на http://www.allbest.ru/

Размещено на http://www.allbest.ru/

Повышение точности определения поправок за влияние рельефа при гравиметрической съемке

ВВЕДЕНИЕ

погрешность высота топографический рельеф

Вычисление поправок за влияние рельефа местности gр является необходимой процедурой при обработке материалов гравиметрических съемок. Величина поправки gр характеризует аномальный эффект, обусловленный отклонениями физической поверхности Земли от плоскопараллельного слоя, использующегося при определении поправки Буге. В случаях, когда вычисление аномалий Буге выполняется с учетом шарообразной формы Земли, плоскопараллельный слой заменяется на плоскосферический [3].

Разработаны различные способы расчета gр, реализующиеся с помощью палеток и номограмм, а также на ЭВМ [2, 3, 5, 6 и др.]. Как правило, область учитываемого влияния рельефа D разбивается на внутреннюю подобласть D1 и внешнюю подобласть D2 [6]. Внутренняя подобласть D1 (центральная зона) охватывает пункт гравиметрических наблюдений и его ближайшие окрестности (минимальный радиус ее на топографической карте не менее 0.4-0.5 см). Внешняя подобласть D2 обычно подразделяется на три зоны (ближнюю, среднюю и дальнюю), характеризующиеся различным шагом сети задания высотных отметок z.

Цифровые модели местности (ЦММ), использующихся в качестве исходных данных при расчете gр, представляют собой матрицы высотных отметок z размером m строк, n столбцов, заданных в узлах равномерной сети точек: z = z(x0+ix, y0+jy); x = y =s; i = 0, 1,.., n; j = 1,2,..,m. Ранее ЦММ формировались вручную: как правило, значения высот снимались по сети 11 см или 0.50.5 см в масштабе используемой топографической карты.

Наиболее рациональной методикой подготовки ЦММ, применяющейся в настоящее время во многих организациях, является сканирование топографических карт и векторизация полученных скан-образов с помощью специализированных программ (EasyTrase, R2V и т.п.). Это позволяет сравнительно быстро создавать «большие ЦММ», сохраняющие в цифровой форме все особенности крупно- и среднемасштабных топографических карт. При этом в результате векторизации карт изогипс рельефа дневной поверхности создаются файлы исходных данных (например: shp- или dat-файлы), содержащие порядка n(105106), n = 1,2,…,10 векторов вида = {x,y,z}, а средняя плотность сети высотных отметок составляет примерно 200300 точек/см2 в масштабе карты (рис. 1).

Совершенно очевидны существенные различия в степени детальности описания особенностей рельефа дневной поверхности при ручном и автоматизированном способах создания ЦММ. Во втором случае можно примерно на порядок уменьшить шаг сети s результативной матрицы z при более высокой точности определения значений высот.

Рис. 1. Результаты векторизации рельефа местности, представленного на топографической карте масштаба 1:25000.

А - исходная топографическая карта и пункты гравиметрических наблюдений; Б - фрагмент векторизованной карты

Необходимо сразу отметить, что серьезных технических проблем, связанных с использованием «больших ЦММ» при вычислении gр не возникает, возможности современных компьютеров позволяют хранить в оперативной памяти и обрабатывать массивы данных указанной размерности. При этом повышается точность определения поправок за рельеф по сравнению с традиционными технологиями, базирующимися на кодировании топографических карт вручную и становится возможным автоматизация расчета gр в пределах внутренней подобласти D1 [5].

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

Приближенная оценка погрешностей интерполяции значений высот в узлы регулярной сети

Важным вопросом, тесно связанным с повышением точности учета влияния рельефа местности на измеренные значения поля силы тяжести, является вопрос оценки различий в значениях картографических и топографических высотных отметок. С этой целью была проведена серия из нескольких вычислительных экспериментов с использованием материалов по Сусанинской площади, расположенной на Западном Урале. Размер площади около 200 км2; диапазон изменения высотных отметок рельефа местности в ее пределах составляет 122 - 260 м, при среднем значении 169.8 м и среднеквадратическом отклонении (СКО) 26.7 м. Большая часть площади покрыта смешанным лесом.

Первоначально была предпринята попытка оценить погрешности, возникающие за счет интерполяции первичных картографических данных в узлы равномерной (х=у) сети при формировании ЦММ. Под первичными картографическими данными в данном случае подразумевается dat-файл, полученный при векторизации топографических карт масштаба 1:25 000 - 1:50 000, содержащий горизонтальные координаты и значения высот, т.е. множество векторов . Первичные картографические данные для рассматриваемой территории (участок гравиметрической съемки плюс зона обрамления - 1123 км2) содержат около 284000 значений , средняя плотность высотных отметок составляет около 254.7 точек/км2.

Известно, что значения высот рельефа в различных точках земной поверхности функционально не связаны между собой. Поэтому теоретически обосновать применение какого-либо способа интерполяции значений высот в узлы регулярной сети при создании матрицы z, как наиболее оптимального, не представляется возможным. Были выбраны 4 широко применяющихся на практике алгоритма интерполяции, реализованные в программе SURFER 8 (Golden Software Inc.):

· №1 - метод взвешенных расстояний (inverse distance to a power);

· №2 - метод крайгинга (kriging);

· №3 - метод триангуляции (triangulation with linear interpolation);

· №4 - метод локальных полиномов (local polynomial)

и с помощью этих методов были созданы четыре различные ЦММ. Шаг всех ЦММ s= 50 м, размеры матриц z составили 307 строк, 285 столбцов. Ниже, в табл. 1, приведены числовые характеристики различий h между построенными ЦММ в метрике Чебышева и в квадратичной метрике .

Таблица 1. Различия в значениях высотных отметок между ЦММ, сформированных с помощью разных методов интерполяции

Примечание: в числителе дроби приведены значения F1, в знаменателе - F2

Как очевидно, значения высот в узлах ЦММ, полученные различными способами интерполяции, существенно отличаются друг от друга. Наиболее близкие результаты дают методы взвешенных расстояний и триангуляции, наиболее различаются между собой - крайгинг и метод локальных полиномов. Точные значения высот рельефа z в узлах регулярной сети нам неизвестны, поэтому судить о преимуществах какого-то одного из методов интерполяции перед остальными не представляется возможным. Однако проведенное сопоставление позволяет получить косвенную, приближенную оценку среднеквадратической погрешности восстановления неизвестных значений высот в узлах равномерной сети zинт: 1 м zинт 3 м.

Оценка различий между картографическими и топографическими высотами

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

Способ 1. В подавляющем большинстве случаев можно предполагать о достаточно слабых изменениях рельефа дневной поверхности исследуемой территории в пределах таксона - отдельной ячейки ЦММ размером 5050 м. Это допускает применение билинейной интерполяции между узлами ЦММ с целью восстановления поверхности рельефа местности z = z(x,y) в любой заданной точке, в т.ч. - непосредственно в пункте гравиметрических наблюдений. При этом значения высот пунктов, заданные в каталоге, никаким преобразованиям не подвергаются. На данном принципе базируется способ оценки различий высот z, представленных в каталоге гравиметрических пунктов и в каждой из ЦММ (табл. 2). Как очевидно, полученные значения z являются весьма высокими, превосходящими на порядок и более точность определения высот гравиметрических пунктов при выполнении полевых топографо-геодезических работ. Все четыре ЦММ характеризуются весьма близкими различиями топографических и картографических высот, о чем свидетельствуют приведенные статистические параметры и гистограммы распределения z (рис. 2).

Таблица 2. Различия в значениях топографических и картографических высот в пунктах гравиметрических наблюдений (оценка способом 1)

Рис. 2. Гистограммы различий z топографических и картографических высот, полученных способом 1

При интерполяции высот в узлы регулярной сети использовались методы: А - взвешенных расстояний; Б - крайгинга; В - триангуляции; Г - локальных полиномов

Способ 2. При небольших расстояниях R между гравиметрическими пунктами и изогипсами рельефа, нанесенными на топографической карте, можно непосредственно по известным значениям координат (х, у) пунктов определять их картографические высоты (по сути - снимать значения высот z с карты). Алгоритмически процесс «снятия высотных высот с карты» сводится к поиску в окружности заданного радиуса (с центром в гравиметрическом пункте), значений высот из файла первичной картографической информации, т.е. к известной в математическом программировании задаче сортировки данных. При наличии в пределах R точек с известными высотами выбирается высота единственной точки, находящейся на минимальном удалении Rmin от пункта.

Рис. 3. Пространственная корреляция высоких значений модуля z и участков повышенной дисперсии высотных отметок рельефа

Растровые карты: А - абсолютных значений разности картографических и топографических высот; Б - дисперсии высот, рассчитанной в скользящем окне

Примечания: черные точки - пункты гравиметрических наблюдений

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

При проведении вычислительного эксперимента были выбраны значения R = 10, 20, 50 м, т.е. максимальный радиус поиска картографической высоты не превышал 2 мм в 1:25000 масштабе. Полученные результаты приведены в табл. 3. Весьма примечательно то, что с ростом величины Rmin не наблюдается существенного увеличения разности высот z (рис. 4).

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

Рис. 4. Кросс-плот значений z и Rmin

Два разных способа оценки расхождений топографических и картографических высот дают очень близкие между собой результаты: (табл. 2, 3): распределение z в пределах Сусанинской площади близко к закону Гаусса при математическом ожидании M 11.5 м и СКО S 6.2 м.

Можно предположить, что выявленные высокие значения z обусловлены двумя факторами: погрешностями плановой привязки гравиметрических пунктов 1 и ошибками 2, связанными с погрешностями значений высот z, представленных на самих топографических картах. Считая оба вышеуказанных фактора независимыми, для определения численных значений 2 можно использовать аппарат дисперсионного анализа: D(z) = D(1) + D(2), где D - символ, обозначающий дисперсию.

С целью определения вклада составляющей 1 в величину z был выполнено имитационное моделирование, позволяющее оценить ее амплитуду. При этом было принято то же допущение, что и при исследовании z способом 1: предполагалось, что при густой сети высотных отметок (5050 м), билинейная интерполяция позволяет достаточно точно провести восстановление рельефа дневной поверхности между узлами ЦММ. По мнению авторов, в геоморфологических условиях Сусанинской площади, характеризующейся сравнительно пологим рельефом и небольшими градиентами высот, данное допущение является вполне оправданным. Моделирование пространственного распределения погрешностей 1, обусловленных отклонениями планового положения пунктов гравиметрических наблюдений, проводилось при использовании ЦММ, созданной методом крайгинга. Процесс имитационного моделирования заключался в следующем:

Картографические высоты z = z(x, y) имеющихся на площади 2247 пунктов гравиметрических наблюдений, полученные при выполнении билинейной интерполяции матрицы z, были приняты за «истинные».

В значения координат (х, у) каждого гравиметрического пункта путем генерации псевдослучайных чисел вносилась помеха xy, равномерно распределенная на интервале [-20м, 20м]. Таким образом, каждый пункт мог произвольно перемещаться в плане в пределах квадрата со стороной 40 м, в центре которого находилось его истинное местоположение (т.е. максимальное смещение пункта в плане могло достигать м относительно первоначального).

Методом билинейной интерполяции данных по матрице z определялись новые значения высот z* = z(x+xy, y+xy)

Для каждого гравиметрического пункта вычислялась разность «истинной» высоты и высоты, полученной при наличии помехи xy: 1 = z-z*.

Осуществлялось построение гистограммы и определялись статистические характеристики разности высот 1.

Установлено, что выборка, включающая в себя 2247 значений 1, характеризуется следующими параметрами: min = -25.55 м; max = 20.69 м; M = -0.08 м; S = 2.10 м. Полученные результаты убедительно свидетельствуют о том, что наблюдаемые в пределах Сусанинской площади различия топографических и картографических высот z не могут полностью объясняться погрешностями планового положения пунктов гравиметрических измерений. Основной вклад в величину z вносит картографическая погрешностью 2 (рис. 5). Погрешность 2>1 характеризуется математическим ожиданием M (11.5) м и СКО S 4.1 м. Следует подчеркнуть, что проведенная оценка характеризует погрешность 2 по минимуму (вероятно, на самом деле она несколько выше, т.к. плановая привязка гравиметрических пунктов проведена более точно).

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

Рис. 5. Гистограммы: А - разности высот z, определенной способом 2 при R = 10 м; Б - погрешности 1, обусловленной отклонениями в плановом положении гравиметрических пунктов

Топографические значения высот из-за низкой плотности сети их определения и ограниченности информации о рельефе размерами участка гравиметрической съемки не могут быть использованы для создания кондиционных ЦММ. В то же время имеющаяся картографическая информация позволяет получать достаточно достоверные сведения о форме земной поверхности (но не об абсолютных значениях высот). В качестве примера можно привести сопоставление аналитического вычисления gр в центральной зоне радиусом 150 м по картографических высотам (при использовании аппроксимации рельефа полиномом 4-ой степени [5]) и определения поправок, рассчитанных по данным нивелирования 8-лучевых «звездочек» в 50 пунктах гравиметрических наблюдений, для одного из горных районов республики Хакасия. Среднеквадратическое расхождение поправок за рельеф составило всего 0.012 мГал. Для условий Западного Урала, где рельеф местности более пологий и характеризуется небольшими градиентами, эти расхождения еще меньше.

Следовательно, необходимо полностью отказаться от участия топографических высот в процессе вычисления поправок за влияние рельефа, считая при этом достоверными относительные изменения картографических высот в пределах учитываемой области D для каждого пункта гравиметрических наблюдений. При таком подходе следует проводить аналитическое определение картографических высот пунктов гравиметрических наблюдений с использованием их горизонтальных координат; где k - количество пунктов. При этом функция (x,y), описывающая рельеф в пределах ЦММ, может быть получена в явном виде, как предлагается В.Н. Страховым [9, 10 и др.], в результате аналитической аппроксимации множества значений {x, y, z} или матрицы z. Для определения дискретных значений (x,y) в заданных точках области D можно также использовать различные приемы, использующиеся в вычислительной геометрии (построение диаграмм Вороного, триангуляцию Делоне и т.д.).

Стохастический метод оценки точности учета влияния рельефа

Достаточно объективную оценку влияния, в том числе и совокупного, различных возмущающих факторов, можно получить с помощью имитационного моделирования вычисления поправки gр для неоднородных по физическим характеристикам геологических сред, при реальных технических условиях выполнения съемок [4]. Область учитываемого влияния рельефа D при этом представляет собой стохастическую модель, для которой "решение прямых задач геофизики имеет вероятностную трактовку и заключается в расчете математических ожиданий, дисперсий или корреляционных функций аномальных эффектов" [1].

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

(1)

где U = {u1, u2, …, um) - m-мерный вектор параметров, характеризующий физические и геометрические параметры отдельной аппроксимационной ячейки и ее местоположение относительно точки расчета gр с координатами (x,y,z); - оператор решения прямой задачи гравиразведки; kk - число элементарных ячеек в пределах области D.

Рассмотрим новый вектор P = {p1, p2,…pm}, структура которого аналогична структуре вектора U, но одна или несколько компонент которого осложнены случайной составляющей : pn = un + , где n = 1, 2, …, n. Погрешность определения поправки g в отдельной точке пространства, обусловленная наличием случайных отклонений в исходных данных, определяется выражением:

(2)

Оценка погрешностей, выраженных приведенной выше формулой, осуществляется с помощью метода Монте-Карло. Моделирование случайной составляющей выполняется путем генерации последовательности псевдослучайных чисел 1, 2, 3, …. . Статистические характеристики данной последовательности, как правило, можно определить, исходя из априорных сведений о геологической среде и технических особенностях съемок. Для оценки одновременного воздействия q возмущающих факторов для единственной точки пространства требуется qkk случайных чисел.

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

Рассмотрим пример оценки точности поправок за рельеф стохастическим методом. Определение g осуществлялось для Верх-Истокской площади (Западный Урал), на которой была проведена гравиметрическая съемка 1:25 000 масштаба. С помощью генерации серий случайных чисел моделировались произвольно ориентированные в пространстве отклонения гравиметрических пунктов в плане от их истинного местоположения. Предполагалось, что отклонения по каждой координате Х и У взаимно независимы, амплитуды смещений распределены по нормальному закону с математическим ожиданием M = 0 и СКО S = 14.14 м. Таким образом, максимальная величина суммарного смещения xy гравиметрического пункта в плане при вероятности 68 % (на уровне M S) не превышала , предусмотренных "Инструкцией по гравиразведке" для съемки масштаба 1:25 000 (рис. 6).

Алгоритмы вычисления поправок, как для внутренней подобласти D1, так и внешней подобласти D2 использовали только картографические высоты пунктов гравиметрических наблюдений. В подобласти D1 высота пункта наблюдений определялась его местоположением на поверхности полинома 4-ой степени, аппроксимирующего рельеф дневной поверхности. Для подобласти D2 высота пункта вычислялась путем билинейной интерполяции высот ближайших четырех узлов ЦММ, а затем выполнялась 3D-интерполяция значений gр из этих узлов в сам пункт.

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

Рис. 6. Полярная диаграмма: модули случайных смещений гравиметрических пунктов по координате Х

Погрешности расчета поправок определялись по разности поправок, вычисленных при первоначальном (указанном в каталоге) и смещенном положении гравиметрических пунктов. Для внутренней подобласти D1 (квадрата со стороной 300 м) эти погрешности изменяются в диапазоне [-0.055, 0.056] мГал, при среднем значении 0.0001 мГал и СКО 0.0075 мГал. Для внешней подобласти D2 (внешний контур которой имеет форму квадрата со стороной 5100 м, внутренний - квадрата со стороной 300 м) эти погрешности изменяются в диапазоне [-0.061, 0.068] мГал, при среднем значении 0.00008 мГал и СКО 0.009 мГал. Статистическое распределение погрешностей подчиняется закону Гаусса. Таким образом, предельная среднеквадратическая погрешность определения суммарных поправок за рельеф местности в радиусе 0 2550 м в условиях исследуемой площади не будет превышать 0.017 мГал (рис. 7). Оценки погрешностей проведены для плотности промежуточного слоя 2.67 г/см3, для плотности 2.30 г/см3 они будут примерно на 15% ниже.

Рис. 7. Гистограммы погрешностей g определения поправок при гравиметрической съемке масштаба 1:25000. Западный Урал

А - внутренняя подобласть D1; Б - внешняя подобласть D2

ЗАКЛЮЧЕНИЕ

Широко применяющиеся на практике способы учета влияния рельефа местности созданы в период "парадигмы ранней компьютерной эпохи (условно 1960 - 1985 гг.)" [8]. Ограниченные вычислительные возможности и высокая стоимость машинного времени использующихся в этот период ЭВМ, а также сложности технологического характера, связанные с формированием ЦММ на машинных носителях наложили свой отпечаток на технологии определения gр. Большинство из этих технологий, перенесенных в начале 90-х годов прошлого века на персональные компьютеры, используют совместно как картографические (ЦММ), так и топографические (каталог гравиметрических пунктов) высоты, что приводит к появлению дополнительных погрешностей при вычислении поправок.

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

Основной причиной возникновения погрешности 2 в условиях залесенной местности является, вероятно, методика построения горизонталей рельефа по данным аэрофотосъемки. Совершенно очевидно, что следует полностью отказаться от использования топографических высот пунктов гравиметрических наблюдений при вычислении поправок за влияние рельефа местности, выполняя замену этих высот на картографические (т.е. при фиксированном плановом положении пунктов осуществлять их «перенос» на поверхность ЦММ). Наиболее рациональным путем исключения влияния 2 является построение метрологических линейных аппроксимаций рельефа местности, предложенное В.Н. Страховым, при использовании в качестве исходных данных первичной картографической информации. Аналитические зависимости вида позволят определять как высоты самих гравиметрических пунктов, так и восстанавливать высотные отметки рельефа в заданных точках пространства, что существенно повышает точность и упрощает технологию расчета gр.

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

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

ЛИТЕРАТУРА

Вахромеев Г.С., Давыденко А.Ю. Моделирование в разведочной геофизике. М.: Недра, 1987. 192 с.

Гордин В.М. Способы учета влияния рельефа дневной поверхности при высокоточных гравитационных измерениях. Обзор ОНТИ ВИЭМС. М.: 1974, сер. IX. 89 с.

Гравиразведка: Справочник геофизика / Под ред. Е.А. Мудрецовой, К.Е. Веселова. 2-е изд. перераб. и доп. М.: Недра, 1990. 607 с.

Долгаль А.С. Моделирование погрешностей учета влияния рельефа при гравиметрической съемке.//Известия РАН. Сер. Физика Земли. № 8. 1997. С. 88-93.

Долгаль А.С. Усовершенствование технологии учета влияния рельефа местности при гравиметрической съемке. // Геофизический журнал. № 2. 1998. С. 51 - 57.

Каленицкий А.И., Смирнов В.П. Методические рекомендации по учету влияния рельефа местности в гравиразведке. Новосибирск, СНИИГГиМС, 1981. 160 с.

Маловичко А.К., Костицин В.И., Тарунина О.Л. Детальная гравиразведка на нефть и газ. 2-е изд., перераб. и доп. М.: Недра, 1989. 224 с.

Страхов В.Н. Что делать ? (о развитии гравиметрии и магнитометрии в России в начале XXI века) М.: ОИФЗ РАН, 1998. 24 с.

Страхов В.Н., Керимов И.А., Страхов А.В. Линейные аналитические аппроксимации рельефа поверхности Земли.//Геофизика и математика: материалы 1-й Всероссийской конференции. М.: ОИФЗ РАН, 1999. С. 199-212.

Размещено на Allbest.ru