york同位素等时线拟合方法的matlab程序设计new

《york同位素等时线拟合方法的matlab程序设计new》由会员上传分享，免费在线阅读，更多相关内容在教育资源-天天文库。

1、万方数据第31卷第5期物探化探计算技术2009年9月文章编号：lool一1749(2009)05珈511---03York同位素等时线拟合方法的Matlab程序设计赵玉岩1，郝立波1，陆继龙1，潘贵2，侯立宾3(1．吉林大学地球探测科学与技术学院，长春130026；2．黑龙江省地质调查总院齐齐哈尔分院，齐齐哈尔161005；3．71217部队，莱阳265200)摘要：等时线拟合方法的选择是提高同住素测年精度的关键因素之一。这里利用Matlab语言编写了基于York双误差回归模型的等时线拟合程序，包括isoacc．in、samw．In、errs．m和resi

2、．m等模块，实现了读入数据，输入衰变常数和斜率迭代初始值，选择权重方法，以及计算截距、年龄、MSWD和绘图的功能。选用Brooks1972的数据，分别作三种权重等时线拟合，结果与Faure1977‘93的计算结果一致。程序简练，易使用，运算速度快，运行稳定。关键词：同位素等时线；程序；双误差回归模型；Matlab中图分类号：TP312文献标识码：A0前言等时线拟合方法的选择，是提高同位素测年龄精度的关键因素之一。国内、外许多学者曾对此问题进行了深人的研究，York、Brooks、Beyer等提出了双误差不相关，双误差相关等一系列等时线拟合方法u娟J。由于这

3、些方法多较为繁琐，编制程序困难，研究人员普遍采用K．R1．aadwig编写的Is)plot软件，但由于其Windows版与中文版Office不兼容。所以国内研究人员目前仍多使用DOS版。该版本数据格式特殊，不同的软件间数据很难进行交流，一般的地质人员较难掌握。在路远发¨1开发的Geokit程序中，包括了基于最／b-乘法的同位素等时线年龄图解，但没有年龄计算程序。因此，作者在本文中用Matlab语言，编写了基于York(1969)双误差回归模型方程的等时线拟合程序。1回归方程的选择等时线性拟合程序使用York‘31给出的方程基金项目：全国危机矿山接替资源找矿

4、专项资助(20089941)收稿日期：2009—02—17改回日期：2009一03—06(6)，来估计相关和不相关数据的加权最小方差直线，使用Wendt和Carl‘81提出的方程(21)，计算平均平方权重偏差MSWD，使用Faure‘91列出的三种权重模式作为权重。2等时线拟合程序的设计程序设计采用MathWorks公司的Matlab语’言。Matlab是适合多学科多种工作平台的功能强劲的大型软件，其程序代码简练，易使用，可绘制各种通用格式的图件n¨15】。2．1程序流程程序流程如下页图l所示。2．2等时线拟合程序基于上述回归方程和程序流程，作者编写了is

5、oacc．m、s≤i／nw．nl、errs．in和resi．In等时线双误差回归拟合程序。其中，isoacc．in是主程序，功能是读入数据，输入衰变常数和斜率迭代初始值，选择权重方法，以及计算截距、年龄、MSWD和绘图。samw．in、errs．In和resi．111分别是三种不同赋权重万方数据512物探化探计算技术3l卷图1程序流程图Fig．1Theflowchartoftheprogram法的斜率计算程序。2．2．1isoace．in主程序％读取数据。％格式：Excel文件，第1列Y数据，第2列Y％分析误差，第3列x数据，第4列X分析％误差，第5列XY

6、误差相关系数data=xlsread(’data．xls’)；Y=data(：，1)；erry2data(：，2)；X=data(：，3)；errx=data(：，4)；r=data(：，5)；％输入衰变常数和斜率迭代初始值disc=inputdlg({’请输入衰变常数(年)’，’请输入斜率迭代初始值’}，’Yorkl969’，l，{’1．39·10‘一11’，’2’})；dis=str2num(disc{l})；b=str2num(disc{2})；％选择权重方法ButtonName=questdlg(’请选择权重方案’，’权重7，’等权重’，’误差平方

7、倒数’，’残差平方倒数’，’误差平方倒数’)；switchButtonName，case’等权重’，samw；etlse’误差平方倒数’，errs；case’残差平方倒数’resi；end％计算年龄和MSWD，绘等时线图nl=log(1+b)／dis；oi=(b'2幸xre．‘2+yre．^2—2宰b事r．木Xl-e．奉yre)．‘(一1)；MSWD=(1ength(oi)-2)‘(一1)·sum(oi．木(Y—b奉X—a)．2)；scatter(X，Y)；holdon；plot(X，a+b}X，’k一’)；2．2．2$an'1w．m程序％利用最／b-"乘

8、法获取斜率估计值，计算权重iso=regress(Y，[ones(