单像空间后方交会程序[JAVA版]

本科实验报告

课

名程：称单空间像后交会编程实现方姓 学专学 名： 恒松李 ：院学学院地业：测 绘工一班程 ：号20312943

指

导师：陈老 强0251年 5 月 20 日

目

录

、 一业作任务................................................................................................................... . 2 -二、计算原 理................................................................................................................... -. 2 、三算 法流程................................................................................................................... . 6 四、- 源程序........................................................................................................................- 7 、 五算计结果....................................................................................................................- 7 六 结、分果析................................................................................................................... . - 7七、心得 体与..会............................................................................................................. -.7 、 附页八........................................................................................................................... - .7 aJva 程.序.................................................................................................................................................................. - 7-

-1-

一

、知条已：

件作业

任务

摄影机主距 f=153

.4mm，20=x0，0y=0 ,片像比例尺 为:41000，0四有对的点像点坐标与相的应面地 坐标下表如 。点 号 21 3 4点坐像标 (xmm )-6.15 8-3.450- 41.78 1.06 y4mm)( -8699 82.21.-7 663. 464.3 (m) X6539.481 7361.0833 190097.40426. 45 地坐面标Y m( 2)52373.2 33241.512493 49.8 3039.81 Z1m( )19521.77 82.692 36.58 7057.31

单像以间后空方交方法会求，该像解片的外位元素方 。、二 算计原

1. 理取已知数据。从航获摄料资查中取均航高与平摄机影主距；获控制点的地面取量测坐并转标换 地面摄测为标。 2. 坐测控量点制的像点坐并作系标误差统改。 3正 确定.知数未的初值始 。竖在直摄且地影面制点控大体对称布分的况情下，按如 下法方定初确值， 即

始0X S

 X Y

,

0 Sn

Y ,

Zn

S

0

mf

 Z n

1 0  0    0

0中： m式 为影摄比尺分母例 ； 为控制点n个 数.4 用三个元素角初始值的下式按算计方向各余值弦，组成旋转阵 矩R

a1 a2R  b 1b 2  c c2 1a 3b3   c 3 

其

中：

2-

-



a 1 cso cos  sni isn isn  a   oc ss in  sni sin  os c 2  a 3  ins cos  b 1 cos sni b2 oc  cso s b  s i  n  3c1 sni cs o  cs o is  sn

i

n  c 2  si n ins  os cs n i co s c 3 co s c os

5

.逐 点计算像点坐标的近似。利用值未数知近似的值控和制点的面地标坐带入以，下共线方程式

， x  f   y  f a 1 ( X AX S ) 1b(YA  YS ) c1 Z(A  ZS a)3( XA X S ) 3 bYA( Y S) 3c( Z A  Z S) 2a (XA X S  )b 2YA( Y S ) c ( Z2A ZS )3 ( X A aX S ) b3(Y  AY ) Sc3 ( Z A Z S)

逐

计点算像坐标的点似值近 x( )、 (y) ， 题将目所给数代入据上式可得点像坐标近的似数据。

逐

点计误算差方程的系数和式常数项组，成误差程式方。

a 11 a 21 a 13 a   41 a5 1  61a a 1 7 81a 1a2 a22a 3 a24 2a52a 2 a72 a628a1 3 a3 23a3a 43a5 3a 6 37a3a8 a13 42a a34 444 aa54 6a4 a47 8a a14 525aa 5 345 a5a5 a5 a67 585 a16 aa26  a3 6 a6 4 a56  a66  7a6 8a 6



.

A68 *6

其中

-3-

f a1 1   Ha1 2 0 xa 1 3 H  a41  f1 (  xy a 15  f  a61  y a21  0  f a 22 H y  a 2  3H  a  xy  2 f4  a2 5 f ( 1  a 26  x

他项类其似得可

。x2

) f

2y

) f2

2

常由数项算计公：式a

(X X S ) 1b( Y  iYS)  c 1 Z i( Z )S  xi  f1i  a 3( X  iX S) b3 Yi( Y S )  c (3Z  Zi S )  xl  iLi    l y  iy  af2( Xi X S )  b2 ( Yi YS ) 2 (cZ i  ZS )  i 3a( X i X S )  b 3(Y iY )  cS 3 (Z i  ZS)  (i  1 ,,32,4)

将题所目数据给入上式代可，得常项数 7. 。计算法方程的是数系矩阵 A A 和T常项数 TA L， 成法组方程。

式V

A X 

L8. 法方解，求得外方位程元的素正数 d改 SX, dY S dZ, S d ,, d, d 。

X  A(TA) 1 TAL

.9用前 次代取迭得的似近，值加次本迭代改正数，计算的方位元素外新值。的

XKS  X K S  1X SdK ,Y SK Y K 1 SdY S ,K ZS K  ZK S  d1Z S

KSK  S 1 Kd SK ,SK  SK 1 dS K  S, K SK 1 d SK

式：中 代K迭表次数代

。4-

-10 将求.得外的位元方素改正数与规的限差比较，定若小限于差则，迭代束结否。则新的近用似重值 复4 9 满~足求要止为 1。.1 误由差方程的计算式：

0vx  a1X s  1a2 1s Y 1aZ s3 a4 1   15a   1a6   x  x  0 vy a21 Xs a22 Ys a 23Z s  24a  a25  a2 6 y y

以及

单位权误差的计中算

式

0 （式中：

r 示多表观余测数 ）及以平差 中6个参数 的因数阵协

VT P V

rQx x

( T AA)

1最

得到终参 数Xi 的 中差误

为 x

 0 Qx x

i

i i

-

5-

一

、法流算

程

输入始原据

数

点像坐标算计系统，差改误正

定外方位确元素始值初

成组旋矩阵 转

R

是

结 束并 提 示 错误信 息

逐点组误成差方程式法并

否化

迭

代次 数小 于 n

否所有

像点否完

是

解法方程

，求方位外元素改正数

计算改

后的正外方位素元

外

方元素改正位是否小数限于

差

是计

中算误差输，出成果，束结

-6

-二、

三、

源 序 计程算果结表 1，最 终算计结

源果程序代请码附见 在经页过三迭次代之后到的最得成终果如表1 所 示 X（s）米

379594.225 5795527

Y

s米）

（72764.64325 2804654

Zs米（

）7527.65987 7780861

（ψ弧）

-度.[1**********] 9597528

8

（弧度ω

）.[1**********] 646751

3

（κ弧度）-

0.[1**********] 971112

、四

果结析

由分计算结可果知拍摄在片照瞬间摄，影中在地心面摄影量坐标测中系的坐标（为3995.7542528， 7472646.325，75872.685988）（ 单：位） 米，航倾向ψ为-角0.009378弧 ，旁向倾度角为ω0 002.114弧 ，像片度角κ旋 0为.607785弧度。 表2精度评，定果

结参数 中误

差.[1**********]697 1.[**************]10.5083 [1**********] 618.[1**********]765E4-41.681 [1**********]5E-34 .[**************]0-E

5

Xs

YsZ s

ψ ω κ、五心 与得体会

过通这大次作业，我的得体会以心几点下：第 ：编一程实确一件是诚非艰的难情事即使整，单像空个后方交间实会的原验理非简常，单晰 明了，清 是要但用程把编它示出来还表一是件常非艰的难事情我遇到的最。的大难是点矩在阵逆的求 骤上步为了，解决这问题。个我了用天的时间，还两请了数教学院学的同学，序程里逆求阵的方法矩是 来自位同学的。这次编程那让也对我线性于代中的数求矩逆阵的问题解理的加更深。刻 第：在发现二问，解决题问的题过程中我，认识到了自也己基础识知不的足对于 Java ，言的 语熟悉让不在我多很常见的b u g前面能无为力。以要后加编强程言语学的习。第 ：很三东多看西似单，简其实节处才体现出细能力。小把的细处理好节，才有能的作大。品 六 附、

页1

.J ava程 序 下第一部以分主程序为i

mportj vaauti..l*;i mpot jrvaama.t.h*

;-7-

uplbi clasc si l

p{ubli csattic voi mdinaSt(irng][a gs){ hren hgnegiFnd=ew nhegn)(; odblue duboe

l[]x[]={-{0.0615,-8.00698}9{,0-0534.0,0.0221}8{-0.,04187,-.00763},6{.001064,.064403}};X[ ][]={{36859.412,275.323,2951.17,{376}31.08,1324351,.27869.}{3,9100.79,4923.49,8382.65 0},{44260.45,30139.1,8775.3}1;}int ,j,niu=m;

0d

obleuX0 [=new]doub le6]; dou[lbe f0=.51324; duboel a=1/04000.;

0doub

l Re[]]=new[ oduleb3[]3];[d oulbe A[]][=nwedoubl e[]8[]6 ;oudbe Ll]=n[ew oubdel8];

[dobuele vaulaetx_[=n]ewdo uleb8][ ;dobule AT_[][=n]wedo ble[6]u[8;] dubloe p=3.i

1

[1**********]932;

d7ublo ATeA_[[]]=ewndo blue6[]6[];

oubld inevesrATe_A][[]n=ew duobel[6]6];[/i/nervesT_A=AehngFnd.inveisr(eT_AA,6); dubole esulrt_[v=n]e woubld[e8;/]/用来放存于用计单算位中误差权v值的do ubl seandtrda][=en doubwe[6l; ]duole rbsuetl=0;/F用来/存v放的置转矩阵do bue rlesut1l][][=ewndoub el[][6];8d ubleo rseult[2]new= odulb[e]6 ;oubldesu mXY[Z=n]ew douleb[]3

System;o.tu.rpintnl"已知(像点坐标：为 )"; or(f=0;ii

jels

Seystmeou.tp.rnilnt"y("(i++1)"+= +""\t"x[+i[]j]);

}

}

System

out.p.rntiln"已(地面四个点的知坐为：标); "ofr(i=0i

-

}

lees

if(

=j1=) }{ }{S ysemtou.tpri.ntl(n"Y+"i(1+)+"= +""\"+tXi][[j);]

lees Ssyet.oumtp.inrtnl"(Z+(i+1")"+ "=+\t""X+[]i[j);

]

}

}

fo

r(=0jj;3;j

+

fo(ri=0i

}{

smXuZY[j+]X=[i]j][;}

froi=0(i;

X[0]=iumsXYZ[]/i;

4X[i0=1]a*/+sfumYZX2][4.0/ d;o{

[5]);

R[0][0]Ma=th.oc(sX[03])M*ah.toc(Xs[50)-Ma]thsi.nX(03[]*)ath.Misn(0X[])4Math.*si(nX

00[5

]);

[0][1]R-=ath.cMo(sX0[])*3Mah.tsi(X0[5])nMath-.ins(0X[]3)*Mah.tisn(X[04)*Math.co](X R[0s[]2]-M=ahtsin(.0[X3]*M)at.chs(oX04][) R;[1[0]=]atMhc.os(X04])[*Mth.ais(X0n[]5; )[1]R[1]=Math.cos(0[4X])M*tha.oscX(05])[ R;1[[]]=2-Math.in(sX[40]);

[5]);

R[2

[0]=Ma]t.hsinX(03][*M)ta.coh(s0X5[]+)Mthac.o(s0X3[)]*Mtha.si(Xn[4])*Mat0h.isn(X0

[5])0;

R

2[[1]]=-athM.ins(0X3[)]M*tahsin(.X0[5])+aMt.hcosX0[(3])Mat*.sih(Xn0[4]*)atM.cho(X sR[2[]2=Math]c.o(X0s3[])*Mtah.co(s0[X])4;

0[X]))/2R(0[]2]*([X0][[0]-X0[0])+[1R[2]*]X[(0[]1-X0[]1]+R)2[[]2*](X0][2][X-0[])2)

;eva

ulae_x[0]=-f*t(R[][0]*0([0X]0]-[0[0]X+R[)1[]]*(X[00[1]-X]0[1)+][2R[0]](X[*]0[2]

X0[-2]))(R/[0][]*(X20[][]-X0[0]0)R[+1]2][(*X[0[1]-]0[X])+1R2][[2*(][X0][2-]X[02))];

e

avuate_x[1]l=f-(R*[0[]]*(X[1][00-]0[X])+0[R]11][*X[0(][]1-X[1]0+)R2[[]1](X[*0]2[]-

veaulteax[2]_=-*(fR0][[0](*[1X]0[-X0][]0)R[+1[]]0*X([]11[-X][01]+R)[2][0*]X[1([]2]--9

X02[))/]R[0][2(*]X[(][10-]X0[])0+R1][[2*]X([][1]1X0-[1)]R[2][+]*(X[12]2[-]0[2]X);) evluaatex_[]=3f*-R(0][[1*(X[1]]0[-]X0[])+0[R1]1][*([X1[]1]-0X[1)]R+[][21](X*[][1]-2

0[X2])/)([0][2]*(R[1X]0]-[0[X0])R+[]1[2](X*1][[1]X-0[])1+[2R[2]]*(X[1[2]-X0]2[))]

;

0X[])2)/(R0[][]*2([2X[]0-]X[00])+R[]12[](*[2X][1-X]01][)+R[2[]2*(X[2]][]2X-0[]));2

veluatae_[x]4=f-*R[0]([]*(X02][0][X0[-]0+)R[1][]0*X[(][21]-X[10)+]R2[][]0*X[2([2]]-

X0[

2]))/(R0[]2[*(]X[2]0[-]0X0][)R[1][2]*(X[+2]1[]-0X[1)]R[2+]2][(X[*2[]2-X0][2]));

eval

auet_[x]5=-*(R[f]0[]*(X12][[]-X00[0]+R)[][11](X[2]*1][-X0[])+1R[2[]1](*[X][22-]

X

[0]2)/)([R0[]2](*X[]30[]X0[-]0+R[)]1[2](X[3*[]]1X-01][+R[2])[2](X[*]3[]2X0-2][);)

eavulae_tx[]=6f*-([R0]0[*(]X[][30]-0X0][)R+[]10[*]([X][3]1X0[1])+R[-]2[0*(][3][X2]-

0X2[)]/)([R0][2](*X[3[]0-X0[]0])+[R1]2]*(X[3[][]-X101])[R+[]2[]2*X[3]([2-]X02[]);)f ori=(;i

{ev

luate_xa7[]-f=*([0R]1[*]X(3][[0-]0X0][+)R1][[]1(X*[3][]1-0X[1]+)R[2[1]*(X[3][2]]

-2

]*([Xi][]-X011]

[)+[R]22]*([[X][2]iX0-[2))]

;/

*11*/Aa2*i[[0]=]R([0[]]0*+f[0R[2]*eval]utea_[2x*])/(iR[0[2]*(X[i]][0-X0[0])]+R[1]

[

2]*([Xi[]]1-0X1][)R[+]22]*(X[[i][2]-0X2[]);)

*a1/*/A[2*2i[1]=(R][]10]*f[+[1R[]]2*veluaae_tx2[i*)/(]R[0[]2](*[iX[]]-00X[]0)R+1[[

]]*(2Xi][1[-]X[10)+][R]22][*X[i(]2][-0[X]2);

)/

*a13/A*[*2i[2]]=R([2[]0*f+R[]2[]]*e2aluavte_x2[*])i(R/[]02[]*X[i](0[]-0[X0)+R[]][1

1[]2[]*([Xi[]]-X1[0]1)R[+]22]*[X[(]i[]2-0[2X]))

/*a2;*/A12*i+[]10]=([[R][10*f]+R0][2]*eva[ulat_ex[*i+2])/(1[R][2]*(0[Xi][]-X000[)+R]

[

]12[*]X(i][[]1-X0[1)]+[2R]2[]*X[i(][]2-0[X]2))

;*a2/2*/[A2*i+1[]1=](R1[[]1*f+R]1[[2]]e*valuatex[_*i+2]1/)([R][2]0*([i][0X]-0[0X)]R

+1][[2*](X[][1i]-X0[1])+[R2]2][*X(i][2][-0[2]X);)

*a/3*2A/[*i2+]12[=(R]2][[]1*+Rf2][[]2*veluataex[_2*i+]1)/([R0[2]*]X[i]([0]-0[X0)+]R

*

i]Mat*.cho(Xs[0])5-ealvateux[_2i+*1]*Mta.sinhX0([]5)+f)M*athcos.X(05[])*Mat).hcosX(0[4 );]

*a/1*/4A2[*i[3]=]valeateux_2*[i1]+*Mat.shni(X[0]4-)ev(auate_x[l2*i/]*f(eavulae_x[t

2

]5)+vaeulat_e[x2i+*1]*aMthc.s(X0[5o)]);

/

*15*aA/[*i]24]=[-*faMht.sn(Xi0[])5-valeateu_[x*2]if*(/evaulate_[x*i2*]Mah.sin(tX[ 0*/a61*/[A2i][5]*=veauale_x[2ti+1]*

;

_xe[*2]i*Mahtco.(sX0[]5-e)alvauet_[x*i21]*M+at.sin(X0h5[)])f-M*thasi.n(X[5])0)*Math.cso-

10 -

/a*4*2A/[2i+1][*]=3-*e1avulte_ax2[i**M]ah.sin(X0[4]t)(-vaelaue_xt2*i[+1]/*(evfalau

t

X0([]);4 /*a52/*A2[i*1][+]=-1*f4M*ta.coh(X0[s5)]-vealuae_x[2*it+]/f*1e(valauetx[2*_i]*Mat.s h/*a62/*[A2i+*1][5]=eva-ualetx[_2i*;

]

n(Xi05])+e[vluataex_2*[+1]iMath.c*s(o0[5X)); }

]//

行进常数项的始初化f roi=(;0i4;i

}

[L2i*+1]x=i][1[]-veluate_ax2*i[1+]

;/

/A的转矩阵 置fo(i=0;i

for(j=0;

j

//实现与AAT相 i乘ntk= 0;

fo

(ir0=i

for+(=0;j

+

fr(o=i;i0

A_ATi][j[]=;

f0rok=0;(

fo

(r=j;0

/

得到A/*A的T矩阵逆放在存invreesTA_A,中此采用面处向象的方法 对niverseT_A=AhegFnid.nivernes(ATA,6);

_

AT_[iA][]+kA_T=[i[j]*A[]j][k]

//实;现矩阵TAA[6]_[]6与T[6][A8的相],乘果存结在res放lut[6][1]8中f ori=(0;

;fr(o=0;i

orf(k=0;

fo+(r=j;0j

//

实现ersutl[16][]与l88[]的乘相得，结到放果re在ults[2]6中for i(0=;

i

esultr[i][k]1+=ivnrseeT_AAi][j]*A[T[j][k_];

for

(j0=j

)un++;m {

erustl2[i+]r=suet1l[][i]*j[jL;]

fo

(r=i0i

-

1 1

-}

S

syte.omut.rintpln("\n+"进"行第+n"m+u"迭代得次X到,Ys,sZ,sψ， ，ω改κ数分别正：")为 ;fro(i=;0i

"sult2e[]*502265.60>)6);

w}hiel(athMabs(resu.l2[3t*2]6260.5)>0||6atMhab.s(rsulet24[*20]66520)>.6|M|th.abasr( ySsemto.t.upint(r\"n"+满足"差条件限束循结环，最结终为果：+""n\") ;ysSet.omt.uriptln(n"Xs+"\"t+"Ys"+""\"t+"sZ"+\"t""+"ρ"+t\+"ω"+""\"t+"k); for("i=;0

Sy

ste.motu.rpitnX0(i]["+t\)"

;////

前之代码所为求方外位素的元个值 以下代码六所为单位

求

权中误差f r(i=o;0

orf(i0;i

re

slu_v[ti+=A[]i[]]*jreulst2[j];

}for(intq 0=q;6;q

{resul

_t[v]q-L[q=];

fro(itnm =0;

resulFt=(re+ultsv[m_*]rsuet_lvm[)]

r;suetFl=aMth.sqr(resultF/t2;) of(rntip=0;p 6

;

Syset.outmp.rntlin("\n""+各求所的中量误差")为; Sytsemo.u.ptinr("tXs"+"t\"+"sY+""\t+"Z"s+"\t"+"ρ"+"\t""+"ω""\+t"+"k+""\n); "ySstmeout.pr.itn(tsnaard[0]d"\t"+stan+drad1]+"[\"ts+tadandr2[+]"\"+tsandtrd[a]+3"t\

}"+sta

dand[4]r"\+"+sttnadrda[]+5"\n); }

"以下

二第分为计部逆算矩阵的向面对程象

p序ubli clcass heng{p ulicbd oulb[]e[ i]vnres(edoble[][u ]ain,tn ) {oudle bb][][n=e wodbue[n]l2[*n];

1- -2

oudlebc[ []=n]ew duboe[l][n]n d;oube led1,tyinhiz,b; bin it,,jk,u ;or(i=0;f

;fo

(rj0=;

}

f

r(i=o;0

{forj(0;=

b[i]j[]=a[]ij[]

;

}

f

roj=0(j;n;

orf(=i;0

{i

f(bi[][]=i0={)

f

orj(=ij

{i

f([bj][i]!0){ }=

htsi.tmp(e[i]bb,[j,n);

}]

}

fo

r(=ik+1k

)yin

hz=i-1)*b([][k]i/[bi[i]] fo;(r=u;0

h

}}

det1=

tih.fus(n,a)n; i(fde1=t=){ }0 ySset.mou.trpitlnn"此(矩不阵可逆");

f(iet1d!=)0

{for

i(=0i

[ib[]]=b[i][j]jb/;b

}

or(fi=-1;i>n;0-i-){ bb=[b][i]k

;

orfk(0;=

b)k[[]]=u[kb[]u-bb]*bi[]u[;

-]13 -

}

}

}

}fr(o=0;ii

f

o(rj0=;

c[

][j]=i[i][b+jn];

}

}r

ternu ;c

pu

lbci vodi etp(dmoble ua[],adoubl eb[b,]intn {) int ; idoulebt mp1;e

for(=i0;

tmp1=ea[a]i;aa[]i=b[b]i;b[i]=tbme1p;

}

pu

lbi doucbe lufndou(le baray[][r,in]t ){n itni ,ij,k,u; jodble uetd11=y,in;

f

or(i=0;iiin;ii++)

i

f(rara[yii]i[i]==0)

{for(j=ji;ij

{

i(afrra[jj]y[ii!]0=) {}

t

empa(ray[iri],raar[yj],jn);

}

}

or(kfii+=1k;n;k++){

n

}

f

r(oii0;=i

{

etd=det1*1array[i]i[ii;

} }

r]eturndet1();

}

实验

程序结如下图：

果 14 -

-

- 15 -

本科实验报告

课

名程：称单空间像后交会编程实现方姓 学专学 名： 恒松李 ：院学学院地业：测 绘工一班程 ：号20312943

指

导师：陈老 强0251年 5 月 20 日

目

录

、 一业作任务................................................................................................................... . 2 -二、计算原 理................................................................................................................... -. 2 、三算 法流程................................................................................................................... . 6 四、- 源程序........................................................................................................................- 7 、 五算计结果....................................................................................................................- 7 六 结、分果析................................................................................................................... . - 7七、心得 体与..会............................................................................................................. -.7 、 附页八........................................................................................................................... - .7 aJva 程.序.................................................................................................................................................................. - 7-

-1-

一

、知条已：

件作业

任务

摄影机主距 f=153

.4mm，20=x0，0y=0 ,片像比例尺 为:41000，0四有对的点像点坐标与相的应面地 坐标下表如 。点 号 21 3 4点坐像标 (xmm )-6.15 8-3.450- 41.78 1.06 y4mm)( -8699 82.21.-7 663. 464.3 (m) X6539.481 7361.0833 190097.40426. 45 地坐面标Y m( 2)52373.2 33241.512493 49.8 3039.81 Z1m( )19521.77 82.692 36.58 7057.31

单像以间后空方交方法会求，该像解片的外位元素方 。、二 算计原

1. 理取已知数据。从航获摄料资查中取均航高与平摄机影主距；获控制点的地面取量测坐并转标换 地面摄测为标。 2. 坐测控量点制的像点坐并作系标误差统改。 3正 确定.知数未的初值始 。竖在直摄且地影面制点控大体对称布分的况情下，按如 下法方定初确值， 即

始0X S

 X Y

,

0 Sn

Y ,

Zn

S

0

mf

 Z n

1 0  0    0

0中： m式 为影摄比尺分母例 ； 为控制点n个 数.4 用三个元素角初始值的下式按算计方向各余值弦，组成旋转阵 矩R

a1 a2R  b 1b 2  c c2 1a 3b3   c 3 

其

中：

2-

-



a 1 cso cos  sni isn isn  a   oc ss in  sni sin  os c 2  a 3  ins cos  b 1 cos sni b2 oc  cso s b  s i  n  3c1 sni cs o  cs o is  sn

i

n  c 2  si n ins  os cs n i co s c 3 co s c os

5

.逐 点计算像点坐标的近似。利用值未数知近似的值控和制点的面地标坐带入以，下共线方程式

， x  f   y  f a 1 ( X AX S ) 1b(YA  YS ) c1 Z(A  ZS a)3( XA X S ) 3 bYA( Y S) 3c( Z A  Z S) 2a (XA X S  )b 2YA( Y S ) c ( Z2A ZS )3 ( X A aX S ) b3(Y  AY ) Sc3 ( Z A Z S)

逐

计点算像坐标的点似值近 x( )、 (y) ， 题将目所给数代入据上式可得点像坐标近的似数据。

逐

点计误算差方程的系数和式常数项组，成误差程式方。

a 11 a 21 a 13 a   41 a5 1  61a a 1 7 81a 1a2 a22a 3 a24 2a52a 2 a72 a628a1 3 a3 23a3a 43a5 3a 6 37a3a8 a13 42a a34 444 aa54 6a4 a47 8a a14 525aa 5 345 a5a5 a5 a67 585 a16 aa26  a3 6 a6 4 a56  a66  7a6 8a 6



.

A68 *6

其中

-3-

f a1 1   Ha1 2 0 xa 1 3 H  a41  f1 (  xy a 15  f  a61  y a21  0  f a 22 H y  a 2  3H  a  xy  2 f4  a2 5 f ( 1  a 26  x

他项类其似得可

。x2

) f

2y

) f2

2

常由数项算计公：式a

(X X S ) 1b( Y  iYS)  c 1 Z i( Z )S  xi  f1i  a 3( X  iX S) b3 Yi( Y S )  c (3Z  Zi S )  xl  iLi    l y  iy  af2( Xi X S )  b2 ( Yi YS ) 2 (cZ i  ZS )  i 3a( X i X S )  b 3(Y iY )  cS 3 (Z i  ZS)  (i  1 ,,32,4)

将题所目数据给入上式代可，得常项数 7. 。计算法方程的是数系矩阵 A A 和T常项数 TA L， 成法组方程。

式V

A X 

L8. 法方解，求得外方位程元的素正数 d改 SX, dY S dZ, S d ,, d, d 。

X  A(TA) 1 TAL

.9用前 次代取迭得的似近，值加次本迭代改正数，计算的方位元素外新值。的

XKS  X K S  1X SdK ,Y SK Y K 1 SdY S ,K ZS K  ZK S  d1Z S

KSK  S 1 Kd SK ,SK  SK 1 dS K  S, K SK 1 d SK

式：中 代K迭表次数代

。4-

-10 将求.得外的位元方素改正数与规的限差比较，定若小限于差则，迭代束结否。则新的近用似重值 复4 9 满~足求要止为 1。.1 误由差方程的计算式：

0vx  a1X s  1a2 1s Y 1aZ s3 a4 1   15a   1a6   x  x  0 vy a21 Xs a22 Ys a 23Z s  24a  a25  a2 6 y y

以及

单位权误差的计中算

式

0 （式中：

r 示多表观余测数 ）及以平差 中6个参数 的因数阵协

VT P V

rQx x

( T AA)

1最

得到终参 数Xi 的 中差误

为 x

 0 Qx x

i

i i

-

5-

一

、法流算

程

输入始原据

数

点像坐标算计系统，差改误正

定外方位确元素始值初

成组旋矩阵 转

R

是

结 束并 提 示 错误信 息

逐点组误成差方程式法并

否化

迭

代次 数小 于 n

否所有

像点否完

是

解法方程

，求方位外元素改正数

计算改

后的正外方位素元

外

方元素改正位是否小数限于

差

是计

中算误差输，出成果，束结

-6

-二、

三、

源 序 计程算果结表 1，最 终算计结

源果程序代请码附见 在经页过三迭次代之后到的最得成终果如表1 所 示 X（s）米

379594.225 5795527

Y

s米）

（72764.64325 2804654

Zs米（

）7527.65987 7780861

（ψ弧）

-度.[1**********] 9597528

8

（弧度ω

）.[1**********] 646751

3

（κ弧度）-

0.[1**********] 971112

、四

果结析

由分计算结可果知拍摄在片照瞬间摄，影中在地心面摄影量坐标测中系的坐标（为3995.7542528， 7472646.325，75872.685988）（ 单：位） 米，航倾向ψ为-角0.009378弧 ，旁向倾度角为ω0 002.114弧 ，像片度角κ旋 0为.607785弧度。 表2精度评，定果

结参数 中误

差.[1**********]697 1.[**************]10.5083 [1**********] 618.[1**********]765E4-41.681 [1**********]5E-34 .[**************]0-E

5

Xs

YsZ s

ψ ω κ、五心 与得体会

过通这大次作业，我的得体会以心几点下：第 ：编一程实确一件是诚非艰的难情事即使整，单像空个后方交间实会的原验理非简常，单晰 明了，清 是要但用程把编它示出来还表一是件常非艰的难事情我遇到的最。的大难是点矩在阵逆的求 骤上步为了，解决这问题。个我了用天的时间，还两请了数教学院学的同学，序程里逆求阵的方法矩是 来自位同学的。这次编程那让也对我线性于代中的数求矩逆阵的问题解理的加更深。刻 第：在发现二问，解决题问的题过程中我，认识到了自也己基础识知不的足对于 Java ，言的 语熟悉让不在我多很常见的b u g前面能无为力。以要后加编强程言语学的习。第 ：很三东多看西似单，简其实节处才体现出细能力。小把的细处理好节，才有能的作大。品 六 附、

页1

.J ava程 序 下第一部以分主程序为i

mportj vaauti..l*;i mpot jrvaama.t.h*

;-7-

uplbi clasc si l

p{ubli csattic voi mdinaSt(irng][a gs){ hren hgnegiFnd=ew nhegn)(; odblue duboe

l[]x[]={-{0.0615,-8.00698}9{,0-0534.0,0.0221}8{-0.,04187,-.00763},6{.001064,.064403}};X[ ][]={{36859.412,275.323,2951.17,{376}31.08,1324351,.27869.}{3,9100.79,4923.49,8382.65 0},{44260.45,30139.1,8775.3}1;}int ,j,niu=m;

0d

obleuX0 [=new]doub le6]; dou[lbe f0=.51324; duboel a=1/04000.;

0doub

l Re[]]=new[ oduleb3[]3];[d oulbe A[]][=nwedoubl e[]8[]6 ;oudbe Ll]=n[ew oubdel8];

[dobuele vaulaetx_[=n]ewdo uleb8][ ;dobule AT_[][=n]wedo ble[6]u[8;] dubloe p=3.i

1

[1**********]932;

d7ublo ATeA_[[]]=ewndo blue6[]6[];

oubld inevesrATe_A][[]n=ew duobel[6]6];[/i/nervesT_A=AehngFnd.inveisr(eT_AA,6); dubole esulrt_[v=n]e woubld[e8;/]/用来放存于用计单算位中误差权v值的do ubl seandtrda][=en doubwe[6l; ]duole rbsuetl=0;/F用来/存v放的置转矩阵do bue rlesut1l][][=ewndoub el[][6];8d ubleo rseult[2]new= odulb[e]6 ;oubldesu mXY[Z=n]ew douleb[]3

System;o.tu.rpintnl"已知(像点坐标：为 )"; or(f=0;ii

jels

Seystmeou.tp.rnilnt"y("(i++1)"+= +""\t"x[+i[]j]);

}

}

System

out.p.rntiln"已(地面四个点的知坐为：标); "ofr(i=0i

-

}

lees

if(

=j1=) }{ }{S ysemtou.tpri.ntl(n"Y+"i(1+)+"= +""\"+tXi][[j);]

lees Ssyet.oumtp.inrtnl"(Z+(i+1")"+ "=+\t""X+[]i[j);

]

}

}

fo

r(=0jj;3;j

+

fo(ri=0i

}{

smXuZY[j+]X=[i]j][;}

froi=0(i;

X[0]=iumsXYZ[]/i;

4X[i0=1]a*/+sfumYZX2][4.0/ d;o{

[5]);

R[0][0]Ma=th.oc(sX[03])M*ah.toc(Xs[50)-Ma]thsi.nX(03[]*)ath.Misn(0X[])4Math.*si(nX

00[5

]);

[0][1]R-=ath.cMo(sX0[])*3Mah.tsi(X0[5])nMath-.ins(0X[]3)*Mah.tisn(X[04)*Math.co](X R[0s[]2]-M=ahtsin(.0[X3]*M)at.chs(oX04][) R;[1[0]=]atMhc.os(X04])[*Mth.ais(X0n[]5; )[1]R[1]=Math.cos(0[4X])M*tha.oscX(05])[ R;1[[]]=2-Math.in(sX[40]);

[5]);

R[2

[0]=Ma]t.hsinX(03][*M)ta.coh(s0X5[]+)Mthac.o(s0X3[)]*Mtha.si(Xn[4])*Mat0h.isn(X0

[5])0;

R

2[[1]]=-athM.ins(0X3[)]M*tahsin(.X0[5])+aMt.hcosX0[(3])Mat*.sih(Xn0[4]*)atM.cho(X sR[2[]2=Math]c.o(X0s3[])*Mtah.co(s0[X])4;

0[X]))/2R(0[]2]*([X0][[0]-X0[0])+[1R[2]*]X[(0[]1-X0[]1]+R)2[[]2*](X0][2][X-0[])2)

;eva

ulae_x[0]=-f*t(R[][0]*0([0X]0]-[0[0]X+R[)1[]]*(X[00[1]-X]0[1)+][2R[0]](X[*]0[2]

X0[-2]))(R/[0][]*(X20[][]-X0[0]0)R[+1]2][(*X[0[1]-]0[X])+1R2][[2*(][X0][2-]X[02))];

e

avuate_x[1]l=f-(R*[0[]]*(X[1][00-]0[X])+0[R]11][*X[0(][]1-X[1]0+)R2[[]1](X[*0]2[]-

veaulteax[2]_=-*(fR0][[0](*[1X]0[-X0][]0)R[+1[]]0*X([]11[-X][01]+R)[2][0*]X[1([]2]--9

X02[))/]R[0][2(*]X[(][10-]X0[])0+R1][[2*]X([][1]1X0-[1)]R[2][+]*(X[12]2[-]0[2]X);) evluaatex_[]=3f*-R(0][[1*(X[1]]0[-]X0[])+0[R1]1][*([X1[]1]-0X[1)]R+[][21](X*[][1]-2

0[X2])/)([0][2]*(R[1X]0]-[0[X0])R+[]1[2](X*1][[1]X-0[])1+[2R[2]]*(X[1[2]-X0]2[))]

;

0X[])2)/(R0[][]*2([2X[]0-]X[00])+R[]12[](*[2X][1-X]01][)+R[2[]2*(X[2]][]2X-0[]));2

veluatae_[x]4=f-*R[0]([]*(X02][0][X0[-]0+)R[1][]0*X[(][21]-X[10)+]R2[][]0*X[2([2]]-

X0[

2]))/(R0[]2[*(]X[2]0[-]0X0][)R[1][2]*(X[+2]1[]-0X[1)]R[2+]2][(X[*2[]2-X0][2]));

eval

auet_[x]5=-*(R[f]0[]*(X12][[]-X00[0]+R)[][11](X[2]*1][-X0[])+1R[2[]1](*[X][22-]

X

[0]2)/)([R0[]2](*X[]30[]X0[-]0+R[)]1[2](X[3*[]]1X-01][+R[2])[2](X[*]3[]2X0-2][);)

eavulae_tx[]=6f*-([R0]0[*(]X[][30]-0X0][)R+[]10[*]([X][3]1X0[1])+R[-]2[0*(][3][X2]-

0X2[)]/)([R0][2](*X[3[]0-X0[]0])+[R1]2]*(X[3[][]-X101])[R+[]2[]2*X[3]([2-]X02[]);)f ori=(;i

{ev

luate_xa7[]-f=*([0R]1[*]X(3][[0-]0X0][+)R1][[]1(X*[3][]1-0X[1]+)R[2[1]*(X[3][2]]

-2

]*([Xi][]-X011]

[)+[R]22]*([[X][2]iX0-[2))]

;/

*11*/Aa2*i[[0]=]R([0[]]0*+f[0R[2]*eval]utea_[2x*])/(iR[0[2]*(X[i]][0-X0[0])]+R[1]

[

2]*([Xi[]]1-0X1][)R[+]22]*(X[[i][2]-0X2[]);)

*a1/*/A[2*2i[1]=(R][]10]*f[+[1R[]]2*veluaae_tx2[i*)/(]R[0[]2](*[iX[]]-00X[]0)R+1[[

]]*(2Xi][1[-]X[10)+][R]22][*X[i(]2][-0[X]2);

)/

*a13/A*[*2i[2]]=R([2[]0*f+R[]2[]]*e2aluavte_x2[*])i(R/[]02[]*X[i](0[]-0[X0)+R[]][1

1[]2[]*([Xi[]]-X1[0]1)R[+]22]*[X[(]i[]2-0[2X]))

/*a2;*/A12*i+[]10]=([[R][10*f]+R0][2]*eva[ulat_ex[*i+2])/(1[R][2]*(0[Xi][]-X000[)+R]

[

]12[*]X(i][[]1-X0[1)]+[2R]2[]*X[i(][]2-0[X]2))

;*a2/2*/[A2*i+1[]1=](R1[[]1*f+R]1[[2]]e*valuatex[_*i+2]1/)([R][2]0*([i][0X]-0[0X)]R

+1][[2*](X[][1i]-X0[1])+[R2]2][*X(i][2][-0[2]X);)

*a/3*2A/[*i2+]12[=(R]2][[]1*+Rf2][[]2*veluataex[_2*i+]1)/([R0[2]*]X[i]([0]-0[X0)+]R

*

i]Mat*.cho(Xs[0])5-ealvateux[_2i+*1]*Mta.sinhX0([]5)+f)M*athcos.X(05[])*Mat).hcosX(0[4 );]

*a/1*/4A2[*i[3]=]valeateux_2*[i1]+*Mat.shni(X[0]4-)ev(auate_x[l2*i/]*f(eavulae_x[t

2

]5)+vaeulat_e[x2i+*1]*aMthc.s(X0[5o)]);

/

*15*aA/[*i]24]=[-*faMht.sn(Xi0[])5-valeateu_[x*2]if*(/evaulate_[x*i2*]Mah.sin(tX[ 0*/a61*/[A2i][5]*=veauale_x[2ti+1]*

;

_xe[*2]i*Mahtco.(sX0[]5-e)alvauet_[x*i21]*M+at.sin(X0h5[)])f-M*thasi.n(X[5])0)*Math.cso-

10 -

/a*4*2A/[2i+1][*]=3-*e1avulte_ax2[i**M]ah.sin(X0[4]t)(-vaelaue_xt2*i[+1]/*(evfalau

t

X0([]);4 /*a52/*A2[i*1][+]=-1*f4M*ta.coh(X0[s5)]-vealuae_x[2*it+]/f*1e(valauetx[2*_i]*Mat.s h/*a62/*[A2i+*1][5]=eva-ualetx[_2i*;

]

n(Xi05])+e[vluataex_2*[+1]iMath.c*s(o0[5X)); }

]//

行进常数项的始初化f roi=(;0i4;i

}

[L2i*+1]x=i][1[]-veluate_ax2*i[1+]

;/

/A的转矩阵 置fo(i=0;i

for(j=0;

j

//实现与AAT相 i乘ntk= 0;

fo

(ir0=i

for+(=0;j

+

fr(o=i;i0

A_ATi][j[]=;

f0rok=0;(

fo

(r=j;0

/

得到A/*A的T矩阵逆放在存invreesTA_A,中此采用面处向象的方法 对niverseT_A=AhegFnid.nivernes(ATA,6);

_

AT_[iA][]+kA_T=[i[j]*A[]j][k]

//实;现矩阵TAA[6]_[]6与T[6][A8的相],乘果存结在res放lut[6][1]8中f ori=(0;

;fr(o=0;i

orf(k=0;

fo+(r=j;0j

//

实现ersutl[16][]与l88[]的乘相得，结到放果re在ults[2]6中for i(0=;

i

esultr[i][k]1+=ivnrseeT_AAi][j]*A[T[j][k_];

for

(j0=j

)un++;m {

erustl2[i+]r=suet1l[][i]*j[jL;]

fo

(r=i0i

-

1 1

-}

S

syte.omut.rintpln("\n+"进"行第+n"m+u"迭代得次X到,Ys,sZ,sψ， ，ω改κ数分别正：")为 ;fro(i=;0i

"sult2e[]*502265.60>)6);

w}hiel(athMabs(resu.l2[3t*2]6260.5)>0||6atMhab.s(rsulet24[*20]66520)>.6|M|th.abasr( ySsemto.t.upint(r\"n"+满足"差条件限束循结环，最结终为果：+""n\") ;ysSet.omt.uriptln(n"Xs+"\"t+"Ys"+""\"t+"sZ"+\"t""+"ρ"+t\+"ω"+""\"t+"k); for("i=;0

Sy

ste.motu.rpitnX0(i]["+t\)"

;////

前之代码所为求方外位素的元个值 以下代码六所为单位

求

权中误差f r(i=o;0

orf(i0;i

re

slu_v[ti+=A[]i[]]*jreulst2[j];

}for(intq 0=q;6;q

{resul

_t[v]q-L[q=];

fro(itnm =0;

resulFt=(re+ultsv[m_*]rsuet_lvm[)]

r;suetFl=aMth.sqr(resultF/t2;) of(rntip=0;p 6

;

Syset.outmp.rntlin("\n""+各求所的中量误差")为; Sytsemo.u.ptinr("tXs"+"t\"+"sY+""\t+"Z"s+"\t"+"ρ"+"\t""+"ω""\+t"+"k+""\n); "ySstmeout.pr.itn(tsnaard[0]d"\t"+stan+drad1]+"[\"ts+tadandr2[+]"\"+tsandtrd[a]+3"t\

}"+sta

dand[4]r"\+"+sttnadrda[]+5"\n); }

"以下

二第分为计部逆算矩阵的向面对程象

p序ubli clcass heng{p ulicbd oulb[]e[ i]vnres(edoble[][u ]ain,tn ) {oudle bb][][n=e wodbue[n]l2[*n];

1- -2

oudlebc[ []=n]ew duboe[l][n]n d;oube led1,tyinhiz,b; bin it,,jk,u ;or(i=0;f

;fo

(rj0=;

}

f

r(i=o;0

{forj(0;=

b[i]j[]=a[]ij[]

;

}

f

roj=0(j;n;

orf(=i;0

{i

f(bi[][]=i0={)

f

orj(=ij

{i

f([bj][i]!0){ }=

htsi.tmp(e[i]bb,[j,n);

}]

}

fo

r(=ik+1k

)yin

hz=i-1)*b([][k]i/[bi[i]] fo;(r=u;0

h

}}

det1=

tih.fus(n,a)n; i(fde1=t=){ }0 ySset.mou.trpitlnn"此(矩不阵可逆");

f(iet1d!=)0

{for

i(=0i

[ib[]]=b[i][j]jb/;b

}

or(fi=-1;i>n;0-i-){ bb=[b][i]k

;

orfk(0;=

b)k[[]]=u[kb[]u-bb]*bi[]u[;

-]13 -

}

}

}

}fr(o=0;ii

f

o(rj0=;

c[

][j]=i[i][b+jn];

}

}r

ternu ;c

pu

lbci vodi etp(dmoble ua[],adoubl eb[b,]intn {) int ; idoulebt mp1;e

for(=i0;

tmp1=ea[a]i;aa[]i=b[b]i;b[i]=tbme1p;

}

pu

lbi doucbe lufndou(le baray[][r,in]t ){n itni ,ij,k,u; jodble uetd11=y,in;

f

or(i=0;iiin;ii++)

i

f(rara[yii]i[i]==0)

{for(j=ji;ij

{

i(afrra[jj]y[ii!]0=) {}

t

empa(ray[iri],raar[yj],jn);

}

}

or(kfii+=1k;n;k++){

n

}

f

r(oii0;=i

{

etd=det1*1array[i]i[ii;

} }

r]eturndet1();

}

实验

程序结如下图：

果 14 -

-

- 15 -