等溫剩磁曲線分解方法及其批量處理工具BatchUnMix

白帆，常燎,2*，薛鵬飛，汪詩(shī)舜

1 北京大學(xué)地球與空間科學(xué)學(xué)院造山帶與地殼演化教育部重點(diǎn)實(shí)驗(yàn)室，北京 100871 2 青島海洋科學(xué)與技術(shù)試點(diǎn)國(guó)家實(shí)驗(yàn)室海洋地質(zhì)過(guò)程與環(huán)境功能實(shí)驗(yàn)室，青島 266071

0 引言

環(huán)境磁學(xué)通過(guò)研究天然樣品的磁學(xué)性質(zhì)，用于示蹤不同環(huán)境、地質(zhì)和地球物理過(guò)程（例如，潘永信和朱日祥，1996；王喜生等，2006；鄧成龍等，2007；Liu et al.,2012）.環(huán)境磁學(xué)研究大多數(shù)是利用全巖磁學(xué)參數(shù)（bulk magnetic parameters）來(lái)表征樣品中所有磁性礦物的綜合響應(yīng)，然而天然樣品中的磁性礦物組分通常高度不均一，含有不同礦物類(lèi)型、形態(tài)粒度、化學(xué)狀態(tài)、以及微觀結(jié)構(gòu)的磁性礦物集合體，使宏觀磁學(xué)參數(shù)的解釋具有不確定性與復(fù)雜性，因此，分離并提取不同磁性組分信息對(duì)于提升環(huán)境磁學(xué)方法的探測(cè)精度和可信度尤為重要（Heslop,2015）.磁性組分是指具有相同礦物學(xué)或者相似物理性質(zhì)（例如顆粒尺寸、形態(tài)、結(jié)晶度等）的磁性顆粒集合（Egli,2003）.目前廣泛使用的分離磁性組分的巖石磁學(xué)方法包括測(cè)量并分析樣品的磁滯回線（Heslop,2015）、磁化率-溫度曲線（Liu et al.,2012）、剩磁曲線（Egli,2003,2004a,b）以及一階反轉(zhuǎn)曲線（Roberts et al.,2000,2014；秦華峰，2008）等.這些方法的廣泛應(yīng)用極大地幫助了人們對(duì)自然環(huán)境中磁性礦物的形成、搬運(yùn)、沉積和轉(zhuǎn)化等過(guò)程的理解.在這些方法中，由于測(cè)量等溫剩磁（isothermal remanent magnetization；IRM）曲線對(duì)樣品無(wú)損且高效簡(jiǎn)便，所以IRM曲線分解成為目前分離磁性礦物組分最常用的方法之一.IRM曲線根據(jù)其獲得方式一般分為IRM獲得曲線和反向場(chǎng)退磁曲線（Egli,2003）.IRM獲得曲線測(cè)量方式為將退磁后的樣品置于逐漸增大的外磁場(chǎng)中，外磁場(chǎng)每增大一步后將其撤去，然后測(cè)量樣品的剩磁.IRM反向場(chǎng)退磁曲線的測(cè)量方式為先將樣品沿著正向外磁場(chǎng)方向飽和，然后施加逐漸增大的反向磁場(chǎng)，測(cè)量樣品在撤去外磁場(chǎng)后的剩磁.

本文首先從數(shù)學(xué)模型的角度綜述了近年來(lái)提出的IRM曲線分解方法，然后重點(diǎn)介紹我們基于MATLAB開(kāi)發(fā)的GUI工具BatchUnMix，并利用實(shí)驗(yàn)數(shù)據(jù)展示和評(píng)估了BatchUnMix處理天然樣品IRM曲線的實(shí)用性，最后討論了IRM曲線分解方法的局限性并提出了相關(guān)方面有待進(jìn)一步研究的科學(xué)問(wèn)題.

1 參數(shù)化方法

1.1 累積對(duì)數(shù)高斯函數(shù)

Robertson和France（1994）通過(guò)分析磁赤鐵礦、赤鐵礦和針鐵礦樣品的IRM獲得曲線發(fā)現(xiàn)，單一磁性礦物的IRM獲得曲線形態(tài)可以用累積對(duì)數(shù)高斯（cumulative log-Gaussian,CLG）分布進(jìn)行擬合：

（1）

其中，B為磁性礦物獲取剩磁時(shí)施加的外磁場(chǎng)大??；B1/2為平均剩磁矯頑力，為磁性礦物獲得半飽和剩磁（SIRM）時(shí)的外磁場(chǎng)大?。籇P（dispersion parameter）為離散參數(shù)，對(duì)應(yīng)于對(duì)數(shù)-高斯函數(shù)的一個(gè)標(biāo)準(zhǔn)差；Mr為該單一磁性組分的飽和剩磁.Robertson和France（1994）進(jìn)一步的研究表明，磁性礦物顆粒間不存在靜磁相互作用的情況下，N種礦物混合后的IRM獲得曲線是每種礦物所對(duì)應(yīng)的IRM獲得曲線的線性疊加，即可以用每種礦物對(duì)應(yīng)的CLG函數(shù)之和擬合磁性礦物混合后的IRM獲得曲線，對(duì)應(yīng)的方程為

（2）

其中i表示第i個(gè)磁性成分.

Eyre（1996）和McIntosh等（1996）在研究中國(guó)黃土?xí)r采用了上述分析方法，利用對(duì)數(shù)高斯概率密度函數(shù)擬合樣品中每個(gè)磁性組分的IRM獲得曲線的一階導(dǎo)數(shù)曲線，從而得到樣品的矯頑力分布.然而在擬合IRM曲線的過(guò)程中，他們對(duì)磁性組分的數(shù)目及擬合效果的好壞是通過(guò)直接觀察判斷的，存在很大程度的主觀因素.此后，Stockhausen（1998）提出使用三個(gè)基于殘差平方和（RSS）的參數(shù)來(lái)判斷IRM曲線擬合優(yōu)度，并用FORTRAN語(yǔ)言編寫(xiě)了擬合程序，但他并沒(méi)有給出何時(shí)使用這三個(gè)參數(shù)的標(biāo)準(zhǔn).

Kruiver等（2001）基于前人提出的累積對(duì)數(shù)高斯分布及對(duì)數(shù)高斯概率密度分布的擬合方法，又提出了用正態(tài)分布的標(biāo)準(zhǔn)分?jǐn)?shù)（Z分?jǐn)?shù)）表示每個(gè)磁性組分的IRM曲線，稱(chēng)其為標(biāo)準(zhǔn)獲得曲線（standardised acquisition plot,SAP；圖1c）.如果樣品只含有一個(gè)磁性組分，那么其IRM曲線的SAP為一條直線.如果樣品含有多個(gè)磁性組分，那么其IRM曲線的SAP為各個(gè)磁性組分的SAP疊加，會(huì)呈現(xiàn)出曲線的形態(tài).例如，當(dāng)樣品中存在赤鐵礦等高矯頑力磁性組分時(shí)，其IRM曲線的SAP高場(chǎng)部分會(huì)有明顯的彎曲（例如，Abrajevitch et al.,2009）.結(jié)合IRM獲得曲線（也被稱(chēng)為線性獲得曲線，LAP；圖1a）、標(biāo)準(zhǔn)獲得曲線（SAP；圖1c）及IRM一階導(dǎo)數(shù)曲線（也被稱(chēng)為梯度曲線，gradient of acquisition plot，GAP；圖1b），可以利用不同數(shù)量的CLG曲線，通過(guò)調(diào)整它們各自的log（B1/2）、DP和Mr參數(shù)對(duì)IRM曲線進(jìn)行組分分析（Kruiver et al.,2001）.對(duì)于擬合過(guò)程中存在的多解性問(wèn)題，例如可以用不同數(shù)目的磁性組分?jǐn)M合剩磁曲線，他們通過(guò)對(duì)不同擬合結(jié)果與原始數(shù)據(jù)的殘差平方和（RSS）進(jìn)行F檢驗(yàn)和t檢驗(yàn)來(lái)判斷兩個(gè)模型是否存在顯著差異，從而得到最優(yōu)擬合參數(shù).基于以上方法開(kāi)發(fā)的EXCEL工具IRM-CLG已在古地磁學(xué)及環(huán)境磁學(xué)研究得到了廣泛使用（Yamazaki and Ikehara,2012；Hu et al.,2013；Chang et al.,2014,2018；Abrajevitch et al.,2015；Wang et al.,2019；Xue et al.,2019）.

Stockhausen（1998）和Kruiver等（2001）提出的方法都需要不斷手動(dòng)調(diào)節(jié)每個(gè)磁性組分IRM曲線的擬合參數(shù)（log（B1/2）、DP和Mr）.當(dāng)樣品的組分較為復(fù)雜時(shí)，例如當(dāng)樣品中包含四個(gè)磁性組分，需要調(diào)節(jié)的擬合參數(shù)多達(dá)12個(gè)，此時(shí)擬合一個(gè)樣品的IRM曲線需要花費(fèi)較長(zhǎng)時(shí)間.Heslop等（2002）提出了基于對(duì)數(shù)似然和期望最大化算法（Dempster et al.,1977）自動(dòng)尋找擬合IRM一階導(dǎo)數(shù)曲線所需最優(yōu)參數(shù)的方法，并用FORTRAN語(yǔ)言編寫(xiě)了相應(yīng)程序.當(dāng)原始數(shù)據(jù)存在較大噪聲時(shí)，首先要對(duì)數(shù)據(jù)進(jìn)行三次樣條平滑，再開(kāi)始進(jìn)行擬合.擬合過(guò)程包括求期望（E步驟）和最大化對(duì)數(shù)似然函數(shù)（M步驟）兩步：E步驟，根據(jù)初始化擬合參數(shù)或M步驟中產(chǎn)生的參數(shù)計(jì)算完整數(shù)據(jù)對(duì)數(shù)似然函數(shù)的期望；M步驟，最大化E步驟中的對(duì)數(shù)似然函數(shù)的期望產(chǎn)生新的擬合參數(shù)值.經(jīng)過(guò)不斷迭代，最終收斂并得到最佳擬合結(jié)果.該方法的最后結(jié)果與初始值無(wú)關(guān)，只需要給定磁組分?jǐn)?shù)目即可自動(dòng)擬合，提高了IRM曲線分解的效率并得到較為客觀的擬合值，有較為普遍的使用（Roberts et al.,2012；Dorfman et al.,2015；Li et al.,2020）.然而，由于該方法無(wú)法對(duì)樣品的先驗(yàn)地質(zhì)背景信息加以約束限定，所以在對(duì)最終擬合結(jié)果進(jìn)行解譯時(shí)需仔細(xì)推敲其物理或地質(zhì)意義.

1.2 偏態(tài)廣義高斯分布函數(shù)

隨著對(duì)數(shù)高斯分布擬合IRM曲線方法的應(yīng)用和發(fā)展，Heslop等（2004）發(fā)現(xiàn)磁性礦物間的靜磁相互作用和熱弛豫會(huì)使得一些磁組分的矯頑力分布呈現(xiàn)出偏離CLG函數(shù)的偏態(tài)分布.鑒于此，Egli（2003）提出使用偏態(tài)廣義高斯分布（skewed generalized Gaussian,SGG）函數(shù)對(duì)磁性組分進(jìn)行擬合.SGG函數(shù)除了包含CLG函數(shù)的B1/2、DP和Mr三個(gè)參數(shù)外，還另外包含了偏度和峰度兩個(gè)參數(shù)，SGG函數(shù)的概率密度函數(shù)為（圖2）：

圖1 IRM獲得曲線在縱坐標(biāo)為線性，梯度及概率尺度的表現(xiàn)形式（a—c） 分別為線性獲得曲線圖，獲得曲線梯度圖，和標(biāo)準(zhǔn)獲得曲線圖.改自Kruiver等（2001）.Fig.1 IRM acquisition curve on linear,gradient,and probability scale（a—c） are linear acquisition plot （LAP）,gradient of acquisition plot （GAP）,and standardised acquisition plot （SAP）,respectively.Modified from Kruiver et al.（2001）.

（3）

Egli（2004a）通過(guò)分析大量天然樣品的IRM曲線分解結(jié)果發(fā)現(xiàn)，SGG函數(shù)中偏度參數(shù)對(duì)擬合結(jié)果的影響遠(yuǎn)比峰度參數(shù)顯著，在p=2的情況下，通過(guò)調(diào)節(jié)方程（3）中的偏度參數(shù)就可以得到理想的擬合結(jié)果.鑒于此，Maxbauer等（2016）提出，可以使用偏斜正態(tài)（skew-normal distribution）分布擬合IRM曲線，偏斜正態(tài)函數(shù)中沒(méi)有引入峰度參數(shù)p，而是在正態(tài)函數(shù)基礎(chǔ)上引入了控制偏度的參數(shù)（s），s=0時(shí)，蛻變?yōu)檎龖B(tài)函數(shù).Maxbauer等（2016）發(fā)布了基于R語(yǔ)言編寫(xiě)的網(wǎng)頁(yè)應(yīng)用程序MAX UnMix，該程序提供了圖形化用戶(hù)界面，操作便捷.使用時(shí)，用戶(hù)需首先確定磁性組分的數(shù)目，并賦予每個(gè)磁性組分B1/2、DP、p和偏度參數(shù)s的初始猜想值（p為磁性組分對(duì)樣品磁性的貢獻(xiàn)），然后應(yīng)用程序基于這些初始值通過(guò)不斷迭代減小擬合模型和原始數(shù)據(jù)的RSS以得到最佳擬合結(jié)果.需要指出的是，MAX UnMix仍對(duì)IRM曲線的一階導(dǎo)數(shù)數(shù)據(jù)進(jìn)行擬合，因此如果原始數(shù)據(jù)信噪比低，可以用該程序提供的樣條平滑功能進(jìn)行降噪.由于MAX UnMix是網(wǎng)頁(yè)應(yīng)用程序，訪問(wèn)便捷，同時(shí)擬合方法相較于Egli（2003）減少了對(duì)峰度參數(shù)的調(diào)節(jié)從而提高了擬合效率，逐漸成為IRM數(shù)據(jù)磁組分分析的熱門(mén)選擇（Maxbauer et al.,2017；Font et al.,2018；Chang et al.,2019；Abrajevitch,2020；He et al.,2020；Usui and Yamazaki,2021；Wang et al.,2021；Xue et al.,2022）.

圖2 偏態(tài)廣義高斯分布示意圖（a） μ=0，σ=1，且q=1時(shí)的SGG分布.p=1為拉普拉斯分布，p=2為高斯分布，p=∞時(shí)為箱形分布；（b） 左偏SGG分布，μ=0，σ=0.5485.改自Egli（2003）.Fig.2 Examples of SGG distributions（a） SGG distribution with μ=0,σ=1,and q=1.p=1 for a Laplace distribution,p=2 for a Gauss distribution;（b） Left-skewed SGG distributions with μ=0 and σ=0.5485.Modified from Egli （2003）.

1.3 Burr Ⅻ 分布

Zhao等（2018）提出使用12種Burr分布（Burr,1942）類(lèi)型之一的Burr Ⅻ 分布的累積分布函數(shù)擬合樣品的矯頑力分布（圖3）.與SGG和偏斜正態(tài)分布相似，Burr XII分布具有豐富的分布形態(tài)；不同的是，其累積分布具有解析的表達(dá)形式，因此可以直接對(duì)IRM原始數(shù)據(jù)進(jìn)行擬合，能避免對(duì)原始數(shù)據(jù)差分而放大噪音.Burr Ⅻ 分布的累計(jì)分布函數(shù)和概率密度函數(shù)分別為

α>0,γ>0,λ>0，

（4）

α>0,γ>0,λ>0，

（5）

其中，B代表磁場(chǎng)強(qiáng)度，λ是與函數(shù)分布寬度有關(guān)的參數(shù)，分布的形狀由參數(shù)α和γ共同決定.當(dāng)γ>1的時(shí)候，分布是單峰的，α或λ的增大會(huì)導(dǎo)致分布變窄.當(dāng)使用該分布擬合具有多個(gè)磁性組分的樣品時(shí)，累積分布函數(shù)的形式為

（6）

每個(gè)磁性組分由四個(gè)參數(shù)加以描述：磁性組分對(duì)樣品磁性的貢獻(xiàn)ci、尺度參數(shù)λi以及形狀參數(shù)（αi和γi）.

圖3 Burr Ⅻ 分布示意圖（改自Zhao et al.，2018）Fig.3 Example of Burr Ⅻ （modified from Zhao et al.,2018）

對(duì)于磁性組分?jǐn)?shù)目的選擇，Zhao等（2018）采用了赤池信息量準(zhǔn)則（Akaike information criterion,AIC；Akaike,1998），提供了權(quán)衡擬合模型復(fù)雜度（即磁性組分?jǐn)?shù)目）和擬合優(yōu)度的標(biāo)準(zhǔn).AIC方程為

=2k+NlnRSS+C，

（7）

相較于前人提出的需要給定參數(shù)初始猜想值的擬合策略（Kruiver et al.,2001；Egli,2003；Maxbauer et al.,2016），Zhao等（2018）根據(jù)磁性組分的數(shù)目設(shè)定參數(shù)范圍，然后在該范圍內(nèi)自動(dòng)隨機(jī)選取參數(shù)值，找到使得擬合數(shù)據(jù)和原始數(shù)據(jù)的RSS最小的參數(shù)值作為最佳擬合參數(shù).該策略提高了擬合效率并且降低了人為給定參數(shù)初始猜想值的主觀影響.He 等（2020）進(jìn)一步優(yōu)化了該方法的擬合策略.他們首先對(duì)原始數(shù)據(jù)點(diǎn)進(jìn)行多次重采樣（重采樣點(diǎn)數(shù)為原始數(shù)據(jù)點(diǎn)的80%），對(duì)每次重采樣的數(shù)據(jù)用Zhao等（2018）提出的方法進(jìn)行處理后得到磁性組分?jǐn)?shù)目和擬合參數(shù)的估計(jì)值，從中選取使得AIC最小的估計(jì)值作為相關(guān)參數(shù)的初始猜想值進(jìn)行優(yōu)化擬合，得到最佳擬合參數(shù)及其置信區(qū)間.該策略對(duì)減少原始數(shù)據(jù)中的噪聲影響有一定幫助.

2 非參數(shù)化方法

Heslop和Dillon（2007）提出使用非負(fù)矩陣分解算法（non-negative matrix factorization,NMF；Lee and Seung,1999）分解一組樣品的IRM獲得曲線以進(jìn)行端元組分分析的方法.運(yùn)用非負(fù)矩陣是因?yàn)閺奈锢斫嵌葋?lái)說(shuō)，IRM獲得曲線隨外磁場(chǎng)變化增大或不變，如果沒(méi)有該物理約束，數(shù)學(xué)上求解時(shí)由于測(cè)量誤差和噪音等原因會(huì)出現(xiàn)與實(shí)際不符的負(fù)端元組分.該方法不依賴(lài)于任何分布函數(shù)，僅使用端元組分的線性疊加表示每個(gè)樣品的IRM獲得曲線.非負(fù)矩陣為

X=AS+ε，

（8）

矩陣A代表貢獻(xiàn)矩陣，大小為n行m列（n為樣品數(shù)目，m為端元組分?jǐn)?shù)目），每行為m個(gè)端元組分對(duì)每個(gè)樣品IRM獲得曲線的貢獻(xiàn)；矩陣S代表端元組分的IRM獲得曲線矩陣，大小為m行l(wèi)列（m為端元組分?jǐn)?shù)目，l為外磁場(chǎng)強(qiáng)度變化步驟數(shù)目），每行為每個(gè)端元組分剩磁隨外磁場(chǎng)變化的大小，即該端元組分歸一化后的IRM獲得曲線；矩陣X為對(duì)n個(gè)樣品測(cè)量得到的IRM獲得曲線矩陣，大小為n行l(wèi)列，每行對(duì)應(yīng)于每個(gè)樣品在每步外磁場(chǎng)下測(cè)量得到的剩磁.ε為測(cè)量誤差或其他情況對(duì)矩陣X的影響，例如靜磁相互作用.

（9）

（10）

其中i、a、μ是矩陣A和S的索引，i從1到n，a從1到m，μ從1到l.

分解矩陣X前如果數(shù)據(jù)存在較大的噪聲，先用三次樣條法平滑數(shù)據(jù)再進(jìn)行求解.求解時(shí)首先需要給定矩陣A和S的初始猜想值以及端元組分的數(shù)目（m）.矩陣S的初始猜想值用對(duì)矩陣X的模糊c-均值聚類(lèi)給出，矩陣A的初始猜想值根據(jù)非負(fù)最小二乘法得到（Lawson and Hanson,1974）.端元組分?jǐn)?shù)目由每步外磁場(chǎng)上的X和AS之間的確定系數(shù)R2的大小來(lái)判斷，隨著端元組分?jǐn)?shù)目的增多，R2隨之增大而后趨于平緩，一般選擇其拐點(diǎn)位置所對(duì)應(yīng)的端元組分?jǐn)?shù)目.然后將這些初始化的值其代入方程（9）和（10）中進(jìn)行迭代即可得到矩陣A和S的解.

Heslop和Dillon（2007）提出的NMF分解法可以快速求解多個(gè)相似樣品IRM獲得曲線的端元組分，是分析復(fù)雜樣品的有力工具（Dekkers,2012；Just et al.,2012）.然而，該方法對(duì)求解過(guò)程沒(méi)有任何約束且解不具有唯一性，所以需要結(jié)合樣品的地質(zhì)背景對(duì)端元組分進(jìn)行謹(jǐn)慎判斷.該方法要求原始數(shù)據(jù)中的每條IRM獲得曲線單調(diào)遞增，每條曲線有相同數(shù)目的數(shù)據(jù)點(diǎn)且外磁場(chǎng)步長(zhǎng)一致，所以有時(shí)需要對(duì)原始數(shù)據(jù)進(jìn)行預(yù)處理，例如插值等.這在一定程度上限制了使用此方法的便捷性.此外，Heslop和Dillon（2007）建議矩陣X最少應(yīng)包含50個(gè)樣品的剩磁曲線以保證結(jié)果的可靠性.

3 BatchUnMix

CLG函數(shù)作為擬合IRM曲線的最常用函數(shù)之一，發(fā)布的EXCEL版（Kruiver et al.,2001）數(shù)據(jù)分析工具為研究人員提供了很大便利.然而，Kruiver等（2001）提出的擬合策略需要提供擬合參數(shù)的初始猜想值，并且需要不斷手動(dòng)調(diào)節(jié)擬合參數(shù)值以找到最佳擬合效果，這使得IRM曲線處理效率不高且結(jié)果受到人為主觀因素的影響，鑒于此，我們開(kāi)發(fā)了基于MATLAB的GUI工具BatchUnMix，以實(shí)現(xiàn)對(duì)IRM數(shù)據(jù)批量擬合，優(yōu)化IRM曲線分解策略的基礎(chǔ)上提高處理效率，以及增強(qiáng)擬合結(jié)果的客觀性.

3.1 方法描述

BatchUnMix采用了CLG函數(shù)模型（Robertson and France,1994），對(duì)歸一化后的IRM一階導(dǎo)數(shù)曲線進(jìn)行分解.擬合參數(shù)為log（B1/2）、DP和C，其中C為磁性組分對(duì)樣品磁性的貢獻(xiàn).擬合使用了MATLAB擬合工具箱（curve fitting toolbox）中的高斯擬合函數(shù)，其表達(dá)式如方程（11）所示：

（11）

其中i表示第i個(gè)高斯組分，參數(shù)與方程（1）中的擬合磁性組分參數(shù)對(duì)應(yīng)關(guān)系為

Bcr=bi，

（12）

（13）

（14）

最佳擬合參數(shù)采用非線性最小二乘法在一定的參數(shù)范圍內(nèi)自動(dòng)尋找得到.我們參考Egli（2004a）對(duì)大量沉積物樣品IRM曲線的擬合結(jié)果以及一些已發(fā)表文獻(xiàn)中樣品的IRM曲線擬合結(jié)果，在BatchUnMix中設(shè)置了默認(rèn)的log（B1/2）和DP的大致范圍，log（B1/2）的范圍為0～3，DP的范圍為0～0.6.但是仍然建議用戶(hù)根據(jù)對(duì)具體樣品中磁性礦物組分的先驗(yàn)地質(zhì)背景知識(shí)設(shè)置參數(shù)范圍，否則有可能得到數(shù)學(xué)上的最優(yōu)解，但卻不符合樣品所含磁性組分的真實(shí)情況.擬合參數(shù)和原始數(shù)據(jù)的95%置信區(qū)間使用自展法（bootstrap）對(duì)原始數(shù)據(jù)進(jìn)行100次重采樣得到.每次采樣中假設(shè)參數(shù)的誤差為正態(tài)分布，置信區(qū)間為2%，以此計(jì)算該次采樣的最佳擬合參數(shù).如果原始數(shù)據(jù)存在較大噪聲，采用滑動(dòng)平均法對(duì)其進(jìn)行降噪處理.滑動(dòng)平均法是將數(shù)據(jù)中某點(diǎn)附近的N個(gè)采樣點(diǎn)的（N為奇數(shù)）算數(shù)平均作為該點(diǎn)平滑后的值.相比于前人采用的三次樣條降噪，滑動(dòng)平均法的數(shù)學(xué)原理較為簡(jiǎn)單，平滑參數(shù)的選取也僅限于整數(shù)，降低了用戶(hù)的使用難度，可以讓用戶(hù)對(duì)曲線平滑有更好的把握.而三次樣條降噪的平滑參數(shù)選取較難把握，比較容易過(guò)度平滑原始數(shù)據(jù)而丟失關(guān)鍵信號(hào)，或者平滑數(shù)據(jù)不到位使噪聲依然對(duì)分解結(jié)果產(chǎn)生影響.需要指出的是，三次樣條法和滑動(dòng)平均法都是常見(jiàn)的降噪方法，這兩種方法對(duì)于有噪音的IRM數(shù)據(jù)都有較好的平滑效果，但都存在當(dāng)噪音較強(qiáng)時(shí)，無(wú)法很好將噪音和真實(shí)信號(hào)分開(kāi)的缺點(diǎn).

BatchUnMix的主要運(yùn)行窗口如圖4所示，處理IRM曲線的主要流程為（圖5）：（1）選擇處理數(shù)據(jù)方式及數(shù)據(jù)格式：?jiǎn)蝹€(gè)文件處理或者批量處理；振動(dòng)樣品磁強(qiáng)計(jì)（VSM）測(cè)量文件格式或包含外磁場(chǎng)和剩磁信息的xlsx.格式文件；剩磁獲得曲線或者反向場(chǎng)退磁曲線；是否對(duì)擬合參數(shù)和原始數(shù)據(jù)求取置信區(qū)間（圖4a）；（2）選擇輸入文件，根據(jù)原始數(shù)據(jù)和數(shù)據(jù)的一階導(dǎo)數(shù)圖判斷是否需要對(duì)數(shù)據(jù)進(jìn)行平滑；（3）選擇擬合所需磁性組分的數(shù)目，然后對(duì)參數(shù)（log（B1/2）、DP和C）的范圍進(jìn)行設(shè)置并開(kāi)始擬合（圖4b）；（4）根據(jù)得到的結(jié)果判斷是否需要返回上步重新調(diào)節(jié)參數(shù)范圍再次擬合（圖4c）；（5）最終結(jié)果自動(dòng)保存為fig.格式圖片和xlsx.格式文件.

3.2 BatchUnMix擬合示例及擬合結(jié)果比較

我們用BatchUnMix對(duì)典型沉積物及火成巖樣品的IRM數(shù)據(jù)分別進(jìn)行了單樣品處理和批量處理，并將擬合結(jié)果和已發(fā)表的使用MAX UnMix（Maxbauer et al.,2016）分析結(jié)果進(jìn)行了比較.分析使用的沉積物樣品來(lái)自國(guó)際大洋發(fā)現(xiàn)計(jì)劃（IODP）位于北大西洋的U1403B鉆孔（Xue et al.,2022），火成巖樣品組來(lái)自中國(guó)大洋航次西南印度洋中脊采集的40塊玄武巖樣品（Wang et al.,2021）.這兩組樣品的磁性組分鑒定有電子顯微直接觀測(cè)結(jié)果的約束——U1403B鉆孔中的磁性礦物組成為碎屑磁鐵礦與軟磁及硬磁性生物磁鐵礦（Xue et al.,2022），西南印度洋中脊玄武巖中的磁性礦物組為微米級(jí)鈦磁鐵礦枝晶與納米級(jí)鈦磁鐵礦包裹體（Wang et al.,2021）.沉積物中的生物磁鐵礦具有指示深海環(huán)境意義（Chang et al.,2018），而洋中脊玄武巖是了解洋脊磁異常機(jī)理和洋殼水巖相互作用等科學(xué)問(wèn)題的重要介質(zhì)（劉隆等，2021；Wang et al.,2021），所以選用這兩類(lèi)樣品展示BatchUnMix擬合效果并與已有結(jié)果對(duì)比.

3.2.1 單樣品處理

對(duì)于U1403B沉積物樣品，參考Xue等（2022）的擬合結(jié)果，首先確定用四個(gè)磁性組分?jǐn)M合IRM一階導(dǎo)數(shù)曲線，然后我們使用了不同的參數(shù)設(shè)置方式以對(duì)比分析：（1）設(shè)置三個(gè)組分的log（B1/2）范圍分別為：1.4～1.5，1.7～1.8和2.3～2.6，其他參數(shù)為默認(rèn)范圍；（2）只設(shè)置高場(chǎng)成分的log（B1/2）范圍為2～3，其他參數(shù)為默認(rèn)范圍.這兩種擬合方式得到了相似的結(jié)果（圖6b,c）：包含兩個(gè)狹窄的矯頑力分布（較小的DP）以及另外兩個(gè)低場(chǎng)和高場(chǎng)組分，顯示了和Xue等（2022）擬合結(jié)果的一致性（圖6a；表1）.Xue等（2022）結(jié)合透射電鏡數(shù)據(jù)，將擬合結(jié)果解釋為生物成因軟磁組分（biogenic soft,BS）和生物成因硬磁組分（biogenic hard,BH），是沉積物中磁鐵礦化石磁小體的典型特征（Egli,2004a），其余兩個(gè)磁性組分分別對(duì)應(yīng)低矯頑力的碎屑磁鐵礦和高矯頑力的赤鐵礦.

圖4 BatchUnMix程序主要運(yùn)行窗口截圖（a） 數(shù)據(jù)選擇窗口；（b） 參數(shù)設(shè)置窗口；（c） 擬合結(jié)果窗口.Fig.4 Screenshots of main user interface of program BatchUnMix（a） Data selection window;（b） Parameter setting window;（c） Fitting result window.

該示例表明，當(dāng)具備對(duì)樣品所含磁性組分的先驗(yàn)知識(shí)時(shí)，可以?xún)H對(duì)較為確定的磁性組分參數(shù)范圍進(jìn)行設(shè)置（圖6c），其他參數(shù)范圍保持默認(rèn)設(shè)置，所得到的擬合結(jié)果和對(duì)參數(shù)范圍進(jìn)行精細(xì)設(shè)置的結(jié)果相差不大（圖6b）.

3.2.2 樣品批量處理

我們對(duì)于40塊西南印度洋中脊玄武巖樣品的IRM一階導(dǎo)數(shù)曲線分為兩種情況進(jìn)行批量處理：（1）同時(shí)設(shè)置log（B1/2）和DP范圍，（2）只設(shè)置log（B1/2）范圍.然后將擬合得到的磁性組分參數(shù)進(jìn)行了100次迭代的K均值聚類(lèi)分析，并與Wang 等（2021）的統(tǒng)計(jì)結(jié)果進(jìn)行了比較.情況（1）得到了和Wang等（2021）相似的統(tǒng)計(jì)結(jié)果（圖7a,b）.情況（2）雖然得到了和圖7a，b相似的聚類(lèi)中心（圖7c），但兩個(gè)聚類(lèi)中心重合度較高.這是由于相比使用MAX UnMix和情況（1）處理的IRM曲線（例如，圖7d,e），當(dāng)參數(shù)設(shè)置為情況（2）時(shí)，部分IRM曲線的擬合結(jié)果并沒(méi)有得到更符合真實(shí)的組分（圖7f）.

圖5 BatchUnMix操作流程圖.虛線框表示該功能僅在處理單個(gè)樣品中具備Fig.5 Procedures of data analysis in BatchUnMix.Dotted boxes indicate that the function is only applicable for single sample processing

表1 使用MAX UnMix和BatchUnMix處理北大西洋沉積物樣品結(jié)果（BatchUnMix采用了兩種參數(shù)設(shè)置）Table 1 Results of decomposing IRM curves for a sediment sample from the North Atlantic using MAX UnMix and BatchUnMix,respectively （BatchUnMix uses two different parameter settings）

圖6 對(duì)比不同分解方法及BatchUnMix不同參數(shù)設(shè)置下處理北大西洋深海沉積物樣品的同一IRM數(shù)據(jù)示例（a） 使用MAX UnMix（Maxbauer et al.,2016）的處理結(jié)果（改自Xue et al.,2022）；（b,c） 使用BatchUnMix在不同參數(shù)設(shè)置情況下的處理結(jié)果.Fig.6 Examples of decomposing IRM curves with different methods and parameter settings in BatchUnMix for the same IRM dataset of a North Atlantic pelagic sediment sample（a） shows result using MAX UnMix modified from Xue et al.（2022）;（b,c） are results using BatchUnMix with two different parameter settings （see text）.

結(jié)合本例及對(duì)其他樣品已發(fā)表數(shù)據(jù)批量處理的測(cè)試，我們發(fā)現(xiàn)批處理過(guò)程對(duì)參數(shù)范圍設(shè)置精細(xì)度要求較高.所以建議在批量處理數(shù)據(jù)前，先隨機(jī)挑選部分樣品對(duì)其剩磁曲線進(jìn)行單獨(dú)擬合，得到相關(guān)磁性組分參數(shù)的大致范圍，由部分估計(jì)整體，再進(jìn)行批量擬合.此外，批量處理一般較適用于所含磁性組分相似的樣品，當(dāng)樣品間磁性特征相差較大時(shí)，可能得不到真實(shí)解.

4 IRM曲線分解方法的局限性

盡管擬合IRM曲線的函數(shù)模型和分解策略經(jīng)過(guò)了幾十年的發(fā)展逐漸完善，但使用IRM曲線分解方法開(kāi)展磁性組分分析仍然存在一定的局限性.

首先，IRM曲線的分解結(jié)果具有多解性.同一樣品的IRM曲線分解時(shí)，在滿(mǎn)足數(shù)學(xué)上的擬合優(yōu)度后，可能有不同的分解結(jié)果.例如，圖8a展示了對(duì)來(lái)自孟加拉灣沉積物樣品的IRM曲線分解后的兩種不同結(jié)果，結(jié)合該樣品的電鏡照片（Xue et al.,2019），這兩種擬合結(jié)果都有一定的合理性，都可以解釋為不同粒度的磁鐵礦導(dǎo)致的兩種不同的磁性組分.SGG函數(shù)和Burr Ⅻ分布具有高度靈活的特點(diǎn)，所以它們相對(duì)于CLG函數(shù)多解性問(wèn)題更為突出（Zhao et al.,2018）.此外，當(dāng)樣品的IRM曲線可以由不同數(shù)目的磁性組分?jǐn)M合時(shí)，增加磁性組分?jǐn)?shù)目通常會(huì)提高擬合優(yōu)度，但是得到的磁性組分可能不具有真實(shí)物理意義.此時(shí)，結(jié)合Heslop和Dillon（2007）提出的非參數(shù)化方法對(duì)樣品的IRM曲線進(jìn)行磁性組分分析可以對(duì)磁性組分?jǐn)?shù)目的選擇起到一定的幫助（Wang and Chang,2022）.

IRM曲線分解結(jié)果對(duì)擬合函數(shù)模型的選取有一定的依賴(lài)性（Egli,2003;Zhao et al.,2018），只有選取與樣品中的真實(shí)組分最接近的函數(shù)模型，才能更好的重建組分信息.圖8b展示了對(duì)合成磁鐵礦樣品IRM曲線擬合時(shí)，用CLG函數(shù)、Burr Ⅻ分布和SGG函數(shù)分別對(duì)其分解后得到了不同的結(jié)果（Zhao et al.,2018）.此外，Egli（2003）中的圖10展示了存在使用一個(gè)SGG組分?jǐn)M合IRM曲線的效果和使用兩個(gè)對(duì)數(shù)正態(tài)分布組分?jǐn)M合IRM曲線的效果相似的情況.所以，當(dāng)使用不同函數(shù)模型擬合IRM曲線后得到差異較大的結(jié)果時(shí)，建議研究人員改變擬合參數(shù)對(duì)樣品的IRM曲線重新進(jìn)行擬合，或者結(jié)合該樣品的其他巖石磁學(xué)數(shù)據(jù)對(duì)擬合函數(shù)模型做出選擇.

此外，樣品中同一物理組分可能需要用多個(gè)數(shù)學(xué)組分描述（He et al.,2020;Bai et al.,2021）.例如，單一尺寸分布的磁鐵礦顆粒集合體，如果顆粒間存在較強(qiáng)的相互作用，那么擬合其IRM曲線需要兩個(gè)數(shù)學(xué)組分（Bai et al.,2021）.噪聲也會(huì)對(duì)IRM曲線的分解結(jié)果有一定的影響（Zhao et al.,2018），CLG函數(shù)相比于SGG、skew-normal、Burr Ⅻ分布有較少的擬合參數(shù)，所以對(duì)噪聲的敏感度更低.對(duì)于有噪聲的IRM數(shù)據(jù)，除了對(duì)其進(jìn)行平滑外，可以使用擬合函數(shù)的cumulative distribution function（CDF）和probability distribution function（PDF）對(duì)IRM曲線分別進(jìn)行擬合并比較結(jié)果，Zhao等（2018）的研究指出，IRM原始數(shù)據(jù)與其一階導(dǎo)數(shù)曲線有不同的噪聲敏感度，所以用擬合函數(shù)模型的CDF擬合IRM原始數(shù)據(jù)和PDF擬合IRM一階導(dǎo)數(shù)曲線就數(shù)據(jù)噪聲而言可以視為兩種獨(dú)立的擬合IRM曲線方法.

圖7 對(duì)比不同IRM曲線分解方法及BatchUnMix不同參數(shù)設(shè)置下處理西南印度洋中脊玄武巖同一組IRM數(shù)據(jù)示例（a—c）對(duì)40塊樣品批量處理后再進(jìn)行K均值聚類(lèi)分析結(jié)果；（a） 使用MAX UnMix（Maxbauer et al.,2016）處理的結(jié)果，改自Wang等（2021）；（b） 在BatchUnMix中同時(shí)設(shè)置log（B1/2）和DP范圍時(shí)的處理結(jié)果；（c） 在BatchUnMix中只設(shè)置log（B1/2）范圍時(shí)的處理結(jié)果；（d—f） 分別對(duì)應(yīng)圖（a—c）批量處理中的單樣品處理結(jié)果示例.Fig.7 Examples of decomposing the same IRM dataset with different methods and parameter settings in BatchUnMix for basalt samples from the Southwest Indian Ridge（a—c） show results of K-means clustering analysis for unmixed IRM components of 40 basalt samples;（a） Result using MAX UnMix,data modified from Wang et al.（2021）;（b） Data processed using BatchUnMix with setting log（B1/2） and DP ranges simultaneously;（c） Data processed using BatchUnMix with only setting log （B1/2） ranges;（d—f） show examples of IRM fitting for single sample from results in （a—c）,respectively.

圖8 IRM曲線分解的多解性及模型依賴(lài)性示例（a） CLG函數(shù)分解IRM曲線多解性示例；樣品為孟加拉灣表層沉積物（Xue et al.,2019）；（b） IRM曲線分解模型依賴(lài)性示例（改繪自 Zhao et al.,2018）.Fig.8 Examples of multiplicity and model dependency for IRM curve decomposition（a） Example of multiplicity for IRM curve decomposition for surface sediment sample from the central Bengal Fan （Xue et al.,2019）;（b） Example of model dependency for IRM curve decomposition.

鑒于上述提到的IRM曲線分解方法存在的局限性，在處理樣品的IRM曲線時(shí)，需要具備對(duì)樣品的先驗(yàn)判斷，并結(jié)合其他磁性參數(shù)，例如磁滯回線、一階反轉(zhuǎn)曲線等，以及運(yùn)用電鏡礦物成像等方法進(jìn)行綜合分析判斷.單一對(duì)IRM曲線進(jìn)行磁組分分析往往不能完全清楚地揭示復(fù)雜樣品的磁性組分信息.

近年來(lái)微磁模擬的快速發(fā)展為建立磁性組分與樣品真實(shí)組分的可靠聯(lián)系提供了一種有效方法（Bai et al.,2021,2022）.微磁模擬可以分別評(píng)估磁性礦物的粒度、形態(tài)及彼此間存在的靜磁相互作用等對(duì)IRM曲線分解的影響，從而為分離的磁組分解釋提供可靠的理論支撐.例如，Bai等（2021）使用微磁模擬的方法系統(tǒng)研究了不同形狀（拉長(zhǎng)和扁平球體、正方體和八面體）和空間分布（不同顆粒間距）的磁鐵礦顆粒集合體對(duì)IRM曲線形狀的影響.研究表明，對(duì)于拉長(zhǎng)顆粒集合體，即使集合體內(nèi)的顆粒粒度分布單一，顆粒間的靜磁相互作用也會(huì)導(dǎo)致其IRM曲線分解得到的高斯組分與富含生物磁小體化石的沉積物磁組分分析結(jié)果相似.

5 小結(jié)與展望

本文系統(tǒng)回顧了IRM曲線分解方法，并基于CLG函數(shù)提出了新的IRM曲線分解MATLAB工具BatchUnMix（下載地址：https:∥doi.org/10.18170/DVN/RVQGLE），該工具加入了批量處理功能，通過(guò)設(shè)定部分?jǐn)M合參數(shù)的范圍自動(dòng)尋找其最優(yōu)解，降低了人為給定擬合參數(shù)初始猜想值的主觀性影響，并顯著提高了IRM數(shù)據(jù)處理效率.

雖然IRM曲線分解是鑒別樣品中磁性礦物組分的有力方法，然而其本身存在的局限性使得研究人員在使用該方法時(shí)需格外謹(jǐn)慎.未來(lái)應(yīng)對(duì)影響IRM曲線分解的因素進(jìn)行更為深入的研究，例如使用微磁模擬的方法，以構(gòu)建起數(shù)學(xué)解和真實(shí)解之間的橋梁，從而更為清楚準(zhǔn)確地解譯樣品中所含磁性組分的信息.

致謝感謝北京大學(xué)的Thomas Berndt博士、澳大利亞國(guó)立大學(xué)的趙翔博士、上海交通大學(xué)的趙翔宇博士的有益討論.感謝中國(guó)海洋大學(xué)的姜兆霞教授和另一位匿名審稿人的建議使本文獲得了提升.