diedai.rar_MATLAB 海浪_参数反演_海流反演_海浪 雷达_海浪反演算法

%---(自适应阈值选取技术迭代法反演流速) %---第一步粗略的流速估计） %-----(较为准确的流速估计:(取总的n%个数据点，模拟的数据在雷达的正北方))------------------------------ clear all; [filename,pathname]=uigetfile({'*.txt','选择文件'}); %选择路径 str=[pathname,filename]; a=load(str); % a=load('f:\image3D-zabo.txt'); a=a'; %(数据是一行一行导入) b=reshape(a,128,128,64); M=128; %采样频率 N=64; %采样频率 g=9.8; dx=7.5; dy=7.5; dt=60/42; T=63*dt; L_x=127*dx; L_y=127*dy; kx=2*pi/L_x; %分辨率 ky=2*pi/L_y; %分辨率 kw=2*pi/T; %分辨率 W_N=pi/dt; %截断频率 kkx=(-M/2:M/2-1).*kx; % k_x频率 kky=((-M/2:(M/2-1)).*ky); % k_y频率 [kkx,kky]=meshgrid(kkx); kky=-kky; kkk=sqrt(kkx.*kkx+kky.*kky); kkw=(0:N/2-1).*kw; %角频率 b1=fftn(b); %三维傅里叶变化 b2=fftshift(b1); %频率移中 IW=abs(b2).*abs(b2); %图像功率谱 I=abs(b2(:,:,33:64)).*abs(b2(:,:,33:64)); %只是考虑正频率部分 I(65,65,:)=0; %去掉（0,0）点 I(:,:,1:3)=0; %去掉0频率点、 % I(19,16,22)=0; %(关键点) I0=I; I1=I; %(迭代时用到) %--------------(选定阈值为N个坐标数)--------- MAX=max(I0(:)); %求出最大值 reshape_I0=reshape(I0,1,128*128*32); sort_I0=sort(reshape_I0); sort_I0=fliplr(sort_I0); value1=sort_I0(20)/MAX; %(所取的点数) % value1=0.6; JL=ones(128,128,32); for i=1:128; for j=1:128; for k=1:32; if (I0(i,j,k)/MAX)>value1; I0(i,j,k)=I(i,j,k); else I0(i,j,k)=0; JL(i,j,k)=0; end end end end %-----------------(第一次估算流速：选择阈值C1=N个数值)--------------------- t1=0.0; t2=0.0; t4=0.0; for i=1:128; for j=1:128; for k=1:32; t1=t1+kkx(i,j)*kkx(i,j)*I0(i,j,k); t2=t2+kkx(i,j)*kky(i,j)*I0(i,j,k); t4=t4+kky(i,j)*kky(i,j)*I0(i,j,k); end end end t3=t2; t5=[t1,t2;t3,t4]; t6=inv(t5); %逆矩阵 t7=0.0; t8=0.0; for i=1:128; for j=1:128; for k=1:32; t7=t7+kkx(i,j)*I0(i,j,k)*(kkw(k)-sqrt(g*kkk(i,j))); t8=t8+kky(i,j)*I0(i,j,k)*(kkw(k)-sqrt(g*kkk(i,j))); end end end t9=[t7;t8]; v=t6*t9; %流速 v_x=v(1,1); v_y=v(2,1); % v_x=10; %设定流速 % v_y=10; P1=0; for i=1:128; for j=1:128; for k=1:32; if I0(i,j,k)~=0; P1=P1+1; end end end end P1; figure(1); plot(v_x,v_y,'b*'); %------------------------(选择阈值C2：比C1小，难点)----------------------------- for num=1:10; JS0=0; %记录0阶数目 JS1=0; %记录1阶数目 I1=I; MAX=max(I1(:)); %求出最大值 reshape_I1=reshape(I1,1,128*128*32); sort_I1=sort(reshape_I1); sort_I1=fliplr(sort_I1); value2=sort_I1(300)/MAX; %(所取的点数) % value2=value1/10; count=1; for k=1:32; for i=1:128; for j=1:128; if (I1(i,j,k)/MAX)>value2; I1(i,j,k)=I(i,j,k); fw(1,count)=kkw(k); bsx(1,count)=kkx(i,j); bsy(1,count)=kky(i,j); xzb(1,count)=i; yzb(1,count)=j; zzb(1,count)=k; count=count+1; else I1(i,j,k)=0; end end end end %--------------（把第一次得到的流速带入计算和一次谐波）------------- %-------------(计算一次谐波，判断一次谐波范围，并移频)------------------------------- len=length(xzb); W_1=zeros(1,len); for i=1:len; W_1(1,i)=sqrt(2)*sqrt(g*kkk(xzb(i),yzb(i)))+bsx(1,i)*v_x+bsy(1,i)*v_y; end figure(2); subplot(221) plot(W_1,'g'); %一阶谐波 hold on plot(ones(1,len)*W_N,'r','linewidth',1); %一倍尼奎斯特频率 hold on plot(ones(1,len)*2*W_N,'r','linewidth',1); %2倍的尼奎斯特频率 for i=1:len; % A=mod(W_1(1,i),W_N); %取余数 % B=round((W_1(1,i)-A)/(W_N)); %取商 A=rem(W_1(1,i),W_N); %取余数 B=floor(W_1(1,i)/W_N); %取商 BB1(1,i)=B; % C=round(mod(B,2)); %判断商是不是偶数 C=rem(B,2); %判断商是不是偶数 CC1(1,i)=C; if (C==0); %是偶数 W_1(1,i)=W_1(1,i)-B*W_N; %波数坐标也要变换 kkx1(1,i)=bsx(1,i); kky1(1,i)=bsy(1,i); elseif (C==1); W_1(1,i)=-(W_1(1,i))+(B+1)*W_N; %(待讨论) % W_1(1,i)=-(sqrt(2)*sqrt(g*kkk(xzb(i),yzb(i)))-bsx(1,i)*v_x-bsy(1,i)*v_y)+(B+1)*W_N; kkx1(1,i)=bsx(1,i); kky1(1,i)=bsy(1,i); end end subplot(222) plot(W_1,'g'); %移频后的频率 hold on plot(kkw(zzb),'r'); %选出的频率 hold on plot(ones(1,len)*W_N,'r','linewidth',1); hold on plot(ones(1,len)*2*W_N,'r','linewidth',1); hold on plot(abs(W_1-kkw(zzb)),'b-','linewidth',1); %移频后与选出值的比较 %------------------(计算基模.判断基模式范围，并移频)--------------------------------------- W_0=zeros(1,len); for i=1:len; W_0(1,i)=sqrt(g*kkk(xzb(i),yzb(i)))+bsx(1,i)*v_x+bsy(1,i)*v_y; end subplot(223) plot(W_0,'g'); hold on plot(ones(1,len)*W_N,'r','linewidth',1); hold on plot(ones(1,len)*2*W_N,'r','linewidth',1); for i=1:len; % A=mod(W_0(1,i),W_N); %取余数 % B=round((W_0(1,i)-A)/(W_N)); %取商 A=rem(W_0(1,i),W_N); %取余数 B=floor(W_0(1,i)/W_N); %取商 BB0(1,i)=B; % C=round(mod(B,2)); %判断商是不是偶数 C=rem(B,2); %判断商是不是偶数 CC0(1,i)=C; if (C==0); %是偶数 W_0(1,i)=W_0(1,i)-B*W_N; kkx0(1,i)=bsx(1,i); kky0(1,i)=bsy(1,i); elseif(C==1); W_0(1,i)=-(W_0(1,i))+(B+1)*W_N; %(待讨论) % W_0(1,i)=-(sqrt(g*kkk(xzb(i),yzb(i)))-bsx(1,i)*v_x-bsy(1,i)*v_y)+(B+1)*W_N; kkx0(1,i)=bsx(1,i); kky0(1,i)=bsy(1,i); end end subplot(224) plot(W_0,'g'); hold on plot(kkw(zzb),'r'); hold on plot(abs(W_0-kkw(zzb)),'m-','linewidth',1); hold on plot(ones(1,len)*W_N,'r','linewidth',1); hold on plot(ones(1,len)*2*W_N,'r','linewidth',1); %--------------（判断选择的属于基模，还是一次谐波）----------------------------------- ww=zeros(1,len); for i=1:len; ww(1,i)=kkw(zzb(i)); %选出的角频率坐标 end chazhi_0=abs(W_0-ww); %判断属于哪一阶次谐波 chazhi_1=abs(W_1-ww); for i=1:len; if (chazhi_0(i))<=(chazhi_1(i)); %---------------------------（属于0阶）------------------------------ JS0=JS0+1; % I1(xzb(i),yzb(i),zzb(i))=0; %只考虑一次谐波 if CC0(1,i)==0; %为截断频率的偶数倍 % I1(xzb(i),yzb(i),zzb(i))=0; %不考虑偶数倍 Dxx(1,i)=kkx0(1,i).^2; Dxy(1,i)=kkx0(1,i)*kky0(1,i); Dyx(1,i)=Dxy(1,i); Dyy(1,i)=kky0(1,i).^2; Bx(1,i)=kkx0(1,i)*(ww(1,i)-sqrt(g*kkk(xzb(i),yzb(i)))+BB0(1,i)*W_N); %(待讨论) By(1,i)=kky0(1,i)*(ww(1,i)-sqrt(g*kkk(xzb(i),yzb(i)))+BB0(1,i)*W_N); %(待讨论) else %为截断频率的奇数倍 % I1(xzb(i),yzb(i),zzb(i))=0; %不考虑奇数倍 Dxx(1,i)