Введение
Процесс вычисления причинных факторов, породивших набор наблюдений, и его применение в геодезии. Инверсная задача в науке – это процесс определения причинных факторов по набору наблюдений, то есть вычисление того, что вызвало наблюдаемые эффекты. Примерами могут служить реконструкция изображения в компьютерной рентгеновской томографии, восстановление источника звука в акустике или определение плотности Земли по измерениям её гравитационного поля. Она называется инверсной, поскольку исходит из следствий для вычисления причин. Это противоположность прямой задаче, которая начинается с причин и вычисляет следствия. Инверсные задачи – одни из важнейших математических задач в науке и математике, поскольку они позволяют узнать о параметрах, которые невозможно измерить напрямую. Они широко используются в идентификации систем, оптике, радиолокации, акустике, теории связи, обработке сигналов, медицинской визуализации, компьютерном зрении, геофизике, океанографии, астрономии, дистанционном зондировании, обработке естественного языка, машинном обучении, неразрушающем контроле, анализе устойчивости склонов и многих других областях.
the use in geodesy
An inverse problem in science is the process of calculating from a set of observations the causal factors that produced them: for example, calculating an image in X ray computed tomography, source reconstruction in acoustics, or calculating the density of the Earth from measurements of its gravity field. It is called an inverse problem because it starts with the effects and then calculates the causes. It is the inverse of a forward problem, which starts with the causes and then calculates the effects. Inverse problems are some of the most important mathematical problems in science and mathematics because they tell us about parameters that we cannot directly observe. They have wide application in system identification, optics, radar, acoustics, communication theory, signal processing, medical imaging, computer vision, geophysics, oceanography, astronomy, remote sensing, natural language processing, machine learning, nondestructive testing, slope stability analysis and many other fields.
История
Начиная с изучения последствий для выявления причин, физиков занимало веками. Историческим примером служат вычисления Адамса и Ле Верье, которые привели к открытию Нептуна на основе возмущенной траектории Урана. Однако формальное исследование обратных задач началось лишь в 20-м веке. Одним из первых примеров решения обратной задачи, обнаруженным Германом Вейлем и опубликованным в 1911 году, было описание асимптотического поведения собственных значений оператора Лапласа — Бельтрами. Сегодня это известно как закон Вейля, и его, пожалуй, легче всего понять как ответ на вопрос, возможно ли определить форму барабана по его звучанию. Вейль предположил, что собственные частоты барабана связаны с его площадью и периметром определенным уравнением, результат, который был усовершенствован более поздними математиками. Позже область обратных задач затронул советский физик армянского происхождения Виктор Амбарцумян. Еще будучи студентом, Амбарцумян глубоко изучал теорию атомной структуры, формирование энергетических уровней, уравнение Шредингера и его свойства, а когда овладел теорией собственных значений дифференциальных уравнений, он указал на явную аналогию между дискретными энергетическими уровнями и собственными значениями этих уравнений. Затем он задался вопросом: зная семейство собственных значений, можно ли определить форму уравнений, для которых эти значения являются собственными? По сути, Амбарцумян исследовал обратную задачу Штурма — Лиувилля, связанную с определением уравнений для вибрирующей струны. Эта статья была опубликована в 1929 году в немецком физическом журнале Zeitschrift für Physik и долгое время оставалась незамеченной. Описывая эту ситуацию спустя десятилетия, Амбарцумян сказал: «Если астроном публикует статью с математическим содержанием в физическом журнале, то, скорее всего, она будет забыта». Тем не менее, к концу Второй мировой войны эту статью, написанную 20-летним Амбарцумяном, обнаружили шведские математики, и она стала отправной точкой для целого направления исследований в области обратных задач, заложив основу для целой дисциплины. Впоследствии значительные усилия были направлены на «прямое решение» обратной задачи рассеяния, особенно Гельфандом и Левитаном в Советском Союзе. Они предложили аналитический конструктивный метод определения решения. С появлением компьютеров некоторые авторы исследовали возможность применения их подхода к аналогичным задачам, таким как обратная задача для 1D волнового уравнения. Однако вскоре стало ясно, что обращение является нестабильным процессом: шум и ошибки могут быть значительно усилены, что делает прямое решение практически невозможным. Затем, примерно в 1970-х годах, появились методы наименьших квадратов и вероятностные подходы, которые оказались очень полезными для определения параметров, участвующих в различных физических системах. Этот подход оказался весьма успешным. В настоящее время обратные задачи исследуются также в областях, выходящих за рамки физики, таких как химия, экономика и информатика. В конечном итоге, по мере распространения численных моделей во многих сферах жизни, можно ожидать, что с каждой из этих моделей будет связана обратная задача.
Элементарный пример: гравитационное поле Земли
Лишь немногие физические системы действительно линейны относительно параметров модели. Одной из таких систем в геофизике является система гравитационного поля Земли. Гравитационное поле Земли определяется распределением плотности Земли в недрах. Поскольку литология Земли меняется весьма значительно, мы можем наблюдать незначительные различия в гравитационном поле Земли на поверхности. Из нашего понимания гравитации (закона всемирного тяготения Ньютона) мы знаем, что математическое выражение для гравитации имеет вид:
где – мера локального гравитационного ускорения, – универсальная гравитационная постоянная, – локальная масса (связанная с плотностью) породы в недрах и – расстояние от массы до точки наблюдения. Дискретизируя это выражение, мы можем соотнести дискретные данные наблюдений на поверхности Земли с дискретными параметрами модели (плотностью) в недрах, которые мы хотим изучить. Например, рассмотрим случай, когда измерения проводятся в 5 точках на поверхности Земли. В этом случае наш вектор данных является столбцом-вектором размерности (5×1): его компонент соответствует точке наблюдения. Мы также знаем, что у нас есть только пять неизвестных масс в недрах (нереалистично, но используется для демонстрации концепции) с известным местоположением: обозначим расстояние между точкой наблюдения и массой. Таким образом, мы можем построить линейную систему, связывающую пять неизвестных масс с пятью точками данных, следующим образом:
Для определения параметров модели, соответствующих нашим данным, мы можем попытаться инвертировать матрицу , чтобы напрямую преобразовать измерения в параметры модели. Например:
Система с пятью уравнениями и пятью неизвестными – это очень специфический случай: наш пример был разработан таким образом, чтобы привести к этой специфике. В общем случае количество данных и неизвестных различно, поэтому матрица не является квадратной. Однако даже квадратная матрица может не иметь обратной: матрица может быть вырожденной (то есть иметь нулевые собственные значения), и решение системы не будет единственным. Тогда решение обратной задачи будет не определено. Это первая трудность. Переопределенные системы (больше уравнений, чем неизвестных) имеют другие проблемы. Кроме того, шум может искажать наши наблюдения, приводя к тому, что может оказаться за пределами пространства возможных ответов на параметры модели, и решение системы может не существовать. Это еще одна трудность.
Инструменты для преодоления первой трудности
Первая трудность отражает важнейшую проблему: наши наблюдения не содержат достаточной информации и требуются дополнительные данные. Дополнительные данные могут поступать из априорной физической информации о значениях параметров, об их пространственном распределении или, в более общем смысле, об их взаимной зависимости. Они также могут быть получены из других экспериментов: например, можно рассмотреть объединение данных, зарегистрированных гравиметрами и сейсмографами, для более точной оценки плотности. Интеграция этой дополнительной информации по сути является задачей статистики. Именно эта дисциплина может ответить на вопрос: как комбинировать величины различной природы? Мы будем более конкретны в разделе "Байесовский подход" ниже. Что касается распределенных параметров, то априорная информация об их пространственном распределении часто включает информацию об их производных. Также, хотя это и несколько искусственно, распространена практика поиска "наиболее простой" модели, которая разумно согласуется с данными. Обычно этого достигают путем штрафования нормы градиента (или полной вариации) параметров (этот подход также называют максимизацией энтропии). Можно также упростить модель, используя параметризацию, которая вводит степени свободы только при необходимости. Дополнительная информация может быть интегрирована также посредством ограничений-неравенств на параметры модели или некоторые их функции. Такие ограничения важны для избежания нереалистичных значений параметров (например, отрицательных). В этом случае пространство, определяемое параметрами модели, перестает быть векторным пространством и становится подмножеством допустимых моделей, которое в дальнейшем будем обозначать как 𝒟.
Инструменты для преодоления второй трудности
Как упоминалось выше, шум может быть настолько сильным, что наши измерения не будут отражать какую-либо модель. В этом случае мы не можем искать модель, которая генерирует данные, а должны искать наилучшую (или оптимальную) модель – то есть, модель, которая наилучшим образом соответствует данным. Это приводит нас к минимизации целевой функции, а именно функционала, который количественно оценивает величину остатков или расстояние между предсказанными и наблюдаемыми данными. Разумеется, при наличии идеальных данных (то есть, без шума) восстановленная модель должна идеально соответствовать наблюдаемым данным. Стандартная целевая функция, , имеет вид:
где – евклидова норма (которая станет -нормой, когда измерениями будут функции, а не отдельные значения) остатков. Этот подход эквивалентен использованию метода наименьших квадратов, широко применяемого в статистике. Однако, евклидова норма известна своей высокой чувствительностью к выбросам. Чтобы избежать этой проблемы, можно рассмотреть использование других метрик расстояния, например, -нормы, вместо евклидовой нормы.
Байесовский подход
Очень похож на метод наименьших квадратов вероятностный подход: если нам известны статистические характеристики шума, искажающего данные, мы можем стремиться к нахождению наиболее вероятной модели m, то есть модели, удовлетворяющей критерию максимального правдоподобия. Если шум имеет гауссовское распределение, критерий максимального правдоподобия сводится к критерию наименьших квадратов, при этом евклидово скалярное произведение в пространстве данных заменяется скалярным произведением с учетом ковариации шума. Более того, если имеется априорная информация о параметрах модели, можно использовать байесовский вывод для формулирования решения обратной задачи. Этот подход подробно описан в книге Тарантолы.
Числовое решение нашего элементарного примера
Здесь мы используем евклидову норму для количественной оценки расхождений между данными и моделью. Поскольку мы имеем дело с линейной обратной задачей, целевая функция является квадратичной. Для её минимизации традиционно вычисляют её градиент, используя ту же логику, что и при минимизации функции одной переменной. В оптимальной модели этот градиент обращается в ноль, что можно записать следующим образом:
где FT обозначает транспонированную матрицу F. Это уравнение упрощается до:
Это выражение известно как нормальное уравнение и предоставляет возможное решение обратной задачи. В нашем случае матрица обычно имеет полный ранг, поэтому указанное выше уравнение имеет смысл и однозначно определяет параметры модели: нам не требуется дополнительная информация для получения единственного решения.
Математические и вычислительные аспекты
Обратные задачи обычно являются плохо поставленными, в отличие от хорошо поставленных задач, с которыми обычно сталкиваются в математическом моделировании. Из трех условий, предложенных Жаком Адамаром для хорошо поставленной задачи (существование, единственность и устойчивость решения или решений), условие устойчивости нарушается чаще всего. В рамках функционального анализа обратная задача представляется отображением между метрическими пространствами. Хотя обратные задачи часто формулируются в бесконечномерных пространствах, ограничения на конечное число измерений и практические соображения, связанные с восстановлением лишь конечного числа неизвестных параметров, могут привести к переформулировке задач в дискретном виде. В этом случае обратная задача, как правило, будет плохо обусловлена. В таких случаях регуляризация может быть использована для введения умеренных ограничений на решение и предотвращения переобучения. Многие примеры регуляризованных обратных задач можно интерпретировать как частные случаи байесовского вывода.
Числовое решение задачи оптимизации
Некоторые обратные задачи имеют очень простое решение, например, когда имеется набор унизольвентных функций, то есть набор из n функций, таких что вычисление их значений в n различных точках дает набор линейно независимых векторов. Это означает, что для любой линейной комбинации этих функций коэффициенты можно вычислить, представив векторы в виде столбцов матрицы и затем инвертировав эту матрицу. Простейшим примером унизольвентных функций являются полиномы, построенные с использованием теоремы об унизольвентности, чтобы обеспечить их унизольвентность. Конкретно, это достигается инвертированием матрицы Вандермонда. Однако это весьма специфическая ситуация. В общем случае решение обратной задачи требует применения сложных алгоритмов оптимизации. Когда модель описывается большим числом параметров (количество неизвестных в некоторых приложениях дифракционной томографии может достигать одного миллиарда), решение линейной системы, связанной с нормальными уравнениями, может оказаться затруднительным. Выбор численного метода для решения задачи оптимизации зависит, в частности, от вычислительных затрат, необходимых для решения прямой задачи. После выбора подходящего алгоритма для решения прямой задачи (прямое умножение матрицы на вектор может быть недостаточно эффективным, если матрица очень большая), подходящий алгоритм для минимизации можно найти в учебниках, посвященных численным методам решения линейных систем и минимизации квадратичных функций (например, в работах Ciarlet или Nocedal). Кроме того, пользователь может захотеть добавить физические ограничения к модели: в этом случае ему необходимо владеть методами оптимизации с ограничениями, которые сами по себе представляют собой отдельную область знаний. Во всех случаях вычисление градиента целевой функции часто является ключевым элементом для решения задачи оптимизации. Как упоминалось выше, информацию о пространственном распределении распределенного параметра можно ввести посредством параметризации. Также можно рассмотреть возможность адаптации этой параметризации в процессе оптимизации. Если целевая функция основана на норме, отличной от евклидовой, необходимо выйти за рамки квадратичной оптимизации. В результате задача оптимизации становится более сложной. В частности, когда норма используется для оценки расхождения между данными и моделью, целевая функция перестает быть дифференцируемой: ее градиент теряет смысл. В этом случае применяются специализированные методы (например, разработанные Lemaréchal) из области недифференцируемой оптимизации. После вычисления оптимальной модели необходимо ответить на вопрос: "Насколько можно доверять этой модели?". Этот вопрос можно сформулировать следующим образом: каков размер множества моделей, которые соответствуют данным "почти так же хорошо", как и эта модель? В случае квадратичных целевых функций это множество содержится в гиперэллипсоиде, являющемся подмножеством (где – число неизвестных), размер которого зависит от того, что подразумевается под "почти так же хорошо", то есть от уровня шума. Направление самой большой оси этого эллипсоида (собственный вектор, соответствующий наименьшему собственному значению матрицы) указывает направление плохо определенных компонент: при движении в этом направлении можно вносить существенные возмущения в модель, не оказывая значительного влияния на значение целевой функции, и, таким образом, получать существенно отличающуюся квазиоптимальную модель. Очевидно, что ответ на вопрос "насколько можно доверять модели" определяется уровнем шума и собственными значениями гессиана целевой функции или, эквивалентно, в случае отсутствия регуляризации, сингулярными значениями матрицы. Разумеется, использование регуляризации (или других видов априорной информации) уменьшает размер множества почти оптимальных решений и, следовательно, повышает уверенность в вычисленном решении.
Некоторые классические линейные инверсные задачи для восстановления распределенных параметров
Упомянутые ниже проблемы соответствуют различным вариантам интегрального уравнения Фредгольма: каждый из них связан с определенным ядром.
Деконволяция
Цель деконволюции — восстановить исходное изображение или сигнал, который выглядит зашумленным и размытым на данных. С математической точки зрения, ядро здесь зависит только от разности между и .
Методы томографии
В этих методах мы пытаемся восстановить распределённый параметр, при этом наблюдение состоит в измерении интегралов этого параметра вдоль семейства линий. Обозначим линией из этого семейства, связанной с точкой измерения . Тогда наблюдение в точке можно записать как:
где – длина дуги вдоль , а – известная функция веса. Сравнивая это уравнение с интегральным уравнением Фредгольма, представленным выше, мы замечаем, что ядро является своего рода дельта-функцией, имеющей максимум на линии . При таком ядре прямое отображение не является компактным.
Компьютерная томография
В рентгеновской компьютерной томографии линии, вдоль которых производится интегрирование параметра, являются прямыми. Томографическая реконструкция распределения параметра основана на обращении преобразования Радона. Несмотря на то, что многие линейные обратные задачи хорошо изучены с теоретической точки зрения, задачи, связанные с преобразованием Радона и его обобщениями, по-прежнему представляют значительные теоретические трудности, а вопросы достаточности данных остаются нерешенными. К таким задачам относятся случаи неполных данных для преобразования рентгеновского излучения в трех измерениях и задачи, связанные с обобщением преобразования рентгеновского излучения на тензорные поля. Среди исследуемых решений – метод алгебраической реконструкции, фильтрованная обратная проекция, а также, с ростом вычислительной мощности, итеративные методы реконструкции, такие как итеративный метод разреженной асимптотической минимальной дисперсии.
Дифракционная томография
Дифракционная томография — классическая линейная обратная задача в сейсморазведке: амплитуда, зарегистрированная в данный момент времени для заданной пары «источник-приёмник», является суммой вкладов от точек, для которых сумма расстояний, измеренных во времени хода волны от источника и приёмника соответственно, равна соответствующему времени регистрации. В трёхмерном случае параметр интегрируется не вдоль линий, а по поверхностям. Если скорость распространения постоянна, такие точки располагаются на эллипсоиде. Обратная задача заключается в восстановлении распределения точек дифракции по сейсмограммам, записанным вдоль профиля, при известном распределении скоростей. Первое прямое решение было предложено Бейлкиным и Lambaré et al.: эти работы стали отправной точкой для методов, известных как миграция с сохранением амплитуды (см. Beylkin и Bleistein). При использовании геометрико-оптических методов (то есть лучей) для решения волнового уравнения, эти методы оказываются тесно связаны с так называемыми методами миграции по методу наименьших квадратов, выведенными из подхода наименьших квадратов (см. Lailly, Tarantola).
Допплеровская томография (астрофизика)
Если мы рассмотрим вращающийся звездный объект, то спектральные линии, которые мы можем наблюдать на спектральном профиле, будут сдвинуты из-за эффекта Доплера. Допплеровская томография ставит своей целью преобразование информации, содержащейся в спектральном мониторинге объекта, в двумерное изображение излучения (как функцию радиальной скорости и фазы периодического вращения) звездной атмосферы. Как объяснил Том Марш, эта линейная обратная задача аналогична томографии: необходимо восстановить распределенный параметр, который был проинтегрирован по линиям для получения наблюдаемого эффекта в записях.
Обратная теплопроводность
Ранние публикации по обратной теплопроводности возникли в связи с определением поверхностного теплового потока при возвращении в атмосферу на основе данных с датчиков температуры, расположенных внутри конструкции. Другие области применения, где требуется знание поверхностного теплового потока, но использование поверхностных датчиков нецелесообразно, включают: внутренние процессы в двигателях внутреннего сгорания, работу ракетных двигателей; и испытания компонентов ядерных реакторов. Для решения проблем не однозначности и чувствительности к погрешностям измерений, вызванным затуханием и запаздыванием в температурном сигнале, разработан ряд численных методов.
Нелинейные инверсные задачи
Нелинейные обратные задачи представляют собой по своей природе более сложный класс обратных задач. В этом случае прямое отображение является нелинейным оператором. Моделирование физических явлений часто основывается на решении частного дифференциального уравнения (см. таблицу выше, за исключением закона всемирного тяготения): хотя сами эти частные дифференциальные уравнения часто линейны, физические параметры, входящие в них, зависят от состояния системы нелинейным образом и, следовательно, от производимых нами наблюдений.
Проблемы обратного рассеяния
В то время как линейные обратные задачи были полностью решены с теоретической точки зрения в конце девятнадцатого века, лишь один класс нелинейных обратных задач был решен до 1970 года – это инверсные спектральные и (в одном пространственном измерении) инверсные задачи рассеяния, после основополагающих работ русской математической школы (Крейн, Гельфанд, Левитан, Марченко). Обширный обзор результатов был дан Чаданом и Сабатье в их книге «Inverse Problems of Quantum Scattering Theory» (два издания на английском языке, одно на русском). В этом типе задач данными являются свойства спектра линейного оператора, описывающие рассеяние. Спектр состоит из собственных значений и собственных функций, вместе образующих «дискретный спектр», а также обобщений, называемых непрерывным спектром. Важный физический момент заключается в том, что эксперименты по рассеянию дают информацию только о непрерывном спектре, и знание его полного спектра является как необходимым, так и достаточным условием для восстановления оператора рассеяния. Таким образом, мы имеем невидимые параметры, гораздо более интересные, чем нулевое пространство, которое обладает аналогичным свойством в линейных обратных задачах. Кроме того, существуют физические движения, при которых спектр такого оператора сохраняется вследствие этого движения. Это явление описывается специальными нелинейными уравнениями эволюции в частных производных, например уравнением Кортевега — де Вриса. Если спектр оператора сводится к единственному собственному значению, соответствующее движение представляет собой распространение одинокого возмущения с постоянной скоростью и без деформации – солитон. Идеальный сигнал и его обобщения для уравнения Кортевега — де Вриса или других интегрируемых нелинейных уравнений в частных производных представляют большой интерес и имеют множество возможных применений. Эта область изучается как раздел математической физики с 1970-х годов. Нелинейные обратные задачи также в настоящее время исследуются во многих областях прикладной науки (акустика, механика, квантовая механика, электромагнитное рассеяние, в частности, радиолокационное и сейсмическое зондирование, а также практически все методы визуализации). Последний пример, связанный с гипотезой Римана, был предложен Ву и Спрунгом: идея заключается в том, что в полуклассической старой квантовой теории обратная величина потенциала внутри гамильтониана пропорциональна половинной производной функции подсчета собственных значений (энергий) n(x).
Соответствие проницаемости в нефтяных и газовых резервуарах
Цель состоит в восстановлении коэффициента диффузии в параболическом уравнении в частных производных, описывающем однофазные течения жидкости в пористых средах. Эта задача привлекала внимание многих исследователей с начала семидесятых годов, со времени пионерской работы в этой области. Что касается двухфазных течений, важной проблемой является оценка относительных проницаемостей и капиллярного давления.