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

白帆，常燎,2*，薛鵬飛，汪詩舜

1 北京大學地球與空間科學學院造山帶與地殼演化教育部重點實驗室，北京 100871 2 青島海洋科學與技術試點國家實驗室海洋地質過程與環境功能實驗室，青島 266071

0 引言

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

本文首先從數學模型的角度綜述了近年來提出的IRM曲線分解方法，然后重點介紹我們基于MATLAB開發的GUI工具BatchUnMix，并利用實驗數據展示和評估了BatchUnMix處理天然樣品IRM曲線的實用性，最后討論了IRM曲線分解方法的局限性并提出了相關方面有待進一步研究的科學問題.

1 參數化方法

1.1 累積對數高斯函數

Robertson和France（1994）通過分析磁赤鐵礦、赤鐵礦和針鐵礦樣品的IRM獲得曲線發現，單一磁性礦物的IRM獲得曲線形態可以用累積對數高斯（cumulative log-Gaussian,CLG）分布進行擬合：

（1）

其中，B為磁性礦物獲取剩磁時施加的外磁場大?。籅1/2為平均剩磁矯頑力，為磁性礦物獲得半飽和剩磁（SIRM）時的外磁場大??；DP（dispersion parameter）為離散參數，對應于對數-高斯函數的一個標準差；Mr為該單一磁性組分的飽和剩磁.Robertson和France（1994）進一步的研究表明，磁性礦物顆粒間不存在靜磁相互作用的情況下，N種礦物混合后的IRM獲得曲線是每種礦物所對應的IRM獲得曲線的線性疊加，即可以用每種礦物對應的CLG函數之和擬合磁性礦物混合后的IRM獲得曲線，對應的方程為

（2）

其中i表示第i個磁性成分.

Eyre（1996）和McIntosh等（1996）在研究中國黃土時采用了上述分析方法，利用對數高斯概率密度函數擬合樣品中每個磁性組分的IRM獲得曲線的一階導數曲線，從而得到樣品的矯頑力分布.然而在擬合IRM曲線的過程中，他們對磁性組分的數目及擬合效果的好壞是通過直接觀察判斷的，存在很大程度的主觀因素.此后，Stockhausen（1998）提出使用三個基于殘差平方和（RSS）的參數來判斷IRM曲線擬合優度，并用FORTRAN語言編寫了擬合程序，但他并沒有給出何時使用這三個參數的標準.

Kruiver等（2001）基于前人提出的累積對數高斯分布及對數高斯概率密度分布的擬合方法，又提出了用正態分布的標準分數（Z分數）表示每個磁性組分的IRM曲線，稱其為標準獲得曲線（standardised acquisition plot,SAP；圖1c）.如果樣品只含有一個磁性組分，那么其IRM曲線的SAP為一條直線.如果樣品含有多個磁性組分，那么其IRM曲線的SAP為各個磁性組分的SAP疊加，會呈現出曲線的形態.例如，當樣品中存在赤鐵礦等高矯頑力磁性組分時，其IRM曲線的SAP高場部分會有明顯的彎曲（例如，Abrajevitch et al.,2009）.結合IRM獲得曲線（也被稱為線性獲得曲線，LAP；圖1a）、標準獲得曲線（SAP；圖1c）及IRM一階導數曲線（也被稱為梯度曲線，gradient of acquisition plot，GAP；圖1b），可以利用不同數量的CLG曲線，通過調整它們各自的log（B1/2）、DP和Mr參數對IRM曲線進行組分分析（Kruiver et al.,2001）.對于擬合過程中存在的多解性問題，例如可以用不同數目的磁性組分擬合剩磁曲線，他們通過對不同擬合結果與原始數據的殘差平方和（RSS）進行F檢驗和t檢驗來判斷兩個模型是否存在顯著差異，從而得到最優擬合參數.基于以上方法開發的EXCEL工具IRM-CLG已在古地磁學及環境磁學研究得到了廣泛使用（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）提出的方法都需要不斷手動調節每個磁性組分IRM曲線的擬合參數（log（B1/2）、DP和Mr）.當樣品的組分較為復雜時，例如當樣品中包含四個磁性組分，需要調節的擬合參數多達12個，此時擬合一個樣品的IRM曲線需要花費較長時間.Heslop等（2002）提出了基于對數似然和期望最大化算法（Dempster et al.,1977）自動尋找擬合IRM一階導數曲線所需最優參數的方法，并用FORTRAN語言編寫了相應程序.當原始數據存在較大噪聲時，首先要對數據進行三次樣條平滑，再開始進行擬合.擬合過程包括求期望（E步驟）和最大化對數似然函數（M步驟）兩步：E步驟，根據初始化擬合參數或M步驟中產生的參數計算完整數據對數似然函數的期望；M步驟，最大化E步驟中的對數似然函數的期望產生新的擬合參數值.經過不斷迭代，最終收斂并得到最佳擬合結果.該方法的最后結果與初始值無關，只需要給定磁組分數目即可自動擬合，提高了IRM曲線分解的效率并得到較為客觀的擬合值，有較為普遍的使用（Roberts et al.,2012；Dorfman et al.,2015；Li et al.,2020）.然而，由于該方法無法對樣品的先驗地質背景信息加以約束限定，所以在對最終擬合結果進行解譯時需仔細推敲其物理或地質意義.

1.2 偏態廣義高斯分布函數

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

圖1 IRM獲得曲線在縱坐標為線性，梯度及概率尺度的表現形式（a—c） 分別為線性獲得曲線圖，獲得曲線梯度圖，和標準獲得曲線圖.改自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）通過分析大量天然樣品的IRM曲線分解結果發現，SGG函數中偏度參數對擬合結果的影響遠比峰度參數顯著，在p=2的情況下，通過調節方程（3）中的偏度參數就可以得到理想的擬合結果.鑒于此，Maxbauer等（2016）提出，可以使用偏斜正態（skew-normal distribution）分布擬合IRM曲線，偏斜正態函數中沒有引入峰度參數p，而是在正態函數基礎上引入了控制偏度的參數（s），s=0時，蛻變為正態函數.Maxbauer等（2016）發布了基于R語言編寫的網頁應用程序MAX UnMix，該程序提供了圖形化用戶界面，操作便捷.使用時，用戶需首先確定磁性組分的數目，并賦予每個磁性組分B1/2、DP、p和偏度參數s的初始猜想值（p為磁性組分對樣品磁性的貢獻），然后應用程序基于這些初始值通過不斷迭代減小擬合模型和原始數據的RSS以得到最佳擬合結果.需要指出的是，MAX UnMix仍對IRM曲線的一階導數數據進行擬合，因此如果原始數據信噪比低，可以用該程序提供的樣條平滑功能進行降噪.由于MAX UnMix是網頁應用程序，訪問便捷，同時擬合方法相較于Egli（2003）減少了對峰度參數的調節從而提高了擬合效率，逐漸成為IRM數據磁組分分析的熱門選擇（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 偏態廣義高斯分布示意圖（a） μ=0，σ=1，且q=1時的SGG分布.p=1為拉普拉斯分布，p=2為高斯分布，p=∞時為箱形分布；（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）類型之一的Burr Ⅻ 分布的累積分布函數擬合樣品的矯頑力分布（圖3）.與SGG和偏斜正態分布相似，Burr XII分布具有豐富的分布形態；不同的是，其累積分布具有解析的表達形式，因此可以直接對IRM原始數據進行擬合，能避免對原始數據差分而放大噪音.Burr Ⅻ 分布的累計分布函數和概率密度函數分別為

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

（4）

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

（5）

其中，B代表磁場強度，λ是與函數分布寬度有關的參數，分布的形狀由參數α和γ共同決定.當γ>1的時候，分布是單峰的，α或λ的增大會導致分布變窄.當使用該分布擬合具有多個磁性組分的樣品時，累積分布函數的形式為

（6）

每個磁性組分由四個參數加以描述：磁性組分對樣品磁性的貢獻ci、尺度參數λi以及形狀參數（αi和γi）.

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

對于磁性組分數目的選擇，Zhao等（2018）采用了赤池信息量準則（Akaike information criterion,AIC；Akaike,1998），提供了權衡擬合模型復雜度（即磁性組分數目）和擬合優度的標準.AIC方程為

=2k+NlnRSS+C，

（7）

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

2 非參數化方法

Heslop和Dillon（2007）提出使用非負矩陣分解算法（non-negative matrix factorization,NMF；Lee and Seung,1999）分解一組樣品的IRM獲得曲線以進行端元組分分析的方法.運用非負矩陣是因為從物理角度來說，IRM獲得曲線隨外磁場變化增大或不變，如果沒有該物理約束，數學上求解時由于測量誤差和噪音等原因會出現與實際不符的負端元組分.該方法不依賴于任何分布函數，僅使用端元組分的線性疊加表示每個樣品的IRM獲得曲線.非負矩陣為

X=AS+ε，

（8）

矩陣A代表貢獻矩陣，大小為n行m列（n為樣品數目，m為端元組分數目），每行為m個端元組分對每個樣品IRM獲得曲線的貢獻；矩陣S代表端元組分的IRM獲得曲線矩陣，大小為m行l列（m為端元組分數目，l為外磁場強度變化步驟數目），每行為每個端元組分剩磁隨外磁場變化的大小，即該端元組分歸一化后的IRM獲得曲線；矩陣X為對n個樣品測量得到的IRM獲得曲線矩陣，大小為n行l列，每行對應于每個樣品在每步外磁場下測量得到的剩磁.ε為測量誤差或其他情況對矩陣X的影響，例如靜磁相互作用.

（9）

（10）

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

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

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

3 BatchUnMix

CLG函數作為擬合IRM曲線的最常用函數之一，發布的EXCEL版（Kruiver et al.,2001）數據分析工具為研究人員提供了很大便利.然而，Kruiver等（2001）提出的擬合策略需要提供擬合參數的初始猜想值，并且需要不斷手動調節擬合參數值以找到最佳擬合效果，這使得IRM曲線處理效率不高且結果受到人為主觀因素的影響，鑒于此，我們開發了基于MATLAB的GUI工具BatchUnMix，以實現對IRM數據批量擬合，優化IRM曲線分解策略的基礎上提高處理效率，以及增強擬合結果的客觀性.

3.1 方法描述

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

（11）

其中i表示第i個高斯組分，參數與方程（1）中的擬合磁性組分參數對應關系為

Bcr=bi，

（12）

（13）

（14）

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

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

3.2 BatchUnMix擬合示例及擬合結果比較

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

3.2.1 單樣品處理

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

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

該示例表明，當具備對樣品所含磁性組分的先驗知識時，可以僅對較為確定的磁性組分參數范圍進行設置（圖6c），其他參數范圍保持默認設置，所得到的擬合結果和對參數范圍進行精細設置的結果相差不大（圖6b）.

3.2.2 樣品批量處理

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

圖5 BatchUnMix操作流程圖.虛線框表示該功能僅在處理單個樣品中具備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處理北大西洋沉積物樣品結果（BatchUnMix采用了兩種參數設置）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 對比不同分解方法及BatchUnMix不同參數設置下處理北大西洋深海沉積物樣品的同一IRM數據示例（a） 使用MAX UnMix（Maxbauer et al.,2016）的處理結果（改自Xue et al.,2022）；（b,c） 使用BatchUnMix在不同參數設置情況下的處理結果.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）.

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

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

盡管擬合IRM曲線的函數模型和分解策略經過了幾十年的發展逐漸完善，但使用IRM曲線分解方法開展磁性組分分析仍然存在一定的局限性.

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

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

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

圖7 對比不同IRM曲線分解方法及BatchUnMix不同參數設置下處理西南印度洋中脊玄武巖同一組IRM數據示例（a—c）對40塊樣品批量處理后再進行K均值聚類分析結果；（a） 使用MAX UnMix（Maxbauer et al.,2016）處理的結果，改自Wang等（2021）；（b） 在BatchUnMix中同時設置log（B1/2）和DP范圍時的處理結果；（c） 在BatchUnMix中只設置log（B1/2）范圍時的處理結果；（d—f） 分別對應圖（a—c）批量處理中的單樣品處理結果示例.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曲線分解的多解性及模型依賴性示例（a） CLG函數分解IRM曲線多解性示例；樣品為孟加拉灣表層沉積物（Xue et al.,2019）；（b） IRM曲線分解模型依賴性示例（改繪自 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曲線時，需要具備對樣品的先驗判斷，并結合其他磁性參數，例如磁滯回線、一階反轉曲線等，以及運用電鏡礦物成像等方法進行綜合分析判斷.單一對IRM曲線進行磁組分分析往往不能完全清楚地揭示復雜樣品的磁性組分信息.

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

5 小結與展望

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

雖然IRM曲線分解是鑒別樣品中磁性礦物組分的有力方法，然而其本身存在的局限性使得研究人員在使用該方法時需格外謹慎.未來應對影響IRM曲線分解的因素進行更為深入的研究，例如使用微磁模擬的方法，以構建起數學解和真實解之間的橋梁，從而更為清楚準確地解譯樣品中所含磁性組分的信息.

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