clear
clc
close?all
%時間間隔delta_t，粗糙度系數(shù)K，10米高度處平均風(fēng)速；
delta_t=0.1;K=0.03;v10=32.6;
%%?輸入數(shù)據(jù)坐標(biāo)；
x=[101010101010108090100];y=[1111111111];z=[102030405060701505050];
%p為AR模型階數(shù)；M為空間點(diǎn)個數(shù)；
p=4;M=7;
%Vz為z高度處的風(fēng)速，Lz為z高度處的湍流積分尺度；
for?i=1:M
Vz（i）=v10*（z（i）/10）^0.16;
Lz（i）=150*（z（i）/30）.^（0.5）;
end

%L為數(shù)據(jù)長度；
L=6000*3;
%產(chǎn)生方差為1的隨機(jī)數(shù)；
n=normrnd（01ML）./5;
%偏度和峰度，當(dāng)高斯時均為0；
pd=0;fd=0;
%求系數(shù)矩陣；
[posiRuR]=posi（Mpxyzdelta_tKv10pdfd）;
%求下三角陣LN；
R0=Ru{11};
RN=R0-posi‘*R;
LT=chol（RN）‘;
%求隨機(jī)數(shù)矩陣N；
N=LT*n;
%求平穩(wěn)高斯脈動風(fēng)速時程u；
u（1:M1:p）=0;
b=（M）*ones（1p）;
posi=mat2cell（posibM）;
for?i=1:p;
????P=zeros（M1）;
????for?k