光束法區(qū)域網(wǎng)空中三角測量作業(yè)報(bào)告

1、 光束法區(qū)域網(wǎng)空中三角測量求像片外方位元素和加密點(diǎn)的地面坐標(biāo) 一、原理：首先給出的原始數(shù)據(jù)像素坐標(biāo)轉(zhuǎn)換成像平面直角坐標(biāo)；然后利用后方交會原理并采用控制點(diǎn)的向平面直角坐標(biāo)和地面坐標(biāo)迭代求解每張像片的外方位元素；再通過前方交會求加密點(diǎn)的地面坐標(biāo)；最后進(jìn)行精度評定。二、計(jì)算流程：1、 求出四個控制點(diǎn)及所有待求點(diǎn)的像點(diǎn)坐標(biāo)（x1,y1）與(x2,y2)；2、 空間后方交會計(jì)算兩像片的外方位元素：根據(jù)編好的計(jì)算機(jī)程序讀取原始文件中的數(shù)據(jù)控制點(diǎn)的地面坐標(biāo)和相應(yīng)的像點(diǎn)像素坐標(biāo)（第一步已經(jīng)轉(zhuǎn)換成像平面直角坐標(biāo)），對兩張像片各自進(jìn)行空間后方交會，計(jì)算各自的6個外方位元素Xs1,Ys1,Zs1,1,w1,k1和X

2、s2,Ys2,Zs2,2,w2,k2；3、 空間前方交會計(jì)算待定點(diǎn)的地面坐標(biāo)：用各片的外方位角元素計(jì)算左右片的方向余弦值，組成旋轉(zhuǎn)矩陣R1與R2，逐點(diǎn)計(jì)算像點(diǎn)的像空間愛你輔助坐標(biāo)(u1,v1,w1)及(u2,v2,w2)；根據(jù)外方位元素計(jì)算基線分量（Bu,Bv,Bw）;計(jì)算投影系數(shù)N1,N2；計(jì)算待定點(diǎn)在各自的像空間輔助坐標(biāo)系中的坐標(biāo)(U1,V1,W1)及(U2,V2,W2)；最后計(jì)算待定點(diǎn)的地面攝影測量坐標(biāo)。4、 進(jìn)行精度評定：利用前方交會求出地面控制點(diǎn)的坐標(biāo)，再與地面實(shí)測坐標(biāo)求差，將兩者差值視為真差值，由這些真差值計(jì)算出點(diǎn)位坐標(biāo)精度。三、基本公式：1、像素坐標(biāo)（I,J）轉(zhuǎn)換成像平面直角坐標(biāo)

3、(x,y)：x = (J- 2136) * 5.19 / 1000 / 1000;y = (1424 - I) * 5.19 / 1000 / 1000;2、后方交會：共線條件方程： 外方位元素系數(shù)矩陣：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、前方交會：計(jì)算旋轉(zhuǎn)矩陣：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;投影系數(shù)：計(jì)算地面點(diǎn)坐標(biāo)：5、 精度評定：四、實(shí)例：使用數(shù)碼相機(jī)拍

6、攝了兩張航測相片Img261和Img262，量測了6對點(diǎn)的影像坐標(biāo)，其中前4個點(diǎn)為地面控制點(diǎn)（GCP）, 后兩個點(diǎn)為加密點(diǎn)。如下表所示。相片參數(shù)：大小為4272（列）*2848（行），像主點(diǎn)位于影像中心，像素大小為5.19micron，相機(jī)的焦距為24mm。.Img261的外方位元素初始值：X=397510.760，Y=3445853.978, Z=1455.153，角元素為0；Img262的外方位元素初始值：X=397513.320，Y=3445979.811, Z=1453.685，角元素為0；請采用光束法區(qū)域網(wǎng)空中三角測量方法，解求每張相片的外方位元素和加密點(diǎn)的地面坐標(biāo)。點(diǎn)名Img261

7、Img262地面坐標(biāo)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？原始數(shù)據(jù)文件：生成結(jié)果文件：程序源代碼：using System.Collections.Generic;using System.Linq;using System.Text;using System.IO;namespace 攝影測量光束法編程作業(yè) class Program static void Main(string args) /原始觀測數(shù)據(jù)：控制點(diǎn)和加密點(diǎn)的像素坐標(biāo)（x1,y1,x2,y2），控制點(diǎn)的地面坐標(biāo)(X,Y,Z) FileStream fs = ne

9、w FileStream("控制點(diǎn)觀測數(shù)據(jù).txt", FileMode.Open, FileAccess.Read); FileStream fs2 = new FileStream("控制點(diǎn)觀測數(shù)據(jù).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); /框標(biāo)像素坐標(biāo)轉(zhuǎn)換成框標(biāo)坐標(biāo)(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; /相機(jī)焦距和相片外方位元素的初始值 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; /計(jì)算旋轉(zhuǎn)矩陣 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、/計(jì)算像點(diǎn)坐標(biāo)近似值，求出常數(shù)項(xiàng)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; /像片的常數(shù)項(xiàng)為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+; /求出外方位元素系數(shù)矩陣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+; /求外方位元素的改正數(shù) 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 ); /相機(jī)焦距和相片外方位元素的初始值

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; /計(jì)算兩張像片的外方位元素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); /計(jì)算旋轉(zhuǎn)矩陣 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); /求待定點(diǎn)的像點(diǎn)坐標(biāo)jx1jj， jy1jj，jx2jj，jy2jj FileStream fs1 = new FileStream("待定點(diǎn)像點(diǎn)坐標(biāo)

35、觀測數(shù)據(jù).txt", FileMode.Open, FileAccess.Read); FileStream fs12 = new FileStream("待定點(diǎn)像點(diǎn)坐標(biāo)觀測數(shù)據(jù).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; /計(jì)算結(jié)果精度評定(實(shí)際精度) /首先用前方交會求地面控制點(diǎn)的坐標(biāo) 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)