Итак, когда в начале учебного года более-менее установилось учебное расписание, я связался со своим научным руководителем Сергеем из [Института солнечно-земной физики]( https://ru.wikipedia.org/wiki/%D0%98%D0%BD%D1%81%D1%82%D0%B8%D1%82%D1%83%D1%82_%D1%81%D0%BE%D0%BB%D0%BD%D0%B5%D1%87%D0%BD%D0%BE-%D0%B7%D0%B5%D0%BC%D0%BD%D0%BE%D0%B9_%D1%84%D0%B8%D0%B7%D0%B8%D0%BA%D0%B8 ), в дальнейшем просто ИСЗФ.
Здесь и началось всё веселье
В прошлом году я писал вместе с ним курсовую по обработке радиоизображений с помощью алгоритма [CLEAN]( https://en.wikipedia.org/wiki/CLEAN_(algorithm) ). Данный алгоритм широко используется для устранения шумов с сырых изображений, полученных с радиотелескопов и, в частности, с антенных решёток. Он, с одной стороны, совсем несложный, с другой - нужен для понимания того, как интерпретировать правильно картинки с радиотелескопов. В курсовой я попробовал с нуля запрограммировать алгоритм, поиграться с параметрами и обработать парочку тестовых картинок с [сибирского радиогелиографа]( https://ru.wikipedia.org/wiki/%D0%A1%D0%B8%D0%B1%D0%B8%D1%80%D1%81%D0%BA%D0%B8%D0%B9_%D1%81%D0%BE%D0%BB%D0%BD%D0%B5%D1%87%D0%BD%D1%8B%D0%B9_%D1%80%D0%B0%D0%B4%D0%B8%D0%BE%D1%82%D0%B5%D0%BB%D0%B5%D1%81%D0%BA%D0%BE%D0%BF ) СРГ-48.
Красивые картиночки и код с моей курсовой можно посмотреть на Гитхабе: https://github.com/vit1-irk/clean_lib
### Вернёмся к настоящему
В ИСЗФ мне понравилось, там достаточно интересно и прикольно (а ещё Сергей линуксоид, прогер, да и в целом с ним приятно было работать), поэтому я захотел большего и решил, что хочется принять участие в каком-нибудь более крутом проекте, уже связанным с реальной физикой. Но при этом включающем в себя немного программирования и обработки данных. В целом говоря, интересна идея изучать Солнце, и рано или поздно, хочу устроиться в ИСЗФ на работу.
Для меня нашлась работёнка в отыскании радиоисточников на Солнце, связанных с явлением гиромагнитного резонанса. Это излучение, порождённое электронами вне атомов, движущимися по замкнутым орбитам. На Солнце оно происходит в на короне и на границе короны с хромосферой, причём сигнал идёт сразу на нескольких гармониках. На частоте 34 ГГц подобные источники практически ни разу не находили, за исключением 2017 года, и целью было проверить, присутствуют ли ещё подобные образования на Солнце в другие года. Проверка не сильно сложная, однако, руки до сих пор ни у кого не дошли.
Использовать надо было данные с [Нобеямской радиообсерватории]( https://ru.wikipedia.org/wiki/%D0%9D%D0%BE%D0%B1%D0%B5%D1%8F%D0%BC%D1%81%D0%BA%D0%B0%D1%8F_%D1%80%D0%B0%D0%B4%D0%B8%D0%BE%D0%BE%D0%B1%D1%81%D0%B5%D1%80%D0%B2%D0%B0%D1%82%D0%BE%D1%80%D0%B8%D1%8F ) в Японии, где были данные аж с 1999 года по 2017. Более поздние тоже имеются, но они порченные, кривые и не подходят для обработки.
Данные Нобеямы включают в себя корреляционные кривые телескопа (показывают общее изменение яркости, которое потом отразится на изображении) и сами картинки, которые надо сначала скачать, а потом синтезировать с помощью специальной программки, которая написана на языке IDL (среда программирования, популярная у солнечников).
## Самый первый поиск источников
Начальная фильтрация - грубый поиск образований высокой яркостной температуры (выше 150 000K), при этом на корреляционной кривой сигнал не должен превышать шумовой порог. Изображения делались только раз в день, поэтому многие из источников, вероятно, были пропущены. Ищутся яркие точки, которые долго остаются на диске Солнца, по несколько часов.
https://ii-net.tk/ii/ii-point.php?q=/f/f/alicorn.blog/dGGy1m9rj59gqEfznqEe
Важно отличать обычные вспышечные события от гирорезонанса, потому что на корреляционных кривых первые дают резкий затухающий выброс, а вторые - нет.
Отфильтровалось несколько десятков дней, для каждого был построен видеоряд (см. в самый конец, где код).
После этого потребовалось выяснить конфигурацию магнитных полей на фотосфере Солнца в нужных активных областях, магнитограммы надо было доставать со спутников [Solar Dynamics Observatory]( https://en.wikipedia.org/wiki/Solar_Dynamics_Observatory ) и [Hinode]( https://en.wikipedia.org/wiki/Hinode_(satellite) ) .
## Магнитограммы
#### Магнитограмма SDO
https://ii-net.tk/ii/ii-point.php?q=/f/f/alicorn.blog/6lUXgHjXzef2rpEbD6eR
#### Увеличенная
https://ii-net.tk/ii/ii-point.php?q=/f/f/alicorn.blog/4QZiaLpBuCdh9qzQjR20
#### Магнитограмма Hinode
https://ii-net.tk/ii/ii-point.php?q=/f/f/alicorn.blog/RuUGVfkzn97VAPY2dtoB
## Фишки Питона
Пишу код и провожу вычисления в интерактивной среде разработки
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()
А вот один из скриншотов из видеоряда. Кому захочется пример полного видоса, могу потом тоже скинуть
https://ii-net.tk/ii/ii-point.php?q=/f/f/alicorn.blog/jPqH5A3gnFAsHzLaOuMN
## Что дальше
Дальше требуется более детально проанализировать магнитограммы, починить некоторые ошибки фильтрации и сменить порог поиска на пониже. А что из этого выйдет, я буду писать тут.
Ссылка в блоге: https://blog.alicorn.tk/posts/the-beginning-sun.html