摄影测量立体相对的前方交会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七、实验心得通过本次前方交会实验，熟悉掌握了前方交会原理，并利用计算机程序对其进行了解决，收益颇多！。
立体像对前方交会绝对定向
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七．⼼得体会经过三周的努⼒，这个当初看来艰巨的任务终于在我的不懈努⼒下圆满的完成了。
- 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。