自己做的後方交會元素解算matlab源程式
廢話不多說，直接上**
**的原始資料都已經直接賦值了，可以再寫一點點**讀入excel或者txt
// a code block
%單幅影像後方交會
%auther:neverforget2020
%time:2020.9.21
close all;clear; clc;
disp('後方交會求解');
x =[-86.15,-53.40,-14.78,10.46];%%input('每次觀測的像平面內的x座標');
y =[-68.99,82.21,-76.63,64.43];%%input('每次觀測的像平面內的y座標');
x =[36589.41,37631.08,39100.97,40426.54];%%input('每次觀測的物空間內的x座標');
y =[25273.32,31324.51,24934.98,30319.81];%% input('每次觀測的物空間內的y座標');
z = [2195.17,728.69,2386.50,757.31];%%input('每次觀測的物空間內的z座標');
x=x/1000;y=y/1000;%修正單位
fk=153.24/1000;
m=50000;
limit1=0.001;limit2=0.0000001;%精度
xs=sum(x)/4;ys=sum(y)/4;zs=fk*m;%待定引數初始值
q=0;w=0;k=0;%外方位角
r = zeros(3,3);%旋轉矩陣
dx=ones(6,1);%改正值
h=fk*m;%航高
round=0;%迭代次數
c=[x(1);y(1);x(2);y(2);x(3);y(3);x(4);y(4)];
while(abs(dx(4,1))>limit2||abs(dx(5,1))>limit2||abs(dx(6,1))>limit2)
r(1,1)=cos(q)*cos(k)-sin(q)*sin(w)*cos(k);
r(1,2)=-cos(q)*sin(k)-sin(q)*sin(w)*cos(k);
r(1,3)=-sin(q)*cos(w);
r(2,1)=cos(w)*sin(k);
r(2,2)=cos(w)*cos(k);
r(2,3)=-sin(w);
r(3,1)=sin(q)*cos(k)+cos(q)*sin(w)*sin(k);
r(3,2)=-sin(q)*sin(k)+cos(q)*sin(w)*cos(k);
r(3,3)=cos(q)*cos(w);
xy_0 = zeros(8,1);%共線方程
a=zeros(8,6);
for i = 1:4
xbar=r(1,1)*(x(i)-xs)+r(2,1)*(y(i)-ys)+r(3,1)*(z(i)-zs);
ybar=r(1,2)*(x(i)-xs)+r(2,2)*(y(i)-ys)+r(3,2)*(z(i)-zs);
zbar=r(1,3)*(x(i)-xs)+r(2,3)*(y(i)-ys)+r(3,3)*(z(i)-zs);
xy_0(2*i-1,1) = -fk*xbar/zbar;
xy_0(2*i,1) = -fk*ybar/zbar;
a(2*i-1,1)=1/ zbar*(r(1,1)*fk+r(1,3)*xy_0(2*i-1,1));
a(2*i-1,2)=1/ zbar*(r(2,1)*fk+r(2,3)*xy_0(2*i-1,1));
a(2*i-1,3)=1/ zbar*(r(3,1)*fk+r(3,3)*xy_0(2*i-1,1));
a(2*i,1)=1/ zbar*(r(1,2)*fk+r(1,3)*xy_0(2*i,1));
a(2*i,2)=1/ zbar*(r(2,2)*fk+r(2,3)*xy_0(2*i,1));
a(2*i,3)=1/ zbar*(r(3,2)*fk+r(3,3)*xy_0(2*i,1));
a(2*i-1,4)=xy_0(2*i,1)*sin(w)-(xy_0(2*i-1,1)/fk*(xy_0(2*i-1,1)*cos(k)-xy_0(2*i,1)*sin(k))+fk*cos(k))*cos(w);
a(2*i-1,5)=-fk*sin(k)-xy_0(2*i-1,1)/fk*(xy_0(2*i-1,1)*sin(k)+xy_0(2*i,1)*cos(k));
a(2*i-1,6)=xy_0(2*i,1);
a(2*i,4)=-xy_0(2*i-1,1)*sin(w)-(xy_0(2*i,1)/fk*(xy_0(2*i-1,1)*cos (k)-xy_0(2*i,1)*sin(k))-fk*sin(k))*cos(w);
a(2*i,5)=-fk*cos(k)-xy_0(2*i,1)/fk*(xy_0(2*i-1,1)*sin(k)+xy_0(2*i,1)*cos (k));
a(2*i,6)=-xy_0 (i*2-1,1);
endround=round+1;
l=c-xy_0;
dx=inv((a')*a)*(a')*l;
xs=xs+dx(1,1);ys=ys+dx(2,1);zs=zs+dx(3,1);
q=q+dx(4,1);w=w+dx(5,1);k=k+dx(6,1);
endformat long g;
disp("r陣")
rdisp("xs")
xs,ys,zs,q,w,k
disp("單位權m0")
m0 = 1000
disp("qxx")
qxx = inv((a')*a)*1000
xs =
39795.472089565627476.50608304877572.69914930955-0.003991094446802740.00210817362230499-0.0675785372578092攝影測量學後方交會
值 r i 1 string u console.readline for int j 0 j n j 2 2.像點y座標的輸入 for int i 0 i n i 值 r i 1 string v console.readline for int j 1 j n j 2 控制點的座標輸入 doub...
德語區國家的攝影測量與遙感
德語區國家 主要包括德國 瑞士和奧地利 在攝影測量與遙感專業方向的研究工作與工業產品是無可爭議的世界最強。這裡列舉幾個比較牛的。工業界 先說公司吧，在德國，有專門飛航片和做航空遙感影像資料處理的hansa luftbild公司，這是一家德國國有大型企業，也是現在的德國漢莎航空公司的祖宗公司。另外還有...
基於halcon的圓環寬度測量
1 建立測量模型 create metrology model metrologyhandle 2 設定測量模型影象尺寸 set metrology model image size metrologyhandle,width,height 3 增加測量模型物件 add metrology obje...