光束法区域网空中三角测量作业报告

1、 光束法区域网空中三角测量求像片外方位元素和加密点的地面坐标 一、原理：首先给出的原始数据像素坐标转换成像平面直角坐标；然后利用后方交会原理并采用控制点的向平面直角坐标和地面坐标迭代求解每张像片的外方位元素；再通过前方交会求加密点的地面坐标；最后进行精度评定。二、计算流程：1、 求出四个控制点及所有待求点的像点坐标（x1,y1）与(x2,y2)；2、 空间后方交会计算两像片的外方位元素：根据编好的计算机程序读取原始文件中的数据控制点的地面坐标和相应的像点像素坐标（第一步已经转换成像平面直角坐标），对两张像片各自进行空间后方交会，计算各自的6个外方位元素Xs1,Ys1,Zs1,1,w1,k1和X

2、s2,Ys2,Zs2,2,w2,k2；3、 空间前方交会计算待定点的地面坐标：用各片的外方位角元素计算左右片的方向余弦值，组成旋转矩阵R1与R2，逐点计算像点的像空间爱你辅助坐标(u1,v1,w1)及(u2,v2,w2)；根据外方位元素计算基线分量（Bu,Bv,Bw）;计算投影系数N1,N2；计算待定点在各自的像空间辅助坐标系中的坐标(U1,V1,W1)及(U2,V2,W2)；最后计算待定点的地面摄影测量坐标。4、 进行精度评定：利用前方交会求出地面控制点的坐标，再与地面实测坐标求差，将两者差值视为真差值，由这些真差值计算出点位坐标精度。三、基本公式：1、像素坐标（I,J）转换成像平面直角坐标

3、(x,y)：x = (J- 2136) * 5.19 / 1000 / 1000;y = (1424 - I) * 5.19 / 1000 / 1000;2、后方交会：共线条件方程： 外方位元素系数矩阵：Zp=a3*(L_X-Xs)+b3 *(L_Y-Ys)+c3 *(L_Z-Zs);a11= 1/Zp*(a1*f+a3*x);a12 =1/Zp*(b1*f+b3*x);a13= 1/Zp*(c1*f+c3*x);a14= y* Sin (w)-(x /f*(x*Cos(k)-y* Sin (k)+f* Cos (k)* Cos (w );a15= -f* Sin (k)-x /f*(x*Si

4、n (k)+y *Cos (k);a16 = y;a21 = 1/Zp*(a2*f+a3*y);a22 = 1/Zp*(b2*f+b3*y); a23 = 1/Zp*(c2*f+c3*y); a24 = -x* Sin (w )-(x/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;3、前方交会：计算旋转矩阵：a = Cos() * Cos(k) - Sin(w) * Sin() * Sin(k); a2 = Cos() * Sin(k) -

5、 Sin(w) *Sin() * Cos(k);a3 = - Sin() *Cos(w);b1 = Cos(w) * Sin(k) ;b2 = Cos(w) * Math.Cos(k);b3 = Sin(w);c1 = =Sin() *.Cos(k) + Cos() * Sin(w) *Sin(k);c2 =- Sin() * Sin(k) + Cos() * Sin(w) * Cos(k);c3 =Cos() * Cos(w);基线分量：Bu = Xs2 - Xs1;Bv = Ys2 - Ys1;Bw = Zs2 - Zs1;投影系数：计算地面点坐标：5、 精度评定：四、实例：使用数码相机拍

6、摄了两张航测相片Img261和Img262，量测了6对点的影像坐标，其中前4个点为地面控制点（GCP）, 后两个点为加密点。如下表所示。相片参数：大小为4272（列）*2848（行），像主点位于影像中心，像素大小为5.19micron，相机的焦距为24mm。.Img261的外方位元素初始值：X=397510.760，Y=3445853.978, Z=1455.153，角元素为0；Img262的外方位元素初始值：X=397513.320，Y=3445979.811, Z=1453.685，角元素为0；请采用光束法区域网空中三角测量方法，解求每张相片的外方位元素和加密点的地面坐标。点名Img261

7、Img262地面坐标x（pixel）y（pixel）x（pixel）y（pixel）X(m)Y(m)H(m)13156.62595.1253064.875903.625397534.5023446082.739655.925123206.857983.3163139.3751792.375397570.8543445925.794659.37273705.375462.6253625.8751251.625397645.9093446035.598653.428161660.8231122.9821598.8751986.375397303.2383445852.908677.94933134

8、.125320.1253047.8751130.625？22188.625624.8752105.6251491.125？原始数据文件：生成结果文件：程序源代码：using System.Collections.Generic;using System.Linq;using System.Text;using System.IO;namespace 摄影测量光束法编程作业 class Program static void Main(string args) /原始观测数据：控制点和加密点的像素坐标（x1,y1,x2,y2），控制点的地面坐标(X,Y,Z) FileStream fs = ne

9、w FileStream("控制点观测数据.txt", FileMode.Open, FileAccess.Read); FileStream fs2 = new FileStream("控制点观测数据.txt", FileMode.Open, FileAccess.Read); StreamReader g_line = new StreamReader(fs); StreamReader r = new StreamReader(fs2); int line = 0; while (g_line.ReadLine() != null) line+;

10、g_line.Close(); double L_X = new doubleline; double L_Y = new doubleline; double L_Z = new doubleline; double I1 = new doubleline; double J1 = new doubleline; double I2 = new doubleline; double J2 = new doubleline; double x1 = new doubleline; double y1 = new doubleline; double x2 = new doubleline; d

11、ouble y2 = new doubleline; int n = 0; for (n = 0; n < line; n+) string asd = r.ReadLine(); string kon; kon = asd.Split(new Char ',', ',', ',' ); I1n = double.Parse(kon0); J1n = double.Parse(kon1); I2n = double.Parse(kon2); J2n = double.Parse(kon3); L_Xn = double.Parse(kon4

12、); L_Yn = double.Parse(kon5); L_Zn = double.Parse(kon6); /框标像素坐标转换成框标坐标(x1line,y1line,x2line,y2line) int j = 0; for (j = 0; j < line; j+) x1j = (J1j - 2136) * 5.19 / 1000 / 1000; y1j = (1424 - I1j) * 5.19 / 1000 / 1000; x2j = (J2j - 2136) * 5.19 / 1000 / 1000; y2j = (1424 - I2j) * 5.19 / 1000 / 1

13、000; /相机焦距和相片外方位元素的初始值 double f = 0.024; double Xs1 = 397510.760; double Ys1 = 3445853.978; double Zs1 = 1455.153; double w1 = 0.0; double 1 = 0.0; double k1 = 0.0; /计算旋转矩阵 Matrix R1 = new Matrix(3, 3); double dX1 =new double 6; double xzhi = 1.0 / 3600 / 180 * Math.PI; int d1 =0; do d1+; double a1

14、= R10, 0 = Math.Cos(1) * Math.Cos(k1) - Math.Sin(w1) * Math.Sin(1) * Math.Sin(k1); double a2 = R10, 1 = -Math.Cos(1) * Math.Sin(k1) - Math.Sin(w1) * Math.Sin(1) * Math.Cos(k1); double a3 = R10, 2 = -Math.Sin(1) * Math.Cos(w1); double b1 = R11, 0 = Math.Cos(w1) * Math.Sin(k1) ; double b2 = R11, 1 = M

15、ath.Cos(w1) * Math.Cos(k1); double b3 = R11, 2 = -Math.Sin(w1); double c1 = R12, 0 = Math.Sin(1) * Math.Cos(k1) + Math.Cos(1) * Math.Sin(w1) * Math.Sin(k1); double c2 = R12, 1 = -Math.Sin(1) * Math.Sin(k1) + Math.Cos(1) * Math.Sin(w1) * Math.Cos(k1); double c3 = R12, 2 = Math.Cos(1) * Math.Cos(w1);

16、/计算像点坐标近似值，求出常数项lx1,ly1; double jx1 = new doubleline; double jy1 = new doubleline; double lx1 = new doubleline; double ly1 = new doubleline; int i = 0; for (i = 0; i < line; i+) jx1i = -f * (a1 * (L_Xi - Xs1) + b1 * (L_Yi - Ys1) + c1 * (L_Zi - Zs1) / (a3 * (L_Xi - Xs1) + b3 * (L_Yi - Ys1) + c3 *

17、(L_Zi - Zs1); jy1i = -f * (a2 * (L_Xi - Xs1) + b2 * (L_Yi - Ys1) + c2 * (L_Zi - Zs1) / (a3 * (L_Xi - Xs1) + b3 * (L_Yi - Ys1) + c3 * (L_Zi - Zs1); lx1i = x1i - jx1i; ly1i = y1i - jy1i; /像片的常数项为l1 double l1 = new doubleline * 2; for (int m2 = 0; m2 < line * 2 - 1; m2+) int m3 = m2 / 2; l1m2 = lx1m

18、3; l1m2 + 1 = ly1m3; m2+; /求出外方位元素系数矩阵HA1,HA2及其中的元素a11,a12,a13,a14,a15. double, HA1 = new doubleline * 2, 6; double a11 = new doubleline; double a12 = new doubleline; double a13 = new doubleline; double a14 = new doubleline; double a15 = new doubleline; double a16 = new doubleline; double a21 = new

19、doubleline; double a22 = new doubleline; double a23 = new doubleline; double a24 = new doubleline; double a25 = new doubleline; double a26 = new doubleline; int n3 = 0; for (n3 = 0; n3 < line * 2 - 1; n3+) int n1 = n3 / 2; double Zp1=a3*(L_Xn1-Xs1)+b3 *(L_Yn1-Ys1)+c3 *(L_Zn1-Zs1 ); a11n1 = HA1n3,

20、 0 = 1/Zp1*(a1*f+a3*x1n1); a12n1 = HA1n3, 1 = 1/Zp1*(b1*f+b3*x1n1); a13n1 = HA1n3, 2 = 1/Zp1*(c1*f+c3*x1n1); a14n1 = HA1n3, 3 = y1n1*Math.Sin (w1 )-(x1 n1/f*(x1 n1*Math .Cos(k1)-y1 n1*Math .Sin (k1)+f*Math .Cos (k1)*Math .Cos (w1 ); a15n1 = HA1n3, 4 = -f*Math .Sin (k1)-x1 n1/f*(x1 n1*Math .Sin (k1)+

21、y1 n1*Math .Cos (k1); a16n1 = HA1n3, 5 = y1n1; a21n1 = HA1n3 + 1, 0 = 1/Zp1*(a2*f+a3*y1n1); a22n1 = HA1n3 + 1, 1 = 1/Zp1*(b2*f+b3*y1n1); a23n1 = HA1n3 + 1, 2 = 1/Zp1*(c2*f+c3*y1n1); a24n1 = HA1n3 + 1, 3 = -x1n1*Math.Sin (w1 )-(x1 n1/f*(x1 n1*Math .Cos(k1)-y1 n1*Math .Sin (k1)-f*Math .Sin (k1)*Math .

22、Cos (w1 ); a25n1 = HA1n3 + 1, 4 = -f*Math .Cos (k1)-y1 n1/f*(x1 n1*Math .Sin (k1)+y1 n1*Math .Cos (k1); a26n1 = HA1n3 + 1, 5 = -x1n1; n3+; /求外方位元素的改正数 Matrix ha1 = new Matrix(HA1); Matrix L1 = new Matrix(line * 2, 1, l1); Matrix HA1T = ha1.Transpose(); Matrix HA1THA1 = HA1T.Multiply(ha1); HA1THA1.In

23、vertGaussJordan(); Matrix z1 = HA1THA1.Multiply(HA1T); Matrix z2 = z1.Multiply(L1); dX1 = z2; Xs1 = Xs1 + dX10; Ys1 = Ys1 + dX11; Zs1 = Zs1 + dX12; 1 = 1 + dX13; w1 = w1 + dX14; k1 = k1 + dX15; while (Math.Abs(dX13) > xzhi | Math.Abs(dX14) > xzhi| Math.Abs(dX15) > xzhi ); /相机焦距和相片外方位元素的初始值

24、double Xs2 = 397513.3200; double Ys2 = 3445979.811; double Zs2 = 1453.685; double w2 = 0.0; double 2 = 0.0; double k2 = 0.0; /计算两张像片的外方位元素Xs1, Ys1, Zs1, 1, w1, k1， Xs2, Ys2, Zs2, 2, w2, k2 Matrix R2 = new Matrix(3, 3); double dX2 = new double 6; int d2 = 0; do d2+; double aa1 = R20, 0 = Math.Cos(2)

25、* Math.Cos(k2) - Math.Sin(w2) * Math.Sin(2) * Math.Sin(k2); double aa2 = R20, 1 = -Math.Cos(2) * Math.Sin(k2) - Math.Sin(w2) * Math.Sin(2) * Math.Cos(k2); double aa3 = R20, 2 = -Math.Sin(2) * Math.Cos(w2); double bb1 = R21, 0 = Math.Cos(w2) * Math.Sin(k2); double bb2 = R21, 1 = Math.Cos(w2) * Math.C

26、os(k2); double bb3 = R21, 2 = -Math.Sin(w2); double cc1 = R22, 0 = Math.Sin(2) * Math.Cos(k2) + Math.Cos(2) * Math.Sin(w2) * Math.Sin(k2); double cc2 = R22, 1 = -Math.Sin(2) * Math.Sin(k2) + Math.Cos(2) * Math.Sin(w2) * Math.Cos(k2); double cc3 = R22, 2 = Math.Cos(2) * Math.Cos(w2); /计算旋转矩阵 double j

27、x2 = new doubleline; double jy2 = new doubleline; double lx2 = new doubleline; double ly2 = new doubleline; int i = 0; for (i = 0; i < line; i+) jx2i = -f * (aa1 * (L_Xi - Xs2) + bb1 * (L_Yi - Ys2) + cc1 * (L_Zi - Zs2) / (aa3 * (L_Xi - Xs2) + bb3 * (L_Yi - Ys2) + cc3 * (L_Zi - Zs2); jy2i = -f * (

28、aa2 * (L_Xi - Xs2) + bb2 * (L_Yi - Ys2) + cc2 * (L_Zi - Zs2) / (aa3 * (L_Xi - Xs2) + bb3 * (L_Yi - Ys2) + cc3 * (L_Zi - Zs2); lx2i = x2i - jx2i; ly2i = y2i - jy2i; double l2 = new doubleline * 2; for (int m2 = 0; m2 < line * 2 - 1; m2+) int m3 = m2 / 2; l2m2 = lx2m3; l2m2 + 1 = ly2m3; m2+; double

29、, HA2 = new doubleline * 2, 6; double aa11 = new doubleline; double aa12 = new doubleline; double aa13 = new doubleline; double aa14 = new doubleline; double aa15 = new doubleline; double aa16 = new doubleline; double aa21 = new doubleline; double aa22 = new doubleline; double aa23 = new doubleline;

30、 double aa24 = new doubleline; double aa25 = new doubleline; double aa26 = new doubleline; for (int n4 = 0; n4 < line * 2 - 1; n4+) int n2 = n4 / 2; double Zp2=aa3*(L_Xn2-Xs2)+bb3 *(L_Yn2-Ys2)+cc3 *(L_Zn2-Zs2 ); aa11n2 = HA2n4, 0 = 1/Zp2*(aa1*f+aa3*x2n2); aa12n2 = HA2n4, 1 = 1/Zp2*(bb1*f+bb3*x2n2

31、); aa13n2 = HA2n4, 2 = 1/Zp2*(cc1*f+cc3*x2n2); aa14n2 = HA2n4, 3 = y2n2*Math.Sin (w2 )-(x2 n2/f*(x2 n2*Math .Cos(k2)-y2 n2*Math .Sin (k2)+f*Math .Cos (k2)*Math .Cos (w2 ); aa15n2 = HA2n4, 4 = -f*Math .Sin (k2)-x2 n2/f*(x2 n2*Math .Sin (k2)+y2 n2*Math .Cos (k2); aa16n2 = HA2n4, 5 = y2n2; aa21n2 = HA2

32、n4 + 1, 0 =1/Zp2*(aa2*f+aa3*y2n2); aa22n2 = HA2n4 + 1, 1 =1/Zp2*(bb2*f+bb3*y2n2); aa23n2 = HA2n4 + 1, 2 = 1/Zp2*(cc2*f+cc3*y2n2); aa24n2 = HA2n4 + 1, 3 = -x2n2*Math.Sin (w2 )-(x2 n2/f*(x2 n2*Math .Cos(k2)-y2n2*Math .Sin (k2)-f*Math .Sin (k2)*Math .Cos (w2 ); aa25n2 = HA2n4 + 1, 4 = -f*Math .Cos (k2)

33、-y2 n2/f*(x2 n2*Math .Sin (k2)+y2 n2*Math .Cos (k2); aa26n2 = HA2n4 + 1, 5 = -x2n2; n4+; Matrix ha2 = new Matrix(HA2); Matrix L2 = new Matrix(line * 2, 1, l2); Matrix HA2T = ha2.Transpose(); Matrix HA2THA2 = HA2T.Multiply(ha2); HA2THA2.InvertGaussJordan(); Matrix z3 = HA2THA2.Multiply(HA2T); Matrix

34、z4 = z3.Multiply(L2); dX2 = z4; Xs2 = Xs2 + dX20; Ys2 = Ys2 + dX21; Zs2 = Zs2 + dX22; 2 = 2 + dX23; w2 = w2 + dX24; k2 = k2 + dX25; while (Math.Abs(dX23) > xzhi | Math.Abs(dX24) > xzhi | Math.Abs(dX25) > xzhi); /求待定点的像点坐标jx1jj， jy1jj，jx2jj，jy2jj FileStream fs1 = new FileStream("待定点像点坐标

35、观测数据.txt", FileMode.Open, FileAccess.Read); FileStream fs12 = new FileStream("待定点像点坐标观测数据.txt", FileMode.Open, FileAccess.Read); StreamReader g_line2 = new StreamReader(fs1); StreamReader r2 = new StreamReader(fs12); int line2 = 0; while (g_line2.ReadLine() != null) line2+; g_line2.Cl

36、ose(); double jI1 = new doubleline2; double jJ1 = new doubleline2; double jI2 = new doubleline2; double jJ2 = new doubleline2; double jjx1 = new doubleline2; double jjy1 = new doubleline2; double jjx2 = new doubleline2; double jjy2 = new doubleline2; int jn = 0; for (jn = 0; jn < line2; jn+) stri

37、ng asd2 = r2.ReadLine(); string han2; han2 = asd2.Split(new Char ',', ',', ',' ); jI1jn = double.Parse(han20); jJ1jn = double.Parse(han21); jI2jn = double.Parse(han22); jJ2jn = double.Parse(han23); int jj = 0; for (jj = 0; jj < line2; jj+) jjx1jj = (jJ1jj - 2136) * 5.19 /

38、1000 / 1000; jjy1jj = (1424 - jI1jj) * 5.19 / 1000 / 1000; jjx2jj = (jJ2jj - 2136) * 5.19 / 1000 / 1000; jjy2jj = (1424 - jI2jj) * 5.19 / 1000 / 1000; /计算结果精度评定(实际精度) /首先用前方交会求地面控制点的坐标 double jX = new doubleline; double jY = new doubleline; double jZ = new doubleline; double Bu = Xs2 - Xs1; double B

39、v = Ys2 - Ys1; double Bw = Zs2 - Zs1; double ux = new doubleline; double uy = new doubleline; double uz = new doubleline; for(int i=0;i<line;i+) double, xyf1 = new double3,line; xyf10,i = x1i; xyf11,i = y1i; xyf12,i = -f; double, xyf2 = new double3,line; xyf20,i = x2i; xyf21,i = y2i; xyf22,i = -f; Matrix _xyf1 = new Matrix(xyf1); Matrix _xyf2 = new Matrix(xyf2)