Размещено на http://allbest.ru

Алгоритмическая и программная реализация прогноза почвенных параметров

Введение

линейный компонента факторный программный

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

Эта связано с тем, что значения почвенных параметров получают на сегодняшний день через большие временные промежутки: 20,30 и т.д. лет.

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

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

Здесь очень важно определение параметров объектов исследования по их спектральным характеристикам.

Это связано с периодическим характером изменения многих почвенных параметров.

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

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

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

Глава I. Спектральный метод

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

Детально описывается алгоритм быстрого преобразования Фурье (БПФ), который аппаратно реализован на ЭВМ высокого класса [1].

1.1 Быстрое преобразование Фурье

Пусть периодический сигнал представлен функцией а(t), кусочной или непрерывной на отрезке [-Т/2, Т/2] с периодом Т.

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

а(t)=a₀/2+с(a_k*cos(2bkt/T)+b_k*sin(2bkt/T)),(1.1.1)

где

a_k=(2/T)*dа(t)*cos(2bkt/T)*dt, (k=0,1,2,3,...)(1.1.2)

b_k=(2/T)dа(t)*sin(2bkt/T)dt, (k=1,2,3,...)

Амплитудный спектр

A_k=fa²_k+b²_k, (k=0,1,2,3...)(1.1.3)

Фазовый спектр сигнала:

t_k=-arctg(b_k/a_k), (k=1,2,...)(1.1.4)

В случае, если а(t) четная функция, разложение будет иметь вид:

а(t)=a₀/2+с(a_k*cos(2bkt/T),(1.1.5)

a_k=(4/T)*dа(t)*cos(2bkt/T)*dt, (k=0,1,2,3,...)(1.1.6)

Для нечетной функции а(t)

а(t)=с(b_k*sin(2bkt/T),(1.1.7)

b_k=(4/T)dа(t)*sin(2bkt/T)dt, (k=1,2,3,...)(1.1.8)

Комплексная форма ряда Фурье имеет вид:

а(t)=с(A_k/2)*EXP[-i2bkt/T],(1.1.9)

где

A_k=2/Тdа(t)*EXP[-i2bkt/T]dt,(1.1.10)

где

A_k=a_k-ib_k, (К = 0, 1, 2,...)(1.1.11)

A_-k=a_k+ib_k

Если а(t) апериодический сигнал, то его изображение в частотной области будет:

F(iz)= dа(t)*EXP[-izt]dt(1.1.12)

Это, как известно, прямое преобразование Фурье.

Обратное преобразование:

а(t)=1/2bd*F[iz]*e^iztdt.(1.1.13)

Условие осуществимости прямого преобразования -- существование интеграла:

d|а(t)|dt(1.1.14)

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

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

Если предположить, что реализация Х(t) обладает периодичностью, и период ее равен Т_р, а основная частота а1=1/Т_р, то реализация может быть представлена рядом Фурье:

X(t)=a₀/2+с(a_q*cos2bqа_?t+b_q*sin2bqа_?t),(1.1.15)

a_q=(2/T_р)*dX(t)*cos2bqа_?t*dt, q=0,1,2,

b_q=(2/T_р)*dX(t)*sin2bqа_?t*dt, q=1,2,

Пусть реализация Х(t) имеет конечную длину Т_r=Т_р, равную основному периоду.

Предположим, что она состоит из четного числа N эквидистантных наблюдений с интервалом дискретности h, который выбран таким образом, что частота среза а=1/2h достаточно высока.

Будем считать, что нулевая ордината реализации равна нулю, и обозначим преобразованную последовательность в виде:

Х_n=Х(nh), n=1,2,... N.(1.1.16)

Вычислим теперь по всем N значениям реализации конечный ряд Фурье.

Для любой точки t, принадлежащей интервалу (0,Т_р),этот ряд имеет вид:

N/2 N/2-1

X(t)=A₀+сA_q*cos(2bqt/t_p)+сB_q*sin(2bqt/t_p),(1.1.17)

q=1 q=1

В точках t=n*h, n=1,2,..,N, где Т_p=N*h будем иметь:

X_n=X(nh)=A₀+сA_q*cos(2bqn/N)+сB_q*sin(2bqn/N),(1.1.18)

Коэффициенты А_q и В_q определяются выражениями:

А_о=1/N сX_n=0

А_q=2/N сX_n*cos(2bqn/N), q=1,2,...,N/2-1

А_N/2=1/N сX_n*cos(bn),

B_q=2/NсX_n*sin(2bqn/N), q=1,2,...,N/2-1(1.1.19)

Программа для расчета величин А_q и B_q должна содержать следующие пункты:

1) определение x=(2pqn)/N при фиксированных значениях q, n;

2) вычисление cosx, sinx;

3) вычисление Х_n*cosx, Х_n*sinx;

4) вычисление суммы для каждого из этих выражений при n=1,2,...,N;

5) приращение аргумента q на единицу и повторение всех перечисленных действий.

Такой способ требует выполнения примерно N² операций умножения и сложения действительных чисел.

Чтобы существенно снизить затраты машинного времени, используется другой способ расчета, получивший название быстрого преобразования Фурье (БПФ).

Преобразование Фурье действительной или комплексной функции Х(t), заданной на бесконечном интервале, представляет собой комплексную величину:

X(а)=dX(t)*e^-2biаtdt.(1.1.20)

Ограничив интервал задания функции Х(t), приняв его равным [0,Т], можно построить конечное преобразование Фурье:

X(а,T)=dX(t)*e^-2biаtdt(1.1.21)

Предположим, что функция Х(t) представлена N эквидистантными наблюдениями с интервалом дискретности h, который выбран таким образом, что частота среза достаточна высока.

Поскольку t₀=0, отсчеты t_n=n*h

Тогда: X_n=X(nh) n=0,1,2,...,N--1.(1.1.22)

Дискретная аппроксимация (1.1.21) при произвольном значении а есть:

X(а,T)=hсX_nEXP[-2biаnh]

Для расчета функции X(а,T) обычно выбираются дискретные значения частоты:

а_k=kа₀=k/(N*h) (k=0,1,2,...,N-1)(1.1.24)

Преобразованная последовательность дает на этих частотах составляющие Фурье:

X_k=X(а_k,T)/h=сX_nEXP[-i(2bkn/N)](1.1.25)

k=0,1,2,...,N-1

Преобразование однозначно только до значений k=N/2, поскольку этой точке соответствует частота Найквиста.

Быстрое преобразование Фурье применяется для вычисления последовательности Х_к, но может быть также использовано для нахождения коэффициентов Фурье А_q и B_q.

Положим:

W(U)=EXP[-i(2bU/N)](1.1.26)

W(N)=1, при всех U и V справедливо равенство:

W(U+V)=W(U)*W(V).(1.1.27)

Положим:

X(k)=X_k, X(n)=X_n (1.1.28)

Тогда:

X_k=сX(n)*W(kn), k=0,1,2,...,N-1(1.1.29)

Для расчета всех значений Х_к необходимо выполнить примерно N² операций умножения и сложения комплексных чисел (одна такая комплексная операция эквивалентна четырем операциям умножения и сложения действительных чисел).

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

Если N может быть представлено в виде произведения целых и больших единицы чисел:

N=mr_i=r₁*r₂*...*r_p,(1.1.30)

то последовательность Х_к может быть найдена итеративно путем расчета суммы p слагаемых. (N/r₁) преобразований Фурье, каждое из которых требует 4r²₁ операций с действительными числами;

Таким образом, общее число операций.

4(N*r₁+N*r₂+N*r_p)=4Nсr_i(1.1.32)

В результате при использовании БПФ, а не стандартного метода получаем коэффициент ускорения вычислений (К.У.В.):

К.У.В.=N²/(4Nсr_i)=N/(4сr_i)(1.1.33)

Для получения результата, постулированного в уравнении (1.1.32), представим индексы k и n в формуле:

X_k=сX(n)*W(kn)

в виде

p-1 v

k=сk_v*mr_i,

v=0 i=0

где k_v=0,1,2,...,r_v+1-1,r₀=1.(1.1.34)

p-1 v

n=сn_v*mr_p+1-i,

где n_v=0,1,2,...,r_p-v-1,r_p+1=1.

v=0 i=0

Последние формулы можно переписать в виде:

k=k₀+k₁*r₁+k₂*r₁*r₂+...+k_p-1(r₁r₂…r_p-1)

n=n₀+n₁*r_p+n₂*r_p*r_p-1+…+n_p-1(r_pr_p-1...r₂),(1.1.35)

где k₀=0,1,2,...,r₁-1, n₀=0,1,2,...,r_p-1

k₁=0,1,2,...,r₂-1, n₁=0,1,2,...,r_р-1-1

k_р-1=0,1,2,...,r_р-1, n_р-1=0,1,2,...,r₁-1

Фиксируя поочередно k_v и n_v, можно непосредственно убедиться в том, что величины k и n принимают значения от 0 до N-1, где N есть произведение всех значений r_i. Перепишем теперь уравнение Х(k):

r_p-1 r_p-1-1 r₁-1

X(k)=X(k₀,k₁,...,k_p-1)= сnnnnnnсnnnbbbnnс

n₀=0 n₁=0 n_p-1=0

X(n₀,n₁,...,n_p-2,n_p-1)*W(kn),(1.1.36)

где

W(kn)=W(k[n₀+n₁*r_p+...+n_p-v(r_p*r_p-1...r_v+1)...+n_p-1(r_p*r_p-1...r₂)]),(1.1.37)

а величина

k=k₀+k₁*r₁+k₂*r₁*r₂+...+k_p-1(r₁*r₂...r_p-1).

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

k=(k₀+k₁*r₁+...+k_v-1*r₁*r₂*...*r_v-1)*n_p-v(r_p*r_p-1...r_v+1)+

N*n_p-v(k_v+k_v+1*r_v+1+...+k_p-1*r_v+1*r_v+2*...*r_p-1).(1.1.38)

Далее при любой целой степени N величина W равна единице, поэтому при V=1,2,...,P:

W(k*n_p-v*r_p*r_p-1*...*r_v+1)=W([k₀+k₁*r₁+...+k_v-1*r₁*r₂*...*r_v-1]*

n_p-v*r_p*r_p-1*...*r_v+1).(1.1.39)

При v=1:

W(k*n_p-1*r_p*r_p-1*...*r₂)=W(k₀*n_p-1*r_p*r_p-1*...*r₂)=

W(k₀*n_p-1*N/r₁)=EXP(-i(2bk₀n_p-1/r₁)). (1.1.40)

Мы получили экспоненциальное выражение для преобразования Фурье функции Х(n_р-1), состоящей из r₁, а не из N членов.

Переменные k₀ и n_р-1 принимают значения от 0 до r₁-1 и поэтому для расчета каждого значения Х(k₀) нужно всего r²₁ операций умножения и сложения.

Нетрудно показать, что при v=2,3,...,p:

W[(k₀+k₁*r₁+...+k_v-2*r₁*r₂*...r_v-2)n_p-v*r_p*r_p-1*...*r_v+1]*

W(k_v-1*r₁*r₂*...*r_v-1*n_p-v*r_p*r_p-1*...*r_v+1)=W(k*n_p-v*r_p*r_p-1*r_v+1). (1.1.41)

Здесь второй сомножитель содержит величину k_v-1.

Он определяет собой выражение

W(k_v-1*n_p-v*N/r_v)=EXP(-i(2bk_v-1*n_p-v/r_v)),(1.1.42)

т. е., экспоненту, необходимую для преобразования Фурье функции

Х(n_р-v), состоящей из r_v членов.

Переменные k_v-1 и n_p-v принимают значения от 0 до r_v-1.

Для вычисления каждого значения Х(k_v-1) нужно всего r²_v операций умножения и сложения.

С учетом приведенных выше преобразований для W(kn) получим:

W(kn)=mT(k₀,k₁,...,k_v-2)*W(k_v-1*n_p-v*N/r_v),(1.1.43)

где

T(k₀,k₁,...,k_v-2)=1 при v=1

T(k₀,k₁,...,k_v-2)=W[(k₀+k₁*r₁*...+k_v-2*r₁*r₂*...*r_v-2)*

(n_p-v*r_p*r_p-1...*r_v+1)]

при v=2,3,...,p.

Подставим (1.1.43) в (1.1.36), получим

rp-1 r_p-1-1

X(k₀,k₁,...,k_p-1)=сT(k₀,k₁,...,k_p-2)*W(k_p-1*n₀*N/r_p)*сT(k₀,k₁,...,k_p-3)*

n₀=0 n₁=0

r₂-1 r₁-1

*сT(k₀)*W(k₁*n_p-2*N/r₂)*сX(n₀,n₁,...,n_p-2,n_p-1)*W(k₀*n_p-1*N/r₁)

n_p-2=0 n_p-1=0(1.1.45)

Согласно формуле (1.1.45), искомое преобразование Фурье может быть построено за p итераций.

1.2 Предварительная обработка данных

Дискретизация процессов.

Пусть в моменты времени, отстоящие друг от друга на интервал времени et=h, производится выборка из отдельной реализации U(t) случайного процесса с непрерывным временем. Величина этого интервала определяет частоту свертывания (или частоту Найквиста) [1]

а_c=1/(2et)=1/2h

Пусть {Un}, n=1,2,...,N.(1.2.1) (1.2.2)

-численные значения реализации в N точках

t_n=t₀+n*h ,n=1,2,...,N.(1.2.3)

Точка t₀ выбирается произвольно и в случае, если принимает значения от 1 до N, а не от 0 до N-1, в дальнейшие формулы не входит.

Рисунок 1.2.1. Дискретное представление непрерывного процесса.

Определение объема выборки.

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

Исходные данные могут быть представлены в виде:

U_n=U(t₀+nh), n=1,2,...,N.(1.2.4)

Вычисление среднего значения.

Выборочное среднее значение находится в виде:

US=1/NсU_n,(1.2.5)

где N -число отсчетов, а U_n -значение отсчетов. Рассчитываемая по этой формуле величина US представляет собой несмещенную оценку истинного среднего значения .

Приведение процесса к нулевому среднему значению.

Для упрощения последующих расчетов необходимо преобразовать процесс таким образом, чтобы среднее его значение было равно нулю. Определим новую реализацию в виде X(t)=U(t)-US.

Тогда последовательность {Хn} значений функции определяется в виде

X_n=X(t₀+nh)=U_n-US, n=1,2,...,N.(1.2.6)

Заметим, что Х=0. Цель представления исходного процесса в виде последовательности {X_n}, а не {U_n} состоит в том, чтобы показать, что среднее значение последовательности {Х_n} равно нулю.

Вычисление стандартного отклонения.

Выборочное стандартное отклонение определяется как

S=[с(X_n)²/(N-1)]^1/2(1.2.7)

где N -- число отсчетов, а X_n -- значение преобразованного процесса с нулевым средним. Рассчитываемые по этой формуле величины S и S² представляют собой несмещенные оценки истинных значений стандартного отклонения s_x и дисперсии s²_x.

Приведение к единичному стандартному отклонению.

При использовании вычислительной машины с фиксированной, а не с плавающей запятой удобно выполнить еще одно преобразование процесса. Умножая преобразованные значения Х_n на 1/S получим последовательность

Z_n=X_n/S, n=1,2,...,N.(1.2.8)

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

Исключение тренда

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

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

Пусть для тренда, содержащегося в анализируемых данных, требуется найти приближение в виде полинома степени К

Un=с^Bk*(nh)^k, n=1,2,...,N.(1.2.9)

Согласно методу наименьших квадратов, последовательность коэффициентов {В_к} выбирается таким образом, чтобы неотрицательная при любых значениях b=(b₀,b₁,...,b_k) величина

Q(b)=с(U_n-U_n)²=с[U_n-сb_k*(nh)^k]²(1.2.10)

была наименьшей. Искомая последовательность коэффициентов отыскивается путем приравнивания нулю частных производных уравнений (1.2.10) по переменной В_e:

rQ/rb_e=с2*[U_n-сb_k*(nh)^k]*[-(nh)^e](1.2.11)

В результате получается система из К+1 уравнений вида:

сb_k*с(nh)^k+e=сU_n*(nh)^e(1.2.12)

решая которую находим искомые значения {b_к}.

Исключение тренда является важным промежуточным этапом цифрового анализа.

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

Применение цифровых фильтров.

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

Общее соотношение между процессами Х(t) на входе и Y(t) на выходе линейного фильтра дается интегралом свертки:

Y(t)=dh(i)*(X(t-i)dio(1.2.13)

где h(i) -- весовая функция фильтра.

Его частотная характеристика представляет собой преобразование Фурье функции h(i):

H(а)=dh(i)*e^-2piаidi(1.2.14)

Эквивалентная уравнению (1.2.13) конечная сумма при t=K*et, К=1,2,...,М, может быть записана в виде симметричного фильтра

Y_n=сh_K*(X_n+K+X_n-K), n=1,2,...,М.(1.2.15)

для которого h_K=h_-K.

В случае симметричного фильтра эквивалентная уравнению (1.2.14) конечная сумма определяет фильтр, фазовая частотная характеристика которого равна нулю:

H(а)=2*сh_K*cos(2bKаet)(1.2.16)

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

Общее число коэффициентов h_K (называемых весами фильтра) в формулах (1.2.15) и (1.2.16) равно М. Эти веса определяются обратным преобразованием Фурье уравнения (1.2.16):

hK=dH(а)*cos(2bKаet)dаb (1.2.17)

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

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

Простейшее выражение рекурсивного фильтра имеет вид:

Y_n=C*X_n+сh_K*Y_n-K(1.2.18)

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

Преобразование Фурье равенства (1.2.18) дает

Y(а)=C*X(а)+Y(а)*сh_K*e^-2biаKet(1.2.19)

причем сумма содержит многочлен по степеням экспоненты:

EXP[-2biаKet].

Обозначая последнее выражение символом Z, можно воспользоваться для изучения свойств цифровых фильтров теорией Z преобразований.

Как следует из формулы (1.2.19), частотная характеристика всей системы имеет вид:

H(а)=y(а)/X(а)=C/(1-сh_K*e^-2biаKet).(1.2.20)

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

1.3 Методы корреляционного и регрессионного анализа

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

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

Существование этих связей и их тесноту можно выразить коэффициентом корреляции g[2].

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

с(X_i-XS)*(Y_i-YS)

где XS,YS-выборочные средние.

Выборочный коэффициент корреляции r_xy лежит в пределах между -1 и +1. Граничные значения достигаются только в том случае, когда наблюдения обнаруживают идеальную линейную зависимость. Если же зависимость отлична от линейной и (или) наблюдается разброс измеренных значений, то независимо от того, обусловлено ли это обстоятельство ошибками измерений или нелинейным характером связи исследуемых величин, коэффициент уменьшается.

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

W=(1/2)ln[(1+r_xy)/(1-r_xy)] (1.3.2)

Как известно, случайная величина W приближенно подчиняется гауссовскому распределению со средним значением:

_W=(1/2)ln[(1+g_xy)/(1-g_xy)](1.3.3)

и дисперсией:

p²_W=1/(N-3)(1.3.4)

Зная оценку r_xy, нетрудно найти на основании этих соотношений доверительные интервалы для коэффициента g_xy.

Область принятия гипотезы о равенстве коэффициента корреляции нулю определяется неравенством:

-Z_l/2sf(N-3)/2*ln[(1+r_xy)/(1-r_xy)]sZ_l/2(1.3.5)

где Z_l/2 -- величина, подчиняющаяся нормированному гауссовскому распределению. Если рассматриваемая величина W лежит вне приведенного интервала, то это означает наличие корреляционной зависимости при уровне значимости l.

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

Методы решения задач подобного рода рассматриваются в регрессионном анализе.

Рассмотрим простой случай двух коррелированных случайных величин Х и Y. Наличие линейной зависимости между двумя этими величинами означает, что У можно вычислить при заданном значении Х, пользуясь уравнением:

Y=A+Bx (1.3.6)

где А и В -- свободный член и тангенс угла наклона прямой линии.

При идеальной линейной корреляции между случайными величинами (r_xy=1) предсказанное значение Y_i при любом Х_i всегда равно измеренному значению Y_i. На практике, однако, данные измерений далеко не всегда обнаруживают идеальную линейную зависимость. Обычно наблюдается некоторый разброс точек, обусловленный наличием посторонних случайных помех, и возможно, искажениями за счет нелинейных эффектов. Тем не менее если допустить существование линейной зависимости между величинами и располагать неограниченно большим количеством измерений, то можно найти такие подходящие значения коэффициентов А и В, которые позволяют предсказать ожидаемые значения Y_i при любых заданных Х_i.

Значение Y_i не обязательно совпадает с измеренной величиной Y_i, соответствующей данному значению Х_i; оно представляет собой некоторое среднее для всех таких измеренных величин.

Принятый на практике способ определения коэффициентов уравнения (1.3.6) состоит в том, что выбираются такие значения А и В, при которых сумма квадратов отклонений измеренных величин от предсказанных была бы минимальной. Этот способ, как известно, называется методом наименьших квадратов. Отклонение измеренной величины от предсказанной

Y_i-Y_i=Y_i-(A+bX_i)(1.3.7)

Поэтому сумма квадратов отклонений выражается в виде:

Q=с(Y_i-A-bX_i)²(1.3.8)

Следовательно, наименьшее значение суммы квадратов отклонений достигается в том случае, когда коэффициенты А и В удовлетворяют условию:

rQ/rA=rQ/rB=0(1.3.9)

Практически имеющиеся данные представляют собой ограниченную выборку, состоящую из N пар измеренных значений Х и Y.

Это означает, что условие (1.3.9) позволяет получить лишь оценки коэффициента А и В, которые обозначим соответственно символами а и b.

Подставляя (1.3.8) в (1.3.9) и решая полученное уравнение относительно оценок коэффициентов А и В, находим:

a=Y-b*X

сX_i*Y_i-N*X*Y

Эти оценки можно использовать с тем, чтобы записать формулу для прогноза величины при заданном Х:

Y=Y+b(X-X)(1.3.11)

Прямая линия, описываемая уравнением (1.3.11), называется линией регрессии Y по Х.

Меняя местами зависимую и независимую переменные в уравнении

X=X+b(Y-Y),(1.3.12)

Перемножая b и b и сравнивая полученное произведение с формулой (1.3.1), нетрудно убедиться, что тангенсы углов наклона линий регрессии Y по Х и Х по Y связаны с выборочным коэффициентом корреляции величин Х и Y соотношением:

r_xy=[b*b]^1/2(1.3.13)

Рассмотрим теперь, какова точность вычисления оценок a и b, если величина Y при заданном Х подчиняется нормальному распределению, то, как известно, величины а и b представляют собой соответственно несмещенные оценки коэффициентов А и В.

Выборочные распределения этих оценок связаны с t-распределением Стьюдента следующим образом:

(b-В)*(с(X_i-X)²)^1/2=S_y/xt_N-2 (1.3.15)

Особый интерес представляет выборочное распределение величины Y, соответствующей некоторому фиксированному значению Х=Х₀.

1.4 Косинор-анализ

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

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

Дадим описание алгоритма косинор-анализа. Пусть известны экспериментально полученные измерения U₁,U₂,...,U_i,...U_m и различные моменты времени t₁,t₂,...,t_i,...t_m (система измерений).

Первый шаг обработки состоит в аппроксимации измерений методом наименьших квадратов косинусоидой:

U(t)=A*cos(W₀t-t)+h=X*cos(W₀t)+Y*sin(W₀t)+h,(1.4.1)

где Х=А*cost, Y=А*sint.

Здесь угловая частота W₀ предполагается заданной, а амплитуда А, акрофаза t и уровень h подлежат определению.

Чтобы найти их, минимизируем выражение:

I=с{U(t_i)-U_i})² или p=fI/m(1.4.2)

откуда приходим к системе трех уравнений с тремя неизвестными -- Х, Y, h:

a₁₁X+a₁₂Y+a₁₃h=a₁₄

a₂₁X+a₂₂Y+a₂₃h=a₂₄(1.4.3)

a₃₁X+a₃₂Y+a₃₃h=a₃₄

Здесь

m

a₁₁=сcos²(w₀t_i),

i=1

m

a₁₂=a₂₁=1/2сsin(2w₀t_i),

i=1

m

a₁₃=a₃₁=сcos(w₀t_i),

i=1

m

a₁₄=сU_icos(w₀t_i),

i=1

m

a₂₂=сsin²(w₀t_i),

i=1

m

a₂₃=a₃₂=сsin(w₀t_i),

i=1

m m

a₂₄=сU_isin(w₀t_i), a₃₃=m, a₃₄=сU_i

i=1 i=1

Решая эту задачу, получаем Х, Y, h, а затем находим А и t по формулам:

A=fX₂+Y₂(1.4.5)

t=arctg(Y/X) если Х>0,X=0.

t=arctg(Y/X)+b если Х<0.

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

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

3.Рассмотрена алгоритмическая сущность метода наименьших квадратов для построения полиномиальных моделей прогноза.

4.Описан косинор-анализ,который используется для построения простейших гармонических трендов почвенных параметров.

Глава II. Применение факторного анализа

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

В настоящей главе изложен один из главных методов факторного анализа, метод главных компонент [4,5]. Кроме этого, приведен разработанный нами алгоритм по автоматическому определению объединений параметров по факторам, который значительно облегчает определение групповых качественных обусловленностей .

2.1 Метод главных компонент

Метод главных компонент совпадает с методом расчленения ковариационной или корреляционной матрицы на совокупность ортогональных векторов (компонент) или направлений по числу рассматриваемых переменных.

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

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

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

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

В факторном анализе основным предположением является равенство:

X_i=сa_ir*F_r+e_i(2.1.1)

где X_i--ая переменная,

F_r- r-ый фактор,

a_ir -- факторная нагрузка,

k -- количество факторов ,

e_i -- остатки, которые представляют источники отклонений, действующие только на X_i.

Эти р случайных величин e_i предполагаются независимыми как между собой, так и с k величинами F_r.

Уравнение (2.1.1) нельзя проверить непосредственно, поскольку р переменных X_i выражены в них через (р+k) ненаблюдаемых переменных.

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

Когда число факторов k>1, то ни факторы, ни нагрузки не определяются однозначно, поскольку в уравнении (2.1.1) факторы F_r могут быть заменены любым ортогональным преобразованием их с соответствующим преобразованием нагрузок. Это свойство использовано для преобразования или вращения факторов, полученных в каком-либо практическом исследовании.

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

2.2 Формализация метода главных компонент

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

Классическая модель факторного анализа имеет вид:

R=A^*C^*A',(2.2.1)

где R -- корреляционная матрица,

А -- матрица факторных нагрузок,

С -- корреляционная матрица, отражающая связи между факторами,

А' -- транспонированная матрица факторных нагрузок.

Если наложить условие некоррелированности факторов, т. е. С=I,

где I -- единичная матрица, то в результате получим:

R=A^*A'(2.2.2)

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

s₁=сa_i1=max(2.1.3)

при m(m-1)/2 независимых друг от друга условиях

r_ik=a_i1*a_k1 (i,k=1,2,...,m,i<k)(2.2.4)

где m -- число переменных в матрице наблюдений.

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

В результате приходим к системе m однородных уравнений с m неизвестными a_i1:

(1-k)a₁₁+r₁₂*a₂₁+...+r_1m*a_m1=0

r₂₁*a₁₁+(1-k)a₂₁+...+r_2m*a_m1=0

r_m1*a₁₁+r_m2*a₂₁+...+(1-k)a_m1=0

Система равенств (2.2.5) составляет проблему собственных значений действительной симметричной матрицы.

В общем она записывается в следующем виде:

Rl_e=k*l_e, (R-kI)l_e=0,(2.2.6)

Где nkn-- собственные значения, они соответствуют собственным векторам l_e матрицы R.

Факторы пропорциональны собственным векторам матрицы R.

Путем нормирования получим искомые значения l_ie матрицы А по компонентам собственных векторов матрицы R:

a_ie=l_iefk_e/fl²_1e+l²_2e+...+l²_me. (2.2.7)

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

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

С учетом характерных факторов уравнение (2.2.2) будет иметь вид:

R=A^*A'+U².(2.2.8)

Для определения максимально правдоподобных оценок l_ie и U²_e максимизируется функция правдоподобия. Это приводит к большому числу возможных результатов. Из них выбирается тот, который удовлетворяет условию, что

I=A'[U²]^-1A(2.2.9)

является диагональной матрицей.

Равенство (2.2.9) приводится к виду [12]:

A=I^-1*A'[U²]^-1(R-U²).(2.2.10)

Кроме этого, имеет место равенство:

U²_i=1-сa²_ie(2.2.11)

Исходя из более или менее произвольно взятых первых приближений для А и U² в правой части (2.2.10), в левой части этого равенства получают новую матрицу А, а из (2.2.11) -- новую матрицу U², которые можно рассматривать как хорошие приближения к истинным матрицам.

Эти полученные матрицы опять подставляются в (2.2.10) и (2.2.11), и итеративная процедура повторяется.

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

Итеративную процедуру решения уравнения (2.2.10) можно рассматривать в качестве особого способа определения собственных значений и собственных векторов матрицы [U²]-¹(R-U²), причем А' содержит собственные векторы, а диагональная матрица I^-1 собственные значения этой матрицы.

Тест проверки значимости числа факторов, найденного методом максимального правдоподобия:

H²=(n-1)ln|R⁺|/|R| (2.2.12)

В этой формуле |R⁺| -- определитель матрицы корреляций, воспроизведенных с помощью выбранной модели; |R| -- определитель исходной корреляционной матрицы, n -- число наблюдений.

Величина H² имеет приблизительно ² распределение с 1/2[(m-r)²-m-r] степенями свободы.

Здесь m -- число переменных, r -- число выделенных факторов.

Если при определенном r вычисленное значение критерия превышает табличное значение H², соответствующее заданному уровню значимости, это указывает на то, что необходимо выделить факторов больше, чем r, по крайней мере r+1.

После определения матрицы факторных нагрузок А для лучшей интерпретации факторов используют вращение A в пространстве общих факторов.

В настоящее время наиболее распространен “метод варимакс” для осуществления вращения матрицы А. Метод варимакс был предложен Кайзером [7].

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

Дисперсия квадратов нагрузок фактора е равна:

S²_e=1/mс(b²_ie)-1/m²*[с(b_2ie)]²(2.1.13)

Просуммируем эту дисперсию по всем факторам. Полученная в результате этого величина будет максимальная в случае, когда дисперсия квадратов нагрузок каждого фактора примет наибольшее значение:

сS²_e=1/mсnсb⁴_ie-1/m²сn[nсb²_ie]²=max(2.1.14)

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

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

Таким образом, при определении положения осей координат имеют дело с нормированными переменными с равными весами.

В отличие от (2.2.14) модифицированный Кайзером варимакс-критерий умножается еще на m:

Mсnnnс(b_ie/h_i)⁴-сn[nсb²_ie/h²_i]²=max(2.1.15)

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

2.3 Алгоритм метода главных компонент

1. Расчет корреляционной матрицы R:

r_jk=S_jk/(fS_jj*fS_kk),(2.3.1)

где

S_jk=с(X_ij-T_j)*(X_ik-T_k)-1/nс(X_ij-T_j)*с(X_ik-T_k), T_j=(сX_ij)/n,

i=1,2,...,n -- наблюдения,

j=1,2,...,m -- переменные.

2. Вычисление собственных значений, собственных векторов корреляционной матрицы.

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

4. Вычисление матрицы факторных нагрузок по собственным значениям и соответствующим собственным векторам корреляционной матрицы.

5. Ортогональное вращение матрицы факторов.

2.4 Алгоритм минимизации количества параметров исследования

Предлагается следующий алгоритм минимизации числа параметров исследования для многопараметрических объектов:

1.Строим матрицу исследования (строчки-наблюдения, столбцы-параметры исследования).

2.Методом главных компонент находим матрицу факторных нагрузок. Осуществляем варимаксное вращение в пространстве факторов(строчки в матрице факторных нагрузок-параметры исследования,столбцы- гипотетические переменные, факторы).

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

4.Определяем по каждому фактору попадание в этот фактор параметров с максимальной по модулю факторной нагрузкой(пункт 3).То есть тем самым определяем объединение параметров по факторам.

5.В объединившихся в каждом факторе параметрах выбираем один параметр с максимальной по модулю факторной нагрузкой.

Число таких выбранных параметров будет равно, очевидно, числу факторов.

6.Строим для всех параметров исследования полиномиальные модели, аргументами в которых будут выбранные в пункте 5 параметры.

7.По построенным моделям для каждого параметра осуществляем определение вкладов параметров-аргументов (оценку количественной обусловленности параметров выбранными параметрами).

8.Сравниваем качественные групповые обусловленности, объединения параметров по факторам, с количественными обусловленностями параметров, полученными в пункте 7.

Если групповые и количественные обусловленности для всех параметров исследования не будут сильно отличаться по числу не совпадений, то выбранные в пункте 5 параметры могут быть приняты за базисные при описании данного многопараметрического объекта, матрица исследования которого была взята за основу в данном алгоритме. То есть тем самым осуществляем минимизацию количества параметров исследования, потому что число факторов меньше числа параметров.

2.5 Программная реализация метода главных компонент

Для программной реализации метода главных компонент был использован алгоритмический язык ФОРТРАН-V и С++(приложение-1).

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

Обращение: CALL CORRE(N,M,ID,X,XBAR,STD,RX,R,B,D,T),

где N -- число наблюдений;

M -- число переменных;

ID -- выбираемый код для входных данных;

О -- данные должны быть прочитаны из устройства ввода специальной подпрограммой DATA;

1 -- все данные в памяти;

Х -- вводимая матрица, содержащая данные, размером NXM;

XBAR -- выходной вектор размера М, содержащий выборочные средние значения;

STD -- выходной вектор размера М, содержащий выборочные стандартные отклонения;

RX -- выходная матрица МxМ, содержащая выборочные ковариации, т. е. S_jk;

R -- выходной массив размера М(М+1)/2, содержащий верхний треугольник симметричной матрицы МхМ выборочных коэффициентов корреляции.

При работе с подпрограммой CORRE использовался режим ID=1

В этом случае во входном потоке должна быть подпрограмма:

SUBROUTINE DATA

RETURN

END

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

Обращение:

CALL EIGEN (А,R,N,MV),

где А -- исходная матрица (симметричная), изменяется в процессе вычислений.

Вычисленные собственные значения располагаются на диагонали матрицы А в порядке убывания.

R -- матрица собственных векторов, расположенных по столбцам соответственно собственным значениям.

N -- порядок матрицы А и R.

MV -- входной параметр, определяющий характер вычислений.

О -- вычисляются собственные значения и собственные векторы.

1 -- вычисляются только собственные значения (массив R в этом случае не используется, но имя его должно быть обязательно в обращении).

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

Обращение:

CALL TRACE (М,R,CON,K,D),

где М -- число переменных;

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

Собственные значения размещены в убывающем порядке.

Размерность матрицы R есть МxМ, но только М(М+1)/2 элемента находятся в памяти.

СОN -- константа, используемая для выбора собственных значений. Вычисляются накопленные отношения собственных значений, которые больше или равны этой константе.

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

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

Последние (М--К) компонент этого вектора нулевые.

Матрица факторных нагрузок определяется с помощью подпрограммы LOAD.

Обращение: CALL LOAD (М,К,R,V),

где М -- число переменных;

К -- число факторов;

R -- входная редуцированная корреляционная матрица (симметричная, верхний треугольник задан по столбцам), содержащая собственные значения на диагонали;

V -- матрица размера МxМ, содержащая нормализованные собственные векторы по столбцам, на выходе -- это матрица факторных нагрузок размера МxК.

Ортогональное вращение факторов осуществляется подпрограммой VARMX.

Обращение: СALL VARMX (M,K,A,NC,TV,H,F,D),

где М -- число переменных и число строк матрицы А;

К -- число факторов;

А -- входная начальная факторная матрица и выходная факторная матрица, полученная в результате вращений;

NC -- выходная переменная, содержащая число выполненных итерационных циклов;

TV -- выходной вектор, содержащий дисперсию факторной матрицы для каждого итерационного цикла, размер вектора TV-51;

Н -- выходной вектор размера М, содержащий начальные суммарные нагрузки простых факторов;

F -- выходной вектор размера М, содержащий конечные суммарные нагрузки простых факторов;

D -- выходной вектор размера М, содержащий разности между начальными и конечными суммарными нагрузками простых факторов.

Выводы

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

2.Рассмотрен алгоритм метода главных компонент для нахождения элементов матрицы факторных нагрузок.

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

4.Дано описание входных и выходных параметров подпрограмм по реализации корреляционного и факторного анализов.

Глава III. Применение регрессионного метода

Для построения нелинейных регрессионных моделей с автоматическим выбором степени аппроксимирующих полиномов в настоящей главе рассматривается ступенчатый регрессионный метод [3].

3.1 Постановка задачи

Первым этапом исследования многопараметрических процессов является отбор параметров, ответственных за процесс. Из полного списка всех возможных параметров ранговыми методами производят их ранжирование и априорное отсеивание [2].

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

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

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

y=b₀+сb_i*X_i+сb_ij*X_i*X_j+…. (3.1.1)

где y -- параметр оптимизации;

b₀,b_i,b_ij,b_ii -- выборочные коэффициенты регрессии, полученные по результатам эксперимента;

X_i,X_iX_j -- параметры и их взаимодействия, i,j=1,2…,

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

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

Упрощенный метод определения коэффициентов уравнения регрессии (3.1.1) предложен в работе Д. Брандона [3].

Этот метод заключается в том, что уравнение (3.1.1) записывается в виде:

Y=a*mа_S(X_S)(3.1.2)

где а_S(X_S) -- любая функция величины Xs.

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