RU2541906C1 - Способ создания многомерных градуировочных моделей аналитического прибора - Google Patents

Abstract

Изобретение относится к средствам создания градуировочных моделей измерительных приборов. Техническим результатом является повышение точности определения анализируемых свойств образца. В способе перед построением градуировочной модели проводят нормировку предобработанных первичных свойств на максимальное значение первичных свойств образцов градуировочного набора с последующим уменьшением размерности матрицы первичных свойств при помощи Фурье-преобразования, построение градуировочной модели осуществляют методами математического программирования с ограничением, после построения градуировочной модели проводят обратное Фурье-преобразование над оптимальной матрицей коэффициентов регрессии. 7 з.п. ф-лы, 8 ил., 5 табл.

Description

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

Важными требованиями при проведении анализа свойств тех или иных образцов являются такие требования, как минимизация времени проведения анализа, а также сохранение состава и свойств образцов в процессе измерений. Это особенно важно для промышленных предприятий, где стремятся проводить быстрый неразрушающий контроль свойств продукции на разных стадиях производства. Однако далеко не всегда можно провести измерения требуемых свойств продукции прямыми методами анализа с учетом указанных требований. В таких ситуациях приходится использовать косвенные методы измерения, где анализируемые свойства образцов определяют путем измерения других свойств образцов, зависящих от анализируемых свойств, которые, однако, в отличие от анализируемых свойств могут быть легко измерены напрямую в течение короткого промежутка времени и без разрушения образца. Например, ближняя инфракрасная (БИК) спектроскопия может рассматриваться как один из методов получения косвенной информации об объекте. На основании полученной информации проводят количественный анализ различных объектов, который включает в себя сочетание различных методов исследования многофакторных зависимостей.

Результатом измерения, выполненного с помощью, например, БИК спектрометра (аналитического прибора), является инфракрасный спектр. Под спектром понимают зависимость интенсивности поглощения (отражения или рассеивания) излучения от длины волны в заданном диапазоне. Спектры поглощения (отражения, рассеивания) молекул являются уникальными для каждого вещества, а интенсивность их связана с содержанием поглощающего компонента в исследуемом объекте. Т.о. величина интенсивности полосы поглощения (отражения, рассеивания) может быть связана градуировочной зависимостью с количественным содержанием интересующего вещества в анализируемой пробе. Однако в БИК-спектроскопии используемые количественные связи зачастую слишком сложны. Вызвано это тем, что спектры в БИК-области содержат множество широких перекрывающихся полос поглощения (отражения, рассеивания), которые накладываются на полосу, интересующего нас компонента. Это не позволяет построить простую градуировочную зависимость между интенсивностями характерной полосы с содержанием интересующего вещества. В таком случае формируют многомерную градуировочную модель, которая связывает весь массив спектральной информации с содержанием интересующего компонента.

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

Величины, измеряемые непосредственно прибором, называются «первичными свойствами», а анализируемые параметры образцов, от которых зависят результаты измерений, но которые не измеряются непосредственно, называют «вторичными свойствами». Для БИК-анализатора "вторичные" свойства образцов (например, количественное содержание протеина, влажности клейковины и т.д.) определяются по результатам измерений их спектральных характеристик, таких как оптические спектры поглощения, отражения или рассеяния, выступающие в качестве "первичных" свойств.

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

Простейшим методом, позволяющим находить оптимальные коэффициенты регрессии для построения многомерной градуировочной модели, является метод множественной линейной регрессии (МЛР) [3]. На фиг.2 представлена упрощенная схема метода МЛР. В данном методе набор переменных в матрице X объединяется в линейную комбинацию, коррелирующую настолько, насколько это возможно с соответствующими значениями матрицы Y. Таким образом, вычисляется регрессионная зависимость между матрицей Y и матрицей X. В блоке 100 представлено основное уравнение МЛР, в котором матрицы Y и X связаны регрессионными коэффициентами b. Задачей данного метода является нахождение такой матрицы $$\stackrel{\wedge }{B}$$

, чтобы величина ошибки f была минимальна. В методе МЛР она решается методом наименьших квадратов (блок 101).

Основным недостатком метода МЛР является необходимость обращения матрицы (X^T X) при вычислении $$\stackrel{\wedge }{B}$$

. Если в матрице X присутствуют коллинеарные данные, то ее обращение оказывается сложной задачей, а в некоторых случаях вообще неразрешимой проблемой. Умножение на (X^T X)^-1 становится неустойчивым и фактически эта операция соответствует делению на ноль. Недостатком метода МЛР является также то, что метод может работать только с ограниченным массивом матрицы первичных свойств X. Если размерность данной матрицы n×p (n - число образцов, p - число переменных в первичных свойствах каждого образца), то существующее ограничение для метода МЛР (p должно быть не больше n-1) существенно ограничивает его применение для градуировки, например, таких приборов, как сканирующие БИК анализаторы, где значение p может быть более 500.

Существенным расширением метода МЛР является метод регрессии на главные компоненты (РГК) [3].

Метод РГК работает в два основных шага.

На первом шаге для сокращения массива данных первичных свойств и исключения мультиколлинеарности их преобразуют, используя метод главных компонент (МГК) так, чтобы вся информация описывалась значительно меньшим числом переменных. МГК позволяет осуществить переход из исходной системы координат в систему координат меньшей размерности. Причем главные компоненты должна быть ортогональны. Уменьшая размерность исходной матрицы первичных свойств, происходит отделение полезной информации от шума. В общем виде разложение матрицы X размерностью n×p записывается в виде: X=TP^T+E,

где T - матрица счетов размерностью n×A, представляющая собой координаты исходных образцов на главных компонентах; P - матрица нагрузок размерностью A×P, представляющая собой набор коэффициентов, связывающих переменные в пространстве главных компонент с переменными в исходном пространстве; A - число главных компонент, выбранных для описания исходных данных. Вычисление матриц T и P осуществляется итерационными методами, например методом NIPALS [4].

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

На фиг.3 представлена упрощенная схема метода РГК. В блоке 110 представлено основное уравнение РГК, в котором матрицы Y и X связаны регрессионными коэффициентами B. Задачей метода РГК, так же как и ранее рассматриваемого метода МЛР, является нахождение такой матрицы $$\stackrel{\wedge }{B}$$

, чтобы величина ошибки F была минимальна. В блоке 111 представлено основное уравнение разложения матрицы первичных свойств X на счета T и нагрузки P. В блоке 117 представлена упрощенная схема МГК по вычислению счетов и нагрузок для каждой главной компоненты. Этапы блока 117 проводятся столько раз, сколько главных компонент определяется для конкретной модели. В блоке 112 в качестве начальных значений вектора счетов i первой компоненты выбирается столбец матрицы X, который имеет наибольшую сумму квадратов отклонений от среднего значения. В блоке 113 представлено решение для уравнения $$X=\stackrel{\wedge }{t}{p}^T+E$$

по нахождению методом МНК нагрузок $$\stackrel{\wedge }{p}$$

. В блоке 114 уточняется вектор счетов $$\stackrel{\wedge }{t}$$

путем решения методом МНК уравнения $$X=t{\stackrel{\wedge }{p}}^T+E$$

. В блоке 115 оценивается значимость расхождения $$\stackrel{\wedge }{t}$$

, определенных в блоках 112 и 114. Если расхождение значимо - то идет возврат в блок 113, в котором для вычисления $$\stackrel{\wedge }{p}$$

используются счета, полученные в блоке 114. Если расхождение не значимо - то повторяются все операции блока 117 заново для следующей главной компоненты, при этом в качестве матрицы первичных свойств берется матрица ошибок, определяемая как $$X-\stackrel{\wedge }{t}{\stackrel{\wedge }{p}}^T$$

. Когда вычислены счета и нагрузки для всех главных компонент - вычисляется матрица регрессионных коэффициентов $$\stackrel{\wedge }{B}$$

методом МНК (блок 116) для уравнения $$Y=\stackrel{\wedge }{T}{B}^T+E$$

.

РГК в отличие от МЛР снимает проблему коллениарности данных (первичных свойств). РГК позволяет работать с существенно большим количеством данных в отличие от МЛР.

Недостатком метода РГК является то, что декомпозиция матрицы X не обязательно приведет к выделению именно той структуры, которая коррелирует с искомыми свойствами Y. Заранее нет никакой уверенности в том, что количество используемых главных компонент будет содержать только ту информацию, которая коррелирует с интересующей нас переменной Y. В выбранных компонентах могут присутствовать и другие источники изменения первичных свойств, а в старших компонентах может остаться часть информации, которая коррелирует с переменной Y, т.е. со вторичными свойствами. Это обстоятельство может снижать точность выполняемых измерений. Также использование градуировочной модели, построенной с применением метода РГК, при выполнении непосредственного анализа неизвестного образца предполагает разложение вектора первичных свойств этого образца в пространстве главных компонент, что требует определенных вычислительных мощностей при рутинном анализе.

Широкое применение нашел метод проекции на латентные структуры (ПЛС), имеющий альтернативное название как метод дробных наименьших квадратов [3]. На фиг.4 представлена упрощенная схема метода ПЛС. В блоке 120 представлено основное уравнение ПЛС, в котором матрицы Y и X связаны регрессионными коэффициентами B. Задачей метода ПЛС, так же как и методов РГК и МЛР, является нахождение такой матрицы $$\stackrel{\wedge }{B}$$

, чтобы величина ошибки F была минимальна. Метод ПЛС, так же как и метод РГК, можно разделить на два шага. На первом шаге методом МГК проводится декомпозиция не только матрицы первичных свойств X на счета T и нагрузки P, но и матрицы вторичных свойств Y на счета T и нагрузки Q (блоки 121 и 122). Но это не два независимых друг от друга разложения, а взаимосвязанные процессы. Особенности вычисления счетов и нагрузок для обоих разложений схематично представлены в блоке 131. Этапы блока 131 проводятся столько раз, сколько главных компонент определяется для конкретной модели. Расчет оптимальных значений счетов и нагрузок для матриц X и Y осуществляется следующим образом. В блоке 123 в качестве начальных значений вектора счетов $$\stackrel{\wedge }{t}$$

первой компоненты выбирается столбец $$\stackrel{\wedge }{u}$$

матрицы вторичных свойств Y, который имеет наибольшую сумму квадратов отклонений от среднего значения. Данный вектор $$\stackrel{\wedge }{u}$$

используется для определения взвешенных нагрузок $$\stackrel{\wedge }{w}$$

матрицы X. Нахождение вектора $$\stackrel{\wedge }{w}$$

(блок 124) осуществляется методом МНК для уравнения $$X=\stackrel{\wedge }{u}{w}^T+E$$

. В блоке 125 представлено решение метода МНК по определению счетов $$\stackrel{\wedge }{t}$$

на основании уравнения $$X=t{\stackrel{\wedge }{w}}^T+E$$

. В блоке 126 вычисляются нагрузки $$\stackrel{\wedge }{p}$$

методом МНК из уравнения $$X=\stackrel{\wedge }{t}{p}^T+E$$

. Далее вычисляются нагрузки для матрицы вторичных свойств Y, увязывая их со счетами, определенными в блоке 125. В блоке 127 представлено решение методом МНК по нахождению нагрузок $$\stackrel{\wedge }{q}$$

из уравнения $$Y=\stackrel{\wedge }{t}{q}^T+E$$

. В блоке 128 уточняется вектор счетов $$\stackrel{\wedge }{u}$$

путем решения методом МНК уравнения $$Y=u{\stackrel{\wedge }{q}}^T+E$$

. В блоке 129 оценивается значимость расхождения $$\stackrel{\wedge }{u}$$

, определенных в блоках 123 и 128. Если расхождение значимо - то идет возврат в блок 124, в котором для вычисления $$\stackrel{\wedge }{w}$$

используются счета, полученные в блоке 128. Если расхождение не значимо - то повторяются все операции блока 131 заново для следующей главной компоненты при этом в качестве матрицы первичных свойств берется матрица ошибок, определяемая как $$X-\stackrel{\wedge }{t}{\stackrel{\wedge }{p}}^T$$

, а в качестве матрицы вторичных свойств берется матрица ошибок, определяемая как: $$Y-\stackrel{\wedge }{t}{\stackrel{\wedge }{p}}^T$$

. Когда вычислены счета и нагрузки для всех главных компонент - вычисляется матрица регрессионных коэффициентов $$\stackrel{\wedge }{B}$$

(блок 130) методом МНК.

Данный метод существенно отличается от методов МЛР и РГК. В ПЛС скрытые переменные связываются как с первичными, так и с вторичными свойствами. Причем эта связь осуществляется так, чтобы как можно большая часть дисперсии первичных свойств была максимально связана со вторичными свойствами. Это может повышать точность градуировочной модели, а также приводит к тому, что при построении модели используется меньшее число главных компонент, чем в РГК. Однако метод ПЛС, так же как и РГК, не идеален. При декомпозиции данных с помощью МГК отбрасывается часть информации, которая может быть полезна при построении модели. Метод ПЛС требует значительной квалификации исследователя, так как сложен в обращении.

Известны [5] также методы математического программирования (МП), которые могут быть использованы, как указано в [6], для поиска регрессионных коэффициентов, связывающих два массива данных, при минимизации значения ошибки регрессионного уравнения.

Известны решения задач квадратичного математического программирования [7], которые можно использовать для расчета градуировочной модели, если целевую функцию представить в виде:

$$g\left(b\right)={\left(Y-Xb\right)}^T\left(Y-Xb\right)\text{, (1)}$$

где g - вектор, содержащий квадраты ошибки градуировочной модели. Матрица Y описывает, как указывалось выше, вторичные свойства, а матрица X - первичные свойства. Вектор b содержит искомые коэффициенты регрессии. Область поиска экстремальных значений целевой функции определяется некими ограничениями. Коэффициенты регрессии b должны находиться внутри пространства гиперкуба. Гиперкуб - фигура в Евклидовом p-мерном пространстве (p - число переменных в первичных свойствах). Например, если p=3, то гиперкуб - это куб в трехмерном евклидовом пространстве, имеющий длину ребра, равную d. Оптимальные размеры гиперкуба могут определяться с помощью процедуры валидации. Т.е. размерность гиперкуба можно считать оптимальной, если критерии оценки качества градуировочной модели удовлетворяют необходимым условиям.

Целевую функцию g(b) необходимо привести к каноническому виду для последующего поиска коэффициентов регрессии. Для этого производится замена переменных в уравнении для целевой функции:

$$C=2\left(X^{T}X\right)\text{ (2)}$$

$$D=-2\left(X^{T}Y\right)\text{ (3)}$$

Тогда целевая функция g может быть представлена в следующем каноническом виде:

$$g\left(b\right)=D^{T}b+b^{T}Cb\to \min \text{ (4)}$$

$$d\ge b,b\ge 0\text{ (5)}$$

Существует много конечных методов решения данной задачи. Например, в качестве эффективного и простого метода решения может быть применен называемый метод решения задач о дополнительности, разработанный Лемке [7].

Упрощенная схема алгоритма построения градуировочной модели с использованием техники математического программирования представлена на фиг.5, в котором задача математического программирования сводится к задаче о дополнительности. В блоке 201 формируется целевая функция и ограничения из уже известного линейного уравнения градуировки, представленного в блоке 200, где вектор b (b₁, b₂, … b_n) содержит указанные выше искомые коэффициенты регрессии, вектор d (d₁, d₂, … d_n) представляет указанное выше пространство ограничений. Систему, сформированную в блоке 201, можно представить как каноническую задачу выпуклого квадратичного программирования (блок 202). Функция Лагранжа для данной задачи записана в блоке 203, где u и ν - множители Лагранжа. Для данной функции Лагранжа, дифференциальные условия Куна-Такера записаны в блоке 204. Произведя замену S=d-b, систему можно представить так, как показано в блоке 205. Представив систему, записанную в блоке 204, в виде набора матриц (блок 205), можно записать эквивалентную задачу о дополнительности (блок 206). Алгоритм решения задачи о дополнительности представлен на фиг.6.

Задача о дополнительности (блок 300) имеет решение, если хотя бы один из элементов q_i<0. Чтобы начать работу алгоритма с начала (w=q, z=0), для получения неотрицательного решения производится преобразование задачи, представленное в блоке 301, где во все уравнения системы вводится искусственная переменная z₀ (при этом e - единичный вектор). Полученную систему можно представить в виде симплекс таблицы (блок 302):

базис	z1	…	zs	…	zk	z0	q
w1	-m11	…	-m1s	…	-m1k	-1	q1
ws	-ms1	…	-mss	…	-msk	-1	qs
wn	-mn1	…	-mns	…	-mnk	-1	qn

В базис вводится z₀, которая заменит базисную переменную. В качестве заменяемой базисной переменной выбирается переменная с наибольшим по абсолютной величине отрицательным значением. Т.е. если $$\underset{i}{q_s=\min q_i<0}$$

, то в этом случае z₀ заменит базисную переменную w_s. Все остальные значения в таблице пересчитываются в блоке 305 относительно выбранной строки следующим образом:

$$q_s^{\text{'}}=-q_s,q_i^{\text{'}}=q_i-q_s,i\ne s\text{ (6)}$$

$$m_{sj}^{\text{'}}=m_{sj},j=1\dots k\text{ (7)}$$

$$m_{ij}^{\text{'}}=-m_{ij}+m_{sj},j=1\dots k,i\ne s\text{ (8)}$$

Таким образом, получается новая симплекс таблица:

базис	z1	…	zs	…	zk	z0	q
w1	$$m_{11}^{\text{'}}$$	…	$$m_{1s}^{\text{'}}$$	…	$$m_{1k}^{\text{'}}$$	0	$$q_1^{\text{'}}$$
z0	$$m_{s1}^{\text{'}}$$	…	$$m_{ss}^{\text{'}}$$	…	$$m_{sk}^{\text{'}}$$	1	$$q_s^{\text{'}}$$
wn	$$m_{n1}^{\text{'}}$$	…	$$m_{ns}^{\text{'}}$$	…	$$m_{nk}^{\text{'}}$$	0	$$q_n^{\text{'}}$$

Далее, в блоке 303 определяется, какая переменная будет выведена из базиса и какая переменная будет введена в базис вместо нее. В базис всегда вводится переменная, дополнительная к базисной переменной, выведенной из базиса на предыдущем шаге. Т.е. если на предыдущем шаге была выведена переменная w_s, то на данном шаге в базис будет вводиться переменная z_s. Переменные w_i и z_i - это пара взаимодополняемых переменных. При определении переменной, которая будет выводиться из базиса, применяется правило минимального соотношения:

$$\underset{m_{is}^{\text{'}}>0}{\min }\left(q_i^{\text{'}}/m_{is}^{\text{'}}\right)=q_k^{\text{'}}/m_{ks}^{\text{'}}\text{ (9)}$$

для всех i=1…n, для которых значения $$m_{is}^i>0$$

.

То есть в блоке 303 в базис вводиться переменная z_s вместо переменной w_k, выводимой из базиса. Далее в блоке 304 проверяется условие, является ли выводимая переменная переменной z₀.

Как только в строке S будет выведена переменная z₀, значит в базисе остались только z_i, где i=1…n, и оптимальное решение достигнуто, z_i=q_i, где i=1…n. Искомый вектор регрессионных коэффициентов (блок 306) есть b_i=z_i, где i=1…n.

Следует отметить, что метод математического программирования не использует нестабильных операций, поэтому коллинеарность анализируемых данных не влияет на полученные результаты в отличии от метода МЛР.

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

Однако применение метода математического программирования, так же как алгоритм МЛР, не может работать с данными, для которых число образцов меньше числа переменных в первичных свойствах, что существенно ограничивает его применение для градуировки, например, таких приборов, как сканирующие БИК анализаторы, где значение p может быть более 500.

Авторы патента [8] показали, что техника математического программирования может эффективно применяться в спектрометрии, при обработке многомерных данных. В патенте представлено устройство для определения наиболее вероятных концентраций компонентов в анализируемом образце. Устройство включает в себя средства для измерения спектральной интенсивности указанной химической смеси, а также средства для вычисления указанных концентраций, средства для хранения данных, для отображения данных и т.д. Для определения концентраций устройство обрабатывает многомерные наборы данных. При этом осуществляется нахождение таких математических соотношений между первичными и вторичными свойствами, чтобы величина ошибки была минимальна. В патенте сказано, что наиболее предпочтительным способом нахождения указанных математических отношений является симплекс метод линейного программирования, так как использование данного метода значительно сокращает время вычисления.

Авторы [9] также предложили использовать метод математического программирования в способе построения многомерной градуировки для БИК анализаторов, применяемых для анализа различных свойств образцов сельскохозяйственной промышленности.

Однако использование метода математического программирования в случаях, предложенных авторами [8] и [9], ограничивает его применение для градуировочных моделей, для которых число образцов меньше числа переменных в их первичных свойствах.

За прототип предлагаемого способа создания многомерных градуировочных моделей аналитического прибора принят способ, изложенный в [1]. Данный способ состоит из поэтапной процедуры расчета градуировочной модели, представлен на фиг 1.

На первом этапе проводится подбор образцов градуировочного и валидационного (проверочного) наборов образцов с известными вторичными свойствами, а также регистрация их первичных свойств (блоки 010 и 011). Вторичные свойства градуировочного и валидационных наборов определяются при помощи независимых референтных методов, имеющих заданную точность и воспроизводимость. Значения вторичных свойств градуировочного и валидационного наборов называются референтными данными. Диапазон изменения вторичных свойств образцов градуировочного и валидационных наборов должен полностью охватывать диапазон возможных изменений этих свойств при последующем анализе неизвестных образцов. Кроме того, для создания наиболее устойчивой модели образцы в градуировочном наборе должны быть однородно распределены в диапазоне изменения анализируемых свойств. Подобные принципы выбора образцов градуировочного и валидационного набора регламентируется, например, в [1]. Таким образом, на первом этапе формируются:

матрица X первичных свойств и матрица Y вторичных свойств (матрица референтных данных) градуировочного набора образцов;

матрица X' первичных свойств и матрица Y' вторичных свойств (матрица референтных данных) валидационного набора образцов.

Валидационный набор обычно отличается от градуировочного меньшим количеством образцов.

На втором этапе для обоих наборов могут быть проведены процедуры предобработки данных. Предобработка представляет собой определенные алгоритмы математического преобразования матрицы X для уменьшения влияния на них различных факторов, не связанных с определяемыми свойствами. Известны, например, следующие виды преобразования данных [2]: выравнивание базовой линии, сглаживание данных (например, с применением фильтра Савицкого-Голея, скользящего среднего и т.п.), взятие производных, центрирование данных, шкалирование данных (нормирование на среднее квадратическое отклонение по каждой переменной), нормирование на среднее квадратическое отклонение по образцу, мультипликативная коррекция и др. Некоторые виды предобработок, например центрирование данных, может применяться, как к матрице X, так и к матрицам X и Y одновременно. Над данными, полученными на этапе 1, может применяться как один определенный алгоритм предобработки, так и несколько алгоритмов в определенной последовательности. Проверка оптимальности выбранных алгоритмов предобработки и последовательности их применения определяется на стадии валидации (этап 4).

Следует отметить, что на градуировочный и валидационный набор накладываются одинаковые алгоритмы предобработок и последовательность их применения (блоки 012 и 013).

На третьем этапе проводится расчет градуировочной модели (блок 014) с использованием предобработанных данных градуировочного набора, заключающийся в нахождении оптимальных соотношений (коэффициентов регрессии) между вторичными (матрица Y) и первичными свойствами (матрица X) градуировочного набора образцов. Этап построения градуировочной модели является наиболее трудоемкой и значимой процедурой в многомерной градуировке и требует от исследователя достаточного опыта и знаний. Для нахождения оптимальных соотношений (коэффициентов регрессии) между вторичными и первичными свойствами образцов применяются различные методы многомерного анализа.

На четвертом этапе проводится проверка созданной модели с применением валидационного набора образцов (блок 015). На данном этапе также могут проводиться и другие виды проверок с использованием, например, перекрестной проверки или проверки на выбросы по оценке расстояния Махаланобиса [2]. По результатам проверки может корректироваться градуировочный набор (отбрасываться аномальные образцы), изменяться предобработки данных, изменяться параметры градуировочной модели.

В качестве методов вычисления регрессионных коэффициентов в [1] предлагается использовать методы МЛР, РГК и ПЛС. Данные методы были описаны выше и, как было показано, имеют определенные недостатки, которые могут влиять на качество градуировочной модели и соответственно точность анализа.

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

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

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

Сущность предлагаемого способа заключается в иной процедуре вычисления регрессионных коэффициентов, связывающих первичные и вторичные свойства градуировочных образцов, которая позволяет ликвидировать недостатки применения методов МЛР, РГК, ПЛС, методов математического программирования, а также повысить точность выполнения анализа по созданным, таким образом, градуировочным моделям.

Предлагаемая процедура вычисления регрессионных коэффициентов работает в четыре последовательных шага. На первом шаге производиться нормировка первичных свойств на модуль максимального значения первичных свойств образцов градуировочного набора. При этом указанной нормировке подвергаются уже предобработанные на этапе 2 общей процедуры градуировки (согласно фиг.1) первичные свойства градуировочных и валидационных образцов. На втором шаге проводится снижение размерности данных первичных свойств с помощью преобразования Фурье. На третьем шаге производится расчет регрессионных коэффициентов с помощью метода квадратичного программирования. На четвертом шаге осуществляется обратное преобразование Фурье над полученными коэффициентами регрессии.

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

Второй шаг - уменьшение размерности матрицы первичных свойств путем применения гармонического преобразования (преобразования Фурье). Это позволяет, с одной стороны, выделить полезную информацию из первичных свойств, отделяя шумовую составляющую, что увеличивает точность получаемых результатов, с другой стороны, позволяет снять ограничение на количество переменных в первичных свойствах (например, количество длин волн в спектральных данных), которое действует при использовании методов математического программирования: число образцов градуировочного набора должно быть не меньше числа переменных первичных свойств. Поскольку следующим шагом является использование метода математического программирования, который не предусматривает не корректных математических операций при коллинеарности данных, Фурье преобразование может выступать как процедура снижения размерности альтернативой методу МГК. Он более прост и требует меньше времени на выполнение соответствующих преобразований. Если f(λ) представляет собой спектр, измеренный при λ длинах волн, где λ=1…N, то эта функция может быть аппроксимирована рядом Фурье:

$$f(\lambda )=\frac{a_0}{2}+{\displaystyle \sum _{k=1}^{M-1}{a}_k}\cos (\frac{\pi k\lambda }{N})+{\displaystyle \sum _{k=1}^{M-1}{b}_k}\sin (\frac{\pi k\lambda }{N}),\text{ (10)}$$

где a _k - коэффициенты ряда Фурье и k=1…M, которые рассчитываются с помощью преобразования Фурье. Таким образом, как уже отмечалось выше, исходный спектр может быть записан в виде пар коэффициентов a _k и b_k. Число коэффициентов M должно быть достаточно большим, чтобы спектр описывался довольно точно, но не настолько чтобы при этом воспроизводился шум, представляющий в исходном спектре относительно высокочастотную составляющую. Оптимальное число коэффициентов определяется с помощью процедуры валидации, т.е. число коэффициентов можно считать оптимальным, если критерии оценки качества градуировочной модели удовлетворяют необходимым условиям. Таким образом, операция на втором шаге позволяет анализировать входные данные, исключая шумовую составляющую, без перехода в новое пространство переменных. Анализ данных в исходном пространстве позволяет интерпретировать данные на любом шаге расчета, так как не надо возвращаться к исходному пространству. Это позволяет снизить количество полезных данных, отброшенных вместе с шумовой составляющей, а также сокращает количество вычислительных операций и повышает качество моделирования. Например, при использовании инфракрасной спектроскопии в исследовании количественного состава зерновых культур применение Фурье преобразования в качестве метода обработки данных позволяет значительно сократить объем исходных данных, не снижая при этом информативности.

Третий шаг - нахождение регрессионных коэффициентов путем решения канонической задачи квадратичного математического программирования для целевой функции, представленной в виде формулы (1), где g - вектор, содержащий квадраты ошибок градуировочной модели. Матрица Y описывает, как указывалось выше, вторичные свойства, а матрица X - первичные свойства. Вектор b содержит искомые коэффициенты регрессии. Область поиска экстремальных значений целевой функции определяется размером гиперкуба, который определяется с помощью процедуры валидации. Решение канонической задачи сводится к решению задачи о дополнительности, например, методом Лемке.

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

Более подробно шаги вычисления регрессионных коэффициентов представлены на фиг.7.

В блоке 400 представлено основное уравнение, в котором матрицы Y и X связаны регрессионными коэффициентами B. Задачей в предлагаемом способе расчета многомерных градуировочных моделей, так же как и ранее рассматриваемых методах, является нахождение такой матрицы $$\stackrel{\wedge }{B}$$

, чтобы величина ошибки F была минимальна. В блоке 401 выполняется операция, описанная для первого шага: нормировка данных первичных свойств на модуль их максимального значения. В блоках 402 и 403 выполняется операция, описанная для второго шага. В блоке 402 формируется матрица Фурье-преобразования, согласно уравнению (10), где $${\nu }_{NM}=\cos (\frac{\pi M\lambda }{N})$$

, $$z_{NM}=\sin (\frac{\pi M\lambda }{N})$$

, N - число спектров. В блоке 403 выполняется преобразование Фурье над матрицей первичных свойств, представляющее собой произведение матрицы полученных коэффициентов Фурье и матрицы первичных свойств. Операции, описанные для третьего шага, выполняются в блоках 404 и 405. В блоке 404 задача математического программирования сводится к задаче о дополнительности с помощью описанного выше алгоритма, представленного на фиг.6. В блоке 405 вычисляются искомые регрессионные коэффициенты с помощью конечных методов решения задачи о дополнительности, например, описанного выше метода Лемке, представленного на фиг.7. Операция заключительного четвертого шага представлена в блоке 406. На данном этапе производится описанная выше операция обратного преобразования Фурье над вычисленными регрессионными коэффициентами.

На фиг.8 представлена поэтапная процедура предлагаемого способа создания градуировочных моделей.

Действия на первом этапе не отличаются от действий первого этапа, представленного в прототипе и изображенного на фиг.1. Аналогичным образом проводится подбор образцов градуировочного и валидационного наборов образцов с известными вторичными свойствами, а также регистрация их первичных свойств (блоки 500 и 501).

На втором этапе для обоих наборов могут быть проведены процедуры предобработки данных. Виды предобработок данных могут применяться аналогичные тем, что и в прототипе (фиг.1). Так же как и в прототипе, на градуировочный и валидационный набор накладываются одинаковые процедуры предобработок и последовательность их применения (блоки 502 и 503). Проверка оптимальности выбранных алгоритмов предобработки и последовательности их применения определяется на стадии валидации (этап 4). Существенной особенностью предлагаемого способа на данном этапе является то, что в перечне применяемых предобработок не должна присутствовать операция нормировки первичных свойств на модуль максимального значения первичных свойств образцов градуировочного набора.

Принципиальные отличия в предлагаемом способе от прототипа определяются в третьем этапе. Для градуировочного и валидационного набора выполняется процедура нормировки первичных свойств на модуль максимального значения первичных свойств образцов градуировочного набора (блоки 504 и 505). Далее, градуировочный набор подвергается процедуре гармонического преобразования с целью уменьшения размерности матрицы первичных свойств (блок 506). Проверка оптимальности выбранного количества коэффициентов гармонического преобразования определяется на стадии валидации (этап 4). Далее вычисляются регрессионные коэффициенты (блок 507) методами математического программирования с ограничениями, описанными выше. Оптимальность параметров ограничений (например, размер гиперкуба) определяется на стадии валидации (этап 4). Следующий шаг (блок 508) - обратное преобразование Фурье регрессионных коэффициентов, чтобы привести их к виду, который бы позволил использовать градуировочную модель на стадии предсказания неизвестных образцов без применения гармонического преобразования к первичным свойствам анализируемых образцов.

На заключительном четвертом этапе проводится проверка созданной модели с применением валидационного набора образцов (блок 509). На данном этапе также могут проводиться и другие виды проверок с использованием, например, перекрестной проверки или проверки на выбросы по оценке расстояния Махаланобиса. По результатам проверки может корректироваться градуировочный набор (отбрасываться аномальные образцы), изменяться предобработки данных, изменяться параметры градуировочной модели. При положительном результате проверки процедура создания градуировочной модели завершается (блок 510).

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

Рассмотрим применение заявляемого способа на примере спектрометров для анализа химического состава образца. Для проверки эффективности предложенного метода использовался прибор ИнфраЛЮМ ФТ-40 компании «Люмэкс» и наборы различных образцов для градуировки анализатора, а также наборы образцов для валидации результатов. Для того чтобы максимально оценить эффективность предложенного способа, проведено сравнение со способом, изложенным в прототипе, в котором применялись классические методы построения многомерной градуировки РГК и ПЛС, а также метод математического программирования (МП) без особенностей, предложенных в предлагаемом способе.

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

$$SEV=\sqrt{\frac{{\displaystyle \sum _{i=1}^n{({\stackrel{⌢}{y}}_i-y_i)}^2}}{d_{\nu }}}\text{ , (11)}$$

где d_ν - общее число образцов дополнительного набора, y_i - референтные значения анализируемого свойства для i-го образца дополнительного набора, $${\stackrel{\wedge }{y}}_i$$

- предсказанные значения анализируемого свойства для i-го образца дополнительного набора.

Расчет и проверка градуировочных моделей проводились с использованием одних и тех же градуировочных и валидационных наборов образцов.

Результаты испытаний представлены ниже в таблицах 1-5.

Таблица 1
Результаты предсказания для образцов валидационного набора по градуировочным моделям, построенным различными способами, для анализа протеина в шроте подсолнечном
Образец	Референтное значение	Прототип			Предлагаемый способ
		РГК	ПЛС	МП
shp-10219	35,57	35,21	35,42	35,62	35,88
shp-10273	37,70	38,46	38,53	38,06	38,05
shp-10229	38,65	39,30	39,05	38,85	39,27
shp-10287	39,10	38,98	39,37	38,97	39,37
shp-10121	39,85	39,75	38,86	39,29	39,42
shp-10127	40,60	41,22	41,28	41,49	41,05
shp-10139	41,68	41,41	41,23	42,06	41,67
shp-10145	42,58	43,31	43,67	42,81	42,63
shp-10169	43,92	42,89	42,15	42,93	42,73
shp-10157	44,37	43,02	42,40	43,72	43,23
Значение SEV		1,03	1,26	1,04	0,97

Таблица 2
Результаты предсказания для образцов валидационного набора по градуировочным моделям, построенным различными способами, для анализа фосфора в подсолнечном масле
Образец	Референтное значение	Прототип			Предлагаемый способ
		РГК	ПЛС	МП
sun_oil-10031	0,05	0,13	0,07	0,08	0,09
sun_oil-10013	0,20	0,32	0,14	0,43	0,31
sun_oil-10053	0,41	0,46	0,44	0,35	0,47
sun_oil-10137	0,58	0,54	0,62	0,51	0,53
sun_oil-10035	0,69	0,64	0,72	0,74	0,69
sun_oil-10120	0,86	0,82	0,81	1,02	0,84
sun_oil-10059	1,87	1,79	1,84	1,78	1,84
sun_oil-10058	1,96	1,88	1,94	1,78	1,93
Значение SEV		0,08	0,06	0,15	0,05

Таблица 3
Результаты предсказания для образцов валидационного набора по градуировочным моделям, построенным различными способами, для анализа влажности в пшеничном шроте
Образец	Референтное значение	Прототип			Предлагаемый способ
		РГК	ПЛС	МП
dps-1065	12,00	12,00	12,01	12,26	11,96
dps-1067	12,20	12,27	12,33	12,36	12,25
dps-1068	13,00	13,09	13,12	13,31	13,19
dps-1004	13,60	13,64	13,71	14,04	13,75
dps-1023	14,40	14,53	14,27	14,37	14,37
dps-1029	15,00	15,27	15,32	14,32	15,00
dps-1050	15,40	15,47	15,41	15,82	15,30
dps-1026	16,00	15,82	15,75	15,62	16,13
dps-1002	16,40	16,29	16,41	16,01	16,37
dps-1042	17,00	16,65	16,61	16,49	16,81
Значение SEV		0,17	0,18	0,49	0,13

Таблица 4
Результаты предсказания для образцов валидационного набора по градуировочным моделям, построенным различными способами, для анализа зольности в шроте пшеничном
Образец	Референтное значение	Прототип			Предлагаемый способ
		РГК	ПЛС	МП
dps-1008	1,44	1,57	1,55	1,57	1,48
dps-1051	1,51	1,57	1,57	1,57	1,54
dps-1026	1,57	1,59	1,59	1,59	1,60
dps-1023	1,61	1,61	1,61	1,61	1,61
dps-1019	1,64	1,69	1,70	1,69	1,69
dps-1009	1,66	1,60	1,61	1,58	1,58
dps-1063	1,67	1,67	1,68	1,67	1,63
dps-1024	1,72	1,68	1,71	1,68	1,72
Значение SEV		0,07	0,06	0,09	0,04

Таблица 5
Результаты предсказания для образцов валидационного набора по градуировочным моделям, построенным различными способами, для анализа содержания белка в дрожжах
Образец	Референтное значение	Прототип			Предлагаемый способ
		РГК	ПЛС	МП
drk-2529	40,26	41,42	41,19	42,13	39,83
drk-2523	42,24	42,49	42,76	43,68	42,86
drk-2532	42,60	42,42	42,48	43,39	42,24
drk-2530	42,88	43,51	43,12	44,17	43,34
drk-2521	43,22	43,19	43,08	43,97	43,82
drk-2528	44,04	42,50	42,78	43,77	43,34
drk-2527	44,64	44,20	43,92	45,69	44,29
Значение SEV		0,87	0,69	1,22	0,52

По результатам, представленным в таблицах 1-5, можно сделать вывод о том, что предлагаемый способ создания градуировочных моделей обеспечивает лучшую точность предсказаний анализируемых свойств ряда продуктов, чем способ, изложенный в прототипе с применением наиболее классических методов нахождения регрессионных коэффициентов, таких как РГК, ПЛС и МП.

Литература

1. ASTM standard, E 1655-00, Practices for Infrared Multivariate Quantitative Analysis.

2. Muhammad A. Sharaf Chemometrics: John Wiley & Sons pp.90-97.

3. Tormod Naes, Tomas Isaksson, Tom Fearn, Tony Davies Multivariate Calibration and Classification.: NIR Publications, 2002. pp 19-55.

4. Kim H. Esbensen Multivariate Data Analysis: Camo ASA, pp.19-35.

5. Masanao Aoki Introduction to optimization techniques.: MacMillan & Co LTD, 1971. 341 p.

6. Harvey M. Wagner Linear Programming Techniques for Regression Analysis // Journal of the American Statistical Association, vol.54 No.285 (Mar., 1959) pp.206-212.

7. Charnes A. Cooper W.W. Ferguson R.O. Optimal Estimation of Executive Compensation by Linear Programming // Management Science, vol.1 No.2 (Jan., 1955), pp.138-151.

8. Патент США №4744028, МПК G06F 1520, опубликована 04.19.1985.

9. Ilya Skutin, Konstantin Zharinov and Maria Sushkova Mathematical Programming Method as a Means of Calibration of NIR Analyzers Journal of Chemistry and Chemical Engineering, vol.6 No 8 (August 2012) pp.692-697.

Claims (8)

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

2. Способ по п.1, отличающийся тем, что нахождение оптимальной матрицы коэффициентов регрессии осуществляют с использованием методов как для линейного, так и нелинейного математического программирования.

3. Способ по п.1, отличающийся тем, что при применении метода математического программирования с ограничением оптимальную матрицу ищут в гиперкубе заданного размера

4. Способ по п.1, отличающийся тем, что при применении метода математического программирования с ограничением используют симплекс таблицы.

5. Способ по п.1, отличающийся тем, что оптимальные параметры Фурье преобразования определяют по процедуре валидации.

6. Способ по п.3, отличающийся тем, что оптимальный размер гиперкуба определяют по процедуре валидации.

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

8. Способ по п.7, отличающийся тем, что используют спектрометры на основе принципа Фурье-спектроскопии.

Images