高教社杯全国大学生数学建模竞赛题目CT系统参数标定及成像之欧阳地创编

1、欧阳地创编2018西交数模第一次模拟赛时间：2021.03. 04创作：欧阳地欧阳地创编欧阳地创编欧阳地创编欧阳地创编数学建模论文首页选题 A队伍编号696班级学号姓名队长钱学森73班2173714071孔令辉队员1钱学森73班2171311098杨峥队员2电气642160400092李佳航018 年230 日摘要本题针对一种二维CT获取样品内部结构信息的工作方式及成像原理，意在通 过借助已知结构样品进行参数标定，消除系统误差，而后对未知结构的样品进行 成像，得出该未知介质的相关信息。问题一中，要求通过标立模板的相关，确立CT系统的旋转中心、探测器单元 之间的距离以及该CT系统使用的X射线的1

2、80个方向。本文利用标立模板几何参 数，通过条数、以及探测器单元间距相等等信息，首先计算出探测器单元之间的 距离为0,2778mm.根据X射线与标左模板的几何特性，以椭圆短轴所在宜线为x 轴，长轴所在直线为y轴建立直角坐标系，得旋转中心在所建立的坐标系中的坐 标为 I 。最后通过近似确立X射线180个方向的旋转角度具有髙度线性相关性，根据部分确定数据拟合出整体旋转角度。前十五组旋转角度为：29.6493 31.0027 31.5572 32.6476 33.6795 34.6478 35.6483 36.6490 37.647538.6486 39.6488 40.6491 41.6487 4

3、2.6497 43.6492问题二中，要求通过已求得的标泄参数，确定未知介质在正方形托盘中的位 置、几何形状和吸收率等信息，并具体给岀所要求十个位宜的吸收率。本文依据 CT层析成像原理，利用逆拉东变换作出重构图像，并利用Excel中数据分布计算出 原位置介质的相关性质。所求十个位置的吸收率为：序号12345吸收率0.01850.9701-0.00011.17031.0287序号678910吸收率1.44251.2815-0.0074-0.00250.0194问题三中，与问题二计算方法相勺以得十个位置的吸收率为：序号12345吸收率0.07362.53057.13270.03320.8015序号

4、678910吸收率3.14875.57140.05308.13620.0272问题四中，要求对问题一q】的标怎模型为上行改进以减、误差并增加稳左性，本文利用在问题一求解过程中遇到的问题进行思考，首先适半增大模板减小偶然误 差，其次做出投射图像后应容易找到极值，并且图像应有一立的对称性；图像扫 描后应尽量少地得出重复数据。关键词：CT层析成像Radon变换与逆变换吸收率MATLAB算法1.问题重述1.1问题背景应用ng yongPAGE2羽构造 ient in parabolic problemsJ.PAGE2PAGE2PAGE2PAGE2PAGE2PAGE2PAGE2PAGE2PAGE2PAG

5、E2PAGE2PAGE2PAGE2PAGE2PAGE2PAGE2PAGE2PAGE2PAGE2PAGE2PAGE2PAGE2PAGE2PAGE2PAGE2PAGE2PAGE2PAGE2PAGE2PAGE2PAGE2PAGE2PAGE2PAGE2PAGE2PAGE2PAGE2PAGE2PAGE2PAGE2PAGE2PAGE2PAGE2PAGE2PAGE2PAGE2PAGE2PAGE2PAGE2PAGE2PAGE2PAGE2PAGE2PAGE2PAGE2PAGE2PAGE2PAGE2PAGE3PAGE3PAGE3PAGE3PAGE3PAGE3PAGE3CTComputedTomography）可以

6、在不破坏样品的情况下，利用样品对射线能量的吸收特性对生物组织和工程材料的样品进行断层成像，由此获取样品内部的结构信息。1.2相关信息本题介绍了一种二维CT系统，平行入射的X射线垂直于探测器平面，每个探 测器单元看成一个接收点，且等距排列。X射线的发射器和探测器相对位置固左 不变，整个发射-接收系统绕某固龙的旋转中心逆时针旋转180次。对每一个X射 线方向，在具有512个等距单元的探测器上测疑经位宜固左不动的二维待检测介 质吸收衰减后的射线能量，并经过增益等处理后得到180组接收信息。并且，为 消除安装误差，需要对安装好的CT系统进行参数标定，即借助于已知结构的样品 （称为模板）标肚CT系统的参

7、数，并据此对未知结构的样品进行成像。题目附件 中提供了标定模板的几何信息，接受信息，待测介质的接收信息以及图3所给位 宜的相应数据。1.3需要解决的问题13.1问题一在正方形托盘上放苣两个均匀固体介质组成的标左模板，模板的几何信息如 图2所示，请根拯这一模板及英接收信息，确泄CT系统旋转中心在正方形托盘中 的位置、探测器单元之间的距离以及该CT系统使用的X射线的180个方向。1.3.2问题二附件3是利用上述CT系统得到的某未知介质的接收信息。利用第一问中得到 的标泄参数，确左该未知介质在正方形托盘中的位置、几何形状和吸收率等信 息。另外，请具体给出图3所给的10个位置处的吸收率。13.3问题三

8、利用给出上述CT系统得到的另一个未知介质的接收信息。利用第一问中得到 的标定参数，给岀该未知介质的相关信息。另外，请具体给出图3所给的10个位 置处的吸收率。1.3.4问题四分析第一问中参数标泄的精度和稳立性。在此基础上自行设汁新模板、建立 对应的标定模型，以改进标立精度和稳定性，并说明理由。图1CT系统示意图图2模板示意图（单位：mm）图310个位置示意图2 问题分析一 4).5 可MTVYV/£卜、IOO21问题一分析结合图2和附件1表中数据，可以首先计算出CT系统探测器个数和模板长度 度量的比值，运用程序2.1可以得岀模板的几何形状如图2.1.1,可大致认为它是 对称的。对于附

9、件2,由于对180个方向尚无淸晰地认识，首先用同样的方式做 出数据分布图2.1.2,观察到图像比较平滑，因此认为按表格的顺序180个方向是 相邻较密、不错位的。（图中有色区域表示该点有吸收率，蓝色部分表示吸收率大于100)*«|V'zii（r 必型 «3qx» wftiw >»“ maw 工j«t）muidi saw m«xh» >4 a b 快W 折 0 ED J ' / a o 、e 心 X 电 Dta 口25020D15010050050100 W 2002303003QC450400350

10、3002tX）150,00K 100 120 140 160 180图2.1.1附件1的数据分布图图2.1.2附件2的数据分布图在图2.1.2中易观查到有一红色条形图案，这是在不同的方向扫描到圆形时留 下的，可以此为突破口首先求出探测器单元之间的距离。求出探测器单元之间的 距离与增益比率之后，可以根拯几何关系，自己设立坐标系并通过数学运算讣算 出旋转中心。最后在建立的坐标系内，将待求解的180个旋转方向转换成X射线 的斜率进行数学运算。2.2问题二分析题目给出了未知介质的接收信息，要求岀介质的相关信息，可以搜索相应的 数学模型，将附件中给岀的按照射线条数与旋转角度列成的表格，一一对应为相 应坐

11、标点的吸收情况，从而根拯各坐标点的不同性质，还原回该未知介质的几何 信息与吸收率等信息。2.3问题三分析问题三与问题二类似，可以大致看出数据分布更具有一般性，不容易描述岀 未知介质的相关信息，可以通过图形大致描绘出介质相关信息。2.4问题四分析问题四要求分析题目所给的二维CT系统，设汁新的标泄模板以提高原系统的 精确度与稳定性。可以搜索相关资料，根据第一问的求解思路与求解过程，以规 避求解过程中因标定模板自身性质而岀现的误差为原则进行思路拓展，以设讣高 精度与稳定性的标定模板。3.模型假设假设在射线经过介质时能量只损失在介质中，及不考虑衍射等现象；假设附件中所给岀数据是正确的、可以直接利用的；

12、假设旋转中心在相邻两条射线的中间直线上：4模型的建立与求解4.1问题一模型的建立与求解4.1.1建立模型从题目中可以得出，由于x射线之间得间距相等，不管x射线怎么旋转，穿 过托盘上圆的x射线条数应该是大致相等的。可将穿过圆形标圮模板的X射线与 模板建立模型示意图如下：图4.i.l.lX射线穿过圆形标泄模板示意图如图4.2.1.2所示可以设穿过圆的在圆心两旁的射线被圆所截的距离分别为 Xi，X2, x3和y】，y2> y3。求解旋转中心时，以椭圆形标泄模板短轴所在直线为X轴，长轴所在直线为 Y轴建立直角坐标系，将托盘进行划分。由附件二的数据分布图的不对称性可 知，旋转中心应在对称轴某一侧，

13、大致确怎旋转中心方位后，根据旋转中心两侧 探测器个数不变且同一探测器接收的与距旋转中心的距离不变具体确左旋转中心 的位置。求解X射线旋转角度时，设穿过椭圆的最边缘的射线到椭圆中心的距离为 R，取椭圆中心为原点，模板的对称中心为x轴建立平面直角坐标系，X射线所在 直线的斜率为ko4.1.2符号说明d相邻探测器之间的距离N对应于圆形模板的射线条数或探测器个数R标定模板的吸收率X处理数据时的增益率4>圆形模板的直径4.1.3模型求解求探测器单元间距离：通过MATLAB编程求岀计算出的条形带平均涉及探测器个数为匚三I,沿直径方向上的平均吸收率巨可:由图二圆形模板直径,计算出探测器之间的平均距离：

14、探测系统平均增益率:求旋转中心： 分析图2.1.2可知，由于蓝色区域仅分布在后部分角度范围内，因此估il旋转 中心的位置在对称轴的某一侧，红色区域分为两部分时表示该角度下由两部分射 线分别照射经过椭圆形和圆形：为了讨论方便，现对正方形托盘做出如下划分：欧阳地创编欧阳地创编简单分析可知，如果旋转点在I区域，则沿180个方向照射后不会出现射线 分成两部分的情况；在II或III或IV区域，当吸收率出现最大值的时候射线也被分 成两部分，而不是像图1.2那样成为一部分，也排除：综合各因素可判断旋转中 心应该在V区域。因为发射-接收系统逆时针旋转，且图2.1.2中两部分红色区 域，间距缩小，说明旋转之后在

15、垂直于发射-接收方向上二者的距离是缩短的，从 而确左旋转中心在短轴的上半侧，先运用Excel对附件二第1列数据作图（即画出 第一个方向上的扫描图像）如图4.1.3.2,发现此时得到两部分图像，说明两模板 之间有一部分射线直接被探测器接收，直到第14列数据两部分图像结合在了一 起，如图4.1.3.3。>3037303 60 03S00D43002003aoa图4.1.3.3第14列数拯分布笫1列第14列图4.1.3.2第一列数拯分布结合托盘的几何特征，在垂直于对称轴方向上应该会岀现最大的吸收率，利 用MATLAB求得出现最大吸收率的方向为第151个方向，在此方向结合增益率得 出的模板长度为

16、匸进一步验证了结果：平行于对称轴方向上最大的吸收率出现在沿对称轴的直线上，计算得出为第61个方向，此方向的模板长度 匸三二1,也验证了结果。根据极近似水平方向为第61方向、最大吸收率出现在235号探测器，极近似竖直方向为151方向、最大吸收率出现在 223号探测器，设旋转中心到竖直轴、水平轴的距离分别为x,y。为求解还需要另一个方向的等疑关系，选取椭圆和圆的一条外公切线的方向,经计算得岀经椭圆与圆公切线所在直线的X射线为第372号射线。则在所建立 坐标系内，切线方程为：进而列出二元方程组:解之得：因此旋转中心在所建立的坐标系中的坐标为匸三Z1°求180个旋转方向：由于x射线的发生装置

17、是连续旋转的，所以在512条射线中穿过模板的射线 条数应该是连续变化的，所以用MATLAB编写程序，计算每次旋穿过模板的x射 线条数，并画岀图像如图4.1.3.4,300穿过模板舱It线条塗与旋芸方向的关丢280260 240 200180 -160 140 -220 -020406080100120140160180淀转次数120图4丄3.4穿过模板射线条数随旋转次数变化曲线图中横坐标为旋转次数，纵坐标为穿过模板的射线条数。从图中可以看出附 表2中的数据是按照射线发射装置旋转的顺序依次给岀的，而且可以看岀，180 次旋转后装宜共旋转了 180度，每次旋转的角度近乎相等。考虑到在前50次旋转中

18、，穿过整个装垃的x射线条数与穿过椭圆的条数相 等。设穿过椭圆的最边缘的射线到椭圆中心的距离为R，取椭圆中心为原点，模板的对称中心为x轴建立平面直角坐标系，由于射线可以近似看成与椭圆和圆X，+ y2 = 2都相切，可以得到以下方程组：x2 + y2 =以y = /<% -f mX2152402欧阳地创编欧阳地创编欧阳地创编欧阳地创编化简以上方程组可得:欧阳地创编欧阳地创编欧阳地创编欧阳地创编利用MATLAB编程得到图像如图4丄35所示90807060504030 20 111111!''05101520253035404550旋法次数图4.1.3.5前50次旋转角度变化曲线

19、用同样方式得出最后而25组数据的图像为：欧阳地创编欧阳地创编3025 20 - 嚴矣角贡与旋無次数陀关系100 11111155160165170175180旋转次数图4.1.3.6后25次旋转角度变化曲线从图像中可以看出每次旋转角度近似为1度。为了更加精确地讣算旋转角度以便据此得到估算苴他角度地依据，本文设讣 了另外一种算法。从附表二中可以看出，前25组数据中穿过椭圆的射线与穿过圆 的射线没有交叉，所以数据中换算岀来的最大吸收距离就是近似穿过椭圆中心的 射线被椭圆所截的距离，弦长公式为：2 2x y= 11524(/,y = kx m且这个距离由直线的斜率k唯一确定,已知k、m时可解得：12

20、00jk2 + 1-.1 = J152/c2 -m2 + 402152/c2 + 402所以用MATLAB编程计算得岀了较为精确的前15次旋转角度的结果如下：29.6493 31.0027 31.5572 32.6476 33.6795 34.6478 35.6483 36.6490 37.6475 38.6486 39.6488 40.6491 41.6487 42.6497 43.6492可以看出，这些结果的线性相关性非常好，据此利用Excel进行拟合，得到如下图表：欧阳地创编欧阳地创编旋转角度与旋转次数的关系O4欧阳地创编欧阳地创编欧阳地创编欧阳地创编246810121416 1S旋转次

21、数图4.13.7旋转角度与旋转次数拟合公式从图表中可以得出拟合公式为：y = °9938 + 28.718,相关系数为：0.9996根据拟合公式利用MATLAB编写程序计算得出所有方向角度。4.2问题二模型的建立与求解4.2.1建立模型通过分析问题与附件数据可以发现，题目所给数据与介质相关性质的二维分 布具有对应关系，X射线将介质一条线上的性质投影为一点。拯此，本文利用 radon变换与radon逆变换进行运算，以通过投影后的讯号重建原始未知介质相关 性质的二维分布。该变换的上义为：令密度函数f(X)=f(x,y)是一个的疋义域为S 的紧致台(compact support)t令R为

22、radon变换的运算子，则Rf(x,y)是一个定义在 回空间中的直线LoRf(Q 二 J f(X)|dX|其基本思想为：radon变换可以理解为图像在空间的投影，空间的每一点对应一条 直线，而radon变换是图像像素点在每一条直线上的积分。因此，图像中高灰度 值的直线会在空间形成亮点，而低灰度值的线段在空间形成暗点。对直线的检测 转化为在变换区域对亮点、暗点的检测。Radon变换是一幅图像在一个特泄的角度下的径向线方向的投影，一幅图像 的radon变换是每一个像素radon变换的集合。对于MATLAB中语句R = radon(l, theta),如果theta是一个标量，R则是一个包含在the

23、ta的列向量。如果theta是 一个向量，R则是一个矩阵，矩阵的每一列是对应其中一个theta的radon变换。 而radon变换的逆运算，就可以将CT系统对于每一直线上的X射线吸收率数据还 原回未知介质的物理性质。radon反变换的公式是：f(x, y) = / g(x cos 0 + y sin®, 8)ddJo该反变换操作比较简单，思路淸晰，可以借助数学运算软件计算。欧阳地创编欧阳地创编4.2.2模型求解 在MATLAB上将附件2和附件3的数据导入，利用radon逆变换将附件中数 据重建未知介质信息，得到图像如图422.1。5010015020025050 1001502002

24、50100 200 300100 200 300100200300100 200 J00100 200 300图4.2.2.1 radon逆变换重建介质图像（1）为使逆拉东变换后得到的图像与原图像大小相等，在逆radon变换公式中取362条射线，每两条射线相距0.2770所以得到的图形长宽也为100,由第一问求 解得知，X射线从约29度位置开始旋转，为消除重建模型与真实模型的旋转角度 差异，将从29度逆radon变换得到的数据导入Excel表格可得变换后图像坐标水平方向平移了 31 X 0. 2770 = 8. 587 ,竖直方向平移了22 X 0. 2770 - 6.094,所以将原图中的坐

25、标按上述数值平移就得到了变换后的坐标。 从模板的逆radon变换产生的矩阵中可以发现，模板中所有点对应的灰度都 近似为0.5,又由于模板的吸收率为1,所以相对比例近似为2。据此利用MATLAB编程可以算出图中对应十个点的吸收率如下表所示： 表422.1所求十点的吸收率（1）序号12345吸收率0.01850.9701-0.00011.17031.0287序号678910吸收率1.44251.2815-0.0074-0.00250.0194然后根据逆radon变换的结果，将得的数据导入Excel表格，在不同范围内的数值 填充成不同颜色对不同吸收率的部分进行色块填充，得到结果如图422.2所示。

26、再利用Excel表格中寻找各椭圆左点的坐标，由于拉东变换中x射线之间的距离都 是0.2770,所以表格中两组数据在实际物体上的距离也是0.2770。由100 - 0.2770 362得知表格取362组数据时，总长度与原图基本相同，此时 可以长度为基础，计算坐标变换公式。根据前文坐标变换的逆变换，用MATLAB编写程序运算各色块（由运算可得 徐色块均为椭圆形）的坐标及长短轴相应数据。将问题二中数据进行radon逆变 换后的图像最低点在Excel行数和列数，将行数和列数乘以倍率0.2770即距离图 像边界的距离，由于第一个图距离两边界的距离已知，可以得到平移的方向和距 离,具体结果如表4222所示

27、。图4.2.2.2不同区域吸收率关系图（"以每个椭圆中心的吸收率代表整个椭圆的吸收率则A, B, C, D, E, F的吸收率分别为：0、1.1870, 1.2914, 0、0.9877、1.0632o表422.2不同吸收率介质在托盘中的位置、几何形状与吸收率椭圆编号ABCDEFx063.017546.259054.707540.580551.245048.3365yo33.101575.067071.189030.331553.184055.5385halfl11.19205.25577.08529.645422.68871.8478half26.11913.465312.8872

28、6.384539.53471.9439吸收率01.18701.291400.98771.0632其中X0、y0、halfl、half2分别表示椭圆中心横纵坐标和两个半轴。数据均以左下 角的点为坐标原点建立坐标系求得。4.3问题三模型的建立与求解4.3.1建立模型与第二问类似，仍利用radon变换的思想建立模型，将空间每一条射线所投影的 点还原回一条直线，将数据合并后重建未知介质的几何性质与吸收率。4.3.2模型求解将附件5的数据导入MATLAB,利用radon逆变换将附件中数据重建未知介质信 息，得到图像如图4.3.2.1。10020030010020030010020030040050010

29、0 200 300 400 500100200300400500100 200 300 400 500图4.3.2.1 radon逆变换重建介质图像（2）其中图三图四显示完全，图二图四经角度修正后得出。 为具体算岀所要求十个点的吸收率，将数据带入MATLAB中运算后结果如下表所 示：表43.2.1所求十点的吸收率（2）序号12345吸收率0.07362.53057.13270.03320.8015序号678910吸收率3.14875.57140.05308.13620.0272根radon逆变换的结果，将得的数据导入Excel表格，在不同范用内的数值 填充成不同颜色对不同吸收率的部分进行色块填

30、充，得到结果如图4.322所示。图4.322不同区域吸收率关系图（2）其中无色处吸收率为0,黄色吸收率为0到2,红色为2到4,浅蓝色为4到 6,绿色为6到8,紫色为8到9。从表格中大致取岀图案中心，调用程序得到图案中心在(51.2450, 47.9210) (若以椭圆中心为原点,则图案中心在(1.2450, -2.0790)附近。4.4问题四的分析与求解：为了便于求相邻探测器之间的距离，考虑最好仍选择圆形模板，但是要适当 增大圆形模板的半径，从而减小偶然误差：为了利用180次旋转得到的投射图像 求出旋转中心的坐标，鉴于在原来的标左模板中椭圆和圆的内公切线的选择有较 大误差，新的模板中应尽量使得

31、切线容易取得，且做出投射图像后容易找到极 值，并且图像应有一泄的对称性；图像扫描后应尽量少地得出重复数据，因此两 个图像差别应较大：考虑到方便地识别投射位置以确左旋转中心，应至少设置两 个模板且相互隔开。基于以上分析，设计出如下新模板：图4.4.1新设计标泄模板标上方法：首先借助于圆形模板很容易求得相邻探测器间距，利用对称性， 更容易求得射线水平、竖直和一条倾斜方向的位程，因此用和第1问相同的思 路，此模板相对来说更能准确地确定CT系统的参数。5. 结果的分析与检验5.1问题一结果分析与检验第一问求得的单元间距，旋转中心与旋转角度与所搜集资料的实际值相差不大且 符合现实认知，在运算旋转角度时，

32、由于线性关系良好，可以印证旋转中心与单 元间距计算误差不大。5.2问题二结果分析第二问十个位置的吸收率计算，由于是直接由题目所给数据汁算得岀，结果较精 确，所得结果与逆radon变换所得图形相对应。在求解具体坐标与几何关系时， 通过Excel表格数据宜接讣算得岀，可能存在误差，但由具体结果运用radon变换 检验后可以看出，基本符合题目所给数据。6. 模型的优缺点分析第一问中求单元间距与解旋转中心时利用了切线的特殊性质，但切线的选择不一 定准确，因为射线宽度远小于探测器宽度且所用射线不一左恰好为标眾模板切线 位置。从题目中可以得出，由于x射线之间得间距相等，不管x射线怎么旋转， 穿过托盘上圆的

33、x射线条数应该是大致相等的。根据附件二可以发现，前13组穿 过圆的x射线条数都是29条，由此可以粗略的计算出x射线之间的距离为 2R厉=0.2758欧阳地创编欧阳地创编但是这种做法并不精确，因为在29条射线中最两边的x射线并不是与圆相切的, 为尽量减小误差，本文利用图4.1.1.1进行误差检验，则如下方程成立：欧阳地创编欧阳地创编欧阳地创编欧阳地创编设模板的吸收率为P，第n个数据为上述方程可以转化为:欧阳地创编欧阳地创编欧阳地创编欧阳地创编利用MATLAB编程带入多组数据求其平均值可得k = 1.7725,进而求得 d = 0.2770.从这里可以看岀在圆两边的射线不管是否与圆完全相切，对结果

34、的影响不是很大，所以该建模方式可以使用7. 模型的改进简单地说第2, 3问根据逆radon矩阵求出坐标值存在一泄的偶然误差，计算 量大。而且本文的模型没有考虑到噪声等其他因素的影响，因此输出图像模糊有 光晕。为了得到清晰的图像，可以进行频域滤波。首先二维傅里叶变换对为：F(Uj v)°°伽決 网）血创oo J-coF(uy 尹齐血 1 lrydudv欧阳地创编欧阳地创编欧阳地创编欧阳地创编引入傅里叶切片立理，其中3是频率分量:/oog(p咖伽叫 QQP OO f OO 广 8II «f3/)da;cos0 + ysiii9 p)e i2nujpdxdydp 8 丿

35、 00 f QQOO POOcos 9 + ysin 9 p)e j2mdpdp dxdy00f g 广8=/ f(x,y)e伽吩g的艸*)必创J_QQ /-OO=F(Q COS 0,3 sill 0)/ f(x.y)QQ J QQf oo这说明一个投影的一维傅里叶变换，是二维投影矩阵的二维傅里叶变换的一个切 片，执行换元操作后，引入窗函数计算积分计算式并滤波，从而得到一个相对较 好的结果。参考文献：工业CT技术刘丰林工业CT系统旋转中心左位方法研究 刘明进附录%程序2.1作附件1的数据分布图axis equalfor i=l:256for j=l:256if(A(i,j)>0)plot

36、(j,257需严)hold onendendend%程序2.2作附件2的数据分布图for i=l:512forj=l:180if(AS(i,j)>0)if(AS(ij)>100)plot(j,513-i,'b)hold onelseplot(j,513，L)hold onendendendendhold off%程序13计算探测系统的增益率及探测器的平均距离ticd=ones(86,l)；%统讣圆模板对应探测器的平均个数 m=zeros(86,l);%统计圆模板对应探测器的最大吸收率平均值 for i=l:14k=0;x=0;forj=374:430if(AS(j,i)&g

37、t;0)k=k+l;endif(AS(j,i)>x)x=AS(jJ);endendd(i)=k;m(i)=x;endfori=109:180k=O;x=O;forj=45:110 if(AS(j,i)>0) k=k+l;endif(AS(j,i)>x)x=AS(j,i);endendd(i-94)=k;m(i-94)=x;enda=mean(d)b=mea n(m)toe%程序1.4求出竖直射线近似方向d=zeros(180/l);k=0 forj=l:180 if(k<max(AS(:J) d(j)=max(AS(:,j); end end d fori=l:180

38、if(k<d(i) k=d(i)Iendend%程序1.5推出水平射线近似方向d=zeros(180,l)；u=1.7721;forj=14:109d(j)=max(AS(:,j);d(j)Jend%程序1.6根据你和公式用matlab编写程序: th=zeros(l/180);fori=l:l:180th(i)=0.9938*i+28.718;endth计算得出所有的方向为：1至15列29.7118 30.7056 31.6994 32.6932 33.6870 34.6808 35.6746 36.668437.6622 38.6560 39.6498 40.6436 41.6374

39、 42.6312 43.625016至30列44.6188 45.6126 46.6064 47.6002 48.5940 49.5878 50.5816 51.575452.5692 53.5630 54.5568 55.5506 56.5444 57.5382 58.532031至45列59.5258 60.5196 61.5134 62.5072 63.5010 64.4948 65.4886 66.482467.4762 68.4700 69.4638 70.4576 71.4514 72.4452 73.439046至60列74.4328 75.4266 76.4204 77.414

40、2 78.4080 79.4018 80.3956 81.389482.3832 83.3770 84.3708 85.3646 86.3584 87.3522 88.346061至75列89.3398 90.3336 913274 92.3212 93.3150 94.3088 95.3026 96.296497.2902 98.2840 99.2778 100.2716 101.2654 102.2592 103.253076至90列104.2468 105.2406 106.2344 107.2282 108.2220109.2158110.2096111.2034 112.1972 1

41、13.1910 114.1848 115.1786 116.1724117.1662118.1600欧阳地创编欧阳地创编91至105列119.1538 120.1476 121.1414 122.1352 123.1290124.1228125.1166126.1104 127.1042 128.0980 129.0918 130.0856 131.0794 106至120列132.0732133.0670134.0608 135.0546 136.0484 137.0422 138.0360139.0298140.0236141.0174 142.0112 143.0050 143.9988

42、 144.9926 145.9864146.9802147.9740148.9678 149.9616 150.9554 151.9492 152.9430153.9368154.9306155.9244 156.9182 157.9120 158.9058 159.8996 160.8934161.8872162.8810欧阳地创编欧阳地创编136至150列163.8748 164.8686 165.8624 166.8562167.8500168.8438169.8376170.8314 171.8252 172.8190 173.8128 174.8066175.8004176.7942177.7880欧阳地创编欧阳地创编151至165列178.7818 179.7756 180.7694 181.7632182.7570183.7508184.7446185.7384 186.7322 187.7260 188.7198 189.7136 190.7074191.7012192.6950欧阳地创编欧阳地创编166至180列193.6888 194.6826 195.6764 196.6702 197.6640 198.6578 199.6