Домашнее задание № 2

Метод наименьших квадратов

В файле ccd_v2.fits приведены данные эксперимента с приёмником изображений на основе прибора с зарядовой связью (детектор Sony ICX424). В указанном файле расположен четырёхмерный массив данных, содержащий 100 пар изображений, каждое размером 659x493 пикселей. Каждая пара представляет из себя изображения, полученные с детектора, равномерно засвеченного с различными интенсивностями (всего 100 вариантов, внутри пары интенсивность одинаковая). Первая по порядку пара — результат мгновенного считывания, и представляет из себя сигнал детектора в отсутствии света. Заметим, что во все остальные изображения аддитивно входит такой же сигнал, который часто называют напряжением смещения или подложкой.

Используя простые рассуждения, можно показать, что из таких данных возможно получить оценки шума считывания (read out noise) и коэффициента усиления (gain). Напомним, что каждый пиксель детектора накапливает заряд, производимый за счёт внутреннего фотоэлектрического эффекта. После оцифровки величина сигнала выражается в условных цифровых единицах (ADU), которые и записаны в файле, и связана с изначальным числом электронов следующим соотношением:

$$x = \frac{N_e}{g}$$

где $g$ — коэффициент усиления измеряемый в $[e^- / ADU]$, а дисперсия величины

$$ \sigma_x^2 = \frac{\sigma_r^2}{g^2}+\frac{\sigma_{N_e}^2}{g^2}=\frac{\sigma_r^2}{g^2} + \frac xg $$

где $\sigma_r$ — шум считывания, для удобства измеряется в количестве $e^-$, а последний переход — следствие того, что для количества фотоотсчетов действует закон распределения Пуассона, и $\sigma^2_{N_e}=N_e$.

Для того, чтобы учесть пространственные неоднородности детектора, на практике удобнее рассматривать поэлементную разность значений пикселей внутри пары кадров, с одинаковым уровнем освещения, тогда

$$2\frac{\sigma^2_r}{g^2}+\frac{2}{g}x = \sigma^2_{\Delta x}$$

где оценку $\sigma^2_{\Delta x}$ можно получить с помощью выборочной дисперсии по всем пикселям поэлементной разности кадров пары, а оценку $x$ можно получить с помощью выборочного среднего по всем пикселям двух кадров пары. При этом для первой пары (в отсутствии света) средний сигнал будет ненулевой, он соответствует напряжению смещения, и не обусловлен световым потоком; поэтому полученную для первой пары величину следует вычесть из всех полученных оценок. Подробности доступны в [1].

Нетрудно вычислить, что точность определения среднего потока $x$ будет на порядок лучше точности оценки $\sigma^2_{\Delta x}$, поэтому последняя стоит в правой части уравнений, а средний поток считается точно известным и стоит в левой части уравнений. Кроме того, можно показать, что дисперсия (точность) оценки $\sigma^2_{\Delta x}$ пропорциональна $4 x^2 + x$.

[1]: Steve B. Howell, "Handbook of CCD Astronomy", Cambridge University Press, 2006

Дедлайн 12 февраля 2026 в 23:55

Вы должны сделать следующее:

• В модуле lsp реализуйте функцию lstsq_ne(a, b) реализующую подход линейных наименьших квадратов с использованием нормальных уравнений. a — матрица задачи размера (N, M), b — вектор задачи размера (N,). Функция должна возвращать кортеж (x, cost, var), где x — решение задачи наименьших квадратов размера (M,), cost — квадрат нормы вектора невязок, var — матрица ковариации оценки параметров размера (M, M).

• В модуле lsp реализуйте функцию lstsq_svd(a, b, rcond=None) реализующую подход линейных наименьших квадратов с использованием вычислений на основ алгоритма сингулярного разложения. Результат работы функции и смысл параметров совпадают с определением функции lsqsq_ne. Параметр rcond принимает необязательный аргумент для регуляризации, если аргумент задан, то все сингулярные числа меньше чем rcond * s_max должны игнорироваться.

• В модуле lsp реализуйте функцию lstsq(a, b, method, **kwargs), которая принимает дополнительный параметр method с возможными значениями ne или svd и возвращает результат работы одной из двух соответствующих функций. С помощью аргумента kwargs должна быть реализована передача специфичных для каждого из алгоритмов параметров.

  Для проверки реализации функций бывает полезно использовать так называемые модульные тесты. В файле test.py приведено несколько примеров, куда можно добавлять свои варианты по образцу. Для запуска тестов можно использовать команду (предварительно нужно установить пакет pytest)

  > pytest test.py

• В файле ccd.py находится заготовка программы, принимающая путь к файлу с данными в качестве аргумента командной строки. Например, загрузите файл ccd_v2.fits:

  > python3 ccd.py ccd_v2.fits
  ccd_v2.fits

  Для задания аргументов командной строки в Spyder используйте меню "Запуск" -> "Настройки для файла..." (Ctrl+F6) -> "Опции командной строки"

  Нужно модифицировать код ccd.py таким образом, чтобы реализовать следующий функционал обработки:

  • загрузите файл с данными и вычислите оценки $\sigma^2_{\Delta x}$ от $x$;
  • точками изобразите на графике ccd.png зависимость $\sigma^2_{\Delta x}$ от $x$;
  • используя метод наименьших квадратов (функцию lstsq, реализованную в предыдущем упражнении), получите оценку двух коэффициентов линейной зависимости;
  • изобразите прямой линией на том же графике ccd.png модельную зависимость $\sigma^2_{\Delta x}$ от $x$;
  • используя полученные оценки коэффициентов линейной зависимости, получите оценки значений величин коэффициента усиления $g$ и шума считывания $\sigma_r$;
  • используя формулу оценки погрешности в косвенных измерениях, оцените ошибки определения коэффициента усиления $g$ и шума считывания $\sigma_r$;
  • используя технику bootstrap, проверьте, насколько хорошо в данном случае работает формула оценки погрешности в косвенных измерениях — для этого сгенерируйте 10000 подвыборок с повторами такого же размера, как исходные данные, для каждой подвыборки оцените коэффициент усиления $g$ и шум считывания $\sigma_r$, а затем рассчитайте их выборочные стандартные отклонения;
  • найдите оценку параметров этой же линейной модели с использованием метода взешенных наименьших квадратов и весами $\left(4x^2 + x + 1\right)^{-1/2}$, нанесите вторую прямую линию на график ccd.png, рассчитайте значения коэффициента усиления $g$ и шума считывания $\sigma_r$ и оцените погрешности их определения;
  • результат сохраните в файле ccd.json в следующем формате:

  [
    {
      "method": "least_squares",
      "ron": 12.34,
      "ron_err": 0.12,
      "gain": 1.234,
      "gain_err": 0.0012
    },
    {
      "method": "weighted_least_squares",
      "ron": 11.89,
      "ron_err": 0.07,
      "gain": 1.011,
      "gain_err": 0.002
    },
    {
      "method": "bootstrap",
      "ron": 12.34,
      "ron_err": 0.11,
      "gain": 1.234,
      "gain_err": 0.001
    }
  ]