梯度下降法在沉积物粒度分布拟合中的应用

梯度下降法在沉积物粒度分布拟合中的应用①①国家863方案〔2007AA09Z118〕、国家973方案〔2005CB422308〕、高等学校博士学科点专项科研基金〔20050423007〕和国家自然科学基金〔40472088〕资助。陈海波②男，1984年生，博士生；研究方向：海洋动力学。吕咸青③通讯作者，E-mail：xqinglv@②男，1984年生，博士生；研究方向：海洋动力学。③通讯作者，E-mail：xqinglv@〔中国海洋大学物理海洋实验室青岛266100〕〔中国地质科学院地质力学研究所北京100081〕摘要本文结合实测数据，以三个对数正态分布函数的和函数为拟合函数，以梯度下降法为主要方法，对沉积物粒度分布进行了数据拟合。通过数值实验我们发现，利用梯度下降法可以有效地优化分布函数的各参数，实现拟合残差的稳步持续减小，具有良好的可操作性；对64个实测样本进行数据拟合，其拟合标准差大局部在0.1%以下，在总体上优于相关文献中的拟合结果，因此利用梯度下降法的拟合效果是令人满意的。梯度下降法为我们进行数据拟合提供了一条新的思路，同时此方法也可以推广到解决其他极值问题。关键词非线性数据拟合，最小二乘，梯度，3个对数正态分布的叠加，激光粒度0引言黄土、冰芯、深海泥土等沉积物被认为是全球环境变化的良好信息载体,反映了第四纪以来的全球气候变化情况。粒度是沉积物沉积的主要特征，是分析沉积环境、沉积及搬运过程和搬运机制的重要参量，同时它也是沉积物一个比拟成熟的古环境指标，因其测定简单、快速、物理意义明确、对气候变化敏感等特点近年来被广泛应用于各种沉积环境的研究中。沉积物粒度的各组分的成因是不相同的，要受到物源区和气候的变化影响。在沉积物的粒度组成中，不同粒级的含量变化具有不同的古气候意义。目前,粒度分析已广泛应用于黄土与古土壤、冰川沉积、海洋、河流、湖泊等多种沉积环境的研究，充分利用各个粒级测量数据的信息对于准确区别不同沉积物的沉积环境类型具有重大意义。自然界大多沉积物均受一种或几种不同的搬运方式、动力类型控制,因而产生了多组分、多模态粒度分布特征,在频率曲线上表现为多峰光滑曲线。从沉积物粒度组分分布特征看,每一组分均属于对数正态分布类型。因对数正态分布的参数均具有明确的物理意义,因此能够合理地表现样品中每一组分的粒度分布特征及其成因意义。因此，本文利用由多个对数正态分布组成的叠加分布对沉积物粒度分布进行数据拟合，同时对各粒度组分进行了数学别离。如何利用多组分分布函数拟合沉积物粒度分布，孙东怀等曾做过大量的研究工作，但由于他们的研究更侧重于分析结果的古环境意义，因此他们只给出了局部拟合结果，而很少涉及到拟合方法的实施过程，并且我们认为拟合精度也是可以提高的。对于由个对数正态分布组成的叠加分布，其概率密度函数应为：其中自变量表示粒径值，和为第个组分的分布函数的位置参数和形状参数，为第个组分的分布函数在总体分布函数密度中的百分比。通过自变量变换，令，我们可以将对数正态分布转化为正态分布的形式：在表达方式上更为简洁。利用多个正态分布函数的和函数对沉积物粒度分布进行拟合，实际上是对参数，和进行估计，最小二乘法是应用最广泛的参数估计方法之一，关于非线性最小二乘问题数值解法也有很多。在应用迭代法寻求使残差的平方和最小的参数之前，我们十分有必要答复下面两个问题：1〕对于特定的问题，参数的最小二乘估计是否存在？2〕怎样选取一个尽可能好的参数初始近似值？对非线性最小二乘问题来说，答复下列问题1〕是极其困难的，但对于拟合函数是一些特殊函数的情况，一些学者给出了相应的最小二乘估计的存在性理论，Jukić讨论了拟合Gaussian型曲线的最小二乘估计的存在性定理，正态分布函数实际上是一种特殊形式的Gaussian函数，参照文献[10]中相关理论的证明方法，我们可以很容易地证明本文中所给定的条件可以充分保证小二乘估计的存在性，证明过程本文不再赘述。对于问题2〕，确定初始近似值的方法有很多，并且根据研究对象的不同以及优化方法的不同，初始近似值的选取方法也不尽相同。关于如何选取非线性最小二乘问题的初始近似值，一些学者提出了有价值的方法：Hölmstrom和Petersson将改良的Prony法应用于初始近似值的选取，他们所应用的拟合函数为假设干个指数函数的叠加函数；对于拟合函数是一些特殊函数的情况，Jukić和Marković等给出了一些易于实现的选取初值的方法；Hasdorff利用梯度推导出了初值的选取公式。本文利用将非线性最小二乘问题线性化的方法，来确定参数的初始近似值。非线性最小二乘问题的数值解法多种多样，究竟选取哪种方法，也要视具体问题而定。Gauss-Newton法是求解非线性最小二乘问题最根本的方法，已被很多学者应用于非线性数据拟合，其中Hölmstrom和Petersson还介绍了拟Newton法和Levenberg-Marquardt法〔LM法〕等其他优化方法的应用，Atieg和Watson对各种数值拟合方法作了总结和比照，其中主要包括Newton法和Gauss-Newton法；Böckmann提出了改良的trust-regionGauss-Newton法，并将其应用于拟合函数是叠加函数的情形；Nyarko和Scitovski将遗传算法〔geneticalgorithms〕应用于参数的辨识；Ahn等给出了应用几何拟合方法拟合圆、球、椭圆、双曲线和抛物线的算例。其中Hölmstrom和Atieg等应用的拟合函数是假设干个函数的线性叠加。梯度下降法是以负梯度方向作为搜索方向的方法，它是所有需要计算导数的无约束优化算法中最简单的方法，将其应用于非线性数据拟合是非常有益的探讨。在本文中，我们结合甘孜黄土激光粒度分析数据，探讨了利用由三个对数正态分布组成的叠加分布进行数据拟合的方法。在第1局部，我们把沉积物粒度分布的数据拟合以非线性最小二乘问题的形式进行描述。在第2局部，我们利用线性化的方法来确定参数的初始近似值。本文介绍线性化方法的目的有两个：一是为梯度下降法提供一个合理的初值来源，二是将其与梯度下降法作比照。在第3局部绍了我们主要的数据拟合方法——梯度下降法。第4局部介绍了我们的数值实验。数值实验说明：梯度下降法是一种有效的数据拟合方法，总体而言，本文拟合结果的精度优于文献[1-3]中的结果。1问题的描述给定组被拟合数据，假设它们可以由一个模型：得到。其中是含有个未知参数的拟合函数，为未知参数向量，为期望为零的随机误差函数，，。在可容许的参数空间内寻求一点，使，〔1〕这种方法称为最小二乘法〔LSmethod〕或最小二乘逼近〔LSapproach〕，而最小参数值称为未知参数关于问题〔1〕的最小二乘估计〔LSestimate〕，在这里我们称为代价函数。我们将拟合函数构造为3个正态分布的叠加函数，即〔2〕，其中，和为第个正态分布函数的位置参数和形状参数，记，为第个正态分布函数在总体分布函数密度中的百分比，且满足：在本问题中，拟合函数共含有个未知参数，而、、中只有两个变量是独立的，因而拟合函数实际上含有8个独立的未知参数。这样，未知参数向量可表示为，其中，，；；。2初始近似值的选择利用数值方法解最小二乘问题，选取一个尽可能好的初始近似值尤为重要。确定初始近似值的方法有很多，即使对于同一个问题，可采用的初值选取方法也不是唯一的。本文利用将非线性最小二乘问题线性化的方法，来确定问题〔1〕的初始近似值，具体步骤如下：1〕作图观察法作图观察法基于人们对数据曲线的观察和目测，通过主观判断来确定各参数值。步骤如下：绘制出被拟合数据，的图像；观察数据曲线的主要特征，根据数据曲线三个峰的主体在横坐标轴上的位置确定前三个参数，和，根据各峰的峰值和峰度来确定剩下的参数，，…，。显然，直接把由作图观察法得到的参数作为初始近似值是不适宜的，必须对其做进一步的处理。2〕线性化方法这里所说的线性化方法，是一种基于对非线性函数或非线性方程〔组〕进行线性Taylor展开的方法。首先定义向量：，我们将残差函数表示为那么代价函数可表示为设为的Jacobi矩阵：，，〔3〕矩阵各元素的具体表达式为：，，那么代价函数的梯度可表示为：〔4〕为使代价函数到达极小值，我们令代价函数的梯度为零，即〔5〕对于非线性方程组〔5〕式，一种典型的线性化方法是Newton法，其迭代公式为：〔6〕其中为在处的梯度；为在处的Hessian矩阵。假设代价函数是正定的二次函数，Newton法一步即可到达最优解，对于非二次函数，Newton法并不能保证经有限次迭代求得最优解，但由于代价函数在极小点附近近似于二次函数，故当初始点靠近极小点时，Newton法的收敛速度一般是很快的，但当初始点远离最优解时，不一定正定，Newton法的收敛性不能保证。但Hessian矩阵中的第二项通常难以计算或者计算量很大，所以利用Newton法求解问题〔1〕是不易于实现的。为此，本文首先将在附近进行线性Taylor展开，得到：〔7〕那么代价函数为：〔8〕其中。再将〔8〕式代入〔5〕式，整理后得到：〔9〕解此线性方程组，并将解记为。从〔9〕式的形式上看，此方法实质上为Gauss-Newton法，它是解非线性最小二乘问题最根本的方法。由此，我们得到了参数的初始近似值。3梯度下降法非线性最小二乘问题的数值解法多种多样，针对拟合函数是叠加函数的情形，本文采用梯度下降法求解问题〔1〕，从而实现对沉积物粒度分布的拟合。梯度下降法以代价函数关于参数梯度的负方向作为极小化算法的搜索方向，是所有需要计算导数的无约束优化算法中最简单的方法。设代价函数在附近连续可微，且梯度。容易证明，负梯度方向〔10〕为代价函数在点处的最速下降方向。梯度下降法的迭代公式为：，〔11〕其中为步长因子；表示单位化后的梯度。对于任意一点，代价函数的梯度是很容易得到的，因而式〔11〕中的也能够比拟容易地得到。梯度下降法的关键在于如何选取适当的步长因子，步长因子既不能太大也不能太小，步长因子太大就有可能导致代价函数增大，使算法失效；而步长因子太小就会导致代价函数减小的速度非常缓慢，降低计算效率。关于步长因子的选取方法可参阅文献[7]和[9]，虽然这些方法在理论上具有很大的优越性，但表达方式比拟复杂，并不适于实际应用。因此，如何选取适当的步长因子，仍然是一个亟待解决的问题。梯度下降法在最优化中具有重要的理论意义，虽然这种方法由来已久，但在实际应用中经常被无视，作为一种经典算法，它的实用价值是值得重视的。对于一些较为复杂的数据拟合问题，梯度下降法作为一种简单的算法，易于实现，具有良好的可操作性和整体收敛性，拟合效果也是令人满意的。4数值实验下面我们应用梯度下降法对甘孜黄土的激光粒度分布进行拟合，同时实现了对3个粒度组分的数学别离。现有一个甘孜黄土激光粒度的实测样本，我们把它作为原始数据，，其中表示第个粒级区间，为粒径值且，为实测值，表示激光粒度在第个粒级区间内的百分含量。对原始数据进行必要的处理后得到被拟合数据，，处理过程如下：对粒径取自然对数：，并对原始数据进行处理：其中为粒级中点的自然对数值，，表示每个粒级区间在对粒径取对数后的宽度。值得一提的是，在本文所拟合的数据中，每个粒级区间在对粒径取对数后的宽度都是相等的，统一记为。对原始数据的拟合残差可由下式：〔12〕得到，其中。数值实验步骤如下：1〕将被拟合数据，绘制在直角坐标系中，其中为横坐标，为纵坐标，应用作图观察法得到各参数的值见表1。在表1中，表示对原始数据进行拟合的拟合标准差。2〕以由作图观察法得到的参数为初值，利用线性化方法，经过4次迭代后得到参数的初始近似值〔表1〕.3〕以为参数的初始近似值，利用梯度下降法对进行拟合。本文针对、、三种参数对步长因子的选取方案如下：，，，迭代10000步后，得到拟合结果见表1和图1。由表1我们可以看出，利用作图观察法估计参数的拟合标准差最大，为；而经过线性化方法得出的初始近似值的拟合标准差显著减小，为，这说明本文使用的初值选取方法是有效的；利用梯度下降法进行数据拟合，拟合标准差进一步减小至，说明我们利用梯度下降法对黄土激光粒度分布进行数据拟合是非常有效的，并且能很好地表达出黄土激光粒度三峰分布的特征〔图1〕。表SEQ表格\*ARABIC1数值实验结果参数估计方法图解法-50.10.43线性化方法-0.32341.70503.57305.43360.81621.94910.00960.6771梯度下降法-0.30931.68443.56853.48400.85731.90910.01840.6495图SEQ图\*ARABIC13个粒度组分及拟合结果1此外，本文应用同样的方法，对另外64个实测样本分别进行了数值实验，均到达了理想的拟合精度〔表2，图2〕。在表2中，绝均差和标准差均是相对于原始数据的拟合误差；在图2中，我们只列出拟合标准差最小的第9号样本的拟合图像〔a〕和拟合标准差最大的第56号样本的拟合图像〔b〕。由表2我们可以看出，对64个实测样本的拟合标准差绝大局部都在以下，最小为，平均为，这再次说明梯度下降法是一种有效的数据拟合方法。表264个样本的数值实验结果样本号绝均差标准差样本号绝均差标准差010.000380660.00054826330.000445680.00058893020.000378050.00052603340.000604520.00073963030.000419950.00063532350.000593980.00074919040.000507760.00064634360.000492280.00062601050.000354670.00050275370.000682870.00086807060.000522020.00069164380.000654660.00081682070.000724330.00090514390.000620560.00077028080.000489590.00066815400.000688480.00082679090.000313310.00043551410.000597180.0007304100.000283450.00044866420.000623190.00078303110.000795410.0010105430.00054030.00068415120.000611020.0008103440.000669610.00085667130.000371370.00051793450.000601210.00074963140.000613540.00080464460.000458670.00059758150.000547660.00073666470.00064320.00082424160.000655750.00088294480.000688950.000924170.000401690.0005404490.000769240.00099171180.000712870.00087265500.000446910.00059015190.000485070.00065388510.000657190.00089951200.000531490.00068144520.000624790.0008126210.000486740.00063836530.000606260.00078845220.000569240.00083757540.000695470.00093944230.000556020.00072327550.000754780.00096269240.000509970.00066889560.000796740.0010187250.000453320.00063725570.000673290.0008643260.000606740.00073622580.000669460.00086196270.000731760.00093603590.000708180.00090343280.000478990.00060638600.000735850.00099672290.000524780.00064484610.000667130.00086381300.000545490.00069155620.000660390.00085916310.000637540.00077518630.000723050.00092306320.000579650.00072029640.000603930.00079226图SEQ图\*ARABIC23个粒径组分及拟合结果25 结论本文将梯度下降法应用于沉积物粒度分布的数据拟合，探讨了利用梯度下降法进行非线性数据拟合的可行性和实施过程，并结合甘孜黄土激光粒度分析数据进行了数值实验。本文第2局部用线性化方法来确定初始近似值的过程，实际上是用Gauss-Newton法进行数据拟合的过程。在数值实验中我们发现，在前4次迭代过程中，每次迭代代价函数都能持续减小，且减小速度较快，而第5步迭代后出现了不稳定现象，我们认为，这是由于拟合函数的非线性程度较大，将其进行线性Taylor展开存在较大的误差所致。因此，在进行非线性数据拟合时，一些在理论上具有很大优越性的数值方法，对于实际应用未必是可行的。梯度下降法为我们进行非线性数据拟合提供了一条新的思路，它有坚实的数学理论根底，是所有需要计算导数的无约束优化算法中最简单的方法。梯度下降法在最优化中具有重要的理论意义，虽然这种方法由来已久，但在实际应用中经常被无视，作为一种经典算法，它的实用价值是值得重视的。数值实验说明：对于一些较为复杂的数据拟合问题，梯度下降法作为一种简单的方法，易于实现，具有良好的可操作性和整体收敛性，拟合效果也是令人满意的。从本文数值实验所得到的结果与文献[1-3]中的结果相比，总体而言，本文数据拟合的误差是最小的。同时此方法也可以推广到解决其他极值问题。参考文献SunDH,BloemendalJ,ReaDK,etal.Grain-sizedistributionfunctionofpolymodalsedimentsinhydraulicandaeolianenvironments,andnumericalpartitioningofthesedimentarycomponents.SedimentaryGeology,2002,152:263–277SunDH,SuRX,BloemendalJ,etal.Grain-sizeandaccumulationraterecordsfromLateCenozoicaeoliansequencesinnorthernChina:ImplicationsforvariationsintheEastAsianwintermonsoonandwesterlyatmosphericcirculation.Palaeogeography,Palaeoclimatology,Palaeoecology,2021,264:39–53SunDH,BloemendalJ,ReaDK,etal.Bimodalgrain-sizedistributionofChineseloess,anditspalaeoclimaticimplications.Catena,2004,55:325–340SunDH,AnZS,SuRX,etal.EoliansedimentaryrecordsfortheevolutionofmonsoonandwesterlycirculationsofnorthernChinainthelast2.6Ma.ScienceinChina(SeriesD),2003,46(10):1049-1059SunDH.MonsoonandwesterlycirculationchangesrecordedinthelateCenozoicaeoliansequencesofNorthernChina.GlobalandPlanetaryChange,2004,41:63–80DearingJA.Sedimentaryindicatorsoflake-levelchangesinthehumidtemperatezone:Acriticalreview.JournalofPaleolimnology,1997,18(1):1～14NocedalJ,WrightSJ.NumericalOptimization.Secondedition.NewYork:Springer,2006FletcherR.PracticalMethodsofOptimization,Vol.1.NewYork:JohnWiley&Sons,1981GillPE,MurrayW,WrightMH.PracticalOptimization.London:AcademicPress,1981.JukićD,ScitovskiR.LeastsquaresfittingGaussiantypecurve.AppliedMathematicsandComputation,2005,167:286–298MarkovićD,JukićD,BenšićM.Nonlinearweightedleastsquaresestimationofathree-parameterWeibulldensitywithanonparametricstart,JournalofComputationalandAppliedMathematics(2021),doi:10.1016/j.cam.2021.09.025JukićD,KralikG,ScitovskiR.LeastsquaresfittingGompertzcurve,JournalofComputationalandAppliedMathematics.2004,169:359–375JukićD.Anecessaryandsufficientcriteriafortheexistenceoftheleastsquaresestimatefora3-parametricexponentialfunction.AppliedMathematicsandComputation,2004,147:1–17JukićD,BenšićM,ScitovskiR.Ontheexistenceofthenonlinearweightedleastsquaresestimateforathree-parameterWeibulldistribution.ComputationalStatisticsandDataAnalysis,2021,52:4502–4511JukićD,ScitovskiR.Existenceofoptimalsolutionforexponentialmodelbyleastsquares,JournalofComputationalandAppliedMathematics.1997,78:317–328JukićD,ScitovskiR,SpäthH.Partiallinearizationofoneclassofthenonlineartotalleastsquaresproblembyusingtheinversemodelfunction,Computing,1999,62:163–178JukićD,SaboK,ScitovskiR.TotalleastsquaresfittingMichaelis–Mentenenzymekineticmodelfunction.JournalofComputationalandAppliedMathematics,2007,201:230–246HölmstromK,PeterssonJ.Areviewoftheparameterestimationproblemoffittingpositiveexponentialsumstoempiricaldata.AppliedMathematicsandComputation,2002,126:31–61HasdorffL.Gradientoptimizationandnonlinearcontrol,NewYork:Wiley,1976AtiegA,WatsonGA.Aclassofmethodsforfittingacurveorsurfacetodatabyminimizingthesumofsquaresoforthogonaldistances.JournalofComputationalandAppliedMathematics,2003,158:277–296BöckmannC.Curvefittingandidentificationofphysicalspectra.JournalofComputationalandAppliedMathematics,1996,70:207–224NyarkoEK,ScitovskiR.Solvingtheparameteridentificationproblemofmathematicalmodelsusinggeneticalgorithms.AppliedMathematicsandComputation,2004,153:651–658AhnSJ,R