Яндекс.Метрика Анатомия появления солнечных пятен и их прогноз
Часть материалов сайта доступна только подписчикам. На период подписки они имеют возможность оперативной консультации по статистическому анализу биомедицинских данных. Запрос на подписку направляйте редактору БИОМЕТРИКИ.

Спектральный анализ численности солнечных пятен

Сарычев В.Т.
Томск, Сибирский физико-технический институт

Введение

    Данные наблюдения солнечных пятен за последние 300 лет с одной стороны указывают на наличие определЈнной периодичности их появления, а с другой стороны иллюстрируют существенную нестационарность этого процесса. Меняются как амплитуды чисел пятен, так и периоды циклов. Следует иметь в виду, что понятие одиннадцатилетнего (двадцатидвухлетнего) цикла лишь приближенно отражает длительность периода вариаций средней за год численности пятен. Наблюдались периоды в 14 и 9 лет. Спектральное оценивание позволяет выявлять не только видимые, но и скрытые периодичности в исследуемых явлениях, поэтому данные наблюдения солнечных пятен неоднократно подвергались спектральному анализу Спектральные оценки чисел солнечных пятен стали традиционными иллюстрациями в работах по спектральному оцениванию. Казалось бы в этих условиях бесперспективно предлагать ещЈ одну работу по анализу динамики числа солнечных пятен, основанную на спектральном оценивании.
    Однако это далеко не так. Спектральные оценки, получаемые большинством современных методов, практически не позволили обнаружить каких либо новых эффектов в динамике изменения их численности. Природу изменение численности пятен до сих пор нельзя считать окончательно выясненной. Считается, что сами пятна образуются в результате сложного турбулентного движения солнечного вещества, которое способно энергию своего движения трансформировать в энергию магнитного поля. Других теорий, кроме теории гидромагнитного динамо, просто не существует. Да, теоретически доказано существование форм трехмерного движения проводящей жидкости, способных увеличивать энергию слабых “затравочных” магнитных полей. Но самосогласованное решение задачи генерации магнитного поля и движения в этом поле проводящей жидкости до сих пор не получено. Т.е. задача решена в кинематической постановке, а не в динамической. Последняя задача слишком сложна. Симптоматично высказывание ведущего специалиста в теории космических магнитных полей Е.Паркера [1]: - “Если магнитное поле Солнца действительно непрерывно воспроизводится динамо-эффектом: неоднородного вращения и циклонической конвекции …”.
    Сами по себе спектральные оценки вариаций солнечных пятен вопрос об их природе, разумеется, не решают, но они могут быть полезны при сопоставлении различных теорий и эксперимента Вне зависимости от природы происхождения солнечных пятен они являются своеобразными легкими Солнца, которые оказывают существенное влияние на процесс обмена веществом и энергией между внутренними и наружными слоями Солнца. Процесс этот происходит хотя и ритмично, но длительности и амплитуды ритмов постоянно меняются. Существуют ли какие-нибудь закономерности в этих ритмах? Можно ли прогнозировать их изменения? Вот те задачи, решению которых посвящена данная работа.

Метод получения спектральных оценок

    Самым важным результатом теории спектрального оценивания за последние три десятилетия следует считать формирование понимания некорректности задачи спектрального оценивания: один и тот же конечный набор данных может быть представлен множеством различных спектров. Иными словами через N точек можно провести сколько угодно кривых, которые отличаются друг от друга не только в пределах окна наблюдения, но и за его пределами, что очень важно для решения задачи прогноза. Какой критерий следует использовать для выбора из этого множества “нужной” кривой?
    Иногда роль критерия может выполнять априорная информация. Но как быть, когда она отсутствует? В 1967 г. Дж. Берг для оценки спектральной плотности мощности (СПМ) предложил использовать принцип максимума энтропии. Очень эмоционально эффект, произведЈнный работой Берга, описан в статье [2]. Когда страсти улеглись, то выяснилось, что оценки СПМ максимума энтропии тождественны авторегрессионым оценкам. Но в данном случае важен был не столько результат, сколько почин. Как указывает Э. Джейнс [3], эти оценки соответствуют автокорреляционной функции, которая может быть реализована наибольшим числом способов.
    К сожалению, СПМ позволяет прогнозировать автокорреляционную функцию процесса, а не сам процесс. Процесс можно прогнозировать, зная его комплексный спектр (КС). Однако не каждая оценка КС годится для прогнозирования. Например, оценкам дискретного преобразования Фурье (ДПФ) соответствует периодическое продолжение сигнала за пределы окна наблюдения. Оценкам конечного преобразования Фурье (КПФ) в пределах полосы частот, равной половине частоты дискретизации сигнала, соответствуют нулевые значения отсчЈтов вне окна наблюдения. Очевидно, что два эти крайние случаи не пригодны для решения задачи прогноза. Кроме того, они дают искаженные представления о спектральном составе сигнала. Можно привести еще много доводов в пользу необходимости поиска новых методов спектрального оценивания, но это бы излишне перегрузило статью. Целесообразней сразу начать с описания метода.
    Наиболее полно метод описан в монографии автора [4] и в статье [5]. В основе метода лежит следующая стохастическая модель сигнала

    Здесь tk – моменты дискретизации сигнала X, Aa- комплексная амплитуда колебания, wa- круговая частота, g – число колебаний.
    Статистические свойства модели (1) описываются функцией плотности распределения колебаний (ФПРК), определЈнной в пространстве частот w и комплексных амплитуд A. В общем случае, при необходимости учитывать корреляцию между различными колебаниями, эта функция многочастичная. В случае, когда основной интерес представляет КС (момент первого порядка ФПРК по амплитуде) можно ограничиться одночастичной ФПРК f(A,w). Эта функция связана с отсчЈтами сигнала X следующим соотношением:

    Интегрирование по A проводится по всей комплексной плоскости, а интервал интегрирования по частоте определяется частотой дискретизации сигнала.
    Для нахождения параметрического вида f(A,w) привлекается принцип максимума энтропии при учете ограничений (2). Однако, чтобы f(A,w) оставалась ограниченной во всей области еЈ определения, необходимо кроме (2) использовать условие конечности энергии сигнала, которое имеет вид

   Удовлетворяющая принципу максимума энтропии функция с учЈтом условий (2) и (3) имеет следующий параметрический вид:


   Здесь l0, l1,..., lM - неопределЈнные множители Лагранжа, обеспечивающие выполнение условий (2), (3), g - среднее число колебаний, S0 - нормирующий множитель, tk – виртуальное время. Выбор значений tk определяется решаемой задачей. Обычно эти значения берутся эквидистантными в пределах окна наблюдения, хотя при решении задачи прогноза некоторые из них могут выходить за эти пределы.
   Три первых момента ФПРК (4) по амплитуде, начиная с нулевого порядка и кончая вторым, представляют следующие функции. r(w) - спектральную плотность колебаний (СПК), X(w) - комплексный спектр, W(w) - спектральную плотность мощности (СПМ). Вид этих функций следующий:

    Здесь, чтобы подчеркнуть, что зависимость от частоты этих функций определяется только значениями множителей l1,..., lM, произведена замена <a(w)>=Ц l0<A(w)>.
    Для разделения переменных l1,..., lMи l0, g подобную замену целесообразно произвести в выражении (2) для компонент модельного сигнала. После интегрирования по комплексной амплитуде это выражение можно представить в следующем виде

    Здесь и далее компоненты сигнала X считаются вещественными. Все множители Лагранжа для таких сигналов вещественны. Вместо более громоздких квадратурных косинусной и синусной компонент используется комплексная экспонента. Значок Re будет далее опускаться.
    Среднее по ФПРК (4) значение квадратов компонент сигнала определяются выражением

    Соответствующее этому выражению средне квадратичное отклонение сигнала определяется формулами

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

    Поиск дна оврага осуществляется покоординатным спуском. Другие способы, например, градиентный, при испытании показали свою бесперспективность. Выбор начальной точки осуществлялся следующим образом. Значения lj оцениваются путЈм решения переопределЈнной системы уравнений

    где B – прямоугольная матрица MґN (N>M) с элементами Bjk = T Sin(p(tj-tk)/T)/p(tj-tk); Nі kі 1, Mі jі 1. При эквидистантном выборе tj значение эффективного интервала дискретизации T=(tM-t1)/(M-1). Значение e выбирается малым настолько, чтобы зависимость КС от l1,...,lMбыла близка к линейной. Для выбранного значения e вычисляются значения l1,...,lM. По этим значениям согласно (7) вычисляются компоненты <xk>, затем вычисляется c по формуле c=XЧ x/|x|2, которая позволяет минимизировать значение невязки. После этого вычисляются значения l0 и g согласно выражениям

    Эти выражения получены на основе (9) и (10) в предположении равенства выборочной невязки F усредненной <F>.
    Значения l1,...,lM , полученные в результате решения системы (11), для большинства сигналов обеспечивают минимальное значение невязки. Почему же нельзя ограничиться этим решением? Дело в том, что КС, соответствующий этому решению, практически совпадает с традиционным конечным преобразованием Фурье, а поэтому наследует все дефекты этого преобразования. Изобилует лишними деталями в виде боковых лепестков, обладает невысокой разрешающей способностью, а главное не годится для решения задачи прогноза. Такая оценка КС обладает малой информацией. Но, чтобы говорить об информации, еЈ нужно уметь оценивать Используемый в данной работе аппарат описания спектральной структуры сигнала позволяет вычислять информацию, соответствующую полученной оценке ФПРК.
    Следуя С. Кульбаку [6], оценку количественной меры информации, приходящейся на одно колебание, предлагается производить на основе выражения

    Где альтернативная ФПРК соответствует белому гауссовскому шуму и определяется выражением


    Параметры этой функции определяются энергией сигнала и условием нормировки.
    Повышение информативности спектральной оценки производится следующим образом. Методом покоординатного спуска находится вектор l, который минимизирует невязку F. При удачном выборе e этот вектор совпадает с решением системы (11). Затем вектор l умножается на некоторый коэффициент b>1, и процедура поиска минимума невязки повторяется. Значение b выбирается таким, чтобы соответствующее новому l значение невязки несущественно отличалось от предыдущего минимального значения. Процесс повторяется многократно до тех пор пока не прекратится рост значения Ik. Опыт показывает, что для большинства сигналов в начале итерационного процесса имеет значение Ik близкое 1. Затем это значение возрастает до 5-7. У очень “хороших” сигналов с низким уровнем шума эта величина может достичь 12-15. Таким образом, Ik является интегральным параметром сигнала, характеризующим как сам сигнал, так и качество оценки его спектра. К сожалению, работы по оценке Ik для реальных сигналов естественного происхождения в научной литературе отсутствуют.

Структурные элементы сигнала

    Большинство методов спектрального оценивания непрерывных сигналов позволяют аппроксимировать эти сигналы и их комплексные спектры суммами

    Качество прогноза в большой степени определяется видом функций Stk, называемых далее структурными элементами, и их КС Swk . Чем лучше адаптированы структурные элементы к анализируемому сигналу, тем качественней будет прогноз.
    Так структурными элементами ДПФ является набор гармонических колебаний. Во временном базисе эти элементы представляются комплексными экспонентами Stk=exp(iwkt), а в частотном Swk=d(w-wk), где d(w-wk) - - дельта-функция Дирака. Структурными элементами КПФ являются функции Stk=Sin(p(t-k))/ p(t-k), спектральный образ которых имеет следующий вид: Swk=exp(-iwk)( q(w)-q(p-w)), где q(w) - функция Хевисайда. Здесь время и частота безразмерные величины, для которых частота дискретизации полагается равной 1.
    Для описанного в предыдущем разделе метода спектрального оценивания в качестве структурных элементов используются следующие функции:

    Здесь r(w) - СПК, определяемая выражением (6).
    Коэффициенты разложения Ck в выражениях (15) равны Фурье коэффициентам в случае ДПФ, равны дискретным отсчЈтам сигнала Xk для КПФ и Ck=lk/Ц l0 для структурных элементов (16). Все структурные элементы ДПФ Stk имеют общий период, равный длительности интервала наблюдений. Структурные элементы КПФ убывают за пределами окна наблюдения по закону Sin(x)/x. Структурные элементы обоих методов совершенно не зависят от сигнала, поэтому эти методы не годятся для решения задачи прогноза.
    В то же время структурные элементы (16) определяются самим сигналом. Во временном базисе структурные элементы Stk и Stj имеют одинаковый вид и лишь смещены друг относительно друга по временной оси на величину k-j. Характерная длительность элементов Ts, а следовательно длительность прогнозируемого интервала, полностью определяются видом функции r(w).
    Между полосой частот Dw, в которой сосредоточена основная энергия исследуемого сигнала, интервалом значений виртуальных отсчЈтов времени Dt=tM-t1, который определяет разрешающую способность спектральных оценок, и M - числом неопределЈнных множителей Лагранжа имеется связь: Dw Dt=p(M-1). Учитывая эту связь, можно сказать, что вещественные сигналы, для которых Dw меньше половины частоты дискретизации, допускают прогнозирование по крайней мере в интервале времени Dt=p(M-1)/Dw. Если спектр сигнала линейчатый, и ширина линий спектра dw достаточно малая (dw<Dw/M), то характерная длительность структурного элемента Ts@ p/dw>Dt. Как будет показано ниже эта ситуация имеет место при прогнозировании численности солнечных пятен.

Примеры прогноза и спектральных оценок динамики численности солнечных пятен

   Для расчЈтов в качестве исходных использовались данные, взятые с сервера SIDC (Sunspot Index Data Center, daily и yearly sunspot number). Ниже представлены результаты расчЈтов прогноза и спектров для четырЈх различных вариантов. В первых трех вариантах в качестве исходных данных использовались ежедневные (daily) данные, в четвЈртом - среднегодовые (yearly).
   Из ежедневных данных специальной фильтрацией, соответствующей частоте дискретизации с периодом 1 месяц, были получены месячные данные, которые далее использовались для расчЈтов спектров и прогноза. Месячные данные SIDC не использовались для расчЈтов, поскольку применяемые при их получении процедуры сглаживания несколько искажали информацию о спектральном составе исходного сигнала.
   1)Целью первого варианта являлось тестирование алгоритма прогноза. В качестве исходных данных использовались 411 отсчЈтов месячных данных численности солнечных пятен, начиная с августа 1943 года. По этим данным оценивался спектр, по которому далее проводилась пролонгация данных на 20 лет вперед за пределы окна в 411 месяцев. Результаты этих расчетов представлены на Рис.1. ЧЈрная линия представляет исходные данные, красная и синяя линии - данные восстановленные по оценкам спектра в различных частотных интервалах.

    2) Второй вариант расчЈтов был полностью аналогичен первому за исключением выбора положения окна наблюдений: 411 отсчЈтов исходных данных брались в интервале от января 1964 года по ноябрь включительно 1998 года По этим данным оценивался спектр, по которому далее восстанавливался сигнал, как в пределах окна наблюдений, так и на 20 лет вперЈд за его пределы. Результаты расчЈтов представлены на Рис.2, где так же как и на предыдущем рисунке, чЈрная линия представляет исходный набор данных, а красная и синяя - модельные.

    3) Вариант 3 отличался от первых двух длительностью интервала определения исходных данных. Этот интервал включал 684 отсчЈтов месячных значений численности солнечных пятен - от января 1942 года по декабрь 1998 года включительно. Интервал прогнозирования составлял 20 лет. Результаты расчета представлены на рис.3, где кроме исходных данных (чЈрная кривая) и модельных (красная и синяя), представлен прогноз NASA[7].

    4) В четвЈртом варианте в качестве исходных использовались среднегодовые значения чисел солнечных пятен, начиная с 1700 г. по 1999 г. Численность за 1999 была взята на основании результатов прогноза по месячным данным. Результаты предварительных расчЈтов показали, что ряд исходных данных из 299 отсчЈтов следует дополнить хотя бы одним отсчЈтом, чтобы произошЈл "захват" следующего 23 цикла. Результаты расчЈта показаны на Рис.4, где также как и на трЈх предыдущих рисунках чЈрная линия представляет исходные данные, а красная и синяя - модельные.

    Как видно из рисунков, модельные (синие линии) с приемлемой точностью совпадают с линиями, отображающими исходные данные, как в пределах окна наблюдений, так и на пролонгируемом участке. Среднее стандартное отклонение в пределах окна наблюдения было ~4%. Отмечается превышение ~10% максимума сглаженных (красная линия) данных для первого прогнозируемого цикла Различие отмечается лишь на уровне мелкомасштабных возмущений длительностью около двух месяцев. Это различие определяется полосой частот, в которой производится оценка спектра. В свою очередь ширина полосы частот определяется частотой дискретизации виртуальных отсчЈтов времени tk. Число неопределенных множителей Лагранжа, а следовательно, и виртуальных отсчЈтов времени бралось 300 в первых трЈх вариантах. Дискрет виртуальных отсчЈтов времени T=(tM-t1)/(M-1) в первых двух вариантах брался равным 59/30 месяца и 923/300 месяца в третьем варианте. В четвЈртом варианте число неопределенных множителей Лагранжа равнялось 290, а дискрет виртуальных отсчЈтов времени T=549/289 года.
    Положение максимумов и минимумов в циклах в первых трЈх вариантах удобно определять по красной линии, представляющей модельные значения сигнала, восстановленные по спектру в полосе частот [0-0.15 г-1]. Ниже в таблице приводятся моменты появления максимумов (tmax) и соответствующие значения Nmax для всех четырЈх вариантов расчЈтов.

Варианты I II III IV
Nцик Tmax Nmax tmax Nmax tmax Nmax tmax Nmax
18 11.1948 168 - - 11.1948 159 1948 148
19 05.1958 218 - - 06.1958 209 1958 186
20 03.1969 122 05.1969 118 05.1969 111 1969 109
21 11.1980 216 прг 12.1980 177 10.1980 185 1980 157
22 06.1990 177 прг 06.1990 169 08.1990 174 1990 156
23 - - 04.2001 206 прг 08.2001 180 прг 2003 148 прг
24 - - 11.2011 157 прг 03.2012 126 прг 2015 96 прг

    При анализе этой таблицы следует учитывать, что среднегодовые данные в максимумах циклов всегда меньше месячных данных. Если посмотреть значение Nmax первого прогнозируемого цикла в I варианте (помечено "прг"), то видно превышение этого значения над соответствующими значениями в двух следующих вариантах на величину ~16%. Этот эффект подобен явлению Гиббса при аппроксимации функции частными суммами рядов Фурье. В данном случае чередование максимумов и минимумов численности солнечных пятен определяется результатом интерференции многих спектральных линий. Причем линии, соответствующие значениям частот большим 0.15 г-1, определяют мелкомасштабные флуктуации на фоне в достаточной степени гладкого одиннадцатилетнего цикла. В свою очередь, число спектральных линий, которые могут быть разрешены, определяется интервалом наблюдений.
    На Рис.5 приведен амплитудные спектры солнечных пятен для первых трЈх вариантов. Интервал наблюдений для первых двух вариантов был 411 месяцев. Число спектральных линий в полосе частот [0-0.15 г-1] для этих вариантов оказалось равным 3 и 4, соответственно. Интервал наблюдений для третьего варианта был больше (684 месяца), и число линий выросло до 7.

    Если предположить, что явление Гиббса в первых 2-х вариантах проявлялось в равной степени, то значение Nmax=206 для 23 цикла варианта II следует уменьшить до 170-175. В третьем варианте в силу большей длительности интервала наблюдений явление Гиббса должно быть слабее (количество основных интерферирующих линий возросло с трех до 7), так что значение Nmax=180 следует уменьшить лишь на величину ~10, полученный результат согласуется с прогнозом 23 цикла второго варианта.
    Как видно из данных таблицы, увеличение интервала наблюдения ведет не только к ослаблению явления Гиббса, но и к смещению положения максимума. Это и следовало ожидать, поскольку при этом происходит изменение количества линий, их положения и формы профилей. ВсЈ это меняет интерференционную картину. Изменяется не только амплитуда циклов, но и их положение на временной оси. Максимальная длина интервала наблюдений для солнечных пятен на сегодняшний день составляет 299 лет. Можно ожидать, что положение интерференционных максимумов, вычисленное по спектральным оценкам для среднегодовых значений численности солнечных пятен, будет отличаться от соответствующих значений, найденных по месячным данным. Это на самом деле имеет место.

    На Рис.6 приведен амплитудный спектр, построенный по среднегодовым данным. Как видно из рисунка, увеличение длительности интервала наблюдений привело к существенному обогащению спектра. В полосе частот [0-0.15 г-1] количество линий возросло от 7 до 34. Толщина линий существенно уменьшилась, соответственно возросла амплитуда. Очень четко проявилась линия, соответствующая периоду колебаний 106 лет. Ширина этой лини оказалась всего 0.0004г-1. Вся спектральная картина распадается на мультиплетные полосы. Одиннадцатилетнему циклу соответствует мультиплет из 9 линий. Следует сказать, что ни один из современных методов не позволяет получать спектральные оценки сравнимые по качеству с оценками, приведенными в данной работе. Достигаемое здесь разрешение спектральных линий позволяет продолжить интерференционную картину далеко за пределы окна наблюдения.
    Особого внимания заслуживает низкочастотные составляющие спектра с периодом более 50 лет (полоса частот [0-0.02 г -1]). Как следует из Рис.4, в настоящее время мы находимся в начале ниспадающего участка восстановленного по этой полосе временного хода численности солнечных пятен (красная линия). Подобная ситуация наблюдалась в окрестности 1800 года, когда в максимумах соседних циклов численность пятен уменьшилась со 132 до 47. Интервал времени между этими максимумами был 17 лет, а между минимумами 14. Согласно прогнозу по месячным данным (3-й вариант) максимум 23 цикла должен отстоять от максимума предыдущего цикла .на 11 лет. Тогда как прогноз по среднегодовым данным этот интервал увеличивает до 13 лет, а прогнозируемая длительность 23 цикла (от минимума до минимума) составит 14 лет. Резкого уменьшения амплитуды цикла не произойдет. Оно наступит лишь в 24 цикле, где среднегодовая численность пятен снизится до 93.
    Как говорилось выше, качество прогноза в большой мере определяется видом структурных элементов, которые используются для аппроксимации сигнала. На Рис.7 представлены структурные элементы для первых трЈх вариантов. Наблюдается чЈтка выраженная структура 11-и летнего цикла. Однако эта структура "зашумлена" мелкомасштабными вариациями, которые приводят к замыванию цикличности с ростом времени. В начальном участке на протяжении ~25 лет структурные элементы всЈх трЈх вариантов практически совпадают.

    На Рис.8 приведены структурные элементы, построенные по усредненным годичным данным численности солнечных пятен. Здесь так же отмечается одиннадцатилетняя цикличность, которая модулирована низкочастотными составляющими спектра и результатом интерференции линий, представляющих тонкую структуру одиннадцатилетнего цикла.
 
 

Выводы

   Полученные в работе результаты позволяют предположить, что динамика численности солнечных пятен в настоящее время претерпевает существенные изменения: на протяжении последующих 30 лет будет наблюдаться сильное монотонное понижение низкочастотной компоненты численности пятен (периоды более 50 лет); в результате интерференции линий мультиплета одиннадцатилетнего цикла следует ожидать заметного увеличения длительности 23 цикла до 14 лет, при этом интенсивность цикла в максимуме будет сравнима с интенсивностью предыдущего 22-го цикла, или несколько ниже его, но интенсивность последующего 24 цикла может резко упасть (почти в два раза).
 
 

Литература

  1. Паркер Е. Космические магнитные поля их образование и появления. М.: Мир 1982.
  2. Робинсон Э.А. История развития теории спектрального оценивания//ТИИЭР. - 1982. - № 9 - С.6-33.
  3. Джейнс Э.Т. О логическом обосновании методов максимальной энтропии// ТИИЭР. - 1982. - № 9 - С.33 - 51.
  4. Сарычев В.Т. Спектральное оценивание методами максимальной энтропии. - Томск: Изд-во ТГУ. - 1994. - 257 с.
  5. Сарычев В.Т. Совместное оценивание спектров нулевого, первого и второго порядков конечных последовательностей//Радиотехника и Электроника. - 1995. - № 9. С.1414-1421.
  6. Кульбак С. Теория информации и статистика. М: Наука. - 1967. - 408 с.
  7. Anderson J., Coffey H, Harvey K, Hathaway D., Heckman G., Hildner E., Mende W., Schatten K., Thompson R., Thompson A.W.P., and White O.R. Solar Cycle 23 Project: Summary of Panel Findings. http://science.msfc.nasa.gov/ssl/pad/solar/predict.htm
Наш адрес:
1997 - 2017.© Василий Леонов

Возврат на главную страницу.

Возврат в "Экспериментальную биомедицинскую гаврилиаду".