Итак, когда в начале учебного года более-менее установилось учебное расписание, я связался со своим научным руководителем Сергеем из Института солнечно-земной физики, в дальнейшем просто ИСЗФ.

Здесь и началось всё веселье

В прошлом году я писал вместе с ним курсовую по обработке радиоизображений с помощью алгоритма CLEAN. Данный алгоритм широко используется для устранения шумов с сырых изображений, полученных с радиотелескопов и, в частности, с антенных решёток. Он, с одной стороны, совсем несложный, с другой - нужен для понимания того, как интерпретировать правильно картинки с радиотелескопов. В курсовой я попробовал с нуля запрограммировать алгоритм, поиграться с параметрами и обработать парочку тестовых картинок с сибирского радиогелиографа СРГ-48.

Красивые картиночки и код с моей курсовой можно посмотреть на Гитхабе: https://github.com/vit1-irk/clean_lib

Вернёмся к настоящему

В ИСЗФ мне понравилось, там достаточно интересно и прикольно (а ещё Сергей линуксоид, прогер, да и в целом с ним приятно было работать), поэтому я захотел большего и решил, что хочется принять участие в каком-нибудь более крутом проекте, уже связанным с реальной физикой. Но при этом включающем в себя немного программирования и обработки данных. В целом говоря, интересна идея изучать Солнце, и рано или поздно, хочу устроиться в ИСЗФ на работу.

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

Использовать надо было данные с Нобеямской радиообсерватории в Японии, где были данные аж с 1999 года по 2017. Более поздние тоже имеются, но они порченные, кривые и не подходят для обработки.

Данные Нобеямы включают в себя корреляционные кривые телескопа (показывают общее изменение яркости, которое потом отразится на изображении) и сами картинки, которые надо сначала скачать, а потом синтезировать с помощью специальной программки, которая написана на языке IDL (среда программирования, популярная у солнечников).

Самый первый поиск источников

Начальная фильтрация - грубый поиск образований высокой яркостной температуры (выше 150 000K), при этом на корреляционной кривой сигнал не должен превышать шумовой порог. Изображения делались только раз в день, поэтому многие из источников, вероятно, были пропущены. Ищутся яркие точки, которые долго остаются на диске Солнца, по несколько часов.

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

Отфильтровалось несколько десятков дней, для каждого был построен видеоряд (см. в самый конец, где код).

После этого потребовалось выяснить конфигурацию магнитных полей на фотосфере Солнца в нужных активных областях, магнитограммы надо было доставать со спутников Solar Dynamics Observatory и Hinode .

Магнитограммы

Магнитограмма SDO

Увеличенная

Магнитограмма Hinode

Фишки Питона

Пишу код и провожу вычисления в интерактивной среде разработки Jupyter Notebook (а точнее - в сборке Jupyter Lab). Очень удобная, позволяет избегать ошибок и строить графики с другими результатами прямо в коде.

Так как данных очень много, приходилось делать долгие промежуточные вычисления, результаты которых для удобства сохранял в бинарных объектах pickle, чтобы не высчитывать каждый раз снова.

# пример записи в файл
with open("150k-curves.obj", "wb") as dump:
    pickle.dump(curves_filtered, dump)

# пример считывания из файла
with open("magnetic-plots.obj", "rb") as dump:
    magnetic_plots = pickle.load(dump)

Построение крутых видосиков

А вот код для построения видеоряда в течение дня. Сохранять картинки в png и склеивать в видео оказалось накладно, поэтому отрисовка идёт напрямую в ffmpeg (через matplotlib FFMpegWriter).

Рендеринг первым способом занимает 4 минуты на видео, вторым - 2 минуты. Надеюсь, кому-нибудь пригодится.

    plt.rc('font', size=12)
    plt.rc('figure', titlesize=16)

    plt.close()
    fig = plt.figure(figsize=(22, 20))
    gs = fig.add_gridspec(nrows=4, ncols=2, width_ratios=[5, 4], hspace=0.4, wspace=0.1)

    magnplot = fig.add_subplot(gs[0, 0])
    ccplot = fig.add_subplot(gs[1, 0])
    maxplot = fig.add_subplot(gs[2, 0])
    maxplot_17 = fig.add_subplot(gs[3, 0])

    picplot = fig.add_subplot(gs[2, 1])
    picplot_17 = fig.add_subplot(gs[3, 1])

    magnplot.set_title("maximum magnetic field by SDO")
    magnplot.plot(magn["times"], magn["maxvals"], ":y")
    magnplot.plot(magn["times"], magn["maxvals"], "o")

    ccplot.set_title("34 GHz correlation curve")
    ccplot.set_ylim(-std * 0.3, std * 1.1)
    ccplot.plot(cc["times"], cc["data"] - np.mean(cc["data"]))
    ccplot.axhline(std)

    maxplot.set_title("34 GHz, maximum value")
    maxplot.set_ylim(-100, cfg.intensity_threshold * 1.5)
    maxplot.plot(pics_times, maxvals, "-o")
    maxplot.axhline(cfg.intensity_threshold, color="green")

    maxplot_17.set_title("17 GHz, maximum value")
    maxplot_17.plot(pics_times_17, maxvals_17, "-o")

    print("making video")
    path = "/mnt/data/filepath/videos/{0}.mp4".format(pic_start)
    writer = FFMpegWriter(fps=3, extra_args=['-vcodec', 'libx264'])

    with writer.saving(fig, path, dpi=110):
        for a in range(0, len(pics)):
            pic = pics[a]
            pic_17 = pics_17[a]
            maxval = maxvals[a]
            maxval_17 = maxvals_17[a]

            tl1 = ccplot.axvline(pic["times"], color="red")
            tl2 = maxplot.axvline(pic["times"], color="red")
            tl3 = maxplot_17.axvline(pic_17["times"], color="red")
            tl4 = magnplot.axvline(pic["times"], color="red")

            picplot.set_title("34 GHz, maxval = " + str(maxval))
            picplot_17.set_title("17 GHz, maxval = " + str(maxval_17))

            picplot.imshow(pic["data"], interpolation=None, origin='low')
            picplot_17.imshow(pic_17["data"], cmap="winter", interpolation=None, origin='low')

            writer.grab_frame()

            tl1.remove()
            tl2.remove()
            tl3.remove()
            tl4.remove()

А вот один из скриншотов из видеоряда. Кому захочется пример полного видоса, могу потом тоже скинуть

Что дальше

Дальше требуется более детально проанализировать магнитограммы, починить некоторые ошибки фильтрации и сменить порог поиска на пониже. А что из этого выйдет, я буду писать тут.