ВведениеВ работе рассматривается подход к решению плоской (двумерной) обратной задачи гравиметрии рудного типа (носителем масс является одиночный изолированный объект с известной эффективной плотностью σ*σ^*) в классе конечноэлементных моделей: исследуемая область среды DD представляется в виде регулярного замощения VV, элементами VjV_j  которого являются горизонтальные четырехугольные призмы бесконечного простирания. Основой рассматриваемого алгоритма является монтажный подход, предложенный в работах [1, 2] и развитый в [3]. Авторами предлагается модификация метода регулируемой направленной кристаллизации (РНК), включающая дополнительный этап в виде решения оптимизационной задачи, направленной на поиск оптимального множества элементов замощения для присоединения к ядру, на определенных этапах итерационного процесса. Рассматриваемая модификация позволяет улучшить содержательность получаемых решений обратной задачи и ускорить вычислительный процесс формирования множества допустимых решений.Монтажный метод решения обратной задачи гравиметрииРНК в классическом виде представляет собой итерационный процесс построения конфигурации Ω*Ω^*, составленной из элементов замощения VjV_j (рис. 1), с соблюдением определенных требований, основными из которых являются условия связности и односвязности носителя аномальных масс. Процесс строится от заданной начальной конфигурации Ω0Ω_0, на каждой итерации i=1,…,ni=1,…,n, к ядру Я(Ωi-1)Я(Ω_{i-1}) добавляется один из элементов оболочки Vj*∈О(Ωi-1)V_j^*∈О(Ω_{i-1}), примыкающий к границе Г(Ωi-1)Г(Ω_{i-1}), образуя новую конфигурацию ΩiΩ_i. Критерием выбора элемента Vj*V_j^* является минимум среднеквадратического отклонения наблюденного и теоретического полей, последнее рассчитывается по элементам Vj*∪Vj∈Я(Ωi-1)V_j^*∪V_j∈Я(Ω_{i-1}). Итерационный процесс завершается при достижении условия σi≤σ*σ_i≤σ^*, где σiσ_i — локально-оптимальное значение эффективной плотности, минимизирующее расхождение наблюденного и теоретического полей на i-й итерации, σ*σ^* — априорно заданная эффективная плотность исследуемого объекта.РНК без ввода ограничений топологического характера на структуру конфигурации ΩΩ имеет недостаток, связанный с формированием в ядре конфигурации малосодержательных фрагментов (отростков) — ветвящихся цепочек элементов замощения. Разработаны алгоритмы [3], позволяющие контролировать степень гладкости получаемого решения. Другая проблема связана с выбором присоединяемого на i-й итерации элемента замощения Vj*V_j^*: выбор конкретного элемента может быть оптимален на данном шаге итерационного процесса, однако при рассмотрении решения задачи в целом этот выбор может оказаться неудачным. Представляет интерес рассмотрение нового варианта РНК с включением в структуру итерации дополнительной оптимизационной задачи, направленной на поиск оптимального подмножества элементов оболочки V*⊂О(Ωi-1)V^*⊂О(Ω_{i-1}) (вместо единственного Vj*∈О(Ωi-1)V_j^*∈О(Ω_{i-1})) для присоединения к ядру Я(Ωi-1)Я(Ω_{i-1}).Рисунок 1. Конфигурация ΩΩ: 1 — элемент замощения VjV_ j, 2 — элементы внутреннего ядра Я(Ω)Я(Ω), 3 — элементы границы Г(Ω)Г(Ω), 4 — элементы оболочки О(Ω)О(Ω). Примечание: красная линия — контур конфигурации ΩМодификация РНК с применением методов эволюционной оптимизацииОдним из вариантов решения вышеперечисленных проблем является использование методов целочисленного программирования. Множеству элементов оболочки Vj∈О(Ωi-1)V_j∈О(Ω_{i-1})  можно поставить в соответствие вектор xi∈{0,1}k,k=|О(Ωi-1)|x_i∈\left\{0,1\right\}^k,k=|О(Ω_{i-1})|, представляющий собой битовую маску и отвечающий за присоединение соответствующего элемента оболочки к ядру Я(Ωi-1)Я(Ω_{i-1}) на i-й итерации. Таким образом, задача поиска оптимального подмножества V*⊂О(Ωi-1)V^*⊂О(Ω_{i-1})  сводится к задаче нахождения битовой маски x=argminx∈{0,1}k⁡F(x)x=\underset{x∈\left\{ 0,1\right\}^k }{argmin}⁡F(x), где F(x)F(x)  — функционал невязки, с последующим определением V*={Vi:xi=1},i=1,…,kV^*=\left\{ V_i:x_i=1\right\},i=1,…,k. Для решения данной задачи используются алгоритм бинарной дифференциальной эволюции (Binary Differential Evolution, BDE) [4] и бинарный генетический алгоритм (Binary Genetic Algorithm, BGA) [5].Алгоритмы BDE и BGA относятся к классу эволюционных стохастических методов многомерной оптимизации. Начальным этапом этих алгоритмов является формирование множества, или популяции, случайных векторов-решений (особей) x10,…,xK0∈{0,1}Nx_1^0,…,x_K^0∈\left\{ 0,1\right\}^N, где KK — размер популяции, NN — размерность задачи. Далее в ходе итерационного процесса происходит изменение состояния исходных векторов с использованием специальных операций мутации, кроссовера и отбора. Результатом работы алгоритма является вектор, обеспечивающий минимум невязки F(x)F(x). Алгоритмы BDE и BGA имеют два общих гиперпараметра C∈[0,1]C∈[0,1] и F∈[0,1]F∈[0,1], регулирующих вероятность возникновения кроссовера и мутации соответственно. Помимо этого, алгоритм BGA имеет еще два дополнительных гиперпараметра: CpC_p — количество точек скрещивания и NtN_t — размер турнирной группы. Для дальнейших расчетов определены оптимальные значения гиперпараметров для поставленной задачи: C=0.4C=0.4 и F=0.025F=0.025 для метода BDE, C=0.2C=0.2, F=0.043F=0.043, Cp=4C_p=4 и Nt=5N_t=5 для метода BGA. Размер популяции для каждого алгоритма K=20K=20, а максимальное количество итераций задано равным 100.Синтетический примерНа рисунке 2 представлены две модели источников поля с постоянной эффективной плотностью σ*=0.3σ^{*}=0.3  г/см3. На профиле длиной 10 км регулярно расположены точки наблюдения с шагом Δx=0.2Δx=0.2 км, для которых рассчитаны значения поля ΔgiΔg_i. Поле ΔgΔg дополнительно осложнено нормальной помехой, среднеквадратическое значение которой составляет приблизительно 4 % и 3 % от максимальной амплитуды аномалии для первой и второй модели соответственно. Для каждого примера заданы ограничения на область поискового пространства DD: 0<Dx<10,0<Dz<40<D_x<10,0<D_z<4, а размер элемента замощения VjV_j квадратного сечения задан равным 0.1 км. Примеры частных решений обратной задачи с использованием модифицированного алгоритма представлены на рисунках 3 и 4. Отметим, что все решения характеризуются значениями невязки F(x)F(x), которые сопоставимы с уровнем помех в наблюдениях.Рисунок 2. Модели № 1 (а) и № 2 (б) источников поля (снизу); отвечающие им гравитационные аномалии ΔgΔg, осложненные помехами (сверху)На рисунках 3а, 4а представлен результат решения обратной задачи монтажным подходом в чистом виде без контроля гладкости границ конфигурации ΩΩ, вследствие чего в структуре Ω*Ω^* присутствуют отростки. Решения, представленные на рисунках 3б, 3в, 3г, 4б, 4в, 4г, получены с использованием вышеописанной процедуры оптимизации. По результатам вычислительных экспериментов установлено, что присоединение подмножества элементов V*⊂О(Ωi-1)V^*⊂О(Ω_{i-1}), найденных в результате оптимизации, вместо единичных элементов Vj∈О(Ωi-1)V_j∈О(Ω_{i-1})  к ядру Я(Ωi-1)Я(Ω_{i-1}) на начальных итерациях или при малой мощности множества О(Ωi-1)О(Ω_{i-1})  приводит к преждевременному завершению итерационного процесса в результате «выхода» эффективной плотности тела на целевое значение σ*σ^*, что связано с избыточным сглаживанием границ конфигурации ΩΩ. В этой связи принято решение включать процедуру оптимизации для определения V*V^* не на каждой итерации (что также не выгодно с точки зрения времени вычислений), а лишь при выполнении условия |О(Ωi-1)|≥m|О(Ω_{i-1})|≥m, где mm — целое число, задаваемое интерпретатором с учетом размеров элементов замощения VjV_j и представлений о структуре или гладкости искомого решения. Решения на рисунках 3б, 3в, 3г и на рисунках 4б, 4в, 4г получены при значении m m = 30, 60 и 90 соответственно. Можно заметить, что степень гладкости границ Ω*Ω^*  зависит от параметра mm: при небольших значениях mm структура Ω*Ω^* характеризуется чрезмерной гладкостью, а при высоких значениях mm наблюдается формирование структурных фрагментов типа отростков. Поэтому для выбора оптимального значения параметра mm необходимо рассмотрение отдельных решений, полученных без включения оптимизационной процедуры в итерационный процесс.Рисунок 3. Частные решения обратной задачи, полученные для модели № 1: классический алгоритм РНК без контроля длины отростков (а); РНК с использованием эволюционных алгоритмов и параметром m, равным 30 (б), 60 (в), 90 (г)Рисунок 4. Частные решения обратной задачи, полученные для модели № 2: классический алгоритм РНК без контроля длины отростков (а); РНК с использованием эволюционных алгоритмов и параметром m, равным 30 (б), 60 (в), 90 (г)Для оценки качества полученных решений обратной задачи рассчитан коэффициент Танимото Ф(S1,S2)Ф(S_1,S_2), отражающий меру сходства двух множеств S1S_1 и S2S_2 и определяемый как Ф(S1,S2)=c(a+b-c)Ф(S_1,S_2)=\frac{c}{(a+b-c)}, где c=|S1⋂S2|c=|S_1⋂S_2 |, a=|S1|a=|S_1 |, b=|S2|b=|S_2 |, S1S_1 — множество элементов Я(Ω*)Я(Ω^*), определенных в процессе решения обратной задачи, S2S_2 — множество элементов замощения VjV_j, принадлежащих истинному объекту.На рисунке 5 представлены графики зависимости коэффициента Танимото от параметра mm конечных решений Ω*Ω^* для каждой модели с использованием рассматриваемых методов BDE и BGA. Как видно, при запуске оптимизационной процедуры при малой мощности множества элементов ядра конфигурации коэффициент Танимото результирующих решений имеет относительно низкие значения, что соответствует малой степени соответствия полученных решений истинной модели. То же самое наблюдается с противоположной стороны: при высоких значениях параметра mm получаемые решения полностью совпадают с конфигурацией, полученной без оптимизационной процедуры, т. е. классическим методом РНК.Рисунок 5. Зависимость коэффициента Танимото Ф(𝑆1,𝑆2)Ф ( 𝑆_ 1 , 𝑆_ 2 ) от параметра mАнализ множества допустимых вариантов решения обратной задачиПостроим множество Q0Q_0 допустимых вариантов решения обратной задачи Ωk*∈Q0Ω_k^*∈Q_0 и исследуем структуру полученного множества путем вычисления значений функции локализации λ(Vj)λ(V_j) [6] на элементах замощения VjV_j. Функция λ(Vj)λ(V_j) характеризует вероятность обнаружения возмущающих масс в отдельной области VjV_j изучаемого пространства DD. Функция λ(Vj)λ(V_j) рассчитана для каждого модельного примера (рис. 2) в трех вариантах: по множеству решений обратной задачи, полученных 1) без включения в итерационный процесс процедуры оптимизации, 2) с включением в итерационный процесс алгоритма BDE с параметром m=60m=60 и 3) с включением в итерационный процесс алгоритма BGA с параметром m=60m=60. Для формирования множества решений Ωk*∈Q0Ω_k^*∈Q_0 использовалась регулярная сетка, элементы которой выступали в качестве центров кристаллизации Я(Ωk0)Я(Ω_{k0}), т. е. начальных приближений. Полученные распределения λ(Vj)λ(V_j) представлены на рисунке 6. Как для первой, так и для второй модели распределение функции λ(Vj)λ(V_j), построенной по результатам решения обратной задачи с использованием модифицированного алгоритма, в большей мере соответствует геометрии источника аномального поля. Без включения в итерационный процесс оптимизационного алгоритма в структуре множества Q0Q_0 проявляются достаточно «резкие» фрагменты, вносящие существенные искажения в общую картину. Можно заметить, что общий вид λ(Vj)λ(V_j) как для первой, так и для второй модели при использовании не модифицированного монтажного подхода во многом соответствует конфигурации Ω*Ω^* единичного решения (рис. 3а, 4а).Рисунок 6. Функция локализации λ(Vj)λ(V_j), построенная на основе множества решений обратной задачи по полю ΔgΔg модели № 1 (слева) и по полю ΔgΔg модели № 2 (справа), полученных классическим алгоритмом РНК (а, г); модифицированным алгоритмом РНК + BDE (б, д); модифицированным алгоритмом РНК + BGA (в, е)На рисунке 7 представлено распределение коэффициента Танимото по множеству допустимых вариантов решения обратной задачи Ωk*∈Q0Ω_k^*∈Q_0: слева для модели № 1, справа для модели № 2. В целом видно, что решения, полученные модифицированным алгоритмом, обладают большим сходством с моделью источника аномального поля.Рисунок 7. Распределение коэффициента Танимото по множеству допустимых решений обратной задачи по полю 𝛥𝑔𝛥 𝑔 модели № 1 (а) и по полю 𝛥𝑔𝛥 𝑔 модели № 2 (б)ВыводыТаким образом, включаемый как дополнительный этап к основному методу, алгоритм оптимизации (поиска подмножества V*⊂О(Ωi-1)V^*⊂О(Ω_{i-1})) выступает в роли регуляризатора, обеспечивая баланс между наблюдениями и представлениями о структурных особенностях изучаемого объекта, и может применяться в качестве альтернативного способа контроля гладкости границ конфигурации ΩΩ. Основным, на наш взгляд, недостатком предложенного подхода является необходимость определения оптимального значения параметра mm, который может резко меняться в зависимости от предполагаемой структуры источника поля и от особенностей используемого замощения. Представляет интерес рассмотрение иных методов оптимизации, помимо рассмотренных алгоритмов BDE и BGA, для решения задачи нахождения подмножества V*⊂О(Ωi-1)V^*⊂О(Ω_{i-1}). Также нужно заметить, что применение данного и ему подобных алгоритмов вполне можно расширить на решение обратных задач для набора аномалиеобразующих тел. Для это следует использовать понятие общей границы и общей оболочки для нескольких моделируемых конфигураций [3]. В целом синтез монтажного подхода и методов эволюционной оптимизации представляется весьма перспективным направлением дальнейших исследований.Исследование выполнено при финансовой поддержке Министерства науки и высшего образования РФ в рамках государственного задания (рег. номер НИОКТР: 126012716041-5).m