Обратные задачи

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

Загвоздка в том, что шум записи превращается в ошибку разреза, и превращается тем сильнее, чем мельче деталь, которую мы пытаемся разглядеть. Это не недостаток приборов. Отклик от тонкого пропластка приходит на поверхность ослабленным в тысячи раз, а восстановление идёт обратным ходом: слабый сигнал приходится умножать на большое число — вместе с ним умножается и шум.

Новосибирская школа Михаила Михайловича Лаврентьева‑младшего выросла именно из этого: обратная задача — не «прямая, решённая наоборот», а задача другого класса, для которой понадобилась отдельная теория. Вывод у теории короткий и практичный: точное решение искать бессмысленно, задачу надо сознательно загрубить — ровно настолько, насколько велик шум в данных. Приём называется регуляризацией, и мера загрубления в нём одна — число $\alpha$. Как его выбрать, разобрано в разделе «Глубже».

Задачу называют корректной по Адамару, если выполнены три условия: решение существует, оно единственно и оно устойчиво — малое изменение данных даёт малое изменение ответа¹. Прямые задачи физики обычно корректны: задали источник — получили поле. Обратные ломают третье условие, и ломают принципиально.

Проще всего это видно на примере самого Адамара — задаче Коши для уравнения Лапласа. Зададим на границе возмущение с амплитудой $1/n$:

$$ u(x,0)=0,\qquad u_y(x,0)=\tfrac{1}{n}\sin(nx) $$

модель

Решение выписывается точно:

$$ u(x,y)=\tfrac{1}{n^2}\sin(nx)\,\operatorname{sh}(ny) $$

модель

Возьмём частоту побольше. Данные при росте $n$ стремятся к нулю — амплитуда $1/n$ исчезает. А решение из-за гиперболического синуса растёт как $e^{ny}$: на любой глубине $y$ множитель $e^{ny}$ обгоняет $1/n^2$. Данные неотличимы от нуля — ответы уходят в бесконечность.

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

Для сейсморазведки отсюда следует практический предел. Отклик от пласта толщиной в метр на глубине двух километров тонет в шуме приборов, и никакое улучшение алгоритма его оттуда не достанет — информации в записи просто нет. Разрешающая способность метода ограничена не программой, а физикой: длиной волны и уровнем шума.

Где работает тот же закон

Всюду, где смотрят внутрь, не вскрывая: компьютерная томография и УЗИ, сейсморазведка, зондирование атмосферы со спутника, астрономия (восстановить источник по излучению), обработка сигналов. Везде малый шум умножается на большое число — и везде спасаются одинаково.

Пусть задача записана как $Ax=y$: оператор $A$ переводит искомую причину $x$ в наблюдаемое следствие $y$. Неустойчивость видна из разложения $A$ по сингулярным числам $\sigma_1\ge\sigma_2\ge\dots\to 0$. Обратный оператор делит на $\sigma_k$, и компоненты с малыми $\sigma_k$ — как раз мелкие детали — усиливаются в $1/\sigma_k$ раз. Отсюда и способ лечения: не дать знаменателю подойти к нулю.

Метод Лаврентьева. Уравнение первого рода заменяют уравнением второго рода, добавив к оператору малую единицу:

$$ (A+\alpha I)\,x_\alpha=y,\qquad \alpha>0 . $$

модель

То есть к оператору прибавляют малую добавку $\alpha$, и деление на почти нулевое сингулярное число заменяется делением на «почти нулевое плюс $\alpha$» — усиление шума перестаёт быть неограниченным.

Знаменатель становится $\sigma_k+\alpha$ и ниже $\alpha$ не опускается. Метод применим, когда $A$ действует в одном пространстве и самосопряжён и неотрицателен (шире — аккретивен): иначе сумма $\sigma_k+\alpha$ не возникает и сдвиг ничего не гарантирует.

Регуляризация Тихонова. Минимизируют невязку со штрафом за величину решения:

$$ |Ax-y|_Y^2+\alpha|x|_X^2\ \to\ \min_x . $$

модель

Словами: ищем решение, которое одновременно не слишком расходится с данными (первое слагаемое) и не слишком дико само по себе (второе). Параметр $\alpha$ задаёт, чем из двух мы больше дорожим.

Условие минимума даёт нормальное уравнение

$$ (A^{}A+\alpha I)\,x_\alpha=A^{}y . $$

модель

Иначе говоря, минимум того выражения находится решением этой системы — и в ней снова стоит та же добавка $\alpha$, что и в первом способе.

Вот и вся разница. Оба метода делают одно и то же: сдвигают спектр на $\alpha$, чтобы деление на малое число заменилось делением на не слишком малое. Тихонов применяет тот же сдвиг не к $A$, а к $A^{}A$ — и потому работает для любого ограниченного оператора, без требований симметрии. Плата за общность видна в спектре: у $A^{}A$ сингулярные числа равны $\sigma_k^2$. Отсюда две оговорки, из-за которых методы всё же не совпадают:

• при равном $\alpha$ Тихонов подавляет мелкие детали сильнее — там, где у Лаврентьева знаменатель $\sigma_k+\alpha$, у Тихонова $\sigma_k^2+\alpha$; — порядок точности разный. На подходящем классе решений ошибка метода Лаврентьева ведёт себя как $O(\delta^{2/3})$, метода Тихонова — как $O(\delta^{1/2})$, где $\delta$ — уровень шума. Поэтому, когда оператор позволяет (самосопряжённый неотрицательный — типично для задач теплопроводности и диффузии), берут Лаврентьева; когда не позволяет — Тихонова¹.

Как выбирать $\alpha$

Ответ проще, чем кажется, и он один на оба метода. Пусть данные измерены с известной погрешностью $\delta$, то есть $|y_\delta-y|\le\delta$. Тогда $\alpha$ берут из условия, чтобы невязка решения сравнялась с уровнем шума, — это принцип невязки Морозова:

$$ |A x_\alpha - y_\delta| = \delta . $$

модель

Читается так: сглаживаем ровно настолько, чтобы расхождение с данными сравнялось с известным уровнем шума. Требовать меньшего расхождения бессмысленно — мы начнём подгонять решение под шум.

Смысл прямой: бессмысленно подгонять модель точнее, чем измерены данные. Невязка меньше $\delta$ — мы объясняем решением собственный шум и получаем ту самую рябь. Больше — недоиспользуем данные и теряем разрешение. Уравнение относительно $\alpha$ монотонно и решается подбором за десяток итераций.

Отсюда же правило на случай, когда приборы улучшили. Чтобы решение сходилось к истинному при $\delta\to 0$, нужно $\alpha(\delta)\to 0$, но медленнее, чем растёт усиление шума: условие $\delta^{2}/\alpha\to 0$. Проще говоря, точность измерений выросла вдвое — загрубление можно ослабить, но не вдвое, а осторожнее.

Когда $\delta$ неизвестна, берут апостериорные правила: L-кривую (излом графика «норма решения против невязки» в логарифмических осях) или перекрёстную проверку. Они слабее принципа невязки и на некоторых классах задач сходимости не дают — поэтому уровень шума стараются измерять, а не угадывать.

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

• Слишком слабое сглаживание — решение пляшет: малое изменение данных даёт совсем другой ответ. Это исходная беда неустойчивой задачи, никуда не девшаяся. — Слишком сильное — решение гладкое, устойчивое и неверное: детали, которые в данных были, оказываются затёртыми. Формально ошибка распадается на два слагаемых, растущих в разные стороны, и минимум суммы достигается в единственной точке. — Выбор параметра — отдельная задача. Способ невязки требует знать уровень шума в данных; L-кривая ищет излом на графике «невязка против нормы решения»; перекрёстная проверка отбрасывает часть измерений и смотрит, предсказывает ли решение отброшенное. Все три дают разные ответы на одних и тех же данных, и универсального критерия нет.

Сколько это в числах Степень неустойчивости измеряют числом обусловленности — отношением наибольшего сингулярного числа к наименьшему, $\sigma_1/\sigma_n$. Словами: во сколько раз задача усиливает ошибку данных. Для задач томографии оно достигает $10^{6}$–$10^{8}$: это значит, что шум в данных на уровне 0,1 % без регуляризации превращается в ошибку решения в сотни процентов — то есть в бессмыслицу. У классической некорректной задачи обратной теплопроводности сингулярные числа убывают экспоненциально, и тогда сколько ни улучшай приборы, вглубь по времени продвинуться почти невозможно: каждая дополнительная цифра точности данных даёт лишь логарифмический выигрыш.

О чём спорят Насколько вообще правомерно вводить в решение внешнее предположение о гладкости — вопрос, по которому единого мнения нет. Классическая регуляризация Тихонова штрафует большие производные, то есть заранее считает решение плавным; но во многих задачах — поиск границы пласта, восстановление изображения с резкими краями — правильный ответ как раз разрывный, и такой штраф систематически его портит. Отсюда конкурирующие подходы: штраф на полную вариацию, разреженные представления, а в последние годы — обучаемые регуляризаторы, где вид штрафа берётся из данных. Спор о том, что считать «разумным ожиданием» о решении, остаётся открытым, и это спор не технический, а о границе между измерением и допущением.

Открытый вопрос

Всё сказанное относится к линейным задачам. В нелинейной постановке — а полноволновая инверсия в сейсморазведке именно такова — функционал невязки перестаёт быть выпуклым: у него множество локальных минимумов, отвечающих разрезам со сдвигом на целую длину волны. Выбор $\alpha$ тогда перестаёт быть отдельной задачей и сплетается с выбором начального приближения и частотного диапазона, а гарантий сходимости к истинному разрезу в общем случае нет.