Способ q томографии

Иллюстрации

Показать все

Изобретение относится к области геофизики и может быть использовано при обработке данных сейсмических исследований. Заявлен способ перестроения моделей (110) Q геологической среды на основании сейсмических данных (10) путем осуществления лучевой Q томографии сдвига центроидных частот. Амплитудный спектр волнового сигнала сейсмического источника аппроксимируют (40) частотно-взвешенной экспоненциальной функцией частоты, имеющей два подбираемых параметра для приведения в соответствие данным о сдвиге частот. В результате чего обеспечивают лучшее соответствие различным асимметричным амплитудным спектрам источника. Боксовые ограничения могут использоваться при выполнении процедуры оптимизации, а многоиндексный способ активных множеств, используемый при томографии скорости, является предпочтительным способом для реализации (100) боксовых ограничений. Технический результат - повышение точности данных сейсмических исследований. 2 н. и 18 з.п. ф-лы, 12 ил.

Реферат

Перекрестная ссылка на родственную заявку

По этой заявке испрашивается преимущество приоритета предварительной заявки № 61/331694 на патент США, поданной 5 мая 2010 года, под названием “Q tomography method”, которая полностью включена в эту заявку путем ссылки.

Область техники, к которой относится изобретение

В общем настоящее изобретение относится к области геофизических исследований, а более конкретно к обработке сейсмических данных. В частности, изобретение относится к области Q томографии.

Предпосылки создания изобретения

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

Затухание сейсмических волн можно численно описать добротностью Q. Простое предположение заключается в том, что затухание сейсмических волн является частотно-зависимым, но добротность является частотно-независимой. Это предположение является справедливым в диапазоне частот применений разведочной геофизики. Q томография является способом оценивания Q. В этом способе стараются перестраивать двумерные или трехмерные модели Q геологической среды на основании сейсмических данных. Обычно алгоритмы Q томографии подразделяют на две основные категории. К одной категории относится лучевая томография (Quan and Harris, 1997; Plessix, 2006; Rossi et al., 2007). К другой категории относится томография на основе решения волнового уравнения (Liao and McMechan, 1996; Hick and Pratt, 2001, Pratt et al., 2003; Watanabe et al., 2004; Gao et al., 2005). Томография на основе решения волнового уравнения физически является более точной, но требующей больших затрат на вычисление и непрактичной в трехмерных случаях. Настоящее изобретение относится к категории лучевой Q томографии.

Одна основная проблема, связанная с Q томографией, заключается в способе установления связи между моделями Q и сейсмическими данными при минимальных приближениях и с максимальной гибкостью. Широко используемый способ основан на зависимости между Q и ослаблением амплитуд сейсмических волн. В другом способе используют сдвиг вниз центроидных частот сейсмических волн для оценивания добротности Q. Этот последний способ считают более робастным, поскольку этот способ не зависит от эффекта геометрического расхождения и потерь на отражение/прохождение. Однако в обычном способе сдвига центроидных частот можно использовать только Гауссову, прямоугольную или треугольную функцию для аппроксимации амплитудного спектра источника, что вносит значительную погрешность, поскольку в большинстве случаев спектр источника нельзя аппроксимировать этими функциями. Настоящее изобретение включает в себя частотно-взвешенную экспоненциальную функцию, которая предназначена для аппроксимации различных асимметричных амплитудных спектров источника, для повышения точности Q томографии при значительно сниженной погрешности аппроксимации амплитудного спектра источника.

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

Томография затухания сейсмических волн (Q томография) исследовалась в течение многих лет и был достигнут значительный прогресс. Двумя основными компонентами алгоритма лучевой Q томографии являются 1) простая, но точная зависимость между сейсмическими данными и моделями Q, предназначенная для построения математической модели Q томографии; 2) надежный и робастный алгоритм оптимизации, предназначенный для решения этой математической задачи. Многочисленные способы были разработаны или предложены для создания этих двух компонентов. Эти способы рассматриваются ниже.

Установление связи между сейсмическими данными и моделями Q

Наиболее простой и прямой способ оценивания Q представляет собой способ спектрального отношения (Spencer et al., 1982; Tonn, 1991), в котором логарифм отношения спектров двух сейсмических волновых сигналов вычисляют в виде функции частоты и эту функцию аппроксимируют линейной функцией частоты, наклон которой трактуют как накопленное затухание сейсмической волны и в конечном итоге связывают со значениями Q вдоль пути пробега волны. В идеальном случае этим способом исключают эффект геометрического расхождения и потери на отражение/прохождение в предположении, что эти эффекты являются частотно-независимыми. При практических применениях этот способ является относительно ненадежным вследствие наложения сейсмических импульсов, неопределенности при линейной аппроксимации и наличия многих других факторов.

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

Основываясь на том, что изменение времени нарастания сейсмического импульса линейно связано с профилем 1/Q вдоль пути пробега волны, Wu и Lees (1996) описали способ томографии на основе затухания сейсмических волн с использованием времени нарастания, предназначенный для сейсмологических исследований. К сожалению, этот способ является непрактичным для разведочной геофизики вследствие неизбежного загрязнения сейсмических импульсов шумом, эффекта рассеяния, наложения и т.д.

Было замечено, что на форму амплитудного спектра сейсмического импульса едва ли не исключительно влияет добротность Q, и был разработан способ изменения пиковой частоты для оценивания Q (Zhang, 2008). Этот способ является привлекательным, но на практике имеются трудности, заключающиеся в неточности измерения изменения пиковой частоты. Кроме того, поскольку используется только информация на индивидуальной частоте, неопределенность оценивания Q может быть большой.

Более робастный способ выдвинули Quan и Harris (1997), в котором информацию в пределах всего диапазона частот сейсмических волновых сигналов используют для вычисления сдвига вниз центроидных частот и затем связывания сдвига центроидных частот с профилем Q вдоль траектории луча с помощью простой формулировки в замкнутой форме. Этот способ по своей сущности невосприимчив к геометрическому расхождению и потерям на отражение/прохождение. Ограничение этого способа заключается в том, что амплитудный спектр источника должен представляться Гауссовой, прямоугольной или треугольной функцией. Хорошо известно, что амплитудный спектр сейсмических волн никогда не представляется прямоугольной или треугольной функцией. Кроме того, обычно он является асимметричным и может очень отличаться от Гауссовой функции. Если этот асимметричный амплитудный спектр аппроксимировать Гауссовой функцией, значительные погрешности будут вноситься при перестраивании моделей Q. Поэтому, если имеется функция, которую можно использовать для более точной аппроксимации различных частотных спектров сейсмических волн без потери простого характера зависимости между центроидной частотой регистрируемых сейсмических данных и профилем Q вдоль траектории луча, то этот робастный способ может быть более точным при практических применениях Q томографии.

Алгоритмы оптимизации с ограничениями для Q томографии

Когда зависимость между сейсмическими данными и моделями Q установлена, задачу лучевой Q томографии можно описывать с помощью линейной задачи оптимизации. В большей части существующих алгоритмов Q томографии эта линейная задача оптимизации решается итеративно с использованием способов подпространства Крылова, таких как способ сопряженных градиентов и способ LSQR, без применения каких-либо ограничений (Quan and Harris, 1997; Plessix, 2006, Rossi et al., 2007). Эти алгоритмы хорошо работают при условии, что сейсмические данные имеют высокое отношение сигнала к шуму. Однако сейсмические данные никогда не являются чистыми; при обработке реальных полевых данных эти алгоритмы оптимизации без ограничений всегда приводят к некоторому количеству отрицательных значений Q или экстремально низких значений Q, которые являются физически нереальными. Кроме того, при некоторых обстоятельствах является известной априорная информация относительно диапазона значений Q. В этих случаях эту информацию необходимо включать в алгоритм томографии через посредство боксовых ограничений, чтобы получать более надежные результаты Q томографии.

Rickett (2006) разработал алгоритм оценивания Q с ограничением. Но в его алгоритме принято неотрицательное ограничение вместо боксового ограничения, что означает исключение отрицательных значений Q, но экстремально низкие положительные значения Q все же могут существовать. Rickett (2006) описал два способа применения неотрицательного ограничения. Первый способ представляет собой нелинейный способ преобразования. В этом нелинейном способе преобразования переменную Q заменяют на e y и находят решение относительно y вместо 1/Q. Выполняя это, принуждают значение Q быть положительным, но вся система может быть сильно нелинейной. Для достижения этой цели результирующую систему оптимизации решают способом Гаусса-Ньютона, что может быть очень затратным. Другой недостаток этого способа заключается в том, что во время оптимизации, когда значения Q являются очень большими, алгоритм оптимизации на градиентной основе является инертным или сходится очень медленно, то есть значения Q остаются неизменными и больше не изменяются. В тяжелом случае, если значения y в исходной модели являются бесконечными, то градиент функции стоимости равен 0 и алгоритм оптимизации не выполняет поиска. Второй способ применения неотрицательного ограничения, который предложил Rickett (2006), заключается в повышении монотонности возрастания характеристики затухания сглаживанием. Этот способ эффективно работает при оценивании Q с использованием данных вертикального сейсмического профилирования, но может перестать действовать при двумерной Q томографии с использованием сейсмических данных об отражениях от поверхности.

Изобретателям настоящего изобретения известно об отсутствии алгоритма Q томографии с боксовыми ограничениями для улучшения оцениваемых значений Q в пределах диапазонов, точно определяемых верхними границами и нижними границами. Однако алгоритмы оптимизации с боксовыми ограничениями используют при некоторых других геофизических применениях, таких как томография скорости. Например, Delbos и соавторы (2006) разработали алгоритм томографии сейсмического отражения с боксовыми ограничениями. В их алгоритме задача оптимизации с ограничениями решается способом Гаусса-Ньютона, дополненным Лагранжем, а связанная задача Лагранжа, еще одна задача оптимизации с ограничениями, решается сочетанием способа проекции градиента, способа активных множеств и способа сопряженных градиентов. Способ активных множеств используется обычным образом и это неэффективно, поскольку в алгоритме обновляется активное множество в соответствии с одним ограничением за один раз (Bierlaire et al., 1991; Löstedt, 1984; Nocedal and Wright, 1999). Когда количество боксовых ограничений очень большое, скорость сходимости алгоритма может быть очень низкой.

В настоящем изобретении последние достижения в области математики, к которым можно отнести многоиндексный способ активных множеств (Morigi et al., 2007), используются для осуществления Q томографии с боксовыми ограничениями, что значительно улучшает характеристики алгоритма Q томографии в части качества перестроения Q и эффективности алгоритма по сравнению с характеристиками алгоритма Q томографии без ограничений и алгоритмом с ограничениями при использовании обычного способа активных множеств.

Краткое изложение изобретения

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

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

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

Краткое описание чертежей

Настоящее изобретение и его преимущества станут более понятными при обращении к нижеследующему подробному описанию и сопровождающим чертежам, на которых:

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

фиг.2 - вид кривых частотно-взвешенной экспоненциальной функции F ( f )   =   A f n   exp ( −   f f 0 ) при n=1, 2, 3, 4 и 5, при этом центроидная частота (n+1)f 0 является фиксированной, 30 Гц;

фиг.3 - блок-схема последовательности действий, показывающая основные этапы осуществления настоящего изобретения при Q томографии;

фиг.4 - 8А-В относятся к применению синтетических примеров, при этом на

фиг.4 - вид амплитудного спектра источника;

фиг.5А-В - иллюстрация скоростной модели, траекторий лучей (5А) и истинной модели 1/Q (5B);

фиг.6А - вид модели 1/Q, перестроенной с использованием обычного способа сдвига центроидных частот Q томографии при аппроксимации Гауссовой функцией и без боксовых ограничений; фиг.6В - иллюстрация различий между перестроенной моделью 1/Q из фиг.6А и истинной моделью 1/Q из фиг.5В;

фиг.7А - вид модели 1/Q, перестроенной с использованием способа сдвига центроидных частот Q томографии при аппроксимации частотно-взвешенной экспоненциальной функцией согласно настоящему изобретению, но без боксовых ограничений; фиг.7В - иллюстрация различий между перестроенной моделью 1/Q из фиг.7А и истинной моделью 1/Q из фиг.5В;

фиг.8А - вид модели 1/Q, перестроенной с использованием способа сдвига центроидных частот Q томографии при аппроксимации частотно-взвешенной экспоненциальной функцией согласно настоящему изобретению и боксовыми ограничениями согласно настоящему изобретению; фиг.8В - иллюстрация различий между перестроенной моделью 1/Q из фиг.8А и истинной моделью 1/Q из фиг.5В.

Вследствие ограничений, накладываемых патентным законом на чертежи, фиг.5А-В - 8А-В являются черно-белыми репродукциями цветных изображений. Копии цветных чертежей могут быть получены при запросе копии экземпляра заявки США от Ведомства США по патентам и товарным знакам и уплате сбора.

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

Подробное описание примеров осуществлений

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

Основными признаками настоящего изобретения согласно по меньшей мере некоторым осуществлениям являются следующие. Амплитудный спектр сейсмического импульса источника анализируют и аппроксимируют специально синтезированной функцией, частотно-взвешенной экспоненциальной функцией. Аппроксимацию амплитудного спектра сейсмического импульса источника осуществляют путем подбора двух параметров частотно-взвешенной экспоненциальной функции. Сдвиги центроидных частот принимаемых сейсмических волновых сигналов относительно центроидной частоты сейсмического импульса источника вычисляют и вводят в алгоритм оптимизации с боксовыми ограничениями для построения модели Q, при этом диапазоны значений Q заранее определяют в соответствии с априорной информацией. В отличие от первоначального способа сдвига центроидных частот (Quan and Harris, 1997) со специально синтезированной функцией аппроксимации амплитудного спектра этим алгоритмом Q томографии можно обрабатывать негауссовские и асимметричные амплитудные спектры источника, чтобы снижать погрешность, вносимую рассогласованием при аппроксимации амплитудного спектра источника. При наличии такой специально синтезированной функции аппроксимации спектра источника задача Q томографии сводится к задаче условной оптимизации с боксовыми ограничениями. Эту задачу оптимизации с ограничениями решают, используя многоиндексный способ активных множеств (Morigi et al., 2007), с помощью которого дополнительно повышают точность и робастность этого алгоритма Q томографии без ухудшения высокоэффективных характеристик. Термин «активное множество» относится к такому подмножеству из множества оптимизируемых неизвестных, которое, как будет показано, не может быть обновлено в конце итерационного цикла, поскольку для него присуще ограничение верхнего предела или нижнего предела.

Далее поясняются несколько основополагающих теорий изобретения.

Если амплитудный спектр сейсмического импульса источника представить как S(f), то амплитудный спектр R(f) принимаемого сейсмического волнового сигнала можно выразить как (Quan and Harris, 1997)

R(f)=GH(f)S(f), (1)

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

H ( f )   =   exp ( −   f ∫ п о   л у ч у π Q v d l ) , (2)

где Q является добротностью и v является скоростью сейсмических (акустических) волн.

В способе сдвига центроидных частот при оценивании Q ключевая часть заключается в использовании аналитической функции для аппроксимации амплитудного спектра S(f) источника и затем в получении явно выраженной зависимости между центроидной частотой принимаемого сейсмического волнового сигнала и профилем Q вдоль пути пробега волны. В первоначальном способе сдвига центроидных частот эту явно выраженную зависимость можно было получать только в предположении, что амплитудный спектр источника является Гауссовой функцией, или прямоугольной функцией (функцией, которая равна нулю на протяжении всей вещественной линии за исключением единственного интервала, где она равна постоянной величине), или треугольной функцией. Это является основным недостатком традиционного способа сдвига центроидных частот, который может приводить к большим погрешностям при оценивании Q (Rickett, 2006; Zhang, 2008).

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

F ( f )   =   A f n   exp ( −   f f 0 ) , (3)

где А является постоянной для масштабирования амплитуды, f 0 называют характеристической частотой и n называют индексом симметрии. Характеристическая частота f 0 является параметром при управлении шириной полосы частот. Если индекс n симметрии является фиксированным, ширина полосы частот этой функции увеличивается с повышением характеристической частоты. Фактически центроидная частота функции F(f) имеет очень простую форму, приведенную ниже:

f F   =   ∫ 0 + ∞ f F ( f ) d f ∫ 0 + ∞ F ( f ) d f   =   ( n   +   1 ) f 0 . (4)

На фиг.1 показано изменение ширины полосы частот функции F(f) при n=2 для пяти различных характеристических частот в пределах от 10 Гц до 50 Гц с размером шага 10 Гц. На фиг.1 при изменении f 0 от 50 Гц до 10 Гц форма F(f) остается относительно неизменной, тогда как ширина полосы частот сокращается, как и ожидалось. С другой стороны, индекс n симметрии используют для управления свойством симметрии F(f). Как показано на фиг.2, чем больше индекс n симметрии, тем более симметрична форма функции F(f), при этом центроидная частота фиксирована на уровне 30 Гц, тогда как индекс n симметрии изменяется от 1 (наименьшая симметрия) до 5 (наибольшая симметрия). Нет необходимости в том, чтобы n было целым числом. Теоретически для точности аппроксимации n может быть любым вещественным числом. Однако на практике обычно n не должно превышать 5. Как упоминалось ранее, поскольку управление центроидной функцией и свойством симметрии этой специально синтезированной функции осуществляют раздельно с помощью двух различных параметров, пользователям очень легко точно аппроксимировать различные асимметричные амплитудные спектры источника этой функцией.

В соответствии с лучевой теорией, если амплитудный спектр источника аппроксимируют частотно-взвешенной экспоненциальной функцией, то есть S(f)=F(f), то после подстановки уравнений (2) и (3) в уравнение (1) амплитудный спектр принимаемого сейсмического волнового сигнала можно записать в форме

R ( f )   =   G H ( f ) S ( f )   =   A G f n exp [ −   f ( ∫ п о   л у ч у π Q v d l   +   1 f 0 ) ] . (5)

Центроидная частота принимаемого сейсмического волнового сигнала может быть вычислена как

f R   =   ∫ 0 + ∞ f R ( f ) d f ∫ 0 + ∞ R ( f ) d f   =   n   +   1 ∫ п о   л у ч у π Q v d l   +   1 f 0 . (6)

Поскольку центроидная частота амплитудного спектра f S источника равна (n+1)f 0, сдвиг центроидной частоты между амплитудным спектром источника и амплитудным спектром принимаемого сигнала можно легко получить в виде

Δ f   =   f S   −   f R   =   ( n   +   1 ) ( f 0 1 ∫ п о   л у ч у π Q v d l   +   1 f 0 ) . (7)

Теперь путем решения уравнения (7) суммарное затухание ∫ п о   л у ч у π Q v d l вдоль траектории луча можно получить на основании сдвига центроидной частоты:

∫ п о   л у ч у π Q v d l   =   Δ f f 0 [ ( n   +   1 ) f 0   −   Δ f ]   =   Δ f f 0 ( f S   −   Δ f )   =   Δ f f 0 f R   =   ( n   +   1 ) Δ f f S f R . (8)

Уравнение (8) показывает, что 1/Q связано со сдвигом центроидной частоты через очень простую линейную зависимость с частотно-взвешенной экспоненциальной функцией.

Дискретная форма уравнения (8) для i-го измерения (i-го луча) имеет вид

∑ j π l j i Q j v j   =   Δ f i f 0 f R i , (9)

где верхний индекс i является индексом измерения и нижний индекс j является индексом сетки, а l j i является длиной луча на j-й сетке при i-ом измерении.

После сбора всех измерений уравнение (9) можно записать в матричной форме

Ax=d. (10)

В уравнении (10) А является матрицей ядра, элементы которой определяются в соответствии с

A i j   =   π l j i v j , (11)

x является вектором неизвестных величин, то есть

xj=1/Qj, (12)

и d является вектором измерения, определяемым в соответствии с

d i   =   Δ f i f 0 f R i . (13)

Оставшейся задачей является решение уравнения (10) относительно Q. Поскольку измеряемые данные неизбежно загрязнены шумом, линейная система уравнений (10) является плохо обусловленной и имеет неоднозначные решения. Поэтому эту систему можно трактовать как задачу наименьших квадратов и находить решение в виде приближенного решения задачи квадратичного программирования

min ‖ A x   −   d ‖ . (14)

где ‖ ... ‖ означает Евклидов вектор.

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

m i n ‖ A x   −   d ‖ при условии l<x<u, (15)

где l и u представляют собой векторы, сохраняющие нижние грани