关于空洞探测模型的报告论文.docx

关于空洞探测模型的报告论文问题重述本题要求我们利用弹性波在介质中和空气中不同的传播速度，来确定矩形平板内的空洞位置。该矩形平板由均匀介质构成，内有一些充满空气的空洞。在平板两个两边分别等距的设置若干波源，在他们的对边对等地设置同样多的接受器，记录弹性波由每个波源到达对边上每个接受器的时间，要求确定平板内位置空洞的位置，并讨论在同样确定空洞位置的前提下，减少波源和接受器的方法。图如下:图一: X轴从左到右分别为P1P7 Y轴从下到上分别为R1R7摘要本模型因引入概率统计求解而可认为是随机模型：我们对平板中分成的小格用概率统计的方法来评判小格成为空洞的可能性。事实上，也就是对所给的数据及其条件建立一个相对科学的处理方法,类似于拟合的方法,从后面的分析和计算可知这种方法是可行和科学的。主要结论: 第一问的结果在模型的求解中已经给出了空洞的位置图(见模型的求解中第一问的求解).第二问的结果: 1如果去掉横向RS 的波源与接收器,可以确定空洞的位置,但是精确度有所降低 2在同样能够确定空洞位置的前提下,可以减少波源:P3,P5,R3,R5接收器Q4,Q6,S4,S6(见模型的求解中第二问的求解). 本模型有效的消除了测量方法带来的系统误差带来的影响，只要波的密度和分布的均匀性达到足够的要求，结果就可以做得很细，很精确。问题分析我们认为该问题是实际应用中的测量问题，主要通过采用适当且有效的方法对已知数据进行分析和处理，来提取所需的信息。在本题中已知波经过的距离及所用的时间，要求木板中空洞的位置。我们可以用数据拟合的方法对其进行处理,但数据拟合由于受变量数的限制不容易做得很细,很精确.因而我们改用统计的方法.首先我们以方板为研究对象，将方板分为尽可能小的细胞，用细胞状态来描述空洞的存在，即0表示空洞，1表示介质。这样可用元素为0和1 的矩阵来表示空洞的分布和形状。由此建立以计算机模拟为主要思想的模型一，用事件步长法穷举求出空洞的位置。但可能的空洞状态数很多,因此用计算机模拟的难度比较大，其优点是在划分足够细的条件下能精确描述空洞的位置和形状。接着我们以波为研究对象出发，可得出波给予我们的两个信息：在这个波方向上出现的空洞距离及交点信息（交点处出现空洞的可能性），得出空洞可能出现的位置与范围，此为模型二。但由于任意两条直线交点的复杂性，得出的范围不够准确。经过对以上两种模型的研究，我们考虑中和方板和波两方面的因素，即采用微元划分与最大可能性判断相结合，用每个微元出现空洞的可能性大小来决定微元的状态（ 0或1），从而得出空洞的分布，由此我们建立模型三。可以发现，该模型具有较好的可行性及较高的可信度。符号系统已知：方板边长L=240米 时间矩阵T1= tij tij:由Pi发出的弹性波到达Qj的时间0.0611 0.0895 0.1996 0.2032 0.4181 0.4923 0.56460.0989 0.0592 0.4413 0.4318 0.4770 0.5242 0.38050.3052 0.4131 0.0598 0.4153 0.4156 0.3563 0.19190.3221 0.4453 0.4040 0.0738 0.1789 0.0740 0.21220.3490 0.4529 0.2263 0.1917 0.0839 0.1768 0.18100.3807 0.3177 0.2364 0.3064 0.2217 0.0939 0.10310.4311 0.3397 0.3566 0.1954 0.0760 0.0688 0.1042 时间矩阵T2=tij tij :由Ri发出的弹性波到达Sj的时间0.0645 0.0602 0.0813 0.3516 0.3867 0.4314 0.57210.0753 0.0700 0.2852 0.4341 0.3491 0.4800 0.49800.3456 0.3205 0.0974 0.4093 0.4240 0.4540 0.31120.3655 0.3289 0.4247 0.1007 0.3249 0.2134 0.10170.3165 0.2409 0.3214 0.3256 0.0904 0.1874 0.21300.2749 0.3891 0.5895 0.3016 0.2058 0.0841 0.07060.4434 0.4919 0.3904 0.0786 0.0709 0.0914 0.0583 波速：在介质中 v1=2880米/秒 在空气中 v2=320米/秒未知：波长矩阵L1=lij lij:由Pi发出到达j 的弹性波的波长.L2=lij lij:由Ri发出到达Si的弹性波的波长同时用lij 表示由Pi或Ri发出到达j或Si的弹性波直线。 空洞距离矩阵X1= xij xij 表示由Pi到j 的弹性波lij上的空洞距离。 X2= xij xij 表示由Ri到Si的弹性波lij上的空洞距离。时间矩阵T3=tij3 tij3:弹性波lij 经过距离全为介质时所需的时间。 划分间距 a (0a240) 则N=L/a为每边被分成的段数。 微元矩阵T= Tkm ( k,m=0,1,2,N-1) Tkm =km ，k表示微元的横坐标，m表示纵坐标。 微元状态矩阵t= tkm ( k,m=0,1,2,N-1) tkm =o.微元为空洞，tkm，微元为介质 概率矩阵PL1= Plij: Plij: : 由Pi到j 的波lij 单位距离内存在空洞的概率。 PL2= Plij: Plij: : 由Ri到Si的波lij 单位距离内存在空洞的概率。P1= Pij: Pij: : 由Pi到j 的波lij 单位距离内存在空洞的加权概率。 P2= Pij: Pij: : 由Ri到Si的波lij 单位距离内存在空洞的加权概率。概率矩阵PT= Ptkm: Ptkm: :微元Tkm:为空洞的概率。 权矩阵W=wij wij :波lij由于波长不同对微元Tkm:作用的权 波的集合Akm= lij: , lij:经过微元Tkm: 元素个数为nkm 微元集合Bij =Tkm ，Tkm在波lij:上元素个数为nij 模型假设我们认为“弹性波”是指在反射、折射及相互干涉中都没有能量损失的波。其波速仅与波在其中传播的介质有关，且在同一均匀介质中波速不变。我们假设有条件使波源发射的弹性波具有很强的单向性，且波源与接受器的仪器精度足够高。因此我们做出以下假设：1. 表中给出的时间数据客观且真实地反映了实际情况，即所记录的数据准确无误。2. 波在传播过程中沿直线单向传播，且不考虑波的反射、折射以及干涉等现象。3. 波的发射与接收为一一对应的，且波的发射与接收均为瞬时。因此我们认为记录的时间即为波沿直线由发射点到接收点的传播时间，且波的传播时间仅由波在该均匀介质中经过的距离和在空气中经过的距离决定。4. 已知弹性波沿板边缘的传播速度与在介质中的传播速度相同，因此我们同等地看待每个波，不因它们的位置不同而有所区别。5. 由实际情况可知，空洞的分布应为分段连续的，因此不能求出每一点出现空洞的概率，而只能用一个区域内出现空洞的概率来逼近每一点的概率.当把板分成微元后,我们可认为一个微元只可能出现全是空洞或全是均匀介质两状态. 模型的建立与求解第一问的求解：1 微元的划分我们选取适当的a为间距,将方板划分成N*N个微元，其中N=L/a. 并对每个微元进行标定，用Tkm 表示微元，k表示微元的横坐标，m表示纵坐标。*注:实际求解中我们发现取a=10，N=24时所得解较优-我们在后面将会证明.2 波上单位距离的空洞概率矩阵PL1,PL2的求解(由程序一完成)(1)求出波长矩阵L1,L2: (2)根据时间矩阵T1,T2,T3 求出空洞距离矩阵X1,X2: xij.=( tij - tij3 )/(1/v2-1/v1 )X1(a=10):-8.0040 1.8075 40.2435 39.6145 114.4660 138.1780 160.8310 5.1915 -8.6880 128.4555 123.8355 138.1825 152.6620 97.9300 78.2595 118.3035 -8.4720 119.0955 118.0035 94.7305 33.0340 82.4185 128.6955 115.0275 -3.4320 33.9915 -4.9725 42.8545 89.5900 129.5065 49.8555 38.5995 0.2040 33.2355 33.5475 98.0020 78.3220 51.5665 78.6915 49.3995 3.8040 6.7035 112.7710 83.2420 92.3260 36.8065 -4.2525 -5.6445 7.5120X2(a=10):-6.7800 -8.7405 -2.3445 93.0385 103.1620 116.2540 163.5310 -3.3045 -4.8000 72.2595 124.6635 92.1385 136.7500 140.2300 92.8035 84.9675 5.0640 116.9355 121.0275 129.9025 75.9820 98.0425 86.7915 122.4795 6.2520 86.5515 45.2115 3.0745 77.8900 53.1865 84.0915 86.8035 2.5440 37.0515 45.0675 59.9140 104.0260 178.6825 76.9635 43.6755 0.2760 -4.9965 117.1990 138.0340 104.4940 -5.2415 -6.0885 2.4915 -9.0120*注: 以上两矩阵中出现负值是由于对测量值没有进行修正,但因为我们采用概率比较的方法对其进行统计,因而不会对最后结果造成影响.对于后面出现的概率值都只是对其大小进行比较,所以其值是否在01之间并无影响.(3)根据L,X求出波上单位距离为空洞的概率矩阵PL1,PL2: Plij= xij./ lij 考虑只有一条波的情况,由于此条波经过的空洞这条波上的可能位置是随机的，因此可假设波lij上单位距离内存在空洞的可能性是相同的，也即概率相同。此分布概率不仅与波上存在的空洞距离xij有关，而且与该波的总长lij有关。且该波的总长越小，所经的空洞距离越大，则单位距离内存在空洞的概率越大。所以我们设每个波出现空洞的概率Plij= xij./ lij3 求经过微元Tkm的波的集合Akm。(由程序二完成)当划分足够细时，认为可以用微元的中心Oij(Xij,Yij)表示微元Tkm.考虑用Oij(Xij,Yij)到波的距离d来决定波lij是否经过微元Tkm ，方法如下:求出波lij的直线方程为Aij*X+Bij*Y+C=0,取r=a/2,用d2=(Aij* Xij +Bij* Yij +C) 2/(Aij2+Bij2)进行判断：(D为微元中心到波线的距离)当d2 r2 时,认为波lij过微元Tkm的有向线段很短或为0，即对微元Tkm影响很小，可以忽略不计。.当d2r2 时,认为波lij过微元Tkm。即lij Akm。4 求微元Tkm为空洞的可能性Ptkm。(由程序二完成) 考察经过Tkm的所有波，则波lij在Tkm内的部分为一有向线段，表示该微元在某一方向上存在空洞的概率，这一概率应与有向线段所在波的Plij 及有向线段的长度均有关 。且线段越长，在Plij相等的情况下对Ptkm: 作用越大，为表现出这一影响，我们引入“权”的概念，定义Plij 的权wij=lij / L 表示波lij 对Ptkm:的影响，同时定义“加权概率”Pij=Plij * wij 。我们假设由已知的波的信息必能求出空洞的位置，在考虑有向线段长度的影响后，用该格内各方向上出现空洞的“加权概率”之和来表征Ptkm。由于经过每个微元的波数nkm不同（如图一所示，在波源处和接收器处以及中间部分波分布较密处的微元，经过它们的波数显著大于其他微元,即线条分布的不均匀），导致在加权概率求和时，对nkm较少的微元有失公平,这是测量方法引起的系统误差。为此我们引入“平均加权概率”表示没经过Tkm的波对Tkm状态的影响，定义=Plij /98.这样可减少nkm不同的影响。Ptkm=Pij+（98-nkm）*,其中加权概率Pij =Plij * wij = Plij *lij / L，=Plij /98 nkm为经过微元Tkm的波数加权概率值如下:P1: -0.0333 0.0075 0.1677 0.1651 0.4769 0.5757 0.6701 0.0216 -0.0362 0.5352 0.5160 0.5758 0.6361 0.4080 0.3261 0.4929 -0.0353 0.4962 0.4917 0.3947 0.1376 0.3434 0.5362 0.4793 -0.0143 0.1416 -0.0207 0.1786 0.3733 0.5396 0.2077 0.1608 0.0009 0.1385 0.1398 0.4083 0.3263 0.2149 0.3279 0.2058 0.0159 0.02790.4699 0.3468 0.3847 0.1534 -0.0177 -0.0235 0.0313P2:-0.0282 -0.0364 -0.0098 0.3877 0.4298 0.4844 0.6814 -0.0138 -0.0200 0.3011 0.5194 0.3839 0.5698 0.5843 0.3867 0.3540 0.0211 0.4872 0.5043 0.5413 0.3166 0.4085 0.3616 0.5103 0.0261 0.3606 0.1884 0.0128 0.3245 0.2216 0.3504 0.3617 0.0106 0.1544 0.1878 0.2496 0.4334 0.7445 0.3207 0.1820 0.0011 -0.0208 0.4883 0.5751 0.4354 -0.0218 -0.0254 0.0104 -0.03755 划分数N的求解讨论:从直观上讲,当N较大时,所要求的数据即所得到的结果都会变精确.但当N太大时会使波经过的微元比例过少,即有很多微元没有波通过而成为盲点.因此,我们引进盲点得置信概率函数(与划分的个数和波的密度有关): PP(N,密度).当PP(N,密度) 5%时,N就由波的密度来决定,即:波的分布越密N就可以分的越大.在波的分布一定时N就有最优值.同时为了使相邻的两个波源间有整数个微元,N必须是6的倍数.在第四步我们求出了N=12和N=24的Ptkm值,下面的第六步中求其解的状况.6 用逼近法求微元Tkm的状态，即是否为空洞。 对波lij进行考虑，可知波lij被它所经过的微元分为nij段，每段为空洞的概率由该段所在微元的Ptkm确定。而每个波实际存在空洞的距离已经唯一确定，因此我们根据每段概率的大小，按由大到小的顺序取出前m段，使得这m段的距离和等于波上实际的空洞长度xij，则认为这m段所在的m个微元为空洞。就可以确定所有n*n个微元的状态，得出平板中空洞的分布。 具体步骤如下:第一步：用最大概率法求得初始解。(由程序三完成)对概率矩阵PT= Ptkm:中的元素按从大到小进行排序，取定值ppt为标准量进行比较，在矩阵t中，当 Ptkm:Pt时，tkm=0; 当 Ptkm:Pt时，tkm=1.由初始解得的图(分成12*12个微元) 由初始解得的图(分成24*24个微元.比左图精确)注: 上图中*号表示其所在的微元整个为空洞.ppt的求法: 我们分别对每条波通过的点的概率进行排序,由于每条波上空洞的距离xij是确定的.即每条波上拥有空洞的微元数是定值Alij.分别求出第Alij个微元的概率,而ppt的值应略大于此.当N=12时,我们取ppt=14,当N=24时,我们取ppt=27得到上图.第二步：用逼近法对初始解进行修正。(由程序四,五,六完成)（1）求出微元集合Bij =Tkm ，Tkm在波lij:上，且按Ptkm:的大小进行排序。 （2） 由于波lij:上的空洞距离xij:唯一确定，即加权概率Pij确定,波lij:被微元分成nij:段，设波lij上满足 tkm=0的微元形成集合Cij, 元素个数为mij ,则应满足mij nij * Pij （3）进行如下判断： 从l00 开始，若mij nij * Pij ,则去掉集合Cij第nij * Pij个元素Tkm，同时令tkm=1,否则不变。 (4) 经过（3）处理后得到新的矩阵t，在重复（2）（3），直到对所有的波lij均满足mij nij * Pij ，得到最后优化解t.，如图所示：优化后的图(分成12*12个微元) 优化后的图(分成24*24个微元比左图精确)注: 上图中*号表示其所在的微元整个为空洞.第二问的求解：1: 当只有PQ方向有接收器和发射源时,我们同样利用第一问的求法得到下图:(由程序七完成)*号表示其所在的微元整个为空洞此图为只有PQ面时空洞的分布(由于波数的减少,使N只能取12)比较此图与第一问中求出的N=12时的最优图(我们有理由认为第一问的解较此处精确),可以发现空洞的分布大致相同,即:利用此模型可以求解只有PQ方向有接收器和发射源时空洞的分布,但很明显,精确度和准确度都降低了.其原因自然是波数的减少而引起波的密度的减少.2:求解减少波源与接收器(由程序八完成)(注:我们理解为减少了波源与接收器后,其余波源与接收器的位置不改变.)由前面的分析可知，当波均匀分布时，通过每个微元的波数nkm相同，则Ptkm仅Pij由决定，而不必引进平均概率的误差。为分析方便，我们用方差D来表示nkm的影响。定义D=(nkm-n0)2 /(N*N) 其中 n0=nkm/(N*N)是对所有的nkm求平均值。可知当nkm的值比较接近时，方差D较小，所得结果的精确度和可信度较高。当减少波源和接收器的数目时,会使方差减少,那么我们的原则是:在逐步减少一个波源或接收器时,选择使方差减少最快的波源或接收器.即以此来判断减少哪一个波源或接收器.在给定一定的准确度(即与第一问中所求得空洞结果相差不大)要求的情况下,减少波源或接收器的数目有个极大值. 即以此来判断减少那几个波源或接收器.我们的目的是:在保证一定的精确度(即:保证平板中有一定的波的密度)的情况下,.使波的分布尽量均匀.结果为:减少波源P3,P5,R3,R5接收器Q4,Q6,S4,S6.以上方法的前提只是波源和接收器的位置，而与数据无关。我们可以利用本题给的数据进行一次检验。计算出空洞位置的结果如下图表示: (由程序九完成) 此图为减少波源与接收器后的示意图 *号表示其所在的微元整个为空洞(12*12格)实线代表波,圆点代表波源与接收器 与第一问中的图例比较是空洞的分布是比较接近的.模型检验对于这类数据处理模型，在“没有真值”进行比较的前提下，所有从“建模理论”出发的评定都只能是一种分析。所以最好检验方法就是类似“用已知长度的标准件检验一个测长仪器的精确度”的方法，给定一个空洞分布已知的板及其测量数据t，然后用我们建立的数据处理模型去求解它，并与真实情况比较，从而检验模型。从这种思想出发，我们就用Monte Carlo方法来建立一个计算机仿真模型来检验前文数据处理模型的优劣。任务即评价其合理性（准确性）和稳定性。用于检验数据处理模型的仿真模型的建立随机产生“标准件”板，即随机产生1个N*N阶01方阵B（0表示空洞，1表示介质）.并给定相应数据（PQ:tij）(RS:tij)注：这时每个待检模型的求解结果都应以这种方式（N*N方阵）表示。 分成n组，每组m次，每组产生一次随时机矩阵，用所建数据处理模型来解“标准件”，得n*m个结果Bij (第i组=1,2,n;第j次=1,2, ,m) 建立统计量表征测量模型的准确性(合理性)和稳定性。a.对B0和Bij的统计量：空洞格总数Nij， 空洞格平均坐Xij =xef/Nij, Yij=yef/Nij其中xefj,yef表示某个中的0元素的坐标.b.每次测量的准确性统计量c.每组测量准确性统计量 d.Monte Cado统计量 综上可得，我们建立了E1,E2, E3 三个统计量来表征测量模型的准确性（合理性）和稳定性。其对数据处理模型的评价原则a. 对某一个数据处理模型的稳定性和合理性的检验原则依测量要求给定一定标准E，EE E表征对空洞大小的判断结论准确E 表征对空洞位置的判断结论准确E 表征判断结论具有稳定性b. 对多个数据处理模型的稳定性和合理性的检验原则这种检验主要目的是通过对这些数据处理模型的稳定性和合理性的检验来选取一种最优模型，故侧重相对性，从而确立如下两个原则：综合因素原则：即以E ，E， E从大到小对各模型排名，选择排名均靠前的模型为最优模型。如：对模型，模型，模型的检验结果为EE E第1名第2名第3名那么模型自然是我们从综合因素出发所要的最优模型。侧重某一因素的原则：即实际测量中我们测量结果所反映的各种信息的关心程度本是有区别的，如有时空洞的大小是我们主要关心的，而有时我们则更关心它的大概位置。那我们就有必要在诸多因素中以某一因素为主要考查对象。如对上表结果，若我主要关心空洞大小则应把模型选取做最优模型。综上，在选定最优模型时我们应把两个原则结合起来进行判断。数据t的确定方法对于随机产生的B方阵，分格方式是确定的，为N*N个。以原题方式分布14对接收器和发射器，则测量波的分布就也定了下来。那么我就可以确定每个波经过了方格总数nl（即“细胞”）及分别是那几个T。即确定了由细胞T波集合L，其中 空洞细胞为mL个.则L波方向上单位长度成为空洞的概率P L= mL/ nl。则L波（其长度为lL）方向上空洞长xL= PL* lL。则由就可给出类似原题的测量数据表。仿真模型对所建数据处理模型的评价结果分为24*24个小方格，板边长240米，每边7对发射器和接收器选定E=48格 E=36格 E=96格模拟分为n=5000组，每组m=5次则=36.3 格； =35.6格E； E=89.7格E上述分析结果说明我们的数据模型是比较准确和稳定的。本模型的优缺点评价 优点：我们在建模过程中充分考虑到各方面的因素，尽量提高了结果的准确度。首先进行了合理的假设，排除了可能存在的随机误差。由于测量时波的分布不均匀，造成了测量方法误差。为了减小系统误差对结果的影响，我们引入了权，加权概率及平均加权概率，剔除了系统误差较大的点，并用方差D来表征测量的精确度。同时我们注意到纵向波和横向波分布的对称性，大大减小了运算量。最终结果用状态矩阵画成了图表，直观形象地反映了空洞的分布。但不能完全消除系统误差，因为我们无法从已知数据准确得出未知波的概率。且在划分很细时，可能有些微元没有直线经过，即出现“盲区”，严重影响我们对微元状态的判断。因此我们引入“置信概率”P，定义为P=1-N0/(N*N) , N0是划分为N*N时出现的盲区数。当N0与N*N相比很小时，可认为它对判断的影响可以忽略不计。为不失一般性，我们取P=95%，即当N0/(N*N)5%时，我们认为划分是合理可信的。同时由于我们进行了一系列理想的假设，而实际过程中必然存在随机误差，也必然降低解的可信度。但由于我们用最大概率法和区间逼近，尽可能如实地反映真实情况。 缺点：我们这种用概率统计方法建立的模型要求足够多的数据量，所以发射器和接收器的个数不能减少太多。而且系统误差的减少有赖于波分布的均匀性。 附录：程序一:pq=240 243.3 252.9 268.3 288.4 312.4 339.4;243.3 240 243.3. 252.9 268.3 288.4 312.4;252.9 243.3 240 243.3 252.9 268.3. 288.4;268.3 252.9 243.3 240 243.3 252.9 268.3;288.4 268.3. 252.9 243.3 240 243.3 252.9;312.4 288.4 268.3 252.9 243.3. 240 243.3;339.4 312.4 288.4 268.3 252.9 243.3 240t1=pq/2880t2=0.0611 0.0895 0.1996 0.2032 0.4181 0.4923 0.5646;0.0989. 0.0592 0.4413 0.4318 0.4770 0.5242 0.3805;0.3052 0.4131. 0.0598 0.4153 0.4156 0.3563 0.1919;0.3221 0.4453 0.4040. 0.0738 0.1789 0.0740 0.2122;0.3490 0.4529 0.2263 0.1917. 0.0839 0.1768 0.1810;0.3807 0.3177 0.2364 0.3064 0.2217. 0.0939 0.1031;0.4311 0.3397 0.3566 0.1954 0.0760 0.0688. 0.1042t3=t2-t1;x1=360*t3;n3=x1./pq;n1=x1/240rs=240 243.3 252.9 268.3 288.4 312.4 339.4;243.3 240 243.3. 252.9 268.3 288.4 312.4;252.9 243.3 240 243.3 252.9 268.3. 288.4;268.3 252.9 243.3 240 243.3 252.9 268.3;288.4 268.3. 252.9 243.3 240 243.3 252.9;312.4 288.4 268.3 252.9 243.3. 240 243.3;339.4 312.4 288.4 268.3 252.9 243.3 240t4=rs/2880t5=0.0645 0.0602 0.0813 0.3516 0.3867 0.4314 0.5721;0.0753. 0.0700 0.2852 0.4341 0.3491 0.4800 0.4980;0.3456 0.3205. 0.0974 0.4093 0.4240 0.4540 0.3112;0.3655 0.3289 0.4247. 0.1007 0.3249 0.2134 0.1017;0.3165 0.2409 0.3214 0.3256. 0.0904 0.1874 0.2130;0.2749 0.3891 0.5895 0.3016 0.2058. 0.0841 0.0706;0.4434 0.4919 0.3904 0.0786 0.0709 0.0914. 0.0583t6=t5-t4;x2=360*t6;n4=x2./rs;n2=x2/240程序二:x1=0 0 0 0 0 0 0 40 40 40 40 40 40 40 80 80 80 80 80 80 80.0. 120 120 120 120 120 120 120 160 160 160 160 160 160 160.0. 200 200 200 200 200 200 200 240 240 240 240 240 240 240;. 0.1 40 80 120 160 200 240 0 40.1 80 120 160 200 240 0.0. 40 80.1 120 160 200 240 0 40 80 120.1 160 200 240 0 40.0. 80 120 160.1 200 240 0 40 80 120 160 200.1 240 0.0. 40 80 120 160 200 240.1;y1=0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0. 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;. 240 240 240 240 240 240 240 240 240 240 240 240 240 240.0. 240 240 240 240 240 240 240 240 240 240 240 240 240.0. 240 240 240 240 240 240 240 240 240 240 240 240 240.0. 240 240 240 240 240 240 240 240 240;k3(7,7)=0;b3(7,7)=0;for i=1:1:49x2=x1(1,i) x1(2,i);y2=y1(1,i) y1(2,i);k=polyfit(x2,y2,1);k3(i)=k(1);b3(i)=k(2);x3=0:.4:240;y3=polyval(k,x3);plot(x2,y2,o,x3,y3)axis(0 240 0 240)set(gca,xtick,0 40 80 120 160 200 240)set(gca,ytick,0 40 80 120 160 200 240)grid on hold onendk1=k3;b1=b3;x1=0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0. 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;. 240 240 240 240 240 240 240 240 240 240 240 240 240 240.0. 240 240 240 240 240 240 240 240 240 240 240 240 240.0. 240 240 240 240 240 240 240 240 240 240 240 240 240.0. 240 240 240 240 240 240 240 240 240;y1=0 0 0 0 0 0 0 40 40 40 40 40 40 40 80 80 80 80 80 80 80.0. 120 120 120 120 120 120 120 160 160 160 160 160 160 160.0. 200 200 200 200 200 200 200 240 240 240 240 240 240 240;. 0 40 80 120 160 200 240 0 40 80 120 160 200 240 0.0. 40 80 120 160 200 240 0 40 80 120 160 200 240 0 40.0. 80 120 160 200 240 0 40 80 120 160 200 240 0.0. 40 80 120 160 200 240;k4(7,7)=0;b4(7,7)=0;for i=1:1:49x2=x1(1,i) x1(2,i);y2=y1(1,i) y1(2,i);k=polyfit(x2,y2,1);k4(i)=k(1);b4(i)=k(2);x3=0:.4:240;y3=polyval(k,x3);plot(x2,y2,o,x3,y3)axis(0 240 0 240)set(gca,xtick,0 40 80 120 160 200 240)set(gca,ytick,0 40 80 120 160 200 240)grid on hold onendk2=k4;b2=b4;hold offw=12;c1(w,w)=0;c2(w,w)=0;f=120/w;pq_c(7,7)=0;pp(2*w,2*w)=0;qq(2*w,2*w)=0;nn1=mean(mean(n1);nn2=mean(mean(n2);for h=1:1:w for j=1:1:w ox=2*f*h-f;oy=2*f*j-f; for u=1:1:7 for v=1:1:7 if u=v d=(40*(u-1)-ox)2 else d=(oy-k1(u,v)*ox-b1(u,v)2/(1+k1(u,v)2) end if d=ppt x=2*f*xx-f;