摄影测量立体相对的前方交会VB程序代码

合集下载

摄影测量-空间前交、后交【精选文档】

空间后交—前交程序设计(实验报告)姓名：班级：学号：时间:空间后交-前交程序设计一、实验目的用 C 、VB或MATLAB语言编写空间后方交会-空间前方交会程序⑴提交实习报告：程序框图、程序源代码、计算结果、体会⑵计算结果：像点坐标、地面坐标、单位权中误差、外方位元素及其精度二、实验数据f=150。

000mm，x0=0，y0=0三、实验思路1。

利用空间后方交会求左右像片的外方位元素(1）.获取m（于像片中选取两点，于地面摄影测量坐标系中选取同点,分别计算距离，距离比值即为m)，x，y，f,X，Y,Z（2).确定未知数初始值Xs,Ys，Zs,q，w,k（3）.计算旋转矩阵R(4）.逐点计算像点坐标的近似值(x），（y）(5）。

组成误差方程式(6)。

组成法方程式（7).解求外方位元素（8）。

检查是否收敛，即将求得的外方位元素的改正数与规定限差比较，小于限差即终止；否则用新的近似值重复步骤(3）-（7）2。

利用求出的外方位元素进行空间前交，求出待定点地面坐标(1）.用各自像片的角元素计算出左、右像片的方向余弦值，组成旋转矩阵R1，R2(2）。

根据左、右像片的外方位元素，计算摄影基线分量Bx，By，Bz(3)。

计算像点的像空间辅助坐标（X1,Y1，Z1）和（X2，Y2,Z2）（4）.计算点投影系数N1和N2(5）。

计算未知点的地面摄影测量坐标四、实验过程⑴程序框图函数AandL%求间接平差时需要的系数％%％已知％a=像点坐标x，b=像点坐标y,f内方位元素主距%φ=q，ψ=w，κ=k%像空间坐标系X,Y，Z%地面摄影测量坐标系Xs,Ys,Zsfunction ［A1，L1，A2，L2］=AandL（a，b，f,q,w,k,X，Y，Z,Xs,Ys,Zs) %％％％％%％％%%％选择矩阵元素a1=cos（q)＊cos(k)—sin(q）＊sin(w）*sin(k）;a2=-cos(q)＊sin(k)—sin(q）＊sin(w）＊cos（k)；a3=-sin(q）*cos(w);b1=cos(w）*sin(k)；b2=cos（w）*cos（k)；b3=—sin(w);c1=sin（q)*cos(k)+cos（q)＊sin(w）*sin（k)；c2=—sin（q）＊sin(k)+cos（q)*sin(w）*cos（k)；c3=cos(q)＊cos(w);％%％%％％％共线方程的分子分母X_=a1*(X—Xs）+b1*（Y-Ys)+c1＊（Z-Zs);Y_=a2*(X-Xs)+b2＊（Y—Ys)+c2＊（Z-Zs);Z_=a3＊(X—Xs)+b3＊（Y—Ys）+c3*(Z-Zs）;％％％％%％%近似值x=-f＊X_/Z_；y=-f＊Y_/Z_;％％％％%％％A组成L组成a11=1/Z_＊（a1＊f+a3＊x)；a12=1/Z_＊(b1＊f+b3＊x）；a13=1/Z_*(c1*f+c3＊x);a21=1/Z_*(a2＊f+a3＊y）;a22=1/Z_*（b2＊f+b3*y);a23=1/Z_＊(c2＊f+c3＊y）;a14=y*sin（w)-（x/f＊(x*cos(k)—y*sin（k））+f＊cos（k）)*cos(w）;a15=-f*sin(k)—x/f＊（x*sin（k）+y*cos(k));a16=y；a24=—x*sin(w）-（y/f＊（x*cos(k)-y*sin(k)）—f＊sin(k））*cos（w）；a25=-f*cos（k）-y/f*(x*sin（k）+y*cos(k));a26=-x;lx=a—x；ly=b-y;％%%%%%％％%组成一个矩阵，并返回A1=[a11,a12，a13,a14,a15，a16];A2=［a21,a22,a23，a24,a25，a26];L1=lx;L2=ly;函数deg2dms％%％%％%％％角度转度分秒function y=deg2dms(x）a=floor(x)；b=floor(（x-a)＊60);c=(x-a—b/60）＊3600；y=a+（b/100）+(c/10000);函数dms2deg%%％％%度分秒转度function y=dms2deg（x）a=floor(x）；b=floor((x-a)＊100）；c=（x-a—b/100)*10000；y=a+b/60+c/3600;函数ok%％％％％%%％%%%％％％目的是为了保证各取的值的有效值%％xy为n*1，a为1＊nfunction result=ok(xy，a)format short gi=size（xy，1)；for n=1:io=xy（n)—floor(xy(n,1));o=round（o＊（10^a（n))）/(10^a(n）);xy(n，1)=floor(xy（n,1）)+o;endformat long gresult=xy;函数rad2dmsxy%%％%求度分秒表现形式的三个外方位元素，三个角度function xydms=rad2dmsxy(xy)［a,b,c，d，e,f]=testvar(xy）;d=deg2dms（rad2deg(d））；e=deg2dms（rad2deg(e))；f=deg2dms(rad2deg(f）);xydms=[a,b，c,d,e，f］';函数spacehoujiao％％%％%%%空间后交％%% f%％输入p（2*n,1）%％像点坐标x,y,X,Y，Z,均为（n,1）function ［xy,m，R］=spacehoujiao(p，x,y，f,X，Y,Z）format long；%%％%％权的矢量化,这是等精度时的，如果非,将函数参数改为PP=diag（p)；％％求nj=size（X，2）；%％初始化Xs=0；Ys=0;Zs=0;for n=1：jXs=Xs+X（n）；Ys=Ys+Y(n）;Zs=Zs+Z(n）;endSx=sqrt((x（2）-x(1）)^2+(y（2)—y(1)）^2);％%％%两像点之间距离Sd=sqrt(（X(2）-X（1))^2+(Y（2）-Y（1))^2);%％%％两地面控制点之间距离m=Sd/Sx; ％％％％图像比例系数Xs=Xs/j;Ys=Ys/j;Zs=m*f+Zs/j；m0=0;q=0;w=0;k=0;i=0;a=rand(2*j,6）；l=rand（2＊j，1);％％％％for n=1:j[a（2＊n—1,:)，l(2*n—1,1）,a(2＊n,：），l(2*n,1）］=AandL(x(n)，y（n),f,q，w，k，X(n)，Y（n），Z(n），Xs,Ys,Zs);enddet=inv（a’＊P*a）*transpose（a)*P*l;％%%%%%％％％循环体while 1％％％％％%%％%％％％％%％％[dXs,dYs,dZs，dq,dw,dk]=testvar(det）;detXs=abs（dXs）；detYs=abs（dYs)；detZs=abs(dZs)；detq=abs(dq)；detw=abs（dw)；detk=abs（dk)；%%％％%%％%%if （（detXs<0。

摄影测量程序汇总(后方交会+前方交会+单模型光束法平差).doc

摄影测量程序汇总(后方交会+前方交会+单模型光束法平差) .程序运行环境为Visual Studio2010.运行前请先将坐标数据放在debug 下。

1. 单像空间后方交会原始数据：内方位元素x0(/mm)y0(/mm)主距f(/mm)比例尺分母00153.2450000像点坐标(/mm)物点坐标(/m)-内方位元素x0(/mm)y0(/mm)主距f(/mm)比例尺分母00153.2450000像点坐标(/mm)物点坐标(/m)：#include #include #include double *readdata();void savedata(int hang,double *data,double *xishuarray,double *faxishu,double *l,int i,double xs,double ys,double zs,double fai,double oumiga,double kapa);void transpose(double *m1,double *m2,int m,int n);void inverse(double *a,int n);void multi(double *mat1,double * mat2,double * result,int a,int b,int c);void inverse(double *a,int n)/*正定矩阵求逆*/{ int i,j,k; for(k=0;k=6.0/206265.0||xcorrect[5][0]=6.0/206265.0||xcorrect[9][0]=6.0/ 206265.0||xcorrect[10][0]=6.0/206265.0||xcorrect[11][0]=6.0/206265.0); savedata(row,controlpoint,A[0],ATA[0],L[0],diedainumber,Xs1,Ys1,Zs1,f ai1,oumiga1,kapa1,Xs2,Ys2,Zs2,fai2,oumiga2,kapa2); printf("迭代次数为：%d\n\n",diedainumber); printf("左像片的外方位元素为：\n"); printf("xs1=%lf,ys1=%lf,zs1=%lf\n\n",Xs1,Ys1,Zs1); printf("fai1=%lf,oumiga1=%lf,kapa1=%lf\n\n",fai1,oumiga1,kapa1); printf("右像片的外方位元素为：\n"); printf("xs2=%lf,ys2=%lf,zs2=%lf\n\n",Xs2,Ys2,Zs2);printf("fai2=%lf,oumiga2=%lf,kapa2=%lf\n",fai2,oumiga2,kapa2); system("pause");}测试结果：如果有正确的数据，请代入进行验证。

数字摄影测量坐标转换(vb代码)

实验报告1 摄影测量坐标变换目的实现摄影测量常用坐标系之间的变换要求用VB进行编程，主要实现像空间坐标系与像辅助坐标系之间的变换方法与详细步骤界面所入空间坐标输入三个角度旋转（度分秒形式）计算旋转变换矩阵得到该点在像空间辅助坐标系中的坐标实验成果1.VB/VC原始代码Const PI As Double = 3.14159265Private Sub Command1_Click()Rem 定义数据类型Dim x1 As Double, x2 As Double, y1 As Double, y2 As Double, z1 As Double, z2 As DoubleDim ω As Single, κ As Single, φ As SingleDim a1 As Single, a2 As Single, a3 As SingleDim b1 As Single, b2 As Single, b3 As SingleDim c1 As Single, c2 As Single, c3 As SingleRem 从text中读取数据x1 = Val(txtx1.Text): y1 = Val(txty1.Text): z1 = Val(txtz1.Text)ω= deg(Val(Txtω.Text)): κ = deg(Val(Txtκ.Text)): φ = deg(Val(Txtφ.Text)) Debug.Print φ, κ, ωRem 求解a1, a2, a3 ，b1, b2, b3，c1, c2, c3a1 = Cos(φ * PI / 180) * Cos(κ * PI / 180) - Sin(ω * PI / 180) * Sin(κ * PI / 180) * Sin(φ * PI / 180)a2 = -Cos(φ * PI / 180) * Sin(κ * PI / 180) - Sin(φ * PI / 180) * Sin(ω * PI / 180) * Cos(κ * PI / 180)a3 = -Sin(φ * PI / 180) * Cos(ω * PI / 180)b1 = Cos(ω * PI / 180) * Sin(κ * PI / 180)b2 = Cos(ω * PI / 180) * Cos(κ * PI / 180)b3 = -Sin(ω * PI / 180)c1 = Sin(φ * PI / 180) * Cos(κ * PI / 180) + Sin(ω * PI / 180) * Cos(φ * PI / 180) * Sin(κ * PI / 180)c2 = -Sin(φ * PI / 180) * Sin(κ * PI / 180) + Sin(ω * PI / 180) * Cos(φ * PI / 180) * Cos(κ * PI / 180)c3 = Cos(ω * PI / 180) * Cos(φ * PI / 180)Debug.Print a1, a2, a3Debug.Print b1, b2, b3Debug.Print b1, b1, b3Rem 定义数组Dim a(3, 3) As Singlea(1, 1) = a1: a(1, 2) = a2: a(1, 3) = a3a(2, 1) = b1: a(2, 2) = b2: a(2, 3) = b3a(3, 1) = c1: a(3, 2) = c2: a(3, 3) = c3Dim x(3, 1) As Doublex(1, 1) = x1x(2, 1) = y1x(3, 1) = z1Rem 求解Dim y(3, 1) As DoubleFor i = 1 To 3y(i, 1) = a(i, 1) * x(1, 1) + a(i, 2) * x(2, 1) + a(i, 3) * x(3, 1)Next ix2 = y(1, 1)y2 = y(2, 1)z2 = y(3, 1)Debug.Print x2, y2, z2Txtx2.Text = Format(x2, "0.000")Txty2.Text = Format(y2, "0.000")Txtz2.Text = Format(z2, "0.000")End SubRem 定义deg函数，即度分秒转换为度Private Function deg(a As Double)sign = Sgn(a)a = Abs(a) + 0.0000000001b = Int(a)c = Int((a - b) * 100)d = a - b - c / 100deg = sign * (b + c / 60 + d / 0.36)End Function➢ 2.程序运行结果图。

前方交会实验报告(含VB程序代码)

立体像对前方交会实验报告一、实验目的在掌握前方交会原理的基础上，自己编写前方交会程序，在计算机上调试，输出计算结果并对计算结果进行检验。

通过上机调试程序加强动手能力的培养，通过对实验过程的掌握以及对实验结果的分析，增强学生综合运用所学知识解决实际问题的能力。

二、实验仪器计算机，VB6.0三、实验数据1.模拟像片一对：左片号23，右片号24；2.航摄机主距：f=150mm；3.左片23号片外方位元素：φ=−0°25′00″ω=−1°00′00″k=−0°10′00″Xs=103007.006117 Ys=139998.994849 Zs=4801.9989994 （m）右片24号片外方为元素：φ=1°39′59″ω=−0°10′00″k=0°40′00″Xs=106002.023762 Ys=140005.002780 Zs=4797.009648 （m）待求像点坐标如下表：四、实验内容利用所给立体像对两张像片的内、外方位元素，编写空间前方交会程序，根据所给像对中若干同名像点在左右像片上的坐标，解求其对应的地面点的物方坐标，实现空间前方交会的过程。

五、实验成果程序流程图：程序设计界面：程序运行界面：运行结果：（注：表上显示地面点坐标依次是：7,9,4,6,5,8）附：excel进行角度转换：六、程序如下：Dim f#, x1#, y1#, x2#, y2#, i%, j%, u1#, u2#, v1#, v2#, w1#, w2#, fai1#, kab1#, omg1#, fai2#, kab2#, omg2# Dim a12#, a13#, b11#, b12#, b13#, c11#, c12#, c13#, a21#, a22#, a23#, b21#, b22#, b23#, c21#, c22#, c23# Dim n1#, n2#, bu#, bv#, bw#Dim xs1#, xs2#, ys1#, ys2#, zs1#, zs2#Private Sub Command1_Click()fai1 = Val(Text1.Text): kab1 = Val(Text3.Text): omg1 = Val(Text2.Text)fai2 = Val(Text11.Text): kab2 = Val(Text10.Text): omg2 = Val(Text7.Text)f = Val(Text13.Text)xs1 = Val(Text4.Text): xs2 = Val(Text9.Text): ys1 = Val(Text5.Text): ys2 = Val(Text8.Text): zs1 = Val(Text6.Text): zs2 = Val(Text12.Text)End SubPrivate Sub Command3_Click()x1 = Val(InputBox("输入23片坐标x的值"))y1 = Val(InputBox("输入23片坐标y的值"))x2 = Val(InputBox("输入24片坐标x的值"))y2 = Val(InputBox("输入24片坐标y的值"))End SubPrivate Sub Command2_Click()Text1.Visible = FalseText2.Visible = FalseText3.Visible = FalseText4.Visible = FalseText5.Visible = FalseText6.Visible = FalseText7.Visible = FalseText8.Visible = FalseText9.Visible = FalseText10.Visible = FalseText11.Visible = FalseText12.Visible = FalseText13.Visible = FalseLabel1.Visible = FalseLabel2.Visible = FalseLabel3.Visible = FalseLabel4.Visible = FalseLabel5.Visible = FalseLabel6.Visible = FalseLabel7.Visible = FalseLabel8.Visible = FalseLabel9.Visible = FalseLabel10.Visible = FalseLabel11.Visible = FalseLabel12.Visible = FalseLabel13.Visible = FalseLabel14.Visible = Falsea11 = Cos(fai1) * Cos(kab1) - Sin(fai1) * Sin(omg1) * Sin(kab1)a12 = -Cos(fai1) * Sin(kab1) - Sin(fai1) * Sin(omg1) * Cos(kab1)a13 = -Sin(fai1) * Cos(omg1)b11 = Cos(omg1) * Sin(kab1)b12 = Cos(omg1) * Cos(kab1)b13 = -Sin(omg1)c11 = Sin(fai1) * Cos(kab1) + Cos(fai1) * Sin(omg1) * Sin(kab1)c12 = -Sin(fai1) * Sin(kab1) + Cos(fai1) * Sin(omg1) * Cos(kab1)c13 = Cos(fai1) * Cos(omg1)a21 = Cos(fai2) * Cos(kab2) - Sin(fai2) * Sin(omg2) * Sin(kab2)a22 = -Cos(fai2) * Sin(kab2) - Sin(fai2) * Sin(omg2) * Cos(kab2)a23 = -Sin(fai2) * Cos(omg2)b21 = Cos(omg2) * Sin(kab2)b22 = Cos(omg2) * Cos(kab2)b23 = -Sin(omg2)c21 = Sin(fai2) * Cos(kab2) + Cos(fai2) * Sin(omg2) * Sin(kab2)c22 = -Sin(fai2) * Sin(kab2) + Cos(fai2) * Sin(omg2) * Cos(kab2)c23 = Cos(fai2) * Cos(omg2)u1 = a11 * x1 + a12 * y1 - a13 * fv1 = b11 * x1 + b12 * y1 - b13 * fw1 = c11 * x1 + c12 * y1 - c13 * fu2 = a21 * x2 + a22 * y2 - a23 * fv2 = b21 * x2 + b22 * y2 - b23 * fw2 = c21 * x2 + c22 * y2 - c23 * fbu = xs2 - xs1bv = ys2 - ys1bw = zs2 - zs1n1 = (bu * w2 - bw * u2) / (u1 * w2 - u2 * w1)n2 = (bu * w1 - bw * u1) / (u1 * w2 - u2 * w1)Dim xx#, yy#, yy1#, yy2#, zz#xx = Fix((xs1 + u1 * n1) * 1000 + 0.5) / 1000yy1 = Val(Text5.Text) + v1 * n1yy2 = Val(Text8.Text) + v2 * n2yy = Fix(((yy1 + yy2) / 2) * 1000 + 0.5) / 1000zz = Fix((Val(Text6.Text) + w1 * n1) * 1000 + 0.5) / 1000Form1.Print "地面坐标为:"Form1.Print "xx="; xx, "yy="; yy; "zz="; zzEnd Sub七、实验心得通过本次前方交会实验，熟悉掌握了前方交会原理，并利用计算机程序对其进行了解决，收益颇多！。

立体像对前方交会绝对定向

间前方 Z 交会方法 Z
1
2
Y2
Y1 s1
Z1 X1
X2
X1
Zs1 Y1
Ztp
Ytp Xs1 M
Ys1 Xtp
摄影基线
s2
B
BZ= Zs2 –Zs1 BY= Ys2 –Ys1
s1
BX= Xs2 –Xs1
同名光线投影
s1
Z1
S1 A S 1 a1

X A X s1 X1
解析绝对定向误差方程设
1,
0

0

0

0
0
则
X 0 Y 0 Y Z 0 l X l X Y 0 l Z
v X v Y vZ
l1 X l 2 Y l 3 Z l x 0 l 4 X l 5Y l 6 Z l y 0
一对同名像点
4个线性方程
3个未知数X、Y、Z 2n个线性方程
若n幅影像中含有同一个空间点
这是一种严格的、不受影像数约束的空间前方交会方法、且不需要空间坐标的初始值
严密解法

Z
p
Z
p
Z
pg
X
tpg

1
n
X n
tpi
X
tp
X
tp
X
•
tpg
Y tpg

1
n
Y tpi n
Y tp Y tp Y tpg Z tp Z tp Z tpg
Z tpg

1
n
Z tpi n
解析绝对定向误差方程设

VB代码

CDg1.Filter = "Text Files(*.TXT)|*.txt|All Files(*.*)|*.*"
CDg1.DialogTitle = "读取已知数据"
CDg1.FileName = "": CDg1.Action = 1
'显示读入的控制点地面坐标
txtShow.Text = txtShow.Text & Xt(1) & " , " & Yt(1) & " , " & Zt(1) & vbCrLf
txtShow.Text = txtShow.Text & Xt(2) & " , " & Yt(2) & " , " & Zt(2) & vbCrLf
Dim fai_R#, omg_R#, kap_R#, XsR#, YsR#, ZsR# '左片外方位元素
Dim Bx#, By#, Bz# '基线分量
Dim R_L#(1 To 3, 1 To 3), R_R#(1 To 3, 1 To 3) '左右像片的旋转矩阵
If CDg1.FileName = "" Then Exit Sub
Open CDg1.FileName For Input As #1
Line Input #1, strTemp '读第一行题头信息
txtShow.Text = txtShow.Text & vbCrLf & strTemp & vbCrLf

06 立体像对前方交会

X a1 Y = a 2 Z a3
b1 b2 b3
c1 X − X s X − X s c2 Y − Ys = R −1 Y − Ys Z −Z c3 Z − Z s s
误差方程
共线条件方程
a1 ( X − X s ) + b1 (Y − Ys ) + c1 (Z − Z s ) X x − x0 = − f =−f a3 ( X − X s ) + b3 (Y − Ys ) + c3 (Z − Z s ) Z a2 ( X − X s ) + b2 (Y − Ys ) + c2 (Z − Z s ) Y y − y0 = − f =−f a3 ( X − X s ) + b3 (Y − Ys ) + c3 (Z − Z s ) Z
A(X,Y,Z) Y X
二、基本公式 1、严密解法、
已知值 x0 , y0 , f , m观测值 x,y 未知数 X, Y, Z 泰勒级数展开
∂x ∂x ∂x vx = ∆X + ∆Y + ∆Z + x0 − x ∂X ∂Y ∂Z ∂y ∂y ∂y vy = ∆X + ∆Y + ∆Z + y 0 − y ∂X ∂Y ∂Z
f sin κ H f a12 = + cos κ H y a13 = + H a 21 = −
Z2
二、基本公式
2、点投影系数法、
Z1
Y2
Y1 s1
Z1 X1
X2
X1
Zs1 Y1 Ztp Ytp Xs1 M
Ys1 Xtp
摄影基线

摄影测量实验报告（空间后方交会—前方交会）

摄影测量实验报告（空间后⽅交会—前⽅交会）空间后⽅交会-空间前⽅交会程序编程实验⼀．实验⽬的要求掌握运⽤空间后⽅交会-空间前⽅交会求解地⾯点的空间位置。

学会运⽤空间后⽅交会的原理，根据所给控制点的地⾯摄影测量坐标系坐标以及相应的像平⾯坐标系中的坐标，利⽤计算机编程语⾔实现空间后⽅交会的过程，完成所给像对中两张像⽚各⾃的外⽅位元素的求解。

然后根据空间后⽅交会所得的两张像⽚的内外⽅位元素，利⽤同名像点在左右像⽚上的坐标，求解其对应的地⾯点在摄影测量坐标系中的坐标，并完成精度评定过程，利⽤计算机编程语⾔实现此过程。

⼆．仪器⽤具计算机、编程软件（MATLAB）三．实验数据实验数据包含四个地⾯控制点（GCP）的地⾯摄影测量坐标及在左右像⽚中的像平⾯坐标。

此四对坐标运⽤最⼩⼆乘法求解左右像⽚的外⽅位元素，即完成了空间后⽅的过程。

另外还给出了5对地⾯点在左右像⽚中的像平⾯坐标和左右像⽚的内⽅位元素。

实验数据如下：内⽅位元素：f=152.000mm，x0=0，y0=0四．实验框图此过程完成空间后⽅交会求解像⽚的外⽅位元素，其中改正数⼩于限差（0.00003，相当于0.1’的⾓度值）为⽌。

在这个过程中采⽤迭代的⽅法，是外⽅位元素逐渐收敛于理论值，每次迭代所得的改正数都应加到上⼀次的初始值之中。

确定Xs，Ys，Zs的初始值时，对于左⽚可取地⾯左边两个GCP的坐标的平均值作为左⽚Xs 和Ys的初始值，取右边两个GCP 的坐标平均值作为右⽚Xs 和Ys的初始值。

Zs可取地⾯所有GCP的Z坐标的平均值再加上航⾼。

空间前⽅交会的数学模型为：五．实验源代码function Main_KJQHFJH()global R g1 g2 m G a c b1 b2;m=10000;a=5;c=4;feval(@shuru); %调⽤shuru（）shurujcp（）函数完成像点及feval(@shurujcp); %CCP有关数据的输⼊XYZ=feval(@MQZqianfangjh); %调⽤MQZqianfangjh（）函数完成空间前⽅、%%%%%% 单位权中误差%%%% %后⽅交会计算解得外⽅位元素global V1 V2; %由于以上三个函数定义在外部⽂件中故需VV=[]; %⽤feval（）完成调⽤过程for i=1:2*cVV(i)=V1(i);VV(2*i+1)=V2(i);endm0=sqrt(VV*(VV')/(2*c-6));输⼊GCP像点坐标及地⾯摄影测量坐标系坐标的函数和输⼊所求点像点坐标函数：function shurujcp()global c m;m=input('摄影⽐例尺：'); %输⼊GCP像点坐标数据函数并分别将其c=input('GCP的总数='); % 存⼊到不同的矩阵之中disp('GCP左⽚像框标坐标：');global g1;g1=zeros(c,2);i=1;while i<=cm=input('x=');n=input('y=');g1(i,1)=m;g1(i,2)=n;i=i+1;enddisp('GCP右⽚像框标坐标：');global g2;g2=zeros(c,2);i=1;while i<=cm=input('x=');n=input('y=');g2(i,1)=m;g2(i,2)=n;i=i+1;end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function shuru()global a;a=input('计算总像对点数='); %完成想计算所需的像平⾯坐标global b1; %坐标输⼊，存⼊不同的矩阵中b1=zeros(a,2); disp('左⽚像点坐标：')i=1;while i<=am=input('x=');n=input('y=');b1(i,1)=m;b1(i,2)=n;i=i+1;end%%b2=zeros(a,2);disp('右⽚像点坐标：')i=1;while i<=am=input('x=');n=input('y=');b2(i,1)=m;b2(i,2)=n;i=i+1;end%%global c;c=input('GCP的总数=');disp('GCP摄影测量系坐标：')global G;G=zeros(3,c);i=1;while i<=cm=input('X=');n=input('Y=');v=input('Z=');G(i,1)=m;G(i,2)=n;G(i,3)=v;i=i+1;end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%空间前⽅交会和后⽅交会函数：function XYZ=MQZqianfangjh()global R1 R2 a f b1 b2 Ra Rb;global X1 X2;R1=Ra;R2=Rb;R1=zeros(3,3);R2=zeros(3,3);global g1 g2 G V1 V2 V WF c QXX QXX1 QXX2;xs0=(G(1,1)+G(3,1))/2;ys0=(G(1,2)+G(3,2))/2;[Xs1,Ys1,Zs1,q1,w1,k1 R]=houfangjh(g1,xs0,ys0); %对左⽚调⽤后⽅交会函数R1=R;V1=V;WF1=WF;save '左⽚外⽅位元素为.txt' WF -ascii %将计算所得的外⽅位元素存⼊到.txt% ⽂件中for i=1:cg1(i,1)=g1(i,1)+V1(2*i-1);g1(i,2)=g1(i,2)+V1(2*i);endsave '左⽚像点坐标.txt' g1 -asciixs0=(G(2,1)+G(4,1))/2;ys0=(G(2,2)+G(4,2))/2;[Xs2,Ys2,Zs2,q2,w2,k2 R]=houfangjh(g2,xs0,ys0); %对右⽚调⽤后⽅交会函数R2=R; V2=V;WF2=WF;QXX2=QXX;save '右⽚外⽅位元素为.txt' WF –ascii %将计算所得的外⽅位元素存⼊到.txt% ⽂件中for i=1:cg2(i,1)=g2(i,1)+V2(2*i-1);g2(i,2)=g2(i,2)+V2(2*i);endsave '右⽚像点坐标.txt' g2 -asciiX1=zeros(a,3);X2=zeros(a,3);xx=zeros(3,1);xxx=zeros(3,1);for i=1:ass=[b1(i,1);b1(i,2);-f];dd=[b2(i,1);b2(i,2);-f];xx=R1*ss;X1(i,:)=xx';xxx=R2*dd;X2(i,:)=xxx';endglobal Xs1 Xs2 Ys1 Ys2 Zs1 Zs2;BX=Xs2-Xs1;BY=Ys2-Ys1;BZ=Zs2-Zs1;global N1 N2;N1=zeros(1,a);N2=zeros(1,a);for i=1:aN1(1,i)=(BX*X2(i,3)-BZ*X2(i,1))/(X1(i,1)*X2(i,3)-X2(i,1)*X1(i,3));N2(1,i)=(BX*X1(i,3)-BZ*X1(i,1))/(X1(i,1)*X2(i,3)-X2(i,1)*X1(i,3));end %计算投影系数，并计算五点的三维坐标global XYZ;XYZ=zeros(a,3);for i=1:aXYZ(i,1)=Xs1+N1(1,i)*X1(i,1);XYZ(i,2)=((Ys1+N1(1,i)*X1(i,2))+(Ys2+N2(1,i)*X2(i,2)))/2;enddisp('左⽚外⽅位元素为:Xs Ys Zs ψωκ');disp(WF1);disp('左⽚外⽅位元素协因素阵为：');disp(QXX1);disp('左⽚像点坐标为：')disp(g1)disp('右⽚外⽅位元素为:Xs Ys Zs ψωκ');disp(WF2);disp('右⽚外⽅位元素协因素阵为：')disp(QXX2)disp('右⽚像点坐标为：')disp(g2)disp('计算所得点摄影测量坐标（X,Y,Z）为：');disp(XYZ);save 'XYZ.txt' XYZ -ascii %将计算所得结果保存到XYZ.txt⽂件中%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [Xs,Ys,Zs,q,w,k R]=houfangjh(g1,Xs0,Ys0) %计算像⽚外⽅位元素%%%%%%%%%%global f G m c b1 b2;f=0.152;Xs=Xs0;Ys=Ys0;Zs=m*f+G(1,3);q=0;w=0;k=0;while 1 %实现⼀个永真循环，是改正数⼩于限差以后跳出循环a1=cos(q)*cos(k)-sin(q)*sin(w)*sin(k);a2=-cos(q)*sin(k)-sin(q)*sin(w)*cos(k);a3=-sin(q)*cos(w);b1_=cos(w)*sin(k);b2_=cos(w)*cos(k);b3=-sin(w);c1=sin(q)*cos(k)+cos(q)*sin(w)*sin(k);c2=-sin(q)*sin(k)+cos(q)*sin(w)*cos(k);c3=cos(q)*cos(w);R=[a1,a2,a3;b1_,b2_,b3;c1,c2,c3];for i=1:caX(i)=a1*(G(i,1)-Xs)+b1_*(G(i,2)-Ys)+c1*(G(i,3)-Zs);aY(i)=a2*(G(i,1)-Xs)+b2_*(G(i,2)-Ys)+c2*(G(i,3)-Zs);aZ(i)=a3*(G(i,1)-Xs)+b3*(G(i,2)-Ys)+c3*(G(i,3)-Zs);endxj=[];yj=[];for i=1:cxj(i)=-f*aX(i)/aZ(i);yj(i)=-f*aY(i)/aZ(i);enda11=[];a12=[];a13=[];a14=[];a15=[];a16=[];a21=[];a22=[];a23=[];a24=[];a25=[];a26=[];for i=1:ca11(i)=(a1*f+a3*g1(i,1))/aZ(i);a12(i)=(b1_*f+b3*g1(i,1))/aZ(i);a13(i)=(c1*f+c3*g1(i,1) )/aZ(i);a21(i)=(a2*f+a3*g1(i,2))/aZ(i);a22(i)=(b2_*f+b3*g1(i,2))/aZ(i);a23(i)=(c2*f+c3*g1(i,2) )/aZ(i);a14(i)=g1(i,2)*sin(w)-(g1(i,1)*(g1(i,1)*cos(k)-g1(i,2)*sin(k))/f+f*cos(k))*cos(w);a15(i)=-f*sin(k)-g1(i,1)*(g1(i,1)*sin(k)+g1(i,2)*cos(k))/f;a16(i)=g1(i,2);a24(i)=-g1(i,1)*sin(w)-(g1(i,2)*(g1(i,1)*cos(k)-g1(i,2)*sin(k))/f-f*sin(k))*cos(w);a25(i)=-f*cos(k)-g1(i,2)*(g1(i,1)*sin(k)+g1(i,2)*cos(k))/f;a26(i)=-g1(i,1);endlx=[];ly=[];for i=1:clx(i)=g1(i,1)-xj(i);ly(i)=g1(i,2)-yj(i);endA=zeros(2*c,6);for i=1:cA(2*i-1,1)=a11(i);A(2*i-1,2)=a12(i);A(2*i-1,3)=a13(i);A(2*i-1,4)=a14(i);A(2*i-1,5)=a15 (i);A(2*i-1,6)=a16(i); A(2*i,1)=a21(i); A(2*i,2)=a22(i); A(2*i,3)=a23(i); A(2*i,4)=a24(i); A(2*i,5)=a25(i); A(2*i,6)=a26(i);endL=zeros(2*c,1);for i=1:cL(2*i-1,1)=lx(i);endX=inv((A')*A)*(A')*L;Xs=Xs+X(1,1);Ys=Ys+X(2,1);Zs=Zs+X(3,1);q=q+X(4,1);w=w+X(5,1);k=k+X(6,1);Xabs=abs(X);aaa=max(Xabs);if aaa<0.00003 %当改正数中绝对值最⼤的改正数⼩于限差0.00003 break; %后跳出循环，计算结果已经收敛endendglobal V;V=L';global WF QXX;WF(1)=Xs;WF(2)=Ys;WF(3)=Zs;WF(4)=q;WF(5)=w;WF(6)=k;QXX=A'*A;六．实验结果左⽚外⽅位元素Xs,Ys,Zs,ψ、ω、κ、为：5.0001950e+003 5.0007250e+003 2.0201583e+003 -7.2888190e-005 2.8193877e-002 9.5130388e-002左⽚外⽅位元素协因素阵为：4.0166895e-008 -3.7263703e-010 1.3218695e-008 7.0720033e-005 1.0001730e-007 -2.5748604e-006-3.7263703e-010 4.0032797e-008 2.6568407e-009 -2.1103715e-007 7.7772275e-005 1.9993587e-0051.3218695e-0082.6568407e-009 1.7931301e-0083.1008915e-005 6.6697659e-006 5.6403374e-0077.0720033e-005 -2.1103715e-007 3.1008915e-005 1.3087511e-001 1.0148977e-003 -1.9981396e-003 1.0001730e-007 7.7772275e-005 6.6697659e-006 1.0148977e-003 1.5539404e-001 3.0264331e-002-2.5748604e-006 1.9993587e-005 5.6403374e-007 -1.9981396e-003 3.0264331e-002 4.0721943e-002左⽚外⽅位元素Xs,Ys,Zs,ψ、ω、κ、为：5.8967023e+003 5.0687355e+003 2.0506347e+003 1.4337709e-002 4.6257617e-0021.1037952e-001右⽚外⽅位元素协因素阵为：3.9305329e-0084.9400147e-010 -1.0339207e-008 6.8065940e-005 -4.2504770e-007 1.8461496e-0064.9400147e-010 3.9051893e-008 3.3958896e-011 -3.9945442e-008 7.6312421e-005 -1.6453951e-005-1.0339207e-008 3.3958896e-011 1.5155886e-008 -2.3705097e-005 3.5940467e-007 -7.3527082e-007 6.8065940e-005 -3.9945442e-008 -2.3705097e-005 1.2229164e-001 -2.3449223e-003 4.8281474e-003-4.2504770e-007 7.6312421e-005 3.5940467e-007 -2.3449223e-003 1.5233230e-001 -2.5374659e-0022.5374659e-0023.6794789e-002GCP在左⽚和右⽚改正后的坐标（x，y）为：1.6019582e-002 7.9954660e-002 -7.3934212e-002 7.8699356e-0028.8559633e-002 8.1141190e-002 -5.2455612e-003 7.8187184e-0021.3352398e-002 -7.9378247e-002 -7.9125440e-002 -7.8877760e-0028.2242309e-002 -8.0017749e-002 -9.8858970e-003 -8.0086832e-002单位权中误差为：±1.515610577029578e-005所求地⾯点的三维坐标（X, Y, Z）为：5.4310348e+003 5.8851463e+003 5.4831646e+0025.1473645e+003 5.0555934e+003 4.8499600e+0025.4957931e+003 5.0826911e+003 5.0668967e+0025.8442434e+003 5.1098033e+003 5.3025650e+0025.5603279e+003 4.2870779e+003 4.6536459e+002七．⼼得体会经过三周的努⼒，这个当初看来艰巨的任务终于在我的不懈努⼒下圆满的完成了。

摄影测量空间后交程序

实验内容:用 C 、VB或MATLAB语言编写空间后方交会-空间前方交会程序程序代码:#include "stdio.h"#include "math.h"#include "iostream.h"#define N 4#define M 2*Nvoid transpose(double *m1,double *m2,int m,int n) //矩阵转置{ int i,j;for(i=0;i<m;i++)for(j=0;j<n;j++)m2[j*m+i]=m1[i*n+j];return;}void inv(double *a,int n)/*正定矩阵求逆*/int i,j,k;for(k=0;k<n;k++){for(i=0;i<n;i++){if(i!=k)*(a+i*n+k)=-*(a+i*n+k)/(*(a+k*n+k)); }*(a+k*n+k)=1/(*(a+k*n+k));for(i=0;i<n;i++){if(i!=k){for(j=0;j<n;j++){if(j!=k)*(a+i*n+j)+=*(a+k*n+j)* *(a+i*n+k); }}}for(j=0;j<n;j++)if(j!=k)*(a+k*n+j)*=*(a+k*n+k);}}}void mult(double *m1,double *m2,double *result,int i_1,int j_12,int j_2)//矩阵相乘{int i,j,k;for(i=0;i<i_1;i++)for(j=0;j<j_2;j++){result[i*j_2+j]=0.0;for(k=0;k<j_12;k++)result[i*j_2+j]+=m1[i*j_12+k]*m2[j+k*j_2];}return;}void main(){int i,m,count;double mm[M+1];double t,w,k,Xs0,Ys0,Zs0,f,S1=0.0,S2=0.0;doubleH[6]={1},a[3],b[3],c[3],Xo[N],Yo[N],Zo[N],A[6*M],B[6*M],l[M],C[36], D[6];double x[N],y[N],X[N],Y[N],Z[N];cout<<"请输入焦距(mm):f=";cin>>f;cout<<"请输入影像比例尺分母:m=";cin>>m;x[0]=-73.93; y[0]=78.706; X[0]=5083.205; Y[0]=5852.099; Z[0]=527.925; x[1]=-5.252; y[1]=78.184; X[1]=5780.02; Y[1]=5906.365; Z[1]=571.549; x[2]=-79.122; y[2]=-78.879; X[2]=5210.879; Y[2]=4258.446; Z[2]=461.81; x[3]=-9.887; y[3]=-80.089;X[3]=5909.264;Y[3]=4314.283;Z[3]=455.484;for(i=0;i<N;i++) {x[i]=x[i]/1000.0; y[i]=y[i]/1000.0; S1 += X[i];S2 += Y[i];}t=w=k=0.0;Xs0=S1/N;Ys0=S2/N;f=f/1000.0;Zs0=m*f;count=0;while(fabs(H[0])>0.00001||fabs(H[1])>0.00001||fabs(H[2])>0.00001||f abs(H[3 ])>0.00001||fabs(H[4])>0.00001||fabs(H[5])>0.00001){count+=1;a[0]=cos(t)*cos(k)-sin(t)*sin(w)*sin(k);a[1]=-cos(t)*sin(k)-sin(t)*sin(w)*cos(k);a[2]=-sin(t)*cos(w);b[0]=cos(w)*sin(k);b[1]=cos(w)*cos(k);b[2]=-sin(w);c[0]=sin(t)*cos(k)+cos(t)*sin(w)*sin(k);c[1]=-sin(t)*sin(k)+cos(t)*sin(w)*cos(k);c[2]=cos(t)*cos(w);for(i=0;i<N;i++){Xo[i]=-f*(a[0]*(X[i]-Xs0)+b[0]*(Y[i]-Ys0)+c[0]*(Z[i]-Zs0))/(a[2]*(X[i]-Xs0)+b[ 2]*(Y[i]-Ys0)+c[2]*(Z[i]-Zs0));Yo[i]=-f*(a[1]*(X[i]-Xs0)+b[1]*(Y[i]-Ys0)+c[1]*(Z[i]-Zs0))/(a[2]*(X[i]-Xs0)+b[2]*(Y [i]-Ys0)+c[2]*(Z[i]-Zs0));Zo[i]=a[2]*(X[i]-Xs0)+b[2]*(Y[i]-Ys0)+c[2]*(Z[i]-Zs0);A[12*i+0]=(a[0]*f+a[2]*x[i])/Zo[i];A[12*i+1]=(b[0]*f+b[2]*x[i])/Zo[i];A[12*i+2]=(c[0]*f+c[2]*x[i])/Zo[i];A[12*i+3]=y[i]*sin(w)-(x[i]*(x[i]*cos(k)-y[i]*sin(k))/f+f*cos(k))*cos(w);A[12*i+4]=-f*sin(k)-x[i]*(x[i]*sin(k)+y[i]*cos(k))/f;A[12*i+5]=y[i];A[12*i+6]=(a[1]*f+a[2]*y[i])/Zo[i];A[12*i+7]=(b[1]*f+b[2]*y[i])/Zo[i];A[12*i+8]=(c[1]*f+c[2]*y[i])/Zo[i];A[12*i+9]=-x[i]*sin(w)-(y[i]*(x[i]*cos(k)-y[i]*sin(k))/f-f*sin(k))*cos(w);A[12*i+10]=-f*cos(k)-y[i]*(x[i]*sin(k)+y[i]*cos(k))/f;A[12*i+11]=-x[i];l[2*i]=x[i]-Xo[i];l[2*i+1]=y[i]-Yo[i]; }transpose(A,B,M,6); mult(B,A,C,6,M,6); mult(B,l,D,6,8,1); inv(C,6);mult(C,D,H,6,6,1); Xs0+=H[0];Ys0+=H[1];Zs0+=H[2];t+=H[3];w+=H[4];k+=H[5];}for(i=0;i<M;i++) {mm[i]=0;}for(i=0;i<M;i++){mm[0]+=l[i]*l[i];}mm[0]=sqrt(mm[0]/(2*N-6)); for(i=0;i<6;i++) {mm[i+1]=sqrt(C[6*i+i])*mm[0];}cout<<"总的叠代次数为:"<<count<<endl;cout<<"像主点的空间坐标为:"<<endl;cout<<"Xs="<<Xs0<<endl;cout<<"Ys="<<Ys0<<endl;cout<<"Zs="<<Zs0<<endl;cout<<"t="<<t<<endl;cout<<"w="<<w<<endl;cout<<"k="<<k<<endl;cout<<"单位权中误差为:"<<mm[0]<<endl;for(i=1;i<7;i++){cout<<"第"<<i<<"个未知数的中误差为:"<<mm[i]<<endl; }}。

前方交会程序

本人系编程爱好者，拿出一些自己编写的程序和大家分享，欢迎大家批评指正，以求共同进步！#include "stdio.h"#include "math.h"void RotationMatrix(double omega, double phi, double kappa, double *matrix){matrix[0]= cos(phi)*cos(kappa)-sin(phi)*sin(omega)*sin(kappa);matrix[1]= -cos(phi)*sin(kappa)-sin(phi)*sin(omega)*cos(kappa);matrix[2]= -sin(phi)*cos(omega);matrix[3]= cos(omega)*sin(kappa);matrix[4]= cos(omega)*cos(kappa);matrix[5]= -sin(omega);matrix[6]= sin(phi)*cos(kappa)+cos(phi)*sin(omega)*sin(kappa);matrix[7]= -sin(phi)*sin(kappa)+cos(phi)*sin(omega)*cos(kappa);matrix[8]= cos(phi)*cos(omega);}void MultMatrix(double *A, double *B, double *Result,int m, int p, int n){int i;int j;int k;for(i=0; i<m; i++){for(j=0; j<n; j++){Result[i*n+j] = 0.0;for(k=0; k<p; k++){Result[i*n+j] += A[i*p+k]*B[k*n+j];}}}return;}//点投影系数的空间前方交会算法void main(){double sx1,sy1,sz1,phi1,omega1,kappa1;double sx2,sy2,sz2,phi2,omega2,kappa2;double N1,N2;//投影系数double Bx,By,Bz;double d1[3],d2[3];//模型点物方坐标double m1[3],m2[3];//同名像点坐标m1[2]=m2[2];double R1[9],R2[9];//旋转矩阵double x,y,z;//模型点物方坐标printf("请输入S1外方位元素：\n");scanf("%lf%lf%lf%lf%lf%lf",&sx1,&sy1,&sz1,&phi1,&omega1,&kappa1);printf("请输入标志点像方坐标：\n");scanf("%lf%lf",&m1[0],&m1[1]);printf("请输入S2外方位元素：\n");scanf("%lf%lf%lf%lf%lf%lf",&sx2,&sy2,&sz2,&phi2,&omega2,&kappa2); printf("请输入标志点像方坐标：\n");scanf("%lf%lf",&m2[0],&m2[1]);//计算Bx，BY，BZBx=sx2-sx1;By=sy2-sy1;Bz=sz2-sz1;//计算旋转矩阵R1和R2RotationMatrix(phi1,omega1,kappa1,R1);RotationMatrix(phi2,omega2,kappa2,R2);//计算d1[3]:dx1,dy1,dz1和d2[3]:dx2,dy2,dz2;MultMatrix(R1,m1,d1,3,3,1);MultMatrix(R2,m2,d2,3,3,1);//计算投影系数N1，N2N1=(Bx*d2[2]-Bz*d2[0])/(d1[0]*d2[2]-d1[2]*d2[0]);N2=(Bx*d1[2]-Bz*d1[0])/(d1[0]*d2[2]-d1[2]*d2[0]);//计算模型点物方坐标x=sx1+Bx+N2*d2[0];y=sy1+By+N2*d2[1];z=sz1+Bz+N2*d2[2];//输出结果printf("模型点坐标为：\n");printf("%lf,%lf,%lf\n",x,y,z);。

摄影测量编程源代码

#include "stdio.h"#include "math.h"#include "Matrixmultiply.c"#include "Matrixtranspose.c"#include "Matrixinverse.c"void main(){int i,j,k,f=0;double x0=0.0, y0=0.0,fk=0.15324; //内方位元素double m=50000; //估算比例尺double R[3][3],XG[6][1],AT[6][8],ATA[6][6],ATL[6][1];double AXG[8][1],V[8][1],VT[1][8],VTV[1][1],m0,D[6][6];double Xs=0.0, Ys=0.0, Zs=0.0,Q=0.0,W=0.0,K=0.0;double X,Y,Z,L[8][1],A[8][6];doubleB[4][5]={-0.08615,-0.06899,36589.41,25273.32,2195.17,-0.05340,0.08221,37631.08,31324.51,72 8.69,-0.01478,-0.07663,39100.97,24934.98,2386.50,0.01046,0.06443,40426.54,30319.81,757.31 };for(i=0;i<4;i++){Xs=Xs+B[i][2];Ys=Ys+B[i][3];Zs=Zs+B[i][4];}Xs=Xs/4; Ys=Ys/4; Zs=m*fk;printf("输出外方位线元素的初始值:\n");printf("%10.7lf\n",Xs);printf("%10.7lf\n",Ys);printf("%10.7lf\n",Zs);do//迭代计算{f++;//组成旋转矩阵R[0][0]=cos(Q)*cos(K)-sin(Q)*sin(W)*sin(K);R[0][1]=-cos(Q)*sin(K)-sin(Q)*sin(W)*cos(K);R[0][2]=-sin(Q)*cos(W);R[1][0]=cos(W)*sin(K);R[1][1]=cos(W)*cos(K);R[1][2]=-sin(W);R[2][0]=sin(Q)*cos(K)+cos(Q)*sin(W)*sin(K);R[2][1]=-sin(Q)*sin(K)+cos(Q)*sin(W)*cos(K);R[2][2]=cos(Q)*cos(W);printf("\n旋转矩阵R为:\n\n");for (i=0;i<3;i++){for(j=0;j<3;j++)printf("%13.5e ",R[i][j]);printf("\n");}//计算系数阵和常数项for( i=0, k=0,j=0;i<=3;i++,k++,j++){X=R[0][0]*(B[i][2]-Xs)+R[1][0]*(B[i][3]-Ys)+R[2][0]*(B[i][4]-Zs);Y=R[0][1]*(B[i][2]-Xs)+R[1][1]*(B[i][3]-Ys)+R[2][1]*(B[i][4]-Zs);Z=R[0][2]*(B[i][2]-Xs)+R[1][2]*(B[i][3]-Ys)+R[2][2]*(B[i][4]-Zs);L[j][0]=B[i][0]-(x0-fk*X/Z);L[j+1][0]=B[i][1]-(y0-fk*Y/Z);printf("输出L矩阵\n");printf("%13.5e\n",L[j][0]);printf("%13.5e\n",L[j+1][0]);j++;A[k][0]=(R[0][0]*fk+R[0][2]*(B[i][0]-x0))/Z;A[k][1]=(R[1][0]*fk+R[1][2]*(B[i][0]-x0))/Z;A[k][2]=(R[2][0]*fk+R[2][2]*(B[i][0]-x0))/Z;A[k][3]=(B[i][1]-y0)*sin(W)-((B[i][0]-x0)*((B[i][0]-x0)*cos(K)-(B[i][1]-y0)*sin(K))/fk+fk*cos(K)) *cos(W);A[k][4]=-fk*sin(K)-(B[i][0]-x0)*((B[i][0]-x0)*sin(K)+(B[i][1]-y0)*cos(K))/fk;A[k][5]=B[i][1]-y0;A[k+1][0]=(R[0][1]*fk+R[0][2]*(B[i][1]-y0))/Z;A[k+1][1]=(R[1][1]*fk+R[1][2]*(B[i][1]-y0))/Z;A[k+1][2]=(R[2][1]*fk+R[2][2]*(B[i][1]-y0))/Z;A[k+1][3]=-(B[i][0]-x0)*sin(W)-((B[i][0]-x0)*((B[i][0]-x0)*cos(K)-(B[i][1]-y0)*sin(K))/fk-fk*sin( K))*cos(W);A[k+1][4]=-fk*cos(K)-(B[i][1]-y0)*((B[i][0]-x0)*sin(K)+(B[i][1]-y0)*cos(K))/fk;A[k+1][5]=-(B[i][0]-x0);k++;}printf("\n\n误差方程系数矩阵A为:\n\n");for (i=0;i<8;i++){for(j=0;j<6;j++)printf("%13.5e ",A[i][j]);printf("\n");}Matrixtranspose(A,AT,8,6);Matrixmultiply(AT,A,ATA,6,8,6);Matrixinverse(ATA,6);Matrixmultiply(AT,L,ATL,6,8,1);Matrixmultiply(ATA,ATL,XG,6,6,1);printf("输出外方位元素近似改正值:\n");printf("%13.5e \n",XG[0][0]);printf("%13.5e \n",XG[1][0]);printf("%13.5e \n",XG[2][0]);printf("%13.5e \n",XG[3][0]);printf("%13.5e \n",XG[4][0]);printf("%13.5e \n",XG[5][0]);Xs=Xs+XG[0][0]; Ys=Ys+XG[1][0]; Zs=Zs+XG[2][0];Q=Q+XG[3][0]; W=W+XG[4][0]; K=K+XG[5][0];printf("%lf",Xs);}while(XG[3][0]>=6.0/206265.0||XG[4][0]>=6.0/206265.0||XG[5][0]>=6.0/206265.0); printf("迭代次数:%d",f);Matrixmultiply(A,XG,AXG,8,6,1);for( i=0;i<8;i++) //计算改正数V[i][0]=AXG[i][0]-L[i][0];Matrixtranspose(V,VT,8,1);Matrixmultiply(VT,V,VTV,1,8,1);m0=VTV[0][0]/2;for(i=0;i<6;i++)for(j=0;j<6;j++)D[i][j]=m0*ATA[i][j];//屏幕输出误差方程系数阵、常数项、改正数printf("\n\n误差方程系数矩阵A为:\n\n");for (i=0;i<6;i++){for(j=0;j<6;j++)printf("%13.5e ",A[i][j]);printf("\n");}printf("\n常数项L为:\n\n");for (i=0;i<8;i++){for(j=0;j<1;j++)printf("%13.5e ",L[i][j]);printf("\n");}printf("\n改正数XG为:\n\n");for (i=0;i<6;i++){for(j=0;j<1;j++)printf("%13.5e ",XG[i][j]);printf("\n");}printf("\n相片的外方位元素为:\n\n");printf(" Xs=%13.7e, Ys=%13.7e, Zs=%13.7e \n\n",Xs,Ys,Zs);printf(" Q=%13.7e, W=%13.7e, K=%13.7e \n",Q,W,K);printf("\n旋转矩阵R为:\n\n");for (i=0;i<3;i++){for(j=0;j<3;j++)printf("%13.5e ",R[i][j]);printf("\n");}printf("\n精度评定结果D为:\n\n");for (i=0;i<6;i++){for(j=0;j<6;j++)printf("%13.5e ",D[i][j]);printf("\n");}}。

摄影测量与遥感：立体像对前方交会

N1
s1
Z1 Y1 X1
点投影系数
S2 A S2a2
X A Xs2 X2
YA Ys2 Y2
ZA Zs2 Z2
N2
A
点投影法前方交会
YA
1 2
[(Ys1
N1Y1
)
(Ys2
N2Y2
)]
X1 x1
Y1
பைடு நூலகம்
R1
y1
,
Z1 f
X 2 x2
Y2
R2
y2
Z2 f
X A X s1 N1X1 X s2 N2 X 2 YA Ys1 N1Y1 Ys2 N Y2 2 Z A Zs1 N1Z1 Zs2 N2Z2
立体像对的前方交会
3.3 立体像对的前方交会
3.3.1 立体像对的前方交会公式
z1
y1
z2
x1 S1
y2
Z
a1(x1,y1)
S2 a2(x2,y2)
x2
• 由立体像对中两张像片的内、外方位元素和像点坐标来确定相应
地面点在物方空
A(X,Y,Z) Y
间坐标系中坐标的方法。
X
Z1
•立体像对空间前方交会
s1
Z1
Zs1 Y1 X1 Ztp
XYst1p Ys1 M
Z2
Y2
Y1
X2
X1
(XA, YA, ZA) Xtp
• 摄影基线
B
s1
BX= Xs2 –Xs1
s2 BZ= Zs2 –Zs1
BY= Ys2 –Ys1
同名光线投影
S1 A S1a1
X A X s1 X1
YA Ys1 Y1

立体像对空间前方交会解算流程

立体像对空间前方交会解算流程下载温馨提示:该文档是我店铺精心编制而成，希望大家下载以后，能够帮助大家解决实际的问题。

文档下载后可定制随意修改，请根据实际需要进行相应的调整和使用，谢谢!并且，本店铺为大家提供各种各样类型的实用资料，如教育随笔、日记赏析、句子摘抄、古诗大全、经典美文、话题作文、工作总结、词语解析、文案摘录、其他资料等等，如想了解不同资料格式和写法，敬请关注!Download tips: This document is carefully compiled by theeditor. I hope that after you download them,they can help yousolve practical problems. The document can be customized andmodified after downloading,please adjust and use it according toactual needs, thank you!In addition, our shop provides you with various types ofpractical materials,such as educational essays, diaryappreciation,sentence excerpts,ancient poems,classicarticles,topic composition,work summary,word parsing,copy excerpts,other materials and so on,want to know different data formats andwriting methods,please pay attention!立体像对空间前方交会解算是摄影测量中的一个重要步骤，用于确定摄影测量中心与地面上某点之间的空间位置关系。

常用的的测量程序vb代码

取一元、二元、五元的硬币共十枚，付给25元钱，有多少种不同的取法?方法一Private Sub Command1_Click()Print "一元", "两元", "五元"For a = 0 To 10For b = 0 To 10For c = 0 To 10If a + 2 * b + 5 * c = 25 And a + b + c = 10 ThenPrint a, b, cEnd IfNext cNext bNext aEnd Sub：方法二Private Sub Command1_Click()Print "一元", "两元", "五元"For a = 0 To 10For b = 0 To 10c = 10 - a - bIf a + 2 * b + 5 * c = 25 And c > 0 ThenPrint a, b, cEnd IfNext bNext aEnd Sub九九乘法表方法一Private Sub Command1_Click() Print Tab(12); "九九乘法表" For i = 1 To 9For j = 1 To iPrint i * j;Next jPrintNext iEnd Sub方法二Private Sub Command2_Click() ShowFontSize = 15Print Tab(12);FontSize = 12PrintFor k = 0 To 9Print Tab(k * 4); k;Next kjiuPrintFor j = 1 To 9Print j;For k = 1 To jPrint Tab(k * 4); j * k;Next kPrintNext jEnd Sub求T = 8! = 1×2×3×…×8 Private Sub Command1_Click()jc = 1n = Val(Text1.Text)For c = 1 To njc = jc * cNext cPrint "jc="; jcEnd Sub用100 元买100 只鸡，母鸡3元1只，小鸡1元3只，问各应买多少只?Private Sub Command1_Click()Dim x As Integer, y As IntegerFor x = 1 To 30y = 100 - xIf 3 * x + y / 3 = 100 ThenPrint "母鸡只数为:"; x,Print "小鸡只数为:"; yEnd IfNext xEnd Sub数组打印数组的上界和下界数值Private Sub Command1_Click() Dim a(1 To 10) As IntegerPrint "下界值", "上界值" Print LBound(a), UBound(a) End Sub数组解决1+2+3+4+5+6+7+8=？Private Sub Command1_Click() Dim a(1 To 10) As IntegerDim sum As IntegerFor b = 1 To 8a(b) = bsum = sum + a(b)Next bText1.Text = sumPrint "1+2+3+4+5+6+7+8=" & sum End Sub任意五个数字之和Private Sub Command1_Click()Dim Data(5) As IntegerDim Sum, I As IntegerFor I = 1 To 5Data(I) = InputBox("输入第" & I & "个数据") Next IFor I = 1 To 5Sum = Sum + Data(I)Next IText1.Text = SumPrint SumEnd Sub连续输入5个数字例如1，2，3，4，51+3+5+7+9=？奇数和Private Sub Command1_Click()Dim a(1 To 5) As IntegerDim sum As IntegerFor x = 1 To 5a(x) = x * 2 - 1sum = sum + a(x)Next xText1.Text = sumPrint sumEnd SubPrivate Function pf(x As Long, y As Long) As Long s = Sqr(x ^ 2 + y ^ 2)pf = sEnd FunctionPrivate Sub Command1_Click()Dim a As LongDim b As LongDim c As Longa = Val(Text1.Text)b = Val(Text2.Text)s = pf(a, b)Print sEnd SubSub过程和Function过程3. 编写过程，求两个数的最大公约数。

立体像对空间前方交会-点投影系数法（python实现）

⽴体像对空间前⽅交会-点投影系数法（python实现）⼀、原理⼆、步骤a.⽤各⾃像⽚的⾓元素计算出左右像⽚的旋转矩阵R1和R2。

b.根据左右像⽚的外⽅位元素计算摄影基线分量Bx,By,Bz。

c.逐点计算像点的空间辅助坐标。

d.计算投影系数。

e.计算未知点的地⾯摄影测量坐标。

f.重复以上步骤完成所有点的地⾯坐标的计算。

三、⽰例代码# -*- coding: utf-8 -*-"""Created on Mon Nov 25 08:18:30 2019@author: L JL"""import numpy as npimport math as mdef r_mat(f,w,k):Rf = np.mat([[m.cos(f), 0, -m.sin(f)],[0, 1, 0],[m.sin(f), 0, m.cos(f)]])Rw = np.mat([[1, 0, 0],[0, m.cos(w), -m.sin(w)],[0, m.sin(w), m.cos(w)]])Rk = np.mat([[m.cos(k), -m.sin(k), 0],[m.sin(k), m.cos(k), 0],[0, 0, 1]])R = Rf*Rw*Rkreturn Rdef SpatialAuxiliaryCoordinate(xy,f,R):coor1 = np.mat([[xy[0]],[xy[1]],[-f]])coor2 = R*coor1return coor2def ProjectionCoefficient(SAC1,SAC2,B):N1 = (B[0,0]*SAC2[2,0]-B[2,0]*SAC2[0,0])/(SAC1[0,0]*SAC2[2,0]-SAC2[0,0]*SAC1[2,0])N2 = (B[0,0]*SAC1[2,0]-B[2,0]*SAC1[0,0])/(SAC1[0,0]*SAC2[2,0]-SAC2[0,0]*SAC1[2,0])return N1,N2#mainleft_HomonymousImagePoints = [0.153,91.798]right_HomonymousImagePoints = [-78.672,89.122]left_In = np.mat([0,0,152.91])left_Ex = np.mat([[970302.448784],[-1138644.971216],[3154.584941],[0.010425],[-0.012437],[0.003380]])right_In = np.mat([0,0,152.91])right_Ex = np.mat([[971265.303768],[-1138634.245942],[3154.784258],[0.008870],[-0.005062],[-0.008703]])R_L = np.mat(np.zeros((3,3)))R_R = np.mat(np.zeros((3,3)))left_SACoordinate = np.mat(np.zeros((3,1)))right_SACoordinate = np.mat(np.zeros((3,1)))baselineComponent = np.mat(np.zeros((3,1)))R_L = r_mat(left_Ex[3,0],left_Ex[4,0],left_Ex[5,0])R_R = r_mat(right_Ex[3,0],right_Ex[4,0],right_Ex[5,0])#left_SpatialAuxiliaryCoordinate = R_L*left_In.T#right_SpatialAuxiliaryCoordinate = R_R*right_In.Tleft_SACoordinate = SpatialAuxiliaryCoordinate(left_HomonymousImagePoints,left_In[0,2],R_L)right_SACoordinate = SpatialAuxiliaryCoordinate(right_HomonymousImagePoints,right_In[0,2],R_R) baselineComponent = right_Ex[0:3,0] - left_Ex[0:3,0]N1,N2 = ProjectionCoefficient(left_SACoordinate,right_SACoordinate,baselineComponent)#GPhotogrammetrCoordinatesGPCoordinates = np.mat(np.zeros((3,1)))GPCoordinates = ((left_Ex[0:3,0]+N1*left_SACoordinate) + (right_Ex[0:3,0]+N2*right_SACoordinate))/2 print("左影像同名点：",left_HomonymousImagePoints)print("左影像同名点：",right_HomonymousImagePoints)print("地⾯点坐标：\n X=%f,\n Y=%f,\n Z=%f"%(GPCoordinates[0,0],GPCoordinates[1,0],GPCoordinates[2,0]))。

1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

Private Sub Command1_Click()Dim zx1 As Single, zy1 As Single, zx2 As Single, zy2 As Single, zx3 As Single, zy3 As Single, zx4 As Single, zy4 As Single, zx5 As Single, zy5 As Single, zx6 As Single, zy6 As SingleDim yx1 As Single, yy1 As Single, yx2 As Single, yy2 As Single, yx3 As Single, yy3 As Single, yx4 As Single, yy4 As Single, yx5 As Single, yy5 As Single, yx6 As Single, yy6 As SingleDim f As SingleDim jd11 As Single, jd12 As Single, jd13 As Single, jd21 As Single, jd22 As Single, jd23 As Single, jd1 As Single, jd2 As SingleDim a1(1 To 3, 1 To 3) As SingleDim a2(1 To 3, 1 To 3) As SingleDim fz1(1 To 6, 1 To 3) As SingleDim fz2(1 To 6, 1 To 3) As SingleDim aa(1 To 6, 1 To 5) As SingleDim p As StringDim bx(1 To 6, 1 To 1) As SingleDim n1(1 To 6, 1 To 1) As Single, n2(1 To 6, 1 To 1) As SingleDim aat(1 To 5, 1 To 6) As SingleDim aataa() As DoubleReDim aataa(1 To 5, 1 To 5)Dim l(1 To 6, 1 To 1) As SingleDim aatl(1 To 5, 1 To 1) As SingleDim atal(1 To 5, 1 To 1) As SingleDim jd11z As Single, jd12z As Single, jd13z As Single, jd21z As Single, jd22z As Single, jd23z As SingleDim jd1z As Single, jd2z As Singlezx1 = Val(Text1(0).Text): zy1 = Val(Text1(1).Text): yx1 = Val(Text1(2).Text): yy1 = Val(Text1(3).Text)zx2 = Val(Text1(4).Text): zy2 = Val(Text1(5).Text): yx2 = Val(Text1(6).Text): yy2 = Val(Text1(7).Text)zx3 = Val(Text1(8).Text): zy3 = Val(Text1(9).Text): yx3 = Val(Text1(10).Text): yy3 = Val(Text1(11).Text)zx4 = Val(Text1(12).Text): zy4 = Val(Text1(13).Text): yx4 = Val(Text1(14).Text): yy4 = Val(Text1(15).Text)zx5 = Val(Text1(16).Text): zy5 = Val(Text1(17).Text): yx5 = Val(Text1(18).Text): yy5 = Val(Text1(19).Text)zx6 = Val(Text1(20).Text): zy6 = Val(Text1(21).Text): yx6 = Val(Text1(22).Text): yy6 = Val(Text1(23).Text)jd11 = Val(Text2(0).Text): jd12 = Val(Text2(1).Text): jd13 = Val(Text2(2).Text)jd21 = Val(Text2(3).Text): jd22 = Val(Text2(4).Text): jd23 = Val(Text2(5).Text)jd1 = Val(Text2(6).Text): jd2 = Val(Text2(7).Text)f = Val(Text3.Text)Do While jd11z < 0.3 * 0.00001 And jd12z < 0.3 * 0.00001 And jd13z < 0.3 * 0.00001 And jd21z < 0.3 * 0.00001 And jd22z < 0.3 * 0.00001 And jd23z < 0.3 * 0.00001 And jd1z < 0.3 * 0.00001 And jd2z < 0.3 * 0.00001jd11z = jd11z + jd11: jd12z = jd12z + jd12: jd13z = jd13z + jd13jd21z = jd21z + jd21: jd22z = jd22z + jd22: jd23z = jd23z + jd23jd1z = jd1z + jd1: jd2z = jd2z + jd2jd11 = jd11z: jd12 = jd12z: jd13 = jd13z: jd21 = jd21z: jd22 = jd22z: jd23 = jd23z: jd1 = jd1z: jd2 = jd2zFor i = 1 To 3For j = 1 To 3a1(1, 1) = Cos(jd11) * Cos(jd13) - Sin(jd11) * Sin(jd12) * Sin(jd13): a1(1, 2) = -Cos(jd11) * Sin(jd13) - Sin(jd11) * Sin(jd12) * Cos(jd13): a1(1, 3) = -Sin(jd11) * Cos(jd12)a1(2, 1) = Cos(jd12) * Sin(jd13): a1(2, 2) = Cos(jd12) * Cos(jd13): a1(2, 3) = -Sin(jd12)a1(3, 1) = Sin(jd11) * Cos(jd13) + Cos(jd11) * Sin(jd12) * Sin(jd13): a1(3, 2) = -Sin(jd11) * Sin(jd13) + Cos(jd11) * Sin(jd12) * Cos(jd13): a1(3, 3) = Cos(jd11) * Cos(jd12)NextNextFor i = 1 To 3p = ""For j = 1 To 3a2(1, 1) = Cos(jd21) * Cn(jd22) * Sin(jd23): a2(1, 2) = -Cos(jd21) * Sin(jd23) - Sin(jd21) * Sin(jd22) * Cos(jd2= Cos(jd22) * Cos(jd23): a2(2, 3) = -Sin(jd22)a2(os(jd21) * Sin(jd22) * Cos(jd23): a2(3, 3) = Cos(jd21) * Cos(jd22)NextNextFor i = 1 To 6For j = 1 To 3fz1(1, 1) = a1(1, 1) * zx1 + a1(1, 2) * zy1 + a1(1, 3) * -f: fz1(1, 2) = a1(2, 1) * zx1 + a1(2, 2) * zy1 + a1(2, 3) * -f: fz1(1, 3) = a1(3, 1) * zx1 + a1(3, 2) * zy1 + a1(3, 3) * -ffz1(2, 1) = a1(1, 1) * zx2 + a1(1, 2) * zy2 + a1(1, 3) * -f: fz1(2, 2) = a1(2, 1) * zx2 + a1(2, 2) * zy2 + a1(2, 3) * -f: fz1(2, 3) = a1(3, 1) * zx2 + a1(3, 2) * zy2 + a1(3, 3) * -ffz1(3, 1) = a1(1, 1) * zx3 +3 + a1(1, 3) * -f: fz1(3, 2) = a1(2, 1) * zx3 + a1(2, 2) * zy3 + a1(2, 3) * -f: fz1(3, 3) =a1(1, 3) * -f: fz1(4, 2) = a1(2, 1) * zx4 + a1(2, 2) * zy4 + a1(2, 3) * -f: fz1(4, 3) = a1(3, 1) * zx4 + a1(3, 2) * zy4 + a1(3, 3) * -ffz1(5, 1) = a1(1, 1) * zx5 + a1(1, 2) * zy5 + a1(1, 3) * -f: fz1(5, 2) = a1(2, 1) * zx5 +a1(2, 2) * zy5 + a1(2, 3) * -f: fz1(5, 3) = a1(3, 1) * zx5 + a1(3, 2) * zy5 + a1(3, 3) * -ffz1(6, 1) = a1(1, 1) * zx6 + a1(1, 2) * zy6 + a1(1, 3) * -f: fz1(6, 2) = a1(2, 1) * zx6 + a1(2, 2) * zy6 + a1(2, 3) * -f: fz1(6, 3) = a1(3, 1) * zx6 + a1(3, 2) * zy6 + a1(3, 3) * -fNextNextFor i = 1 To 6For j = 1 To 3fz2(1, 1) = a2(1, 1) * yx1 + a2(1, 2) * yy1 + a2(1, 3) * -f: fz2(1, 2) = a2(2, 1) * yx1 + a2(2, 2) * yy1 + a2(2, 3) * -f: fz2(1, 3) = a2(3, 1) * yx1 + a2(3, 2) * yy1 + a2(3, 3) * -ffz2(2, 1) = a2(1, 1) * yx2 + a2(1, 2) * yy2 + a2(1, 3) * -f: fz2(2, 2) = a2(2, 1) * yx2 + a2(2, 2) * yy2 + a2(2, 3) * -f: fz2(2, 3) = a2(3, 1) * yx2 + a2(3, 2) * yy2 + a2(3, 3) * -ffz2(3, 1) = a2(1, 1) * yx3 2(3, 1) * yx3 + a2(3, 2) * yy3 + a2(3, 3) * -ffz2(4, 1) = a2(1, 1) * yx4 + a2(1, 2) * yy4 + a2(1, 3) * -f: fz2(4, 2) = a2(2, 1) * yx4 + a2(2, 2) * yy4 + a2(2, 3) * -f: fz2(4, 3) = a2(3, 1) * yx4 + a2(3, 2) * yy4 + a2(3, 3) * -ffz2(5, 1) = a2(1, 1) * yx5 + a2(1, 2) * yy5 + a2(1, 3) * -f: fz2(5, 2) = a2(2, 1) * yx5 + a2(2, 2) * yy5 + a2(2, 3) * -f: fz2(5, 3) = a2(3, 1) * yx5 + a2(3, 2) * yy5 + a2(3, 3) * -ffz2(6, 1) = a2(1, 1) * yx6 +2(1, 3) * -f: fz2(6, 2) = a2(2, 1) * yx6 + a2(2, 2) * yy6 + a2(2, 3) * -f: fz2(6, 3) = a2(3, 1) * yx6 + a2(3, 2) * yy6 + a2(3, 3) * -fNextNextbx(1, 1) = zx1 - yx1: bx(2, 1) = zx2 - yx2: bx(3, 1) = zx3 - yx3: bx(4, 1) = zx4 - yx4: bx(5, 1) = zx5 - yx5: bx(6, 1) = zx6 - yx6For i = 1 To 6For j = 1 To 1n1(1, 1) = (bx(1, 1) * fz2(1, 3) - bx(1, 1) * jd2 * fz2(1, 1)) / (fz1(1, 1) * fz2(1, 3) - fz2(1, 1) * fz1(1, 3))n1(2, 1) = (bx(2, 1) * fz2(2, 3) - bx(2, 1) * jd2 * fz2(2, 1)) / (fz1(2, 1) * fz2(2, 3) - fz2(2, 1) * fz1(2, 3))n1(3, 1) = (bx(3, 1) * fz2(3, 3) - bx(3, 1) * jd2 * fz2(3, 1)) / (fz1(3, 1) * fz2(3, 3) - fz2(3, 1) * fz1(3, 3))n1(4, 1) = (bx(4, 1) * fz2(4, 3) - bx(4, 1) * jd2 * fz2(4, 1)) / (fz1(4, 1) * fz2(4, 3) - fz2(4, 1) * fz1(4, 3))n1(5, 1) = (bx(5, 1) * fz2(5 1) * jd2 * fz2(5, 1)) / (fz1(5, 1) * fz2(5, 3) - fz2(5, 1) * fz1(5, 3))n1(6, 1) = (bx(6, 1) * fz2(6, 3) - bx(6, 1) * jd2 * fz2(6, 1)) / (fz1(6, 1) * fz2(6, 3) - fz2(6, 1) * fz1(6, 3))NextNextFor i = 1 To 6For j = 1 To 1n2(1, 1) = (bx(1, 1) * fz1(1, 3) - bx(1, 1) * jd2 * fz1(1, 1)) / (fz1(1, 1) * fz2(1, 3) - fz2(1, 1) * fz1(1, 3))n2(2, 1) = (bx(2, 1) * fz1(2, 3) - bx(2, 1) * jd2 * fz1(2, 1)) / (fz1(2, 1) * fz2(2, 3) - fz2(2, 1) * fz1(2, 3))n2(3, 1) = (bxn2(5, 1) = (bx(5, 1) * fz1(5, 3) - bx(5, 1) * jd2 * fz1(5, 1)) / (fz1(5, 1) * fz2(5, 3) - fz2(5, 1) * fz1(5, 3))n2(6, 1) = (bx(6, 1) * fz1(6, 3) - bx(6, 1) * jd2 * fz1(6, 1)) / (fz1(6, 1) * fz2(6, 3) - fz2(6, 1) * fz1(6, 3))NextNextFor i = 1 To 6For j = 1 To 5aa(1, 1) = bx(1, 1): aa(1, 2) = -(fz2(1, 2) / fz2(1, 3)) * bx(1, 1): aa(1, 3) = -((fz2(1, 1) * fz2(1, 2)) / fz2(1, 3)) * n2(1, 1): aa(1, 4) = -(fz2(1, 3) + (fz2(1, 2) * fz2(1, 2)) / fz2(1, 3)) */ fz2(2, 3)) * bx(2, 1): aa(2, 3) = -((fz2(2, 1) * fz2(2, 2)) / fz2(2, 3)) * n2(2, 1): aa(2, 4) = -(fz2(2, 3) + (fz2(2, 2) * fz2(2, 2)) / fz2(2, 3)) * n2(2, 1): aa(2, 5) = fz2(2, 1) * n2(2, 1)aa(3, 1) = bx(3, 1): aa(3, 2) = -(fz2(3, 2) / fz2(3, 3)) * bx(3, 1): aa(3, 3) = -((fz2(3, 1) * fz2(3, 2)) / fz2(3, 3)) * n2(3, 1): aa(3, 4) = -(fz2(3, 3) + (fz2(3, 2) * fz2(3, 2)) / fz2(3, 3))4, 2) / fz2(4, 3)) * bx(4, 1): aa(4, 3) = -((fz2(4, 1) * fz2(4, 2)) / fz2(4, 3)) * n2(4, 1): aa(4, 4) = -(fz2(4, 3) + (fz2(4, 2) * fz2(4, 2)) /-(fz2(5, 2) / fz2(5, 3)) * bx(5, 1): aa(5, 3) = -((fz2(5, 1) * fz2(5, 2)) / fz2(5, 3)) * n2(5, 1): aa(5, 4) = -(fz2(5, 3) + (fz2(5, 2) * fz2(5, 2)) / fz2(5, 3)) * n2(5, 1): aa(5, 5) = fz2(5, 1) * n2(5, 1)aa(6, 1) = bx(6, 1): aa(6, 2) = -(fz2(6, 2) / fz2(6, 3)) * bx(6, 1): aa(6, 3) = -((fz2(6, 1) * fz2(6, 2)) / fz2(6, 3)) * n2(6, 1): aa(6, 4) = -(fz2(6, 3) + (fz2(6, 2) * fz2(6, 2)) / fz2(6, 3)) * n2(6, 1): aa(6, 5) = fz2(6, 1) * n2(6, 1)NextNextFor i = 1 To 5For j = 1 To 6aat(i, j) = aa(j, i)NextNextFor i = 1 To 5For j = 1 To 5For t = 1 To 6aataa(i, j) = aataa(i, j) + aat(i, t) * aa(t, j)NextNextFor o = 1 To 5For m = 1 To 1NextNextNextFor i = 1 To 6For j = 1 To 1l(1, 1) = n1(1, 1) * fz1(1, 2) - n2(1, 1) * fz2(1, 2) - bx(1, 1) * jd1 l(2, 1) = n1(2, 1) * fz1(2, 2) - n2(2, 1) * fz2(2, 2) - bx(2, 1) * jd1 l(3, 1) = n1(3, 1) * fz1(3, 2) - n2(3, 1) * fz2(3, 2) - bx(3, 1) * jd1 l(4, 1) = n1(4, 1) * fz1(4, 2) - n2(4, 1) * fz2(4, 2) - bx(4, 1) * jd1 l(5, 1) = n1(5, 1) * fz1(5, 2) - n2(5, 1) * fz2(5, 2) - bx(5, 1) * jd1 l(6, 1) = n1(6, 1) * fz1(6, 2) - n2(6, 1) * fz2(6, 2) - bx(6, 1) * jd1 NextNextFor i = 1 To 5For j = 1 To 1For t = 1 To 6aatl(i, j) = aatl(i, j) + aat(i, t) * l(t, j)NextNextFor o = 1 To 5For m = 1 To 1NextNextNextFor i = 1 To 5For j = 1 To 1For t = 1 To 5atal(i, j) = atal(i, j) + aataa(i, t) * aatl(t, j)NextNextFor o = 1 To 5For m = 1 To 1NextNextNextLoopEnd SubPrivate Function MRinv(n As Integer, mtxA() As Double) As Boolean ReDim nIs(n) As Integer, nJs(n) As IntegerDim i As Integer, j As Integer, k As IntegerDim D As Double, p As DoubleFor k = 1 To nD = 0#For i = k To nFor j = k To np = Abs(mtxA(i, j))If (p > D) ThenD = pnIs(k) = inJs(k) = jEnd IfNext jNext iIf (D + 1# = 1#) ThenMRinv = FalseExit FunctionEnd IfIf (nIs(k) <> k) ThenFor j = 1 To np = mtxA(k, j)mtxA(k, j) = mtxA(nIs(k), j)mtxA(nIs(k), j) = pNext jEnd IfIf (nJs(k) <> k) ThenFor i = 1 To np = mtxA(i, k)mtxA(i, k) = mtxA(i, nJs(k))mtxA(i, nJs(k)) = pNext iEnd IfmtxA(k, k) = 1# / mtxA(k, k)For j = 1 To nIf (j <> k) Then mtxA(k, j) = mtxA(k, j) * mtxA(k, k)Next jFor i = 1 To nIf (i <> k) ThenFor j = 1 To nIf (j <> k) Then mtxA(i, j) = mtxA(i, j) - mtxA(i, k) * mtxA(k, j)Next jEnd IfNext iFor i = 1 To nIf (i <> k) Then mtxA(i, k) = -mtxA(i, k) * mtxA(k, k)Next iNext kFor k = n To 1 Step -1If (nJs(k) <> k) ThenFor j = 1 To np = mtxA(k, j)mtxA(k, j) = mtxA(nJs(k), j)mtxA(nJs(k), j) = pNext jEnd IfIf (nIs(k) <> k) ThenFor i = 1 To np = mtxA(i, k)mtxA(i, k) = mtxA(i, nIs(k))mtxA(i, nIs(k)) = pNext iEnd IfNext kMRinv = TrueEnd Function。