MATLAB 利用SVD分解擬合三維平面

方法簡述：

利用svd分解求擬合平面 擬合方程為ax+by+cz=d，約束條件為a2+b2+c2=1，目標是使得盡可能多的點在平面上。構建矩陣為ax=0，將a矩陣進行奇異值分解之後，最小奇異值對應的特徵值向量即為擬合平面的係數向量。

詳細的原理可以檢視：svd分解平面擬合原理

在這裡我舉的例子是利用餘震的三維空間座標擬合餘震的可能發震平面。

as=load('data.txt');

%根據需要選擇合適經緯度範圍的資料點
as1 = ;
for ir = 1:m
if (as(ir,1)>= xx &&as(ir,1)<= xx)
as1 = [as1;as(ir,:)];
endendas=as1;
lon=as(:,2);
lat=as(:,1);
depth=as(:,3);
mag=as(:,4);
% 將經緯度轉化為utm座標（m）
% function [east, north, z] = ll2utm(lat,lon, hgt,z)
[east, north, z]=ll2utm(lat,lon); %z=43s
x=east/1000; %轉換單位為km
y=north/1000; %轉換單位為km
z=depth;
planedata=[x,y,z];
% 協方差矩陣的svd變換中，最小奇異值對應的奇異向量就是平面的方向
xyz0=mean(planedata,1);
centeredplane=bsxfun(@minus,planedata,xyz0);
[u,s,v]=svd(centeredplane);
% 對a進行奇異值分解，最小奇異值對應的特徵向量即為擬合平面的係數向量 
a=v(1,3);
b=v(2,3);
c=v(3,3);
d=-dot([a b c],xyz0);
% 圖形繪製
figure;
% subplot(1,2,1)
scatter3(x,y,depth,40,'filled','blue');
set(gca,'zdir','reverse')
grid on;
title(strcat(num2str(a),'x+',num2str(b),'y+',num2str(c),'z=',num2str(d)));
xlabel('east (km)');
ylabel('north (km)');
zlabel('depth (km)');
set(gca,'xlim',[650 740]);%x軸的資料顯示範圍
hold on;
xfit = min(x):0.1:max(x);
yfit = min(y):0.1:max(y);
[xfit,yfit]= meshgrid (xfit,yfit);
zfit = -(d + a * xfit + b * yfit)/c;
mesh(xfit,yfit,zfit,'edgecolor',[0.9290 0.6940 0.1250]); %平面的顏色可以自主選擇，用rgb三元組表示即可

% view([-90,0]);%左檢視

% view([0,90]);%下檢視

p.s., 如果想將所有點繪製乙個曲面，可以檢視我往期的博文：

matlab中如何實現mesh三維圖

基於matlab的矩陣奇異值（SVD）分解

目錄 1 計算原理 1.1求解v 1.2求解d 1.3求解u 2 matlab程式 2.1 注意 3 設矩陣a的大小m n m n a ud 首先求出的特徵值及特徵值的特徵值從大到小按順序排列，對應的正交特徵向量也要按特徵值的順序排列。排列好的正交的特徵向量即為v v1,v2 其中v1為非零特徵值的...

SVD分解原理詳解

在介紹svd之前，先補充一些基礎知識 1.酉矩陣 2.正規 正定 矩陣 3.譜分解 4.svd分解 作為譜定理的泛化，svd 分解對於原矩陣的要求就要弱得多。4.手動svd分解的乙個例項 svd的分解實際可以將矩陣 m寫成乙個求和形式 5.svd分解的應用 1 分析了解原矩陣的主要特徵和攜帶的資訊 ...

機器學習 SVD分解

3 svd的舉例 4 svd的應用 5 svd的優缺點 參考經常看到svd奇異值分解，但一直沒有去了解它講的什麼，剛好在李航老師統計學習方法第二版上是單獨的一章，下面看了一些部落格總結一下 任意乙個m n m nm n的矩陣 m mm階正交陣 降序排列的非負對角線元素組成的m n m nm n對角陣...