范文博，王文龍 ，陳岳龍，李大鵬

1) 中國地質(zhì)大學(北京)地球科學與資源學院，北京，100083; 2)北京理工大學數(shù)學學院，北京，100081

內(nèi)容提要：最小二乘法(least squares method)是直線擬合常用的方法，在自然科學和社會科學內(nèi)被廣泛應用，尤其在同位素地質(zhì)年代學領域更是必不可少。僅考慮x或y誤差的普通最小二乘法(normal/ordinary least squares)廣為人知，但事實上最小二乘法并非如此簡單，尤其是在同時考慮x、y誤差(乃至誤差相關)并采用加權(quán)處理時，其數(shù)學處理方法會變得十分復雜，而此時普通最小二乘法顯得極不合理。本文在前人研究的基礎上，結(jié)合同位素地質(zhì)年齡計算的需要，對直線擬合的最小二乘法進行了系統(tǒng)地總結(jié)研究，詳細介紹了普通最小二乘法、加權(quán)普通最小二乘法(normal/ordinary weighting least squares)、標準加權(quán)最小二乘法(standard weighting model)、標準獨立加權(quán)最小二乘法(standard independent weighting model)、獨立加權(quán)最小二乘法(independent weighting model)及誤差相關最小二乘法(error correlated independent weighting model)的數(shù)學原理及相關變量的計算過程。在此基礎上，進一步闡述了MSWD(mean squared weighted deviation)這個同位素地質(zhì)年代學中經(jīng)常使用的參數(shù)的數(shù)學意義，以及MSWD對計算結(jié)果的評判意義。準確理解這些數(shù)學方法，對于我們合理選擇同位素地質(zhì)年齡相關參數(shù)的計算方法，科學評價直線擬合結(jié)果并探討其地質(zhì)意義至關重要。我們的研究同時有助于拓展和完善數(shù)學領域最小二乘法的基本理論，并可用于其他領域相似的研究之中。

同位素地質(zhì)年代學研究是確定地質(zhì)體絕對年齡的主要手段，包括Rb-Sr、Sm-Nd、Re-Os、U-Th-Pb、K-Ar(Ar-Ar)等測年方法，在當前地質(zhì)學領域內(nèi)得到了廣泛應用。在同位素地質(zhì)年齡的計算中，常常需要根據(jù)一組樣品的同位素測試結(jié)果進行直線擬合，最常見的是等時線擬合，主要有Rb-Sr、Sm-Nd、K-Ar、Re-Os、Pb-Pb等等時線。除此之外，U-Pb不一致年齡的計算，以及利用26Al-26Mg、129I-129Xe等絕滅了的放射性(extinct radioactivity)同位素方法研究隕石的相對年齡時，也要用到直線擬合。事實上，在地球科學內(nèi)的諸多領域，都用到了直線擬合。

直線擬合常用的是最小二乘法，其在自然科學和社會科學領域內(nèi)得到了廣泛的應用。其所要解決的問題是：根據(jù)兩個理論上存在線性關系的變量的幾組測試數(shù)據(jù)點(多數(shù)情況下，這些數(shù)據(jù)點不見得恰好位于理論直線之上)，尋找這兩個變量之間的線性函數(shù)關系。當這些測試數(shù)據(jù)點與擬合直線的“偏差平方和”最小時，得到最佳擬合直線，即為所求。通常情況下，直線擬合多數(shù)采用最簡單的普通最小二乘法，這也是諸多數(shù)學書籍中介紹最多、并被多數(shù)人所熟知的直線擬合方法(例如中國科學院數(shù)學研究所數(shù)理統(tǒng)計組，1974；張啟銳，1988；周復恭和黃運成，1989；張小蒂，1991；周紀薌，1993；吳翊等，1995；李慶陽，2001；同濟大學應用數(shù)學系，2002；盛聚等，2008；周惟公，2008；王黎明等，2008；褚寶增和王翠香，2010；何曉群和劉文卿，2011)。這種方法隱含地假設變量x不存在誤差，并將所有誤差都歸于變量y(類似，可將所有誤差都歸于變量x)，即便在x和y同為測量值的情況下，此時普通最小二乘法在確定未知參數(shù)及估計其真實誤差方面，失去了準確性(Borcherds and Sheth, 1995; Sheth et al., 1996)。

僅考慮x或y誤差的普通最小二乘法廣為人知，但實際上最小二乘法并非如此簡單，尤其是在同時考慮x、y誤差(乃至誤差相關)并采用加權(quán)處理時，其數(shù)學處理方法會變得十分復雜，而這對于許多應用同位素地質(zhì)年代學方法進行地質(zhì)學研究的人來說，都十分陌生。之所以會如此，主要有兩方面的原因：① 數(shù)學處理方法的復雜性以及相關資料的相對缺乏；② 采用已有的軟件(Isoplot)而“無需”思考其數(shù)學原理。

雖然同位素地質(zhì)年代學方法被廣泛應用，但地質(zhì)學研究者多熱衷于年齡結(jié)果本身的地質(zhì)意義，同時由于同時跨越了數(shù)學和地質(zhì)學兩個學科，因此很少有人同時結(jié)合地質(zhì)需要、對最小二乘法的數(shù)學原理作系統(tǒng)深入的總結(jié)介紹。如前所述，許多數(shù)學書籍也只涉及最簡單的普通最小二乘法，何曉群和劉文卿(2011)、費業(yè)泰(2010)對考慮加權(quán)的最小二乘法也只是做了簡單介紹。李志昌等(2004)介紹了York(1966,1969)、McIntyre 等(1966)算法的一部分，而更多的同位素書籍(韓吟文和馬振東，2003；陳駿和王鶴年，2004；陳岳龍等，2005；陳道公等，2009)，則忽略具體算法而只強調(diào)MSWD(見下文介紹)對等時線擬合結(jié)果評價的用處。事實上，想要真正了解MSWD對擬合結(jié)果的評判意義，就必須對最小二乘法本身有所把握。國外一些學者在前人[如Deming，1943；據(jù)York(1966)，Borcherds and Sheth(1995)所述]研究的基礎上，對最小二乘法的數(shù)學算法及相關問題作了進一步的研究和探索(如McIntyre et al., 1966; York, 1966, 1967, 1969; Brooks et al., 1972; Titterington and Halliday, 1979; Borcherds and Sheth, 1995; Sheth et al., 1996; Mahon, 1996; Bruzzone and Moreno, 1998)，取得了突破性進展。但這些研究各有新的嘗試和側(cè)重，同時難免存在著一些分歧、甚至錯誤，因此想要系統(tǒng)了解最小二乘法并非易事。Ludwig(2008)開發(fā)的Excel加載宏Isoplot，提供了同位素年齡計算的一個很好的平臺，其以York(1969)、Titterington和Halliday(1979)、Brooks 等(1972)等的算法(Ludwig，2008)為基礎，但未對具體數(shù)學方法作詳細介紹。這樣的話，雖然Isoplot被廣泛使用，但相關同位素年齡的計算方法，對多數(shù)地質(zhì)學研究者來說仍然是一個“謎”。正如Allegre(2008)所述，如今相關測試方法已經(jīng)有了很大的進步，而理論模型——或許是統(tǒng)計模型應該跟上分析方法前進的步伐，這點十分重要，但事實并非如此，這使得我們不得不思考如何對待通過自動化的分析方法已經(jīng)積累的大量數(shù)據(jù)。可見，將最小二乘法這個在同位素地質(zhì)年代學內(nèi)普遍應用的統(tǒng)計方法“顯性化”，是非常有意義的，這也有助于解決實際研究中存在的一些困擾，例如如何理解MSWD對年齡結(jié)果的評價及參考意義。

鑒于此，以下將在前人研究的基礎上，結(jié)合同位素地質(zhì)年齡計算的相關知識，對直線擬合的最小二乘法作詳細總結(jié)和研究，并進一步在此基礎上討論同位素地質(zhì)年代學應用的相關問題。

1 相關測年原理簡介

首先以Rb-Sr等時線法為例，介紹應用最小二乘法所要解決的地質(zhì)問題。核素87Rb通過一次β-衰變變成穩(wěn)定同位素87Sr，根據(jù)放射性衰變定律，進一步推演得到

87Sr=87Sr0+87Rb(eλt-1)(Faure，1986)

式中，λ為87Rb的衰變常數(shù)，t為從同位素體系封閉到現(xiàn)在的時間，亦即所求的地質(zhì)年齡。而87Sr、87Sr0和87Rb均表示物質(zhì)的量，即是以mol為單位，本當寫為n(87Sr)、n(87Sr)0和n(87Rb), 但本文中太多，為簡便計寫為現(xiàn)式，后邊其他同位素亦仿此處理。

兩邊同時除以穩(wěn)定核素86Sr的含量，得：

此式為Rb-Sr同位素定年的基礎，其中87Rb/86Sr、87Sr/86Sr為現(xiàn)今樣品的實測同位素含量比值，(87Sr/86Sr)0為距今t時刻時這兩種同位素的初始比值。由于 (87Sr/86Sr)0未知，因此不能直接計算年齡。對于一套同源、同時、封閉的樣品，其87Sr/86Sr 初始值均一，也就是不同樣品間的(87Sr/86Sr)0一致，而其距今t時刻時的87Rb含量可能不同，從而導致放射性成因累積的87Sr含量不同。這樣的話，這一套樣品內(nèi)的不同樣品之間，87Rb/86Sr、87Sr/86Sr就會有所差異。這時，t、λ為常數(shù)，上式便可看作直線方程

y=a+xb(Faure，1986)

這里b=eλt-1，y=87Rb/86Sr，x=87Sr/86Sr。

對于上述一組樣品，根據(jù)其(87Rb/86Sr，87Sr/86Sr)測試值及其誤差，運用最小二乘法擬合直線，就能夠得到其斜率與截距。根據(jù)斜率b=eλt-1，就可得到所關注的年齡t，而其截距則為(87Sr/86Sr)0，可用于同位素地球化學示蹤。

除Rb-Sr同位素體系之外，Sm-Nd、K-Ar、Re-Os同位素體系定年也常用到等時線法，其同位素定年方程分別為(Faure，1986) ：

對于U-Pb同位素體系，最小二乘法直線擬合主要用于Pb-Pb等時線擬合以及U-Pb諧和圖中的不一致線擬合。前者的斜(率)截(距)式直線方程為(Rosholt et al., 1973; Faure, 1986; 陳道公，2009)：

后者的點斜式直線方程為:

(Wetherill, 1956; Allegre, 2008)

在應用絕滅了的放射性同位素對隕石進行相對年齡測定時，同樣需要首先用到直線擬合，如26Al—26Mg、129I—129Xe隕石測年，其對應直線方程分別為:

(Allegre, 2008)

綜上可知，這里所要解決的問題是，根據(jù)樣品的一組同位素比值測試數(shù)據(jù)(往往帶有實驗誤差)(xi,yi)，i=1,2,3，…，n，選擇最優(yōu)的最小二乘法處理模型，擬合最佳直線并計算斜率、截距及其誤差，并進一步根據(jù)擬合結(jié)果，盡可能地甄別、評判其是否有真正的地質(zhì)年代意義。

2 直線的最小二乘法擬合

現(xiàn)給出樣品的一組同位素比值測試數(shù)據(jù)點(xi、yi)，i=1,2,3，…，n(n>2)，據(jù)此求自變量x和因變量y之間的函數(shù)關系L:y=a+xb。其中a、b為待定參數(shù)，當擬合直線與數(shù)據(jù)點之間的(加權(quán))偏差平方和S最小時，得到最佳擬合直線。關于偏差平方和，根據(jù)York(1969)，有兩種等價的最一般形式:

(1)

(2)

這里(xi,yi)是第i個測試值，(Xi,Yi)是對應的調(diào)整值(圖1)，即最佳擬合直線上對應的點。

(3)

ρi為第i個點x和y變量的誤差相關系數(shù)(error correlation，詳見誤差相關最小二乘法部分)，對于U-Th-Pb測年，由于縱橫坐標變量xi、yi是同時同步測試的，往往需要考慮誤差相關系數(shù)，而對于其它測年方法，在實際應用中多不考慮這個參數(shù)。

圖1 最小二乘擬合中，相關量之間的一般幾何關系Fig. 1 Illustration of relationship between various variables in least squares fitting of a straight lineL為最佳擬合直線，L′為調(diào)整線，(xi,yi)為測試點，(Xi,Yi)為對應調(diào)整點，yri表示y方向的殘差，xri表示x方向的殘差L is the best straight line，L′ is the adjustment line，(xi,yi)is the data point while(Xi,Yi)is the adjustment point corresponding. yri denotes residual in y direction，and xri in x direction

當ρi=0時，x、y的誤差不相關，這時得到Y(jié)ork(1966)、Borcherds和Sheth(1995)中所提到的Deming(1943)所定義的偏差平方和的公式，即與(1)式對應的(4)式，以及McIntyre 等(1966)、Wendt(1969)[轉(zhuǎn)引自Brooks, 1972]、Borcherds和Sheth(1995)、Sheth等(1996)、何曉群和劉文卿(2011)所用的與(2)對應的(5)式，

(4)

(5)

(6)

圖1給出了擬合直線與相關量之間的一般幾何關系。L為最佳擬合直線y=a+xb，連結(jié)測試點(xi,yi)與(Xi,Yi)的直線Li′為調(diào)整線，其斜率為每個點的y方向殘差(yi-residual，以下用yri表示)與x方向殘差(xi-residual，以下用xri表示)之比，即

(7)

可以看出，隨著擬合直線(即參數(shù)a、b)的變化，S也會相應地發(fā)生變化，當S取最小值時，所有點沿各自調(diào)整線回歸到擬合直線上的距離之和最短，這時就得到了最佳擬合直線。方便起見，這里定義最佳擬合直線所對應的S的最小值為S*。以下從最簡單的普通最小二乘法開始，逐步深入，分別介紹基于不同假設所得到的最小二乘直線擬合方法。實際應用中，則需要根據(jù)數(shù)據(jù)本身特征，選擇最適合的最小二乘法進行直線擬合。

2.1 普通最小二乘法(normal/ordinary least squares)

顧名思義，這是介紹并使用最多的最小二乘擬合方法(例如吳翊等，1995；同濟大學應用數(shù)學系，2002；李慶陽，2001；盛聚等，2008；王黎明等，2008；褚寶增和王翠香，2010；何曉群和劉文卿，2011)，在Isoplot被廣泛使用之前，等時線擬合多用這種方法[如申浩澈(1983)所述]。假設xi不存在誤差(xi=Xi)，yi存在誤差(yi≠Yi)，令wi=w(yi)=1(可認為w(xi)=)，這時所有數(shù)據(jù)點沿與y軸平行的調(diào)整線調(diào)整到最佳擬合直線(圖2a)，

(8)

對于固定的一組數(shù)據(jù)，(8)式是S關于a、b的二元函數(shù)，根據(jù)多元函數(shù)的極值條件(同濟大學應用數(shù)學系，2002)，分別對S求關于a、b的偏導數(shù)，并令其為0，得：

進一步變換得到所謂的“法方程”或“正規(guī)方程”，

(9)

(10)

據(jù)(9)、(10)，解得(吳翊等，1995；王黎明等，2008；盛聚等，2008；何曉群和劉文卿，2011)：

(11)

(12)

與之類似，也可將所有的誤差都歸于xi(xi≠Xi，yi=Yi)，令w(xi)=1，wi=w(xi)/b2=1/b2(可認為w(yi)=)，進行處理，如圖2b。

對于一些研究對象，普通最小二乘法假設一個變量不存在誤差，這是成立的。而對于同位素地質(zhì)年代學，由于所有xi、yi都是測試值，這點顯然不成立。筆者認為，對于Rb-Sr等時線法，考慮到87Rb/86Sr的測試誤差往往比87Sr/86Sr的誤差大許多，因此采用將所有誤差歸于87Rb/86Sr的(加權(quán))普通最小二乘法，或許是可以接受的。Sm-Nd等時線法類似。

2.2 加權(quán)普通最小二乘法(normal/ordinary weighting least squares)

這是對普通最小二乘法的改進，與前者不同的是其給予每個點不同的權(quán)重。假設xi不存在誤差(xi=Xi)，yi存在誤差(yi≠Yi)，這時wi=w(yi)=1/σ2(yi)(可認為w(xi)=)，所有調(diào)整線都平行于y軸(圖2a)，

(13)

同理，可以根據(jù)?S/?a=0，?S/?b=0，解得(Bruzzone and Moreno, 1998)：

(14)

對(14)式中b的表達式進一步變換可得：

(15)

類似，也可將所有的誤差都歸于xi(xi≠Xi，yi=Yi)并采用加權(quán)處理，令w(xi)=1/σ2(xi)，wi=w(xi)/b2(可認為w(yi)=)，進而確定最佳擬合直線(如圖2b)。

2.3 標準加權(quán)最小二乘法(standard weighting model)

普通最小二乘法以及加權(quán)普通最小二乘法的一個顯著缺陷是只考慮一個變量的誤差，在同位素年代學研究中，兩個變量都是測試值，從而不可避免都存在誤差，這顯然不夠精確。

(16)

同上，令?S/?a=0，得

(17)

令?S/?b=0，得：

(18)

將(17)式代入(18)，化簡會得到關于b的一元二次方程，解之并取其中一個解，得到:

(19)

a的解法同上，由最佳擬合直線過重心計算得到[(17)式，與(12)式本質(zhì)相同]。

2.4 標準獨立加權(quán)最小二乘法(standard independent weighting model)

圖 2 不同類型的最小二乘擬合方法中，數(shù)據(jù)點的調(diào)整方式：(a)(加權(quán)的)普通最小二乘模型，只考慮縱坐標的誤差，所有調(diào)整線都平行于y軸，需要說明的是，對于同一組數(shù)據(jù)點，加權(quán)的普通最小二乘擬合與普通最小二乘擬合所擬合直線可能會不同；(b)(加權(quán)的)普通最小二乘模型，只考慮橫坐標的誤差；(c) 標準加權(quán)模型，所有調(diào)整線都互相平行且垂直于最佳擬合直線；(d)標準獨立加權(quán)模型，所有調(diào)整線都互相平行；(e)獨立加權(quán)模型，一般情況下，不同數(shù)據(jù)點的調(diào)整線不平行，但限制在三角陰影區(qū)內(nèi)；(f)誤差相關的加權(quán)最小二乘擬合模型，同樣不同數(shù)據(jù)點的調(diào)整線不再平行，但調(diào)整線不再限制在三角陰影區(qū)內(nèi)(參考York, 1966, 1967, 1969; Borcherds and Sheth, 1995; Sheth et al., 1996; Mahon, 1996; Bruzzone and Moreno, 1998，有增改)Fig. 2 Various types of adjustments in different models of least squares fitting: (a) normal/ordinary least squares and normal/ordinary weighting least squares, only the errors in y direction are considered. It should be noted that weighting least squares may have different results from the not-weighting in this situation; (b) normal/ordinary least squares and normal/ordinary weighting least squares, only the errors in x direction are considered; (c) standard weighting model, all the adjusting lines are parallel to each other and are perpendicular to the best fitting line at the same time; (d) standard independent weighting model, the adjusting lines are parallel to each other; (e) independent weighting model, the adjusting lines are not necessary to be parallel but should be situated within the shadow area; (f) error correlated independent weighting model, the adjusting lines are not necessary to be parallel or be restricted within the shadow area(after York, 1966, 1967, 1969; Borcherds and Sheth, 1995; Sheth et al., 1996; Mahon, 1996; Bruzzone and Moreno, 1998)

(20)

令?S/?a=0，?S/?b=0，同上解得：

(21)

(22)

S的最小化給出關于b的二次方程，解之并選取其中一個合適的根作為擬合直線的斜率，得到(Bruzzone and Moreno，1998):

(23)

這里b*指(14)中的b，即對于此組數(shù)據(jù)，假若xi沒有誤差時計算所得的斜率。S、Sx、Sy、Sxy、Sxx、Syy如2.2節(jié)中所述。

2.5 獨立加權(quán)最小二乘法(independent weighting model)

這是雙變量最小二乘擬合的更一般情形，其同時考慮每個點xi、yi的不確定度，并認為σ(xi)、σ(yi)二者之間不存在確定的數(shù)學關系。與上述2.1～2.4節(jié)的情形不同，此時調(diào)整線一般不再平行，并被限制在假設x、y都不存在誤差時的豎直和水平調(diào)整線以及最佳直線所包圍的三角形內(nèi)(圖2e)，而且得不到斜率b和截距a的解析解，必須使用迭代法，尋找認為達到我們所需要精度的近似解。下面給出兩種處理辦法。

2.5.1 獨立加權(quán)模型解法一

York(1966)從表達式(4)著手解決這個問題，延續(xù)其思路，以下進一步更詳細地解析相關過程：

2.5.1.1 步驟 1

S是關于Xi、Yi的泛函，根據(jù)泛函極值條件及相關運算法則，當S取得極值時，其一階變分為0，即有：

(24)

(25)

這里δs、δ(xi)、δ(yi)分別表示S、xi、yi的一階變分[此概念來自“變分法”這一數(shù)學分支學科，而變分法則是研究求泛函極大值和極小值的方法，簡單理解，“一階變分”類似于高等數(shù)學中的微分，有關其詳細的數(shù)學定義可參考老大中(2004)、邵克勇等(2011)]。

由于(Xi,Yi)位于最佳擬合直線上，因此有：

Yi=a+bXi

(26)

兩邊同時求變分，得到：

δa+Xiδb+bδ(Xi)-δ(Yi)=0

(27)

對于不同的i，(27)式乘以各自的一個未定乘數(shù)λi，并求和，得到：

(28)

(25)+(28)，得：

(29)

令δ(Xi)、δ(Yi)、δa、δb的系數(shù)為0，得到：

(30)

(31)

(32)

(33)

將(30)、(31)分別代入(26)中，分別替換Xi、Yi，得到：

(34)

由(34)式解出λi，得到：

λi=wi(yi-a-bxi)

(35)

wi由(6)式給出。

2.5.1.2 步驟2

現(xiàn)在，(32)、(33)、(35)一共有(n+2)個等式，(n+2)個未知數(shù)：a、b和λi。根據(jù)(32)和(35)，有：

(36)

根據(jù)(33)和(35)，有：

(37)

將(30)代入(37)，替換Xi，得到：

(38)

進一步化簡，有：

(39)

將(35)式代入(39)式，替換λi得到：

(40)

這里可以通過(36)、(40)同時解出a、b。

2.5.1.3 步驟 3

現(xiàn)在，從等式(36)中解出a，得到[與(12)式形式相同]：

(41)

這里

(42)

將(41)式代入(40)式，化簡即得到Y(jié)ork(1966)的“最小二乘立方方程”。事實上，正如Mahon(1996)所認識到的，由于wi中也含有b，因此這也只是一個形式上的立方方程：

(43)

(43)式中，左邊第一項系數(shù)化為1，得：

b3-3αb2+3βb-γ=0

(44)

其中

(45)

(46)

(47)

將(44)式看作一元三次方程，得到其三個實根為：

這里

York(1966)指出b3是我們所需的根，因此：

(48)

由于α、β、γ都含有wi，而wi中也含有b，因此并不能直接從(48)式中解出所需的斜率b。這時就需要用到方程求根的迭代法(參見李慶揚，2001)，用計算機編程計算(參見彭國倫，2002)。簡要說明其原理：

① 首先給出一個b的近似值b(1)，代入(6)計算wi并將所得代入(45)、(46)、(47)計算α、β、γ，進一步將所得代入(48)得到b(2)；

② 將b(2)代入(6)計算wi并將所得代入(45)、(46)、(47)計算α、β、γ，進一步將所得代入(48)得到b(3)；

③ 依次重復以上過程，直到b(k+1)-b(k)<ε。這里ε是一個我們可以接受的、使得b(k+1)與bk接近幾乎完全相等的數(shù)值，取決于我們所要求的斜率b的精度。b(k+1)即為所求斜率b。

這時，a同樣由(12)式(也就是(41)式)解出。

將(30)、(31)代入(7)式，得到一個極有意義的表達式(York，1967)，

(49)

據(jù)此可以確定調(diào)整線的方向，可以看出，調(diào)整線的斜率與最佳擬合直線的斜率b符號相反，亦即調(diào)整線位于圖2 所示的三角陰影區(qū)內(nèi)。此等式對于2.1節(jié)至2.5節(jié)都成立。

2.5.2 獨立加權(quán)模型解法二

其它學者(McIntyre et al, 1966; Wendt, 1969，轉(zhuǎn)引自Brooks, 1972; Borcherds and Sheth, 1995等)則從公式(5)著手，根據(jù)S取最小值時的條件，得出斜率和截距的計算辦法。下面我們嘗試推導，

令?S/?a=0，得：

(50)

令?S/?b=0，得：

(51)

將(50)代入(51)，化簡得到：

(52)

同樣可將(52)式看作關于b的三次方程，用上述類似的迭代法處理即可。

Borcherds和Sheth(1995)并未先將(50)式代入(51)，而是將解三次方程的復雜性轉(zhuǎn)移到迭代過程中去。根據(jù)(51)化簡得到一個二次方程：

Ab2+Bb+C=0

(53)

其中

A=a∑w2(yi)w(xi)xi-∑w2(yi)w(xi)xiyi

B=a2∑w2(yi)w(xi)-2a∑w2(yi)w(xi)yi

C=∑w2(xi)w(yi)xiyi-a∑w2(xi)w(yi)xi

這時開始用迭代法：

① 用普通最小二乘法得到b的初始值b0，代入(50)計算a，用這個a代入(51)式，計算得到b的兩個根。② 將b的兩個根及①中a代入(5)分別求S。選擇較小的S對應的b作為新的b值，記為b1。③b1再代入(51)式，計算得到b的兩個根，與②類似處理，得到b2。依次重復，得到b(j)，直到b(k+1)-b(k)<ε。ε與2.5.1節(jié)中相同。這時的b(k+1)與對應的a即為所求。

2.5.3 獨立加權(quán)模型與其它最小二乘模型的關系

可以認為2.1～2.4節(jié)所述的四種模型是2.5節(jié)中的獨立加權(quán)最小二乘模型的特殊情形，雖然如此，但由于前四種情形下，能夠得到斜率b的解析解，而且實踐中會有與之類似的問題，因此詳細討論是有必要的。而獨立加權(quán)最小二乘模型則無解析解，需用迭代法進行數(shù)值計算。

2.1～2.4節(jié)以及2.5.2節(jié)都是用偏導數(shù)為0作為取得極值的條件進行推導的，因此其解是相通的；而2.5.1節(jié)則采用變分法得到，給定2.5.1節(jié)中(43)式不同的條件，進一步演算，則可得到公式(11)、(15)、(20)、(22)。

2.6 誤差相關最小二乘法(error correlated independent weighting model)

York(1969)將獨立加權(quán)模型進一步復雜化，提出了(1)(2)兩種計算公式，并引進誤差相關(系數(shù))ρi誤差相關這個參數(shù)，即σ(xi)和σ(yi)之間的線性相關系數(shù)，ρi1，當其取0時，σ(xi)和σ(yi)相互獨立；當其絕對值為1時，說明其完全線性相關，+1意味著正相關，-1則負相關(Ludwig，2008；梁晉文，2001；費業(yè)泰，2010；褚寶增和王翠香，2010)。通常所見的U-Pb諧和圖中，誤差橢圓的方位也就是誤差相關系數(shù)的反映。對于此種最小二乘法，不僅考慮每個點的誤差，而且考慮每個點縱橫坐標誤差之間的相關系數(shù)，這時所有測試點的調(diào)整線一般不平行，且不會限制在上述2.5節(jié)中所述的陰影三角形內(nèi)(圖2f)。

與2.5節(jié)類似，也可以從(1)式著手，令S的一階變分為0，或從(2)式著手，分別對a、b求偏導數(shù)并令其結(jié)果為0，推演b的計算過程。這里不再做詳細推導，其具體原理同上。York(1969)從(2)式著手，得到的最小二乘平方方程為：

(54)

據(jù)此的最佳擬合直線斜率為：

(55)

(56)

其中：

wi由(3)式給出。

x、y的殘差分別為：

(57)

(58)

3 相關參數(shù)的誤差估計

關于最佳擬合直線的斜率b和截距a的誤差(標準差)計算，較多采用的是誤差傳遞公式(McIntyre et al., 1966; York, 1966, 1969; Brooks et al., 1972; Kullerud, 1991; Mahon, 1996)，而York(1969)所給的誤差公式則存在錯誤(Titterington and Halliday, 1979; Mahon, 1996)。Titterington and Halliday(1979)用最大似然法所討論的結(jié)果與誤差傳遞公式的計算結(jié)果相同。Sheth等(1996)則提出了另一種標準差的計算方法。無論哪種算法，都是對標準差的近似估計，這里主要介紹前者。而另一個至關重要的問題是在得到標準差后，如何確定置信區(qū)間。

3.1 標準差的計算

這里以誤差相關的最小二乘模型為例介紹，其它模型可以通過相關參數(shù)取特殊值得到，也可用類似的方法進行推導。記(56)式為F(b,xi,yi)=0，即：

(61)

根據(jù)誤差傳遞公式(省略泰勒序列的高次項)，有(Titterington and Halliday, 1979；Mahon，1996；費業(yè)泰，2010)：

(62)

根據(jù)隱函數(shù)求導公式(同濟大學應用數(shù)學系，2002)，

(63)

以上兩式聯(lián)立，就可得到Mahon(1996)所用的公式：

(64)

或

(65)

對于截距a的誤差，也可采用類似的方法，將(12)式看作關于xi、yi的函數(shù)，有：

(66)

將?F/?xi、?F/?yi、?F/?b、?a/?xi、?a/?yi代入(64)(66)式，得到以下兩式(Mahon, 1996; Titterington and Halliday, 1979)：

(67)

(68)

York(1969)中認為可以將點關于最佳擬合直線的分散程度考慮在內(nèi)，進行誤差估計，得到擴展的不確定度(York, 1966，1969; Kullerud, 1991)，即：

(69)

有關MSWD的介紹見下。Mahon(1996)認為這種做法是從數(shù)學角度來說是不嚴格的，如果MSWD小于1的話，就會降低誤差。因此在實踐中如若用到，一定要指明。

3.2 確定置信區(qū)間

通常給標準差(σb和σa；代表68%的置信區(qū)間)乘以2(更準確地說是1.96，表1)，得到斜率b和截距a95%的置信區(qū)間。事實上，這樣做的前提是用于進行等時線擬合的樣品數(shù)據(jù)點足夠多(即n較大)，且斜率和截距的計算結(jié)果服從正態(tài)分布。這時，平均值附近的2σ(1.96σ)區(qū)間，代表了概率密度曲線中95%的面積區(qū)域。當擬合數(shù)據(jù)點較少時，就需要乘以t-分布的α(這里取0.025，對應于95%的置信區(qū)間)分位點tn-2(α=0.025)(參見褚寶增和王翠香等，2010)(Mahon, 1996; Ludwig, 2008)。 這時：

b=b*±tn-2(α=0.025)·σb，a=a*±tn-2(α=0.025)·σa(b*表示最佳斜率，a*同)。

當數(shù)據(jù)點較少時，斜率的誤差會非常大，正因為如此，Ludwig(1998, 2000, 2003, 2008等)一直強調(diào)避免3點或4點等時線。

4 MSWD

根據(jù)Wendt和Carl(1991)，MSWD(mean squared weighted deviation)計算方法為

(70)

其中f=n-k，為自由度，對于加權(quán)且固定截距的線性回歸，k取1，而對于無限制條件的線性回歸，k取2(Mahon, 1996)，n表示擬合所用的數(shù)據(jù)點數(shù)目。wi如(3)式或(6)式所示。據(jù)此可知，這里：

(71)

即最佳擬合直線所對應的偏差平方和(S*)的平均值。根據(jù)以上3中的討論可知，S*的大小反映了所有數(shù)據(jù)點偏離最佳等時線的總距離，當所有數(shù)據(jù)點嚴格位于最佳等時線時，S*取值為0，當然MSWD也為0，當然這種情況是絕對理想化的。

接下來一個重要的問題是在同位素地質(zhì)年代學中，如何應用MSWD對擬合結(jié)果進行解釋。通常可能會這樣認為，當測試所得的MSWD小于對應自由度的HMSWD時，說明擬合效果較好，具有地質(zhì)意義；反之高于時，則需要重新考慮是否樣品完全符合測年條件。是否如陳駿和王鶴年(2004)所述“當MSWD>2.5時，很可能是一條誤差等時線”，或如韓吟文和馬振東(2003)所述，“該值越小，等時線質(zhì)量越好。當存在地球化學誤差時，MSWD>1；當不存在地球化學誤差時，MSWD≤1”。問題似乎并沒有這么簡單，Kalsbeek (1992)、Wendt (1992)對此已經(jīng)做了一些討論。

首先，需要明確的是，MSWD值明顯受對實驗分析誤差估計準確性的影響(陳道公，2009；這點可從公式(70)與(1)至(6)明顯地看出)，只有在使用最小二乘法(或最大似然法)擬合等時線，且對數(shù)據(jù)點的誤差估計準確時，MSWD才可作為檢驗擬合結(jié)果的參數(shù)(Mahon, 1996)。但實際研究工作中，這點并沒有得到重視，筆者查閱了一些近年來用Rb-Sr等時線法測年的文章發(fā)現(xiàn)，許多文章并未交代等時線擬合中相關變量的(分析)誤差估計，卻使用MSWD對擬合結(jié)果進行評價，如李真真等(2009)、李光來等(2011)、張臣(1998)、藍廷廣等(2011)、祝禧艷等(2010)、周雄等(2010)、陳江峰等(2003)、魏俊浩等(2003)、毛光周等(2008)。而給出了誤差的文獻中，對分析誤差估計的方式依然存在著不明確性。李秋立等(2006)明確指出擬合時，87Rb/86Sr比值采用2%誤差、Sr同位素比值采用保守外精度誤差0.05%；姚海濤等(2001)選用87Rb/86Sr 、87Sr/86Sr 的分析誤差分別為1.5%、0.02%，并提及這是“一般的全巖Rb-Sr等時線法定年的分析誤差”；楊進輝和周新華(2000)則只提及87Rb/86Sr的輸入誤差為0.5%～2%。

其次，根據(jù)上述3、4部分的詳細介紹可知，這里僅考慮了實驗分析誤差，且給出了僅由分析過程中的隨機誤差引起數(shù)據(jù)點對最佳擬合直線偏離時，MSWD的概率分布情況。但事實上，數(shù)據(jù)點對最佳擬合直線的偏離(甚至靠近)，不僅可能來自實驗分析誤差，而且還有可能來自地質(zhì)背景、每個樣品本身、以及測年方法的選擇，但對于除實驗分析誤差之外的其它誤差，目前似乎并沒有合理的方法將其定量化(Allegre, 2008)，因此在擬合時，也僅考慮了實驗測試誤差。實際研究中所用樣品，更多地是接近而非完全滿足同位素測年所需的條件，這就必然降低MSWD對擬合結(jié)果的評判意義。以等時線法為例，測年樣品不一定真正具有完全相同的初始同位素比值(見Zheng, 1989)，同位素體系能夠在所研究對象內(nèi)不一定能夠完全封閉。這樣的話，一方面，95%的置信邊界只是表明當數(shù)據(jù)點的分散僅是由測試過程中的隨機誤差所引起時，MSWD會有95%的可能性落在這個區(qū)間內(nèi)，但并不排除存在其它因素導致MSWD落入這個區(qū)域的可能性，并且這種可能性有多大是未知的；另一方面，對于MSWD落入其95%置信區(qū)間外的情況，也有可能是可以接受的。

基于此，在應用同位素地質(zhì)年代學方法對地質(zhì)體進行年代測定時，需盡量選擇符合所用方法基本假設的樣品，并且盡可能地準確估計分析測試誤差。在此基礎上，當MSWD落入95%的置信區(qū)間內(nèi)時，說明測試結(jié)果線性相關程度良好，在盡可能地降低其他因素影響的情況下，所獲年齡具有較大的參考意義。而當MSWD落入95%的置信區(qū)間外時，且測試誤差估計準確、MSWD不至于過大時，擬合結(jié)果仍有一定參考意義。當然MSWD十分大時，則需要慎重對待測年結(jié)果。除此之外，更要結(jié)合地質(zhì)背景，對所得結(jié)果進行檢驗。

5 討論

在同位素地質(zhì)年代學研究中，最小二乘法直線擬合是一種常用的數(shù)學處理方法。雖然普通最小二乘法廣為人知，但由于擬合所用數(shù)據(jù)點均為測試所得同位素比值，縱橫坐標變量均存在誤差甚至誤差相關，而普通最小二乘法假設一個變量不存在誤差的前提條件已經(jīng)不能成立。因此，需要建立更加符合數(shù)據(jù)本身特點的擬合方法。

基于此，我們在前人研究的基礎上，對最小二乘法進行了系統(tǒng)性地總結(jié)與研究。詳細介紹了普通最小二乘法、加權(quán)普通最小二乘法、標準加權(quán)最小二乘法、標準獨立加權(quán)最小二乘法、獨立加權(quán)最小二乘法以及誤差相關最小二乘法，給出了各自的“最小二乘”原理，以及擬合直線斜率、截距及其誤差的計算方法。在將相關數(shù)學計算方法進一步完善并系統(tǒng)化的同時，也闡明了各種最小二乘模型之間的相互關系。誤差相關最小二乘法可以看作是最一般的情形，而獨立加權(quán)、標準獨立加權(quán)、標準加權(quán)、普通最小二乘法則可看作是其一步步“特殊化”的情形。誤差相關、獨立加權(quán)這兩種最小二乘法沒有數(shù)值解，需要借助計算機編程后用迭代法得到其近似解，文中給出了迭代的具體算法；而前四種最小二乘法則具有數(shù)值解。

可以看出，不同的擬合模型對應于不同的假設條件，這允許根據(jù)數(shù)據(jù)本身的特點，選擇合適的模型進行擬合。如對于U-Pb不一致曲線，只能選擇誤差相關最小二乘法，而對于Rb-Sr、Sm-Nd、K-Ar、Re-Os等時線法等，則可選用最一般形式的獨立加權(quán)最小二乘法。當然，如果數(shù)據(jù)有符合標準獨立加權(quán)、標準加權(quán)、(普通)最小二乘法的假設條件的，即可采用這些簡化的模型，從而簡化計算。這里，筆者發(fā)現(xiàn)對于Rb-Sr、Sm-Nd等時線法，由于87Rb/86Sr(147Sm/144Nd)的實驗誤差比87Sr/86Sr(143Nd/144Nd)要大許多，因此采用普通最小二乘法或加權(quán)普通最小二乘法，將數(shù)據(jù)點偏離擬合直線的原因全部歸于橫坐標x(87Rb/86Sr，147Sm/144Nd)，或許是可以接受的。

我們進一步澄清了參數(shù)MSWD對于同位素地質(zhì)年代測定結(jié)果的參考意義。理論上，在擬合直線時，應該考慮包括每個樣品的采樣、地質(zhì)背景、方法應用及實驗分析所帶來的所有誤差，但實際上只有實驗分析誤差是可定量估計的。因此，擬合過程中我們只考慮了實驗分析誤差。鑒于此，在所選樣品盡量符合所用測年方法原理的基礎上，如果實驗分析誤差估計準確，那么當MSWD落入其95%的概率區(qū)間內(nèi)或不至于太大時，擬合結(jié)果具有較大的參考意義，但也不排除錯誤的測年結(jié)果存在的可能性。當MSWD過于大時，則需慎重。當然，由于地質(zhì)“未知”世界的復雜性，在評價擬合結(jié)果時，也需要結(jié)合已有的地質(zhì)知識進行判別。實際研究中，對于每一個數(shù)據(jù)點，更需將其看作“一個個同位素比值”，而非“一個個年齡”。

同位素年代學研究方法是建立在一定的假設基礎之上的，根據(jù)地質(zhì)樣品本身的規(guī)律，對這些假設予以統(tǒng)計學修正，應該是今后發(fā)展的一個方向。這也有利于以地質(zhì)學內(nèi)容為素材，開拓新的方法。我們的研究對于其他自然科學及社會科學領域內(nèi)，相似問題的解決也具有一定的參考意義。

6 結(jié)論

本文對最小二乘直線擬合這一同位素地質(zhì)年代學中常用的統(tǒng)計方法進行了總結(jié)，詳解了其數(shù)學原理與實現(xiàn)過程，實現(xiàn)了從普通最小二乘模型到誤差相關最小二乘模型的統(tǒng)一。不同類型的最小二乘模型對應于不同的假設條件，研究者可據(jù)其數(shù)據(jù)特點進行選擇。對于U-Pb不一致曲線，只能選擇誤差相關最小二乘法，而對于Rb-Sr、Sm-Nd、K-Ar、Re-Os及其它等時線，除可選用最一般形式的獨立加權(quán)最小二乘法，也可據(jù)數(shù)據(jù)特點，選擇標準獨立加權(quán)、標準加權(quán)、普通最小二乘法等簡化模型計算。在此基礎上，我們進一步澄清了MSWD這一參數(shù)對擬合結(jié)果的評判意義。準確理解這些，對于我們更好地應用相關測年方法，評價測年結(jié)果的可靠性至關重要。我們的研究也可用于其他領域內(nèi)相似問題的研究中。

致謝：本文完成過程中得到了北京離子探針中心宋彪研究員，中國科學院地質(zhì)與地球物理研究所姜能研究員、儲著銀副研究員，中國地質(zhì)大學(北京)諸慧燕老師、蘇文博教授、高世臣教授、王翠香副教授，中國軍事醫(yī)學科學院張春茂、上海財經(jīng)大學馬馨瑞及中國地質(zhì)大學(北京)甘彩虹、薛玉山、楊寬、張志強等的幫助，匿名評審人及章雨旭編輯的認真工作使本文進一步完善，在此一并表示衷心感謝。