matlab 四向剪切最小二乘相位解缠

点云侠

于 2024-08-19 08:28:23 发布

阅读量181

点赞数 4

分类专栏： InSAR学习文章标签： matlab 开发语言算法计算机视觉人工智能

本文链接：https://blog.csdn.net/qq_36686437/article/details/141277268

版权

InSAR学习专栏收录该内容

29 篇文章 9 订阅

订阅专栏

一、算法原理

1、四向最小二乘

在这里插入图片描述

图1 干涉图在x,y,p,q四个方向的相位差分

在这里插入图片描述

2、参考文献

[1]郭媛,吴全,毛琦,等. 基于横向剪切的四向最小二乘相位解包裹算法 [J]. 光电子·激光, 2015, 26 (10): 1953-1959.

二、代码实现

clc;
clear;
close all;

%% --------------------------------初始相位--------------------------------
phi0 = peaks(200)*3;
figure(1);
surf(phi0);
title('山峰3D图')
%% -------------------------------包裹相位---------------------------------
phi = angle(exp(1i*phi0));         
figure(2);
%imshow(phi,[]);
imagesc(phi)
xlabel('X/Pixels','FontSize',14);
ylabel('Y/Pixels','FontSize',14);
%title('Wrapped Phase','FontSize',14)
set(figure(2),'name','包裹相位','Numbertitle','off');
axis on

%% ------------------------------相位解包裹--------------------------------
% 对包裹相位求一阶偏微分
[M,N] = size(phi);
deltX = zeros(M,N);
deltY = zeros(M,N);
deltP = zeros(M,N);
deltQ = zeros(M,N);

% deltX(1:M-1,:) = phi(2:M,:)-phi(1:M-1,:);
% deltX = deltX - pi*round(deltX/pi);           % 去除差分中的跳跃
deltX(1:M-1,:) = angle(exp(1i*(phi(2:M,:)-phi(1:M-1,:))));           % 公式（2）
deltY(:,1:N-1) = angle(exp(1i*(phi(:,2:N)-phi(:,1:N-1))));           % 公式（3）
deltP(1:M-1,1:N-1) = angle(exp(1i*(phi(2:M,2:N)-phi(1:M-1,1:N-1)))); % 公式（4）
deltQ(1:M-1,1:N-1) = angle(exp(1i*(phi(2:M,1:N-1)-phi(1:M-1,2:N)))); % 公式（5）

% 对包裹相位求二阶偏微分
DdeltX = zeros(M,N);
DdeltY = zeros(M,N);
DdeltU = zeros(M,N);
DdeltV = zeros(M,N);
DdeltX(1:M-1,:) = angle(exp(1i*(deltX(2:M,:)-deltX(1:M-1,:))));      % 公式（10）
DdeltY(:,1:N-1) = angle(exp(1i*(deltY(:,2:N)-deltY(:,1:N-1))));
DdeltU(1:M-1,1:N-1) = angle(exp(1i*(deltP(2:M,2:N)-deltP(1:M-1,1:N-1))));
DdeltV(1:M-1,1:N-1) = angle(exp(1i*(deltQ(2:M,1:N-1)-deltQ(1:M-1,2:N))));

Rou = DdeltX + DdeltY + DdeltU + DdeltV;                             % 公式（8）
Wij = sqrt(DdeltX.^2 + DdeltY.^2 + DdeltU.^2 + DdeltV.^2);           % 公式（9）
k = 0; % 用于调节权重强度，其值由实际情况确定。
Rou_ij = (1+k*Wij).*Rou;                                                % 公式（11）

%% --------------------------DCT求解泊松方程-------------------------------
% 基于离散余弦变换DCT计算ρ_ij
tic
pp = dct2(Rou)+eps; 
% 计算Φ_ij在DCT域的精确解
fi = zeros(M,N);
for m = 1:M                                         
   for n = 1:N
       k1 = 2*cos(pi*(m-1)/M); 
       k2 = 2*cos(pi*(n-1)/N); 
       kk = k1 + k2 - 4;
      fi(m,n) = pp(m,n)/(kk+eps);
   end
end
fi(1,1) = pp(1,1);                   % 赋值DCT域的Φ11
phs = idct2(fi);                     % 用iDCT计算解缠相位在空域中的值
toc
phi333 = phs;

figure(3);
surf(Rou,'FaceColor','interp', 'EdgeColor','none','FaceLighting','phong');
camlight left, axis tight
xlabel('X/Pixels','FontSize',14);
ylabel('Y/Pixels','FontSize',14);
zlabel('Phase/Radians','FontSize',14);
%title('lou3','FontSize',14)
set(figure(3),'name','R(x,y) 3D','Numbertitle','off');
figure(4);
imagesc(Rou)
xlabel('X/Pixels','FontSize',14);
ylabel('Y/Pixels','FontSize',14);
set(figure(4),'name','R(x,y) 2D','Numbertitle','off');
figure(5);
surf(phi333,'FaceColor','interp', 'EdgeColor','none','FaceLighting','phong');
camlight left, axis tight
xlabel('X/Pixels','FontSize',14);
ylabel('Y/Pixels','FontSize',14);
zlabel('Phase/Radians','FontSize',14);
%title('Four LSBLS Phase Unwrapping','FontSize',14)
set(figure(5),'name','Four LSBLS 相位解缠结果','Numbertitle','off');