clc; clear; close all;
% 生成数据
% 原代码段：
% [x, y, z] = peaks(200);
% 修改为马鞍面函数：
x_range = linspace(-3, 3, 200); % 创建-3到3的200点均匀分布
y_range = linspace(-3, 3, 200);
[x, y] = meshgrid(x_range, y_range); % 生成网格矩阵
z = x.^2 - y.^2; % 马鞍面方程 z=x²-y²
v = z/max(abs(z(:))) * 0.04; % 标准化到[-0.04,0.04]范围 % 使用peaks函数生成基础数据
% 绘制等值线填充图
contourf(v, 40, 'LineStyle', 'none') % 40级等值线，取消线条显示
% 设置颜色映射和颜色条
colormap(jet) % 使用jet颜色映射（红-黄-绿-蓝）
caxis([-0.04 0.04]) % 固定颜色范围
cb = colorbar;
cb.Label.String = '沉降速度 (m/a)';
cb.Ticks = -0.04:0.02:0.04; % 设置刻度位置
% 添加标题和修饰
title('图8-2 模拟产生的沉降速度场')
axis equal tight % 等比例坐标轴
set(gca, 'Visible', 'off') % 隐藏坐标轴
% 生成200个随机PS点坐标（在1到200范围内）
num_ps = 200;
ps_x = rand(num_ps, 1) * 199 + 1; % 生成1-200之间的x坐标
ps_y = rand(num_ps, 1) * 199 + 1; % 生成1-200之间的y坐标
hold on;
% 绘制PS点（黑色圆点）
plot(ps_x, ps_y, 'ko', 'MarkerSize', 4, 'MarkerFaceColor', 'k');
% 构建PS网络（距离限制连接）
threshold = 40; % 连接距离阈值
distances = pdist2([ps_x, ps_y], [ps_x, ps_y]); % 计算欧氏距离矩阵
[row, col] = find(triu(distances < threshold, 1)); % 获取满足条件的点对索引
% 绘制连接线（浅灰色细线）
for i = 1:length(row)
line([ps_x(row(i)), ps_x(col(i))],...
[ps_y(row(i)), ps_y(col(i))],...
'Color', [0.6 0.6 0.6], 'LineWidth', 0.3);
end
% 保持图形修饰一致性
set(gcf, 'Color', 'w'); % 设置白色背景
hold off;
%% 高程修正模拟（新增代码段）
% 生成高斯分布的高程修正量（均值为0，方差为3）
sigma = sqrt(3); % 标准差=√方差
ps_z_correction = sigma * randn(num_ps, 1); % 生成200个修正值
% 可选：显示统计特征验证
fprintf('高程修正统计量:\n均值=%.4fm\n方差=%.4fm²\n标准差=%.4fm\n',...
mean(ps_z_correction), var(ps_z_correction), std(ps_z_correction));
% 将修正量叠加到原始高程（假设原始高程为peaks生成的z值）
% 获取PS点所在位置的高程索引
[~,idx] = min(pdist2([x(:),y(:)], [ps_x,ps_y]), [], 1);
original_z = z(idx); % 提取原始地形高程
corrected_z = original_z + ps_z_correction; % 应用高斯修正
% 可视化高程修正分布（新增图形窗口）
figure;
histfit(ps_z_correction, 20); % 带拟合曲线的直方图
title('PS点高程修正量分布');
xlabel('高程修正量(m)'); ylabel('频数');
grid on;
%% 干涉相位计算（续写部分）
% 从表8-2获取雷达系统参数
lambda = 0.056; % 波长 (m)
R_near = 5536.35965; % 近地点斜距 (m) 修正单位问题
theta = deg2rad(23); % 假设入射角23度（需根据实际情况调整）
% 从表1获取基线参数（示例数据，请核对完整参数）
time_baselines = [-2159,-2054,-1844,-1739,-1634,-1121,-946,-806,-771,...
-770,-701,-700,-175,-70,-35,350,560,595,735,770,875,1295,1505,1540,1575]; % 时间基线(天)
perp_baselines = [504,146,-36,274,-639,-207,178,505,-1144,-1000,...
-1253,-1104,-762,-1239,-487,247,-348,-141,303,-158,290,-198,1048,144,-1021]; % 垂直基线(m)
% 生成高程修正场（基于PS点插值）
F = scatteredInterpolant(ps_x, ps_y, ps_z_correction, 'natural', 'none');
[x, y] = meshgrid(1:200, 1:200);
delta_h = F(x,y); % 插值得到全场高程修正量
% 相位计算核心循环
for i = 1:length(time_baselines)
% 时间参数转换
dt_year = (time_baselines(i)+2159)/365.25; % 转换到年单位
% 形变相位分量（全场计算）
phi_def = (4*pi/lambda) * v * dt_year;
% 高程相位分量（考虑基线空间去相关）
B_perp = perp_baselines(i)-504;
phi_top = (4*pi/lambda) * (B_perp.*delta_h) ./ (R_near*sin(theta));
% 合成干涉相位（考虑相位缠绕）
total_phase = mod(phi_def + phi_top, 2*pi);
% 可视化设置
figure('Name',sprintf('Interferogram %02d',i),'Color','w');
imagesc(total_phase)
colormap(jet)
caxis([0 2*pi])
title(sprintf('干涉对%02d 相位图\n相对时间基线:%.1f年 相对垂直基线:%dm',...
i,abs(dt_year),B_perp))
colorbar('Ticks',0:pi/2:2*pi,...
'TickLabels',{'0','π/2','π','3π/2','2π'})
axis equal off
set(gca,'Position',[0.1 0.1 0.8 0.8])
% 添加PS网络叠加（可选）
hold on
plot(ps_x, ps_y, 'k.','MarkerSize',4)
for k = 1:length(row)
line([ps_x(row(k)),ps_x(col(k))],...
[ps_y(row(k)),ps_y(col(k))],...
'Color',[0.6 0.6 0.6],'LineWidth',0.3);
end
hold off
% 保存结果（可选）
% print(sprintf('Intf_%02d.png',i),'-dpng','-r300')
end
%% 续写：施加不同噪声级别的干涉相位图
% 设定目标基线参数
target_dt_year = 3.7; % 时间基线3.7年
target_B_perp = 1; % 垂直基线1m
% 计算相位分量（复用原始参数）
phi_def = (4*pi/lambda) * v * target_dt_year;
phi_top = (4*pi/lambda) * (target_B_perp * delta_h) ./ (R_near * sin(theta));
total_phase = mod(phi_def + phi_top, 2*pi); % 初始无噪相位
% 生成四幅噪声图
noise_sigmas = [0.00, 0.20, 0.50, 0.65]; % 噪声标准差序列
figure('Name','Noisy Interferograms','Color','w', 'Position', [100 100 1200 800])
for k = 1:4
% 添加相位噪声
sigma = noise_sigmas(k);
phase_noise = sigma * randn(size(total_phase)); % 生成高斯白噪声
noisy_phase = mod(total_phase + phase_noise, 2*pi);% 带噪相位
% 绘制子图
subplot(2,2,k);
imagesc(noisy_phase)
colormap(jet)
caxis([0 2*pi])
title(sprintf('σ=%.2f rad', sigma), 'FontWeight','bold')
colorbar('Ticks',0:pi/2:2*pi,...
'TickLabels',{'0','π/2','π','3π/2','2π'})
axis equal off
% 叠加PS网络
hold on
plot(ps_x, ps_y, 'k.', 'MarkerSize', 4)
for i = 1:length(row)
line([ps_x(row(i)),ps_x(col(i))],...
[ps_y(row(i)),ps_y(col(i))],...
'Color',[0.6 0.6 0.6],'LineWidth',0.3);
end
hold off
end
% 设置主标题
sgtitle(sprintf('相对时间基线%.1f年 相对垂直基线%.1fm 干涉相位噪声对比',...
target_dt_year, target_B_perp), 'FontSize',14, 'FontWeight','bold')续写代码，要求在各级噪声水平下0.00, 0.20, 0.50, 0.65，0.80，分别从时序干涉图中提取出PS点上的干涉相位。基于PS差分相位观测值，采用VLT算法建立网络连接边，得到各边上的复相关系数，其绝对大小可作为网络基线边稳健性的评价指标。 相关系数越大，基线边的观测相位与PS差分模型吻合程度越高，得到的PS沉降速度差和高程修正差越可靠。反之，基线估计精度越差。在不 同噪声级别下，对网络中所有基线的相关系数进行统计，得到其最大值 、 最小值 、 均值 、标准差，以此作为对不同噪声水平下网络基线连接稳健性的衡量指标。
最新发布
扫一扫
3万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
