黄金分割法在煤热动力学方程中求古地质时间的运用

Application of Gold Section Method in Calculating the Paleogeologic Time of the Coal Thermal Kinetic Equation

摘要: 知道古地温和镜质体反射率，可以利用迭代算法求解煤热动力学方程中的古地质时间。选择所要求的古地质时间为迭代变量，建立一个能够进行迭代的关系式，即将一个有古地质时间数学计算放在方程等式的右边成为自变量，而另一个有古地质时间数学计算放在方程等式的左边成为因变量。用前后所迭代值的差值控制迭代过程以免迭代过程无休止地执行下去。用黄金分割值0.618优化新迭代值与前一个迭代值的关系。用5个分别在石炭纪、侏罗纪、和白垩纪的例子说明计算的简洁和有效。

Abstract: Knowing the vitrinite reflectance and the paleogeologic temperature, the iterative algorithms can be used to solve the paleogeologic time in the coal thermal dynamics equation. The required paleogeologic time is selected as the iterative variable to establish a relationship that can be iterative. The iterative process is controlled with the difference of the iterative values to limit the loop times. The relationship between the new iteration value and the previous iteration value was optimized with the gold section value of 0.618. The simplicity and effectiveness of the calculations is illustrated with 5 examples in the Carboniferous, Jurassic, and Cretaceous, respectively.

1. 引言

煤化作用是泥炭的物理化学结构在温度、压力及时间等因素的作用下经历褐煤、亚烟煤、烟煤最终成为无烟煤的物理化学性质变化，也称为煤变质。二个外在影响因素(温度和时间)对煤化作用的影响和控制是综合性的，也是至关重要的。在地质学的研究中，常用镜质体反射率来表示煤变质程度。吴冲龙和杨起根据地热学与化学动力学的一般原理，参照Bostick曲线和我国中、新生界实测数据，采用双重回归的方法来建立温度–时间–反射率的三变量经验公式 [1] [2] [3] [4] [5]，就是定量表示煤变质热动力学。根据温度–时间–反射率的三变量方程 [6] [7] [8]，有文章讨论自然变质煤化与人工提质煤化热动力学的差异。还有文章计算已知绝对年龄和其镜质体反射率的煤层最低变质温度。也有文章计算已知绝对年龄和变质温度的最大变质极限，即煤层镜质体反射率。还有文章定量说明，提高温度对提高煤变质程度的影响大于延长时间因素。这些通过时间和温度求反射率的例子、通过时间和反射率求温度的例子、以及计算褐煤或长焰煤临界温度的例子虽然证明温度–时间–反射率经验公式能在一定程度上解释煤变质作用，但从数学的角度上，属于有精确解的方程。即所有自变量(已知量)都在方程等式的右边，且右边不包含因变量。而因变量(待解量)只在方程等式的左边，而且只有一个数学运算。因此本文要探讨在已知镜质组平均反射率和古地温的条件下，求古地质时间。

2. 吴冲龙煤热动力学方程

当知道古地质时间和变质程度，求古地质温度时：

$\ln \left(T-273\right)=\frac{646.32}{\ln t+111.85}-\frac{0.492t^{0.093}}{R_0}$ (1)

可以看出，自变量都在方程等式的右边，而因变量则在等式的左边。直接代入自变量，可以求得因变量。

当知道古地质时间和古地质温度，求变质程度时：

$R_0=\frac{0.492t^{0.093}}{\frac{646.32}{\ln t+111.85}-\ln \left(T-273\right)}$ (2)

式中：t：作者先用“岩层绝对年龄”后用“有效受热时间”两种说法，百万年；R₀：镜质体反射率，%，T：古地温，K。

可以看出，自变量都在方程等式的右边，而因变量则在等式的左边。直接代入自变量，可以求得因变量。

但当知道古地质温度T和变质程度R₀，求古地质时间t时，就会出现不管如何移项，除了自变量都在方程等式的右边，而因变量则同时在等式的左边和右边。如：

$\ln t=646.32\left(\frac{R_0}{R_0\ln \left(T-273\right)+0.492t^{0.093}}\right)-111.85$ (3)

方程(3)不存在求根公式，因此求精确根非常困难，从而寻找方程的近似根就显得特别重要。

3. 成煤期和镜质体反射率

从中国六个成煤期的研究报告中，选择三个均能产生从低阶煤(长焰煤)到高阶煤(无烟煤)地质时期作为相互验证的古地质期：① 古生代石炭纪晚期至二叠纪早期：距今约3.20亿年~2.78亿年。时间跨度为四千二百万年，简称石炭纪。石炭纪是最早的煤炭资源形成期，我国的石炭纪煤基本上分布在黄河流域，煤种范围从长焰煤到无烟煤。② 中生代侏罗纪早中期：距今约2.05亿年~1.59亿年。时间跨度为四千六百万年，简称侏罗纪。侏罗纪煤主要集中在内蒙古，陕西，甘肃，宁夏四省区交界地带和新疆北部，煤种范围从褐煤、长焰煤到无烟煤。③ 中生代白垩纪早期：距今约1.42亿年~0.99亿年。时间跨度为四千三百万年，简称白垩纪。白垩纪煤分布于内蒙古东部和东北三省。东北三省的白垩纪煤种范围从长焰煤到无烟煤。因为所要求的古地质时间应该是在确定的地质期里，前后相差一百万年，而决不可能在不同的地质期里。

所有烟煤煤种(气煤、肥煤、焦煤、瘦煤、贫煤、弱粘煤、和不粘煤)的镜质体反射率参考李文华“烟煤镜质组平均最大反射率与煤种之间的关系”一文 [9]。而煤化程度分级(低煤级煤、7级的中煤级煤和3级的高煤级煤)可以参照中国煤炭行业标准 [10]。

4. 求古地质时间的运算

利用迭代算法解决问题。迭代，顾名思义就是不停的代换的意思。主要步骤涉及选择迭代变量，建立迭代关系式，和控制迭代过程。

4.1. 选择迭代变量和建立迭代关系式

所谓迭代关系式，指如何从变量的前一个值推出其下一个值的公式(或关系)。迭代关系式的建立是解决迭代问题的关键，通常可以使用递推或倒推的方法来完成。在可以用迭代算法解决的问题中，至少存在一个可直接或间接地不断由旧值递推出新值的变量，这个变量就是迭代变量。将方程(3)设计成为能够进行迭代的关系式。

$\ln t_2=646.32\left(\frac{R_0}{R_0\ln \left(T-273\right)+0.492t_1^{0.093}}\right)-111.85$ (4)

很显然，所要求的古地质时间t就是迭代变量。

4.2. 黄金分割法

具体操作是将已知的古地温T、镜质体反射率R₀和假设的t₁代入方程(4)，得t₂。一般通过①设定迭代次数或②设定确定前后所迭代值的差值两种方法来控制迭代过程以免迭代过程无休止地执行下去。求得t₂后，按下列方程进行后续计算：

$dt=t_2-t_1$ (5)

$\Delta =0.618\ast dt$ (6)

从方程(6)中可以看到用了黄金分割值0.618。黄金分割 [11] [12] 是指将整体一分为二，较大部分与整体部分的比值等于较小部分与较大部分的比值，其比值约为0.618。这个比例被公认为是最能引起美感的比例，因此被称为黄金分割。

而重新代入方程(4)的新t₁与前一个t₁的关系是：

$新t_1=老t_1+\Delta$ (7)

本文采用确定前后所迭代值的差值方法来控制迭代过程，即

$\left|dt\right|\le 1.0$ (8)

4.3. 计算实例

已知对于来自石炭纪的煤样，古地温T = 373 K、镜质体反射率R₀ = 0.95，迭代结果如表1，图1为10次迭代值的趋势图。

Table 1. Iterative results of Carboniferous coal samples with paleogeothermal temperature T = 373 K and vitrinite reflectance R0 = 0.95

表1. 古地温T = 373 K、镜质体反射率R₀ = 0.95石炭纪煤样的迭代结果

Figure 1. Trend chart of 10 iterations with paleogeothermal temperature T = 373 K and vitrinite reflectance R₀ = 0.95

图1. 古地温T = 373 K、镜质体反射率R₀ = 0.95石炭纪煤样10次迭代的趋势图

通过10次迭代，求得古地质时间为三亿三千一百万年。但这个结果超出石炭纪这个地质期。这说明所给定的古地温偏低。

已知对于来自石炭纪的煤样，古地温T = 433 K、镜质体反射率R₀ = 1.7，迭代结果如表2，图2为迭代的趋势图。

Table 2. Iterative results of Carboniferous coal samples with paleogeothermal temperature T = 433 K and vitrinite reflectance R0 = 1.7

表2. 古地温T = 433 K、镜质体反射率R₀ = 1.7石炭纪煤样的迭代结果

Figure 2. Trend chart of iteration with Carboniferous coal samples with paleogeothermal temperature T = 433 K and vitrinite reflectance R₀ = 1.7

图2. 古地温T = 433 K、镜质体反射率R₀ = 1.7石炭纪煤样迭代的趋势图

通过5次迭代，求得古地质时间为一亿四千一百万年。但这个结果未到石炭纪这个地质期。这说明所给定的古地温偏高。

已知对于来自侏罗纪的煤样，古地温T = 414 K、镜质体反射率R₀ = 1.4，迭代结果如表3，迭代趋势图如图3所示。

Table 3. Iterative results of Jurassic coal samples with paleogeothermal temperature T = 414 K and vitrinite reflectance R0 = 1.4

表3. 古地温T = 414 K、镜质体反射率R₀ = 1.4侏罗纪煤样的迭代结果

Figure 3. Trend chart of iteration with paleogeothermal temperature T = 414 K and vitrinite reflectance R₀ = 1.4

图3. 古地温T = 414 K、镜质体反射率R₀ = 1.4侏罗纪煤样迭代的趋势图

通过6次迭代，求得古地质时间为一亿八千七百万年。但这个结果在侏罗纪这个地质期内。

已知对于来自白垩纪的煤样，古地温T = 440 K、镜质体反射率R₀ = 1.84，迭代结果如表4，迭代趋势图如图4。

Table 4. Iterative results of Cretaceous coal samples with paleogeothermal temperature T = 440 K and vitrinite reflectance R0 = 1.84

表4. 古地温T = 440 K、镜质体反射率R₀ = 1.84白垩纪煤样的迭代结果

Figure 4. Trend chart of iteration with paleogeothermal temperature T = 440 K and vitrinite reflectance R₀ = 1.84

图4. 古地温T = 440 K、镜质体反射率R₀ = 1.84白垩纪煤样迭代的趋势图

通过3次迭代，求得古地质时间为一亿二千九百万年。但这个结果在白垩纪这个地质期内。

已知对于来自白垩纪的煤样，古地温T = 353 K、镜质体反射率R₀ = 0.52，迭代结果如表5，图5为迭代趋势图。

Table 5. Iterative results of Cretaceous coal samples with paleogeothermal temperature T = 353 K and vitrinite reflectance R0 = 0.52

表5. 古地温T = 353 K、镜质体反射率R₀ = 0.52白垩纪煤样的迭代结果

Figure 5. Trend chart of iteration with paleogeothermal temperature T = 353 K and vitrinite reflectance R₀ = 0.52

图5. 古地温T = 353 K、镜质体反射率R₀ = 0.52白垩纪煤样迭代的趋势图

通过5次迭代，求得古地质时间为二千万年。但这个结果未到白垩这个地质期。这说明所给定的古地温偏高。

同时，可以利用计算机编程解决，在知道古地温和镜质体反射率的前提下，进行迭代计算求解煤热动力学方程中的古地质时间。

5. 结语

因为知道古地温和镜质体反射率，直接求解煤热动力学方程中的古地质时间是非常困难的，所以利用迭代算法来解决这一问题。

1) 选择所要求的古地质时间为迭代变量，建立一个能够进行迭代的关系式，即将一个有古地质时间的数学计算放在方程等式的右边成为自变量，而另一个有古地质时间的数学计算放在方程等式的左边成为因变量。用前后迭代值的差值控制迭代过程，以免迭代过程无休止地执行下去。

2) 用黄金分割值0.618优化新迭代值与前一个迭代值的关系。用5个分别在石炭纪、侏罗纪、和白垩纪的例子说明计算的简洁和有效。