Физика Земли, 2021, № 1, стр. 133-143

2.1. Натурные данные

Мировая база данных по напряжениям WSMP 2016 содержит 42 870 ориентаций главных напряжений, полученных различными методами, в большинстве по механизмам очагов землетрясений. Качество данных различно: так данные качества А указывают ориентации с ошибкой ±15°; качество B – ±20°; качество C – ±25°. Таких данных в базе более 32000, другие данные для моделирования обычно не используются. На рис. 1 показаны ориентации главных сжимающих горизонтальных напряжений в области, ограниченной географическими координатами 36:46 с.ш., 56:76 в.д., т.е. той же, что рассмотрена и в статье И.У. Атабекова. Следует отметить, что число данных высшего качества А всего 4, данных качества В – 16, остальные 381 ориентации имеют качество С. Также следует отметить отсутствие данных в большой части данного региона, что вызывает необходимость рассмотрения нескольких областей меньшего размера. Далее рассматриваются две: первая область ограничена географическими координатами 36:39 с.ш., 56:60 в.д., где находятся 43 данных (показаны в левом нижнем углу на рис. 1); вторая область ограничена географическими координатами 36:44 с.ш., 66:76 в.д., где находятся 334 данных. Следует отметить, что во второй области расположено около ста геологических разломов, в частности глубинные Северо-Тянь-Шаньский, Нарынский, Таласо-Ферганский, Южно-Тянь-Шаньский, Памиро-Гиндукушский и второстепенные разломы [Уломов, 1974], которые могут рассматриваться как активные. Учет всех разломов в моделировании нецелесообразен, поскольку алгоритмы для многосвязных областей используют непрерывность векторов напряжений на разломах в некотором (достаточно большом) числе точек. Если число данных по ориентациям относительно невелико, например, как в данном случае 3–4 данных (в среднем) на каждую подобласть, ограниченную разломами, то количество уравнений из-за условий непрерывности значительно больше, чем число уравнений, порождаемых данными, что приводит к практически непрерывным полям напряжений. Однако непрерывные поля могут быть получены, если пренебречь влиянием разломов, что значительно упрощает расчеты. Ранее проведенные реконструкции полей напряжений в некоторых регионах (см. [Мухамедиев, Галыбин, 2006; Irša, Galybin, 2011]), показали, что существенное различие в результатах реконструкции напряжений возникает, если данных в каждой подобласти десятки, поэтому далее вторая область рассматривается в виде двух связных подобластей с границей вдоль Таласо-Ферганского разлома, как показано на рис. 2. Выбор данного разлома объясняется тем, что он полностью пересекает вторую область, что облегчает проведение моделирования. Несмотря на то, что в отличие от других крупных разломов сейсмическая активность вдоль этого разлома ниже, он, как будет показано в дальнейшем, оказывает существенное геомеханическое влияние на поле напряжений рассматриваемой области.

2.2. Описание алгоритма

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

Предполагается, что вертикальное сжимающее напряжение ${{S}_{V}}$ является главным и определяется весом вышележащих пород. Главные горизонтальные напряжения ${{S}_{H}}$ и ${{S}_{h}}$, действующие в плоскости пласта, являются сжимающими и при этом $\left| {{{S}_{H}}} \right| \geqslant \left| {{{S}_{h}}} \right|$. В соответствии с представлениями Колосова–Мусхелешвили вводятся гармоническая функция $P = \frac{1}{2}\left( {{{S}_{H}} + {{S}_{h}}} \right)$ – среднее давление и комплекснозначная функция девиатора напряжений $D = \frac{1}{2}\left( {{{S}_{H}} - {{S}_{h}}} \right){{e}^{{i{\alpha }}}},$ где аргумент $\alpha = \pi - 2\theta $ выражается через главные направления плоского тензора напряжений – функцию θ. Далее используются выражения напряжений через два комплексных потенциала $\Phi (z)$ и $\Psi (z)$ в виде $P(z,\bar {z}) = 2\operatorname{Re} \left( {\Phi (z)} \right)$, которые подлежат нахождению из следующих условий.

Пусть Ω – односвязная упругая область (конечная или бесконечная), ограниченная несамопересекающимся замкнутым контуром Γ. Эта область может быть подобластью большей области, которая необязательно односвязная или ограниченная. Пусть данные о главных направлениях θ_j тензора напряжений известны в точках z_j (j = 1…N), которые расположены как внутри области, так и могут быть на ее границе.

Требуется определить двумерное поле напряжений в Ω, которое:

− обеспечивает наилучшее согласие данных θ_j и расчетных главных направлений в заданных точках z_j;

− удовлетворяет всем уравнениям плоской теории упругости, сформулированным в напряжениях.

Поскольку главные направления определяются только аргументом комплекснозначной биголоморфной функции девиатора напряжений $D = D(z,\bar {z})$, то задача может быть рассмотрена как задача определения биголоморфной функции D с последующим нахождением гармонической функции $P = P(z,\bar {z})$ из уравнений равновесия. Таким образом, исходная задача эквивалентна нахождению комплексных потенциалов $\Phi (z)$ и $\Psi (z)$, которые обеспечивают минимизацию суммы квадратов отклонений аргумента комплекснозначной функции девиатора напряжений от углов ${{\alpha }_{j}} = \pi - 2{{\theta }_{j}}$, т.е.:

Отметим, что если считать, как это принято в геомеханике, что максимальное главное напряжение есть наибольшее сжимающее напряжение (геомеханическое соглашение), то надо положить ${{\alpha }_{j}} = - 2{{\theta }_{j}}$, что не приводит к изменению в виде (1).

Комплексные потенциалы ищутся в виде многочленов:

Комплексные константы c_k неизвестны и должны определяться из решения задачи (1), после чего аргумент биголоморфной функции D становится известным, как и сама эта функция. Количество комплексных неизвестных, 2n + 1, не должно превосходить половины числа данных N, в противном случае разрешающая система линейных уравнений будет недоопределенной.

Отметим, что можно использовать разное количество членов для комбинаций голоморфных функций в (2). Здесь представление для Φ′(z) имеет на один член меньше, чем для Ψ(z), что не является обязательным, но обеспечивает одинаковое количество членов в функциях напряжений Р и D.

Использование (2) приводит к следующему выражению для функции девиатора напряжений:

Поскольку ${\text{Im}}({{e}^{{ - i{\alpha }}}}D) = 0,$ то подстановка из (3) приводит к функциональному уравнению:

Уравнение (4) далее сводится к системе линейных алгебраических уравнений относительно неизвестных констант c_k. Поскольку функция аргумента α известна в точках z_j, то методом коллокаций получаем следующую систему N реальных уравнений для определения 2n + 1 комплексных констант:

Уравнение (5) является однородным, таким образом, по крайней мере одно дополнительное условие необходимо для того, чтобы найти нетривиальное решение системы. Как видно, умножение системы (5) на любую положительную постоянную не нарушает этой системы, что означает, что любое решение, удовлетворяющее (5), можно нормировать на константу. Эта константа может быть выбрана из условия, что средний модуль D по области равен единице. Так как |D| = e^-iαD, то дополнительное уравнение принимает вид:

Удобно переписать систему (5)–(6) в вещественной матричной форме, записывая комплексные константы как c_k = X_k + iX_{2n +k+ 1} (следовательно, вводя вектор С из 4n + 2 вещественных неизвестных), что приводит к выражению:

где: A – (N + 1)(4n + 2) матрица системы с вещественными коэффициентами A_k,j, которые определены ниже; С – (4n + 2) вектор вещественных неизвестных X_k; b – (N + 1) вектор, все компоненты которого равны нулю за исключением последнего, равного N, как следует из (6):

Переопределенная система (8) решается далее методом наименьших квадратов, при этом для ее разрешимости необходимо, чтобы N ≥ 4n + 2.

В случае, когда область Ω состоит из М подобластей, необходимо записать уравнения (5) для каждой подобласти и добавить к уравнениям (6) условия непрерывности вектора напряжений на всех внутренних границах Γ_j между соседними подобластями, j = 1…N_int. Формально задача сводится к оптимизации целевой функции:

где weight > 0 – некоторый вес, который регулирует взаимный вклад данных и условий непрерывности, а угловые скобки означают скачок величины при переходе через контур. Задача сводится к решению системы вида (7), в которой число уравнений равно N + 1 + KN_int, где K – число точек коллокаций на одном интерфейсе, а число неизвестных М(4n + 2).

Приближенное решение системы (7) по методу наименьших квадратов есть:

Однако этот метод сам по себе не в состоянии обеспечить устойчивость полученного решения, неустойчивость может возникнуть из-за некорректности задачи. Для контроля точности и устойчивости решения используются соответственно невязка R и число обусловленностей CN. При высоком числе обусловленностей используется SVD-регуляризация. Невязка подсчитывается по формуле:

где: D_rec – реконструированная функция девиатора напряжений, вычисляемая по формуле (3) после нахождения коэффициентов c_k.

Как уже отмечалось, D_rec из-за нормировки определяется с точностью по произвольной положительной мультипликативной постоянной. Поскольку ее аргумент определяется однозначно, то эта постоянная относится к модулю комплексного девиатора, который представляет собой функцию максимального касательного напряжения ${{\tau }_{{\max }}}(z,\bar {z})$ в горизонтальной плоскости. Следовательно, среднее давление (которое получается интегрированием уравнений равновесия) зависит от двух параметров, которые не влияют на траектории поля напряжений. Соответственно общее решение в терминах главных напряжений представляется в следующем виде:

где: $\tau (z,\bar {z}) = \left| {{{D}_{{rec}}}(z,\bar {z})} \right|,$ $p(z,\bar {z}) = 2\operatorname{Re} {{\Phi }_{{rec}}}(z)$ являются частными решениями; a и b > 0 – вещественные константы.

Следует отметить, что как только определена функция комплексного девиатора, то по ее аргументу α$(z,\bar {z})$ определяются главные направления, после чего поле траекторий наибольших главных напряжений однозначно находится путем интегрирования системы дифференциальных уравнений первого порядка:

второе семейство траекторий ортогонально к первому. В случае использования “геомеханического соглашения” эти траектории меняются местами. Интегрирование (13) не встречает значительных трудностей и везде далее выполняется численно. Поэтому, когда аргумент D известен в любой точке области, поле траекторий можно считать известным везде, за исключением особых (изотропных) точек, которые могут быть определены в процессе решения.

3.1. Область 1

Расчеты для первой области проводились для n = 0…4, при этом обнаружилось, что невязки при n = 1–3 приблизительно одинаковы и округленно равны 4, при n = 2 невязка минимальна, при n = 4 R = 7, числа обусловленностей невысоки (n = 1: CN = 30; n = 2: CN = 70; n = 3: CN = 166; n = 4: CN = 540). Поэтому все иллюстрации приведены для n = 2. На рис. 3 показаны траектории главных напряжений, а на рис. 4 – изолинии максимальных касательных напряжений в первой области. Средние углы наклона траекторий получаются при использовании n = 0 и составляют 64^о.

Как видно из рис. 3, реконструированные траектории удовлетворительно коррелируют с данными WSMP. Рис. 4 показывает, что ${{\tau }_{{\max }}}$ минимально в точке с координатами 36° с.ш., 58° в.д. и по мере удаления возрастает в любом направлении, достигая максимума в правой верхней части области. Следует отметить, что здесь и далее напряжения на всех графиках приводятся в нормализованной (безразмерной) форме, т.е. в формулах (12) a = 0 и b = 1, при этом среднее значение по области всегда ${{\tau }_{{\max }}}$ = 1.

3.2. Область 2

3.2.1. Без учета разлома. Невязки при n = 1–4 различаются незначительно, с округлением R = 17, числа обусловленности также низкие (n = 1: CN = 46; n = 2: CN = 98; n = 3: CN = 181; n = 4: CN = 319). Ниже приведены результаты для n = 2. На рис. 5 показаны траектории главных напряжений в области 2 без учета разлома и данные WSMP. Несмотря на то, что данные имеют существенный разброс, траектории в целом соответствуют дискретным направлениям. На рис. 6 показана карта изолиний максимальных касательных напряжений в горизонтальной плоскости, из которой видно, что интенсивность ${{\tau }_{{\max }}}$ возрастает в направлении с запада на восток.

3.2.2. С учетом разлома. При расчетах с учетом разлома число точек коллокаций на разломе, в которых выполнялись условия непрерывности вектора напряжений, выбиралось равным половине числа данных, а именно 167. Таким образом, число уравнений, ответственных за непрерывность, совпадало с числом уравнений, использующих данные, т.е. суммарное число уравнений было равным 669, из которых только одно неоднородное, которое обеспечивает условие равенства единице среднего девиатора напряжений (функции D) по всем точкам данных. Число неизвестных коэффициентов в каждой подобласти варьировалось в диапазоне n = 1–4, при этом было установлено, что невязки изменяются незначительно и составляют порядка 17, в то время как числа возрастают или число возрастает от 318 при n = 1, далее 1535 при n = 2, 5687 при n = 3 и до 15 510 при n = 4. Оптимальное число коэффициентов достигается при n = 2, что обеспечивает наилучшую устойчивости к разбросу данных. При этом система вещественных уравнений имеет матрицу размерности 669 × 21.

Результаты расчетов для n = 2 показаны на рис. 7 и рис. 8. На рис. 7 представлены поля траекторий главных напряжений. Как видно из рисунка, траектории напряжений в нижней подобласти качественно совпадают с траекториями, показанными на рис. 5, где разлом не учитывался. Однако они существенно отличаются в верхней подобласти, где имеют другую кривизну, хотя и сохраняют тенденцию загибания на север. На рис. 8 приведены изолинии максимальных горизонтальных касательных напряжений ${{\tau }_{{\max }}}$. Интенсивность ${{\tau }_{{\max }}}$ в нижней подобласти меньше, чем в верхней, но сохраняет тенденцию возрастания в направлении с запада на восток, как и в случае, когда разлом не учитывался (см. рис. 6). В верхней подобласти изолинии практически параллельны разлому, около которого они имеют минимальное значение, а затем быстро возрастают в направлении на северо-восток, где остаются похожими на изолинии, приведенные на рис. 6.

Таким образом, существенное изменение напряжений наблюдается только в окрестности разлома. При переходе через него нормальная и касательная компоненты напряжения остаются непрерывными, а латеральная (тангенциальная) компонента ${{\sigma }_{t}}$ терпит разрыв, также как и максимальное касательное напряжение и оба главных напряжения ${{S}_{H}},{{S}_{h}}$. Разрыв тангенциальной компоненты совпадает с разрывом удвоенных средних напряжений (поскольку нормальная компонента непрерывна), соответственно получаем $\left\langle {{{\sigma }_{t}}} \right\rangle = \left\langle {{{S}_{H}} + {{S}_{h}}} \right\rangle $. На рис. 9 показаны эпюры разрывов латеральных и максимальных касательных напряжений вдоль разлома в зависимости от долготы. Как отмечалось выше оба напряжения даны в безразмерной форме, что соответствует случаю a = 0 и b = 1.