反演计算

2024-06-20

反演计算（共9篇）

反演计算篇1

0引言

为提高沉陷计算的准确性,一些学者对开采沉陷计算参数与地质采矿条件(如开采方法、覆岩岩性、结构、矿层倾角、采动程度、冲积层厚度等)之间的定性、定量关系进行了大量研究[1,2,3,4,5,6,7,8]。然而受设备、成本等方面的制约,对开采沉陷计算参数进行准确测量十分困难,使得计算参数值往往与实际值相差很大。因此,利用观测站实测资料,对开采沉陷计算参数进行反演,已越来越受到研究者重视。目前, 开采沉陷计算参数反演方法主要有特征点求参、曲线拟合求参、空间问题求参、正交试验设计法求参、最小二乘法求参[1]、模矢法求参[1]和遗传算法求参[8]等,但这些方法一般存在早熟收敛现象,容易陷入局部最优解。蚁群算法[9]属于随机类搜索算法, 是一种具有正反馈、自组织、较强鲁棒性和本质上并行等特点的算法。本文将蚁群算法引入开采沉陷计算参数反演领域,以矿区开采沉陷计算方法中应用较为广泛的概率积分法[10,11]为例,给出了基于蚁群算法的概率积分法计算参数反演的具体实现,并通过应用实例验证了其可行性和有效性。

1基于蚁群算法的概率积分法求参

1.1问题描述

常用的开采沉陷计算参数位移反分析计算模型为

式中:[K]为与开采沉陷计算参数相关的刚度矩阵; {A}为测点下沉值或水平移动值列阵;{F}为测点等效下沉值或水平移动值列阵;T(P)为目标函数,P为开采沉陷计算参数向量;M为测点总数;Wd(P) 为第d个测点基于P的计算下沉值或水平移动值; W′d为第d个测点的实测下沉值或水平移动值。

概率积分法计算参数有8个:下沉系数q、主要影响角正切tanβ、水平移动系数b、开采影响传播角 θ0、左拐点偏移距SL、右拐点偏移距SR、上拐点偏移距SU、下拐点偏移距SD。设P=(q,tanβ,b,θ0,SL, SR,SU,SD),B为P的搜索空间,则式(1)可描述为在8维空间B中找一个向量P0,使目标函数T(P0) 的值最小。

1.2搜索空间离散

蚁群算法适合于求解组合优化问题,因此,需将式(1)的连续变量优化问题转化为组合优化问题。离散搜索空间B的方法[12]:将P中的第n(n=1, 2,…,8)个计算参数Pn的取值范围分为Mn段,用每一段中点的值代替该段的值。搜索空间B被离散成了含有有限个向量的离散空间。参数Pn的离散程度越高,计算精度也越高。

1.3蚁群算法实现

设m为蚂蚁总数;(i,j)为从元素i到元素j的连线,称为边;τij(t)为t时刻边(i,j)上信息素浓度; t时刻包含边(i,j)的路径共有gij(t)条;路径w对应的目标函数为Tw(P);若路径w包含边(i,j),记成w(i,j);在t时刻,边(i,j)对应的Tw(P)的标准差为σij(t)。m只蚂蚁同时完成1步(从一个元素到达另外一个元素)时,时间步自动加1。m只蚂蚁完成对所有8个元素的访问后,算法完成1个循环。

初始时刻,将m只蚂蚁随机放在离散后的搜索空间B中。蚂蚁k(k=1,2,…,m)在运动过程中, 根据τij(t)与ηij(t)决定其下一步运动方向。其中 ηij(t)为元素i,j的能见度,用于防止算法早熟收敛,取ηij(t)=σij(t)。在t时刻,蚂蚁k由元素i转移到元素j的概率为

式中:Dk为当前蚂蚁k可以选择向其移动的计算参数编号,Dk={N-Tk},其中N为8维向量,元素为8个计算参数编号,Tk为8维动态向量,记录到当前为止蚂蚁k已经访问过的计算参数编号;h为小的常数,避免蚂蚁以零概率选择对应σij(t)较小的边; α,β为控制信息素浓度与转移期望相对重要性的参量。

在蚁群完成1个循环后,对每条边上的信息素浓度进行更新:

式中:ρ为外激素残留系数;f(sbest)为迭代最优解。

参数α,β,ρ的最佳组合可由试验确定[13]。同时,为避免搜索过程中发生停滞状态,参照最大-最小蚁群系统[13],对每条边上的信息素浓度τij(t)做如下限制:

式中:τmax为在每个选择点上,其中一个解元素上的信息素浓度;τmin为所有其他可选择的解元素上的信息素浓度,τmin=0.02τmax。

τmax按式(6)选取:

式中:f(sgb)为渐进的最大值估计。

基于最大-最小蚁群算法的概率积分法计算参数反演流程如图1所示。其中,C为迭代次数,Cmax为最大迭代次数。

2应用实例

处于平顶山矿区的某矿22071工作面走向长1 760m,斜长176m,煤层倾角21°,采高3.6m,采深717~780m,第四系黄土夹卵石层厚度20.86~ 41.23m,平均厚度35m,其下为泥岩、砂质泥岩、细砂岩、煤层等岩层,属中硬岩层。开采时间为2005年12月至2007年12月,推进速度约61m/月,采用全陷法管理顶板。

某矿22071工作面观测站布置如图2所示。观测线总长2 205 m,共布置控制点3个,工作测点73个。其中,走向观测线长1 470 m,工作测点48个;在工作面西侧设半条(下山)倾向观测线,倾向观测线长735 m,工作测点25个。观测时间为2005年11月至2008年8月,共观测了35次,在观测期间,观测线上测点缺失比较严重,走向观测线缺失22个测点,倾向观测线缺失6个测点,但观测数据仍能反映22071工作面地表沉陷状况。

倾向观测线和走向观测线的部分实测下沉曲线分别如图3、图4所示。倾向观测线实测最大下沉值为815mm(Q12测点),走向观测线实测最大下沉值为510mm(Z43测点)。通过分析所得的实测数据,该观测站地表沉陷规律可采用概率积分法计算模型来描述。

将获得的45个下沉值和45个水平移动值作为观测值。概率积分法计算参数取值范围和参数离散数目见表1。蚁群规模m=100,α=1,β=1,ρ=0.8。

利用蚁群算法得到的倾向观测线和走向观测线的下沉拟合曲线分别如图5、图6所示,可看出拟合效果较好。需要指出的是,图5中倾向观测线上山一侧测点拟合值与实测值最大相差142 mm,这是由于反演中未考虑临近采空区的残余沉陷变形影响所致。

为考察蚁群算法的有效性,与最小二乘法、模矢法和遗传算法反演计算结果进行对比,见表2。采用蚁群算法时,倾向观测线拟合中误差为51mm,相对中误差 (相对于实测最大下沉值)为6.1%;走向观测线拟合中误差为23mm,相对中误差为5.3%。从整条观测线的拟合误差上来看,蚁群算法拟合效果优于最小二乘法和模矢法,与遗传算法相当,蚁群算法能有效解决开采沉陷计算参数求参问题。

3结语

提出了基于蚁群算法的概率积分法计算参数反演的求参方法。应用实例表明,采用蚁群算法能较好地进行概率积分法模型计算参数求取,可有效避免算法早熟收敛现象,对观测站测点缺失具有较强的抗干扰能力。本文对基于蚁群算法的开采沉陷计算参数反演仅是一种初步研究,该方法求参的准确性、稳定性、反演结果精度与参数取值范围和离散数目的影响关系,以及蚁群算法与其他算法(如遗传算法)结合用于开采沉陷计算,将是进一步研究的方向。

反演计算篇2

储层反演技术是利用地震和测井信息,在地质理论指导下对储层的空间展布和几何形态进行宏观描述.然而,某些储层反演模型重构方法由于对测井曲线录取时的影响因素、地层埋深对测井曲线的影响、复杂地质条件下不同岩性在测井曲线上响应的差异等问题缺乏深入分析,在一定程度上影响了利用储层反演结果区分地层岩性的能力,同时也制约了储层预测技术的广泛应用.研究以二连断陷盆地为例,在分析自然伽马曲线的各种不同影响因素的基础上,总结出一套自然伽马岩性反演模型的重构方法.认为由于地层受不同沉积时期等因素的影响,造成了自然伽马测井曲线上泥岩基线的差异和砂岩判别标准的`不统一,同时由于钙质泥岩等特殊岩性的存在而影响了利用自然伽马测井曲线对砂岩和泥岩的识别,对这些影响因素做校正处理,可提高储层反演结果的岩性识别能力.

作者：石兰亭郑荣才 SHI Lan-ting ZHENG Rong-cai 作者单位：石兰亭,SHI Lan-ting(成都理工大学“油气藏地质及开发工程”国家重点实验室,四川,成都,610059;中国石油勘探开发研究院西北分院,甘肃,兰州,730020)

郑荣才,ZHENG Rong-cai(成都理工大学“油气藏地质及开发工程”国家重点实验室,四川,成都,610059)

反演计算篇3

岩土工程计算中, 岩体初始地应力和弹性模量是计算中的基本参数, 但是往往难于准确获得。然而现场量测洞周位移却是比较容易的, 因而目前计算中, 以岩体周围位移作为计算参数, 以取代地应力, 或由周边位移反算得地应力和岩体弹性模量。位移反分析法目前已用于岩土工程的监控设计中, 然而应用弹塑性边界元法进行位移反分析目前国内外还不多见。本文是在文献[1]的基础上, 提出以应变空间弹塑性理论为基础的反分析计算的边界元法。文中还提出了求解非线性方程的近似中点刚度法, 它能加速收敛。计算表明, 当荷载的分级较密时, 几乎不用迭代即能获得予期的精度。

1 应变空间中表述的弹塑性本构方程

弹塑性本构关系的统一表达式

(1) 式中α (g) 为加, 卸载因子

式 (1) 可用矩阵形式表示

其中

分别称为弹性矩阵、塑性矩阵和弹塑性矩阵, 矩阵中的元素可由将指标ij和kl分别作为行和列的位置, 按11, 22, 33, 23, 31, 12的顺序取值而得。

2 近似中点刚度法

弹塑性问题的应力增量与应变增量是由系数矩数[Dep]联系的

(2) 式中[Dep]是应变εij的函数。对于第i+1级加荷, 如果i与i+1点之间某点t的εij值能使[Dep]正好等于i与i+1点连线的斜率 (图1) , 则由式 (2) 计算的将是精确值。然而t是未知的, 只能用已知的近似点来代替。显然, 用中点的值来代替, 所得计算结果是比较理想的, 因而有人提出了予测中点刚度法, 不过这种方法计算较为复杂。本文提出了较为简易的用泰勒级数展开式求解中点刚度的近似方法。

[Dep]在i点展成泰勒级数, 取一次项作为近似, 则有

$[D_{e p}]_{t} = [D_{e p}]_{i} + \frac{\partial [D_{e p}]_{i}}{\partial ε} [δ ε] = [D_{e p}]_{i} + [D_{e p}] ´_{i} (3)$

其中

$[D_{e p}]^{'}_{i} = \frac{\partial}{\partial ε} [D_{e p}]_{i} [δ ε], [δ ε] = \frac{1}{2} [\dot{ε}] 。$

将式 (3) 代入 (2) 式得

$[\dot{σ}] = [\dot{σ}]_{i} + [\dot{σ}]^{'}_{i} (4)$

(4) 式中 $[\dot{σ}]_{i} = [D_{e p}]_{i} [\dot{ε}]$ 是以i点的刚度计算的应力增量, 应力修正值 $[\dot{σ}]^{'}_{i} = [D_{e p}]^{'}_{i} [\dot{ε}]$ 写成分量形式则有

${\dot{σ}}^{'}_{i j} = \frac{1}{2} \frac{\partial D_{i j k e}^{e p}}{\partial ε_{r s}} \dot{ε}_{r s} \dot{ε}_{k l} (5)$

(5) 式中

$D_{i j k l}^{e p} = C_{i j k l} - \frac{α (\hat{g})}{A} \frac{\partial Q}{\partial ε_{i j}} \frac{\partial Q}{\partial ε_{k l}} 。$

可见, 只要知道屈服函数Q, 即可求得应力修正值σ′ij。一般情况下, 式 (5) 求导十分复杂。但当Q为Jle及 $\sqrt{J^{'}_{2 e}}$ 的线性函数时, 则因Cijkl, β1, β2以及A都是常数, 求导十分容易, 因此可以推得

$\begin{array}{l} \frac{\partial D_{i j k l}^{e p}}{\partial ε_{r s}} = - \frac{β_{1} β_{2}}{2 A \sqrt{J^{'}_{2 e}}} (δ_{k l} δ_{i r} δ_{j s} + δ_{i j} δ_{r k} δ_{s l} - \frac{2}{3} δ_{i j} δ_{k l} δ_{r s}) + \\ \frac{β_{1} β_{2} e_{r s}^{e}}{4 A (J^{'}_{2 e})^{3 / 2}} (e_{i j}^{e} δ_{k l} + e_{k e}^{e} δ_{i j}) + \frac{β_{2}^{2} β_{i j}^{e} e_{r l}^{e} e_{r s}^{e}}{4 A (J^{'}_{2 e})} - \frac{β_{2}^{2}}{4 A J^{'}_{2 e}} \times \\ [(δ_{i r} δ_{j s} - \frac{1}{3} δ_{i j} δ_{r s}) e_{k l}^{e} + (δ_{r k} δ_{s l} - \frac{1}{3} δ_{r s} δ_{k l}] e_{i j}^{e}) 。 \end{array}$

将上式代入 (5) 式可得

$\begin{array}{l} σ^{'}_{i j} = \frac{x (\hat{g}) G β_{2}}{4 A} [\frac{β_{2} W_{2}}{Ι^{'}_{2}^{2}} τ_{i j} + \frac{β_{1} \dot{w}}{(Ι^{'}_{2}^{})^{3 / 2}} (\dot{w} δ_{i j} + \dot{J}_{1} Ι_{i j}) - \\ \frac{β_{2}}{Ι^{'}_{2}^{}} (\dot{w} \dot{e}_{i j} + 2 {\dot{J}}^{'}_{2} τ_{i j}) - \frac{2 β_{1}}{\sqrt{Ι^{'}_{2}}} (2 {\dot{J}}^{'}_{2} δ_{i j} + \dot{J}_{1} \dot{e}_{i j})] (6) \end{array}$

(6) 式中

$\dot{J}_{1} = \dot{ε}_{i j} {\dot{J}}^{'}_{2} = \frac{1}{2} \dot{e}_{k r} \dot{e}_{k s},$

有了式 (6) , 即可由式 (4) 进行增量计算。

3 弹塑性反演计算的边界元法

用矩阵表示的弹塑性边界元法的基本公式为[1]:

$\begin{array}{l} [Η] [\dot{U}] = [G] [\dot{Τ}] + [D] [\dot{σ}^{p}], \\ [U \dot{Τ}] = [G Τ] [\dot{Τ}] + [Η Τ] [\dot{U}] + [D Τ] [\dot{σ}^{p}] (7) \end{array}$

(7) 式中[ $\dot{U}$ ], [ $\dot{Τ}$ ]为边界上结点的位移、面力增量组成的列阵, [ $\dot{U}$ T]为内部网格结点的位移增量列阵, 而[ $\dot{σ}$ p]为内域各网格的塑性应力 (初应力) 列阵。[H]·[G]·[GT]和[HT]为大家熟知的影响系数矩阵。[D]和[DT]为非线性系数矩阵, 笔者在文献[1]中给出了在常应力单元中的解析表达式。

岩土工程在开挖前, 岩体在地应力场 (σ $_{i j}^{0}$ ) 作用下处于静力平衡状态。由于开挖, 围岩各点将产生位移U′ $_{i}^{}$ , 此位移可看成是由释放应力σ′ $_{i j}^{}$ 引起的。当洞室开挖, 围岩应力重分布达到新的平衡后, 各点的位移和应力为

$\begin{array}{l} U_{i} = U_{i}^{o} + U^{'}_{i} ‚ \\ σ_{i j} = σ_{i j}^{o} + σ^{'}_{i j} 。 \end{array}$

其中U $_{i}^{o}$ 为开挖已有的位移, 计算中不必算入, 一般可取为零。位移反分析计算是根据洞周边界上测得的释放位移U′i推算出地应力σ $_{i j}^{o}$ 和弹性模量E。

锚喷支护洞室中, 边界上的释放分布力可写成

$[Τ] = [F] - [Ρ] (8)$

其中喷层对围岩的分布力[F]由下式求得 (喷层按有限元处理) :

$- [Μ Τ] [F] = [Κ Τ] [U B] (9)$

[MT]为分布力与等效结点力的转换矩阵, [KT]为喷层的刚度矩阵, [UB]为设置喷层后, 喷层与围岩共同作用所产生的洞周位移。

(8) 式中的[P]为对应于地应力σ $_{i j}^{0}$ 的边界表面力;其矩阵形式为

$\begin{array}{l} [Ρ] = [Μ] [σ^{0}], \\ [Μ] = [\begin{matrix} [n (q^{1})] \\ \cdot \\ \cdot \\ [n (q^{Ν})] \end{matrix}] (10) \end{array}$

(10) 式中N为边界点的个数, q为边界结点的位置, n是方向余弦, $[σ^{0}]$ 为地应力列阵。

对二维情况有

$\begin{array}{l} [σ^{o}] = [σ_{11}^{0}, σ_{22}^{0}, σ_{12}^{0}] ‚ \\ [n (q^{i})] = [\begin{matrix} n_{1} (q^{i}) & 0 & n_{2} (q^{i}) \\ 0 & n_{2} (q^{i}) & n_{1} (q^{i}) \end{matrix}] 。 \end{array}$

式中n1 (qi) 为边界外法线方向余弦的第j个分量在第q个边界点处的值。若用[UO]表示加喷层前洞周已产生的位移, 则总位移用增量表示为

$[\dot{U}] = [\dot{U} B] + [\dot{U} Ο] 。$

由式 (10) , 式 (9) , 式 (8) 及式 (7) 联立求解可得

$\begin{array}{l} [\dot{Ρ}] = [X] [\dot{σ}^{p}] + [\dot{Y}], \\ [U \dot{Τ}] = [X Τ] [\dot{σ}^{p}] + [\dot{Y} Τ] (11) \end{array}$

其中

$\begin{array}{l} [X] = [G]^{- 1} [D], [X Τ] = [D Τ] - [G Τ] [G]^{- 1} [D] \\ [\dot{Y}] = - [Μ Τ]^{- 1} [Κ Τ] [\dot{U} B] - [G]^{- 1} [Η] [\dot{U}] \\ [\dot{Y} Τ] = ([Η Τ] + [G Τ] [G]^{- 1} [Η]) [\dot{U}] (12) \end{array}$

注意到[G], [GT], [D]及[DT]每项中都含有1/E (E为弹性模量) , 因而可将式 (11) 系数一致化, 并写为

$\begin{array}{l} [\dot{Ρ}^{*}] = [X] [\dot{δ}^{* p}] + [\dot{y}^{*}], \\ [\dot{U} Τ] = [X Τ^{*}] [\dot{σ}^{* p}] + [Y Τ] | (Δ - 11) (13) \end{array}$

其中

$\begin{array}{l} [\dot{p}^{*}] = \frac{1}{E} [\dot{p}], [\dot{y} *] = \frac{1}{E} [\dot{Y}] ‚ \\ [X Τ^{*}] = E [X Τ], [σ^{* j}] = \frac{1}{E} [\dot{σ} p] 。 \end{array}$

相应地式 (10) 可写为

$[\dot{Ρ}^{*}] = [Μ] [\dot{σ}^{* o}], \begin{matrix} \end{matrix} [\dot{σ}^{* o}] = \frac{1}{E} [\dot{σ}^{o}] (14)$

4 优化计算与增量迭代

由式 (12) 可见, 若知[UO]和[UB], 则[ $\dot{y}$ ]及[ $\dot{y}$ T]便可由每一级荷载的 $[\dot{U} Ο]$ 和 $[\dot{U} B]$ 求得, 由式 (13) 经过迭代运算即可求出[P*]和[UT], 然后由式 (14) 求得标准地应力[σ*0]。但因矩阵[P*]的元素个数一般大于地应力的个数, 所以必须进行优化运算。文献[1]用函数极值原理导出了如下的优化计算式

$[σ^{* o}] = ([Μ]^{Τ} [Μ])^{- 1} ([Μ]^{Τ} [Ρ^{*}]) (15)$

这样, 如果已知E值, 则可由 (14) 式求得地应力[σo]。反之, 如果假定地应力重直分量σ $_{22}^{0}$ 为上部复盖层重量, 即σ $_{22}^{0}$ =γH (γ为地层容重, H为洞室埋深) , 则由式 (14) 可求得弹性模量

$E = γ Η / σ_{22}^{0} 。$

一旦求出E值, 即可由标准地应力求得其他方向的地应力。

迭代中各网格的应变增量可由各结点的位移增量 $[\dot{U} Τ]$ , 根据有限元公式求得

$[\dot{ε}] = [B] [\dot{U} Τ] (16)$

(16) 式中[B]为单元应变矩阵。

塑性应力增量由下式求得

$[\dot{σ}_{i j}^{p}] = \dot{σ}_{i j}^{e} - \dot{σ}_{i}_{j}$

其中 $\dot{σ}$ $_{i j}^{e}$ 和 $\dot{σ}$ ij是同一应变增量 $\dot{ε}$ ij分别按弹性和弹塑性本构方程求得, 具体迭代步骤如下:

(1) 第一次迭代令 $[\dot{σ}^{* p}] = 0$ , 否则令[ $\dot{y}$ ]和[ $\dot{y}$ T]为零。由式 (12) 和式 (15) 求得初地应力增量[ $\dot{σ}$ o]和各节点的位移 $[\dot{U} Τ]$ 。

(2) 按式 (16) 由 $[\dot{U} Τ]$ 求单元应变增量 $\dot{ε}$ ij, 再由广义虎克定律求弹性应力增量 $\dot{σ}$ $_{i j}^{e}$ , 简写为 $\dot{σ}$ e。

(3) 由 $\dot{ε}$ ij求真实应力增量 $\dot{σ}$ ij。设初始屈服时的K (wp) 值为K, 该级荷载施加以前的应力状态为σ, 并做

令f1=g (σ) , f2=g (σ′) , 分下述三种情况进行讨论:

① 若f1<K0, f2<K0, 则 $\dot{σ}_{i j} = \dot{σ}_{i j}^{e}$ ;

② 若f1≥K0, f2≥K0, 则按应力状态由式 (4) 求得 $\dot{σ}$ ij;

③ 若f1<K0, f2≥K0, 则该单元为过渡单元, 需进行如下计算

按应力状态σep由式 (4) 求得应力 $\dot{σ}$ ep, 则有

(4) 计算初应力增量 $\dot{σ}$ $_{i j}^{p}$ , 检查是否收敛, 即检查 $\dot{σ}$ $_{i j}^{p}$ 是否小于规定值;

(5) 计算地应力, 真实应力和各结点的位移值。

(6) 从第二步起继续计算下一个网格, 直至所有网格都算完为止。

(7) 回到第一步作新的迭代。

迭代直到每个网格都有收敛为止, 然后再加荷载, 重复上述7个步骤, 直到荷载加完为止。如果是已知垂直地应力σ $_{22}^{0}$ 的问题, 则得到E值可作下一次的初始值, 然后重新加载, 重复上述步骤直到两次E值的相对误差在允许范围内为止。

5 算例

现对某地下洞室进行了计算, 其几何尺寸如图2所示。岩体内摩擦角φ=31°, 粘聚力C=0.104MPa, 泊松比γ=0.34, 容重γ=0.0027kg·m-2, 洞室埋深为42.5m。

为了验算计算的正确性, 先按下述数据作了正分析计算,

令 σ $_{x x}^{0}$ =0.95, σ $_{y y}^{0}$ =-0.15σ $_{x y}^{0}$ =0, E=980

计算结果如下 (算例中单位为MPa) :

按上表中位移值, 分两种情况再进行反分析计算。

(11) 已知E=980 MPa, 反算地应力场。计算结果列表如下:

(22) 取σ $_{22}^{0}$ =1.15, 反算E及其他地应力

当E的初始值取1 000 MPa的情况下, 计算结果如下

当E初始值为900 MPa的情况下, 计算结果如下表。

由上可见, 反算的精度是较高的, 而且不管E的初始值如何取, 都可收敛于恒定值, 表明计算结果是稳定的。

参考文献

[1]郑颖人, 张德澎, 高效伟.弹塑性问题反演计算的边界元法.中国土木工程学会第三届年会议论文集, 同济大学出版社, 1986

[2]郑颖人, 徐建刚.岩土工程中的三维弹塑性边界元法及程序.第四届全国岩土力学数值分析与解析方法讨论会论文集, 北京测绘出版社, 1991

[3]Shimigee N, Sakurai S.Application of boundary element method for back analysis associoted with tamelling problems.Boundary Element, Hroshima, Japan, 2003

济南市城市热岛遥感反演篇4

摘要：随着城市建设的高速发展,热岛效应越来越被人们所重视.本文采用的Landsat ETM+遥感影像,阐述了地表温度反演的原理和计算方法;用ERDAS IMAGINE软件构建的`模型反演济南地表温度得到地表温度图;分别从较低温区域、中温区、较高温区、特高温区、热核区,研究了济南市热岛效应的空间分布特征及其成因,并提出相应的措施.作者：史同广王丽娟孟飞 SHI Tong-guang WANG Li-juan MENG Fei 作者单位：史同广,孟飞,SHI Tong-guang,MENG Fei(山东建筑大学土木工程学院,山东,济南,250101)

王丽娟,WANG Li-juan(山东师范大学人口・资源与环境学院,山东,济南,250014)

反演计算篇5

关键词：叠前反演,Zoeppritz方程,偏导计算,反演精度,方法验证

叠前地震反演技术能够利用不同偏移距地震变化信息提取地层弹性参数, 对储层岩性、物性以及含流体性研究具有十分重要意义。不断发展的石油勘探对叠前反演储层预测精度要求越来越高, 尽管精度影响因素众多, 但叠前反演方法是影响叠前反演技术精度的重要原因之一。

叠前反演技术的理论基础平面波非垂直入射理论, 在同一界面中地震波从不同角度入射会产生地震波的反射和折射。并且反射系数和折射系数随着入射角度变化而变化。Zoeppritz方程给出了其解析数学表达式, 精确描述了反射系数与入射角度、速度以及密度间的关系。叠前反演的核心就是计算反射系数与弹性参数的梯度变化关系;但由于Zoeppritz方程解析表达式非常复杂, 在叠前反演技术中直接构建介质弹性参数与反射系数的梯度变化关系难度大, 地球物理学者基于弱反射界面以及小角度入射的近似假设推导出了该方程的简化表达式[1—3], 目前大多数叠前反演方法通常采用简化表达式推导构建反射系数梯度方程;但由于简化式的近似假设会存产生两方面的误差:一方面上下地层弹性参数变化大的情况下, 简化式反演会产生一定的计算误差;另一方面在大入射角度地震反演计算时, 误差也会比较大。随着油气勘探深入, 地质目标越来越复杂, 对勘探精细程度要求越来越高, 需要提高叠前地层弹性参数反演的精度;同时由于地震采集技术的发展, 大偏移距宽入射角度的地震资料也越来越多, 基于简化式的叠前反演方法在这些地质目标和地震资料应用中会产生较大的误差, 反演准确性差和精度低。

针对目前叠前反演方法存在的问题, 直接由Zoeppritz方程计算反射系数与弹性参数的关系, 推导构建了基于Zoeppritz方程的高法精度和高效率的叠前高精度反演新方法。该方法没有近似条件限制, 对复杂地质条件和大入射角度地震资料都有较好适应性, 可为储层预测研究提供更可靠的反演成果。利用二维模地震型数据和实际地震资料对研究的高精度叠前反演方法进行了验证, 并与目前常用的基于简化式叠前反演方法进行了对比, 分析认为新方法对地震资料适应性更强, 储层表征和分辨能力更好。

1 叠前反演目标函数构建

根据广义线性反演的思路需要建立以下叠前反演目标函数

式中Vp表示纵波速度, Vs表示横波速度, ρ为密度, D为实际角道集记录, S为期望地震模型响应。R (Vp, Vs, ρ) 为纵波反射系数, W为地震子波。

将R (Vp, Vs, ρ) iΔ在初始模型响应处R (Vp, Vs, ρ) i泰勒展开, 为使运算简单通常将泰勒展开式中的二阶和二阶以上的省略, 得到

对两端求取ΔVp、ΔVs和Δρ的一阶偏导并令其为零从而达到目标函数最小, 令Δdi=D (Vp, Vs, ρ) iS (Vp, Vs, ρ) i, 则推出

通过这样一个方程来求解三个弹性参数的摄动量ΔVp、ΔVs和Δρ, 该方程涉及三个未知量, 为欠定方程求解时必然有很强的多解性。为减少其反演产生的多解性, 引入多个角度的地震资料来联合反演这三个参数。

若有三个角度的角道集数据, 分别为小角度D1、中角度D2和大角度D3, 以G (vp) 代替, 以G (vp) 代替, 以G (ρ) 代替, 则式 (3) 可化为:

将以上三式写为矩阵形式有

根据式 (7) 采用最小二乘法可得摄动量的表达式为

由于矩阵中元素G (vp) 为为, 而, 所以通过求解便可求的G (vp) , 即可以得到矩阵元素。求解方程 (8) 便可得到方程中的参数即反演结果, 而求解方程的核心是地震反射系数对介质弹性参数偏导关系的构建。

2 基于Zoeppritz方程的精确偏导计算

Zoeppritz方程给出了反射系数和入射角度、纵横波速度及密度等的复杂关系, 是叠前反演理论基础, 如式 (9) 。

式 (9) 中:Rp、Rs、Tp和Ts分别是纵波反射系数、转换横波反射系数、透射纵波、透射横波反射系数, α为纵波反射角, β转换横波反射角, α'透射纵波折射角, β'透射横波折射角, vp1、vs1、ρ1、vp2、vs2、ρ2分别为波在上下介质纵横波速度以及密度。

可令

则Zoeppritz方程可简写化为

式中:

对Zoeppritz方程进行简化得出的纵波反射系数与纵横波速度及密度的近似关系式, 相比复杂解析式描述准确性差, 基于简化式偏导方程的叠前反演运算精度低, 因此开展了基于Zoeppritz方程直接进行偏导方程计算构建研究。

利用利用斯奈尔定律和三角关系式用入射角度替代折射角和反射角, 减少Zoeppritz方程中的未知量个数。再在Zoeppritz方程中分别对上下地层纵横波速度以及密度6个变量求偏导得到6个偏导数方程组,

对所有6个变量可以统一用下式来表示:

式 (10) 中:R为反射系数, m=[vp1, vp2, vs1, vs2, ρ1, ρ2]。

通过计算A的逆以及A和B对弹性参数的偏导数就可以得到纵横波反射系数以及纵横波透射系数对上下地层纵横波速度及密度的偏导数, 式 (11) 为A对上介质纵波速度偏导。

对偏导为零的值进行了单独的分离, 提高计算精度的同时也提高了计算效率。这样就导出Jacobi偏导矩阵方程, 建立了各种反射系数对纵、横波速度、介质密度关系, 其中纵波反射系数对各个弹性参数的偏导是叠前高精度反演主要应用方程, 实现了基于Zoeppritz方程反射系数梯度矩阵的精确计算[4—6]。该偏导矩阵直接从Zoeppritz方程精确计算所得, 相比用其近似式构建的偏导方程更准确, 能够更高精度实现地层弹性参数的反演迭代运算。

3 二维模型验证

模型正演具有弹性参数明确、正演地震噪音少、反演子波确定等优势, 利用模型验证能够排除以上非方法因素对反演结果的影响, 因此首先基于二维正演数据开展高精度反演方法验证。设计二维理论模型弹性参数如图1所示, 弹性参数包含纵波速度、横波速度, 考虑到实际的地质弹性参数变化情况, 地层速度向深部逐渐增大。二维理论模型包含断层、楔状体、透镜体、薄互层等常见地质构造, 这些特殊构造的大小规模不等, 厚度不一, 模拟了地下复杂的地下情况。

考虑实际地震采集情况, 二维模型设计最大偏移距6 000 m, 最大深度2 300 m左右, 炮间距10 m, 估算最大入射角度达50°左右, 正演地震大偏移距信息多, 入射角度宽。采用声学波动方程正演得到了三个角度范围 (0°~15°, 16°~25°, 26°~30°) 地震剖面, 图2, 地震正演子波主频40 Hz左右, 频宽5~90 Hz, 与实际地震资料频谱相当。

基于正演分角度叠加地震分别开展了近似式和高精度反演处理, 并对两种反演结果开展了对比。图3为两种反演方法纵波速度反演结果对比剖面, 高精度反演对于层间差异表现更好, 无论是对薄层还是厚层, 高精度反演结果精度更高, 地层参数差异更加明显, 与模型吻合程度更好, 同时高精度反演结果对于断层面具有很好的表现。

4 实际资料验证

在模型验证基础上又开展了实际资料的验证, 选择目标为胜利油田KD地区浅层河流相储层, 目的层埋藏较浅 (1 100~1 500) m, 地震资料品质相对较好, 频带宽, 分辨率能力高, 能够避免一些资料因素带来的影响。地质目标单砂体薄 (2~12) m, 横向沉积变化快, 砂体纵向上相互叠置, 砂体关系复杂, 同时由于埋藏浅大偏移距信息相对丰富, 因此储层预测对叠前反演方法提出了更高的要求, 以往开展了很多叠前反演方法的应用未能满足储层预测高精度要求。

根据叠前反演对地震资料需要, 首先开展叠前精细速度分析和叠前偏移处理, 同时开展了能量均衡以及去噪处理, 保证了地震资料具有高信躁比和高保真性, 最后对分角度叠加地震资料分别开展了叠前高精度方法以及简化式方法反演处理, 可以得到纵横波速度、密度, 进一步可以得到提取储层敏感参数纵横波速度比、泊松比、拉梅系数乘密度, 利用地层敏感弹性参数对储层岩性流体识别精度更高。

图5为过KD地区KD104井的纵横波速度比反演剖面, 高精度反演结果可以清楚的表现三个河道砂体的展布, 砂体的形态和叠置关系清晰, KD104井钻遇了左边和中间的两个砂体, 反演结果与测井吻合较好。而简化式反演结果对砂体形态和叠置关系反映不清, 不能很好表现地震中储层的变化特征。实际资料表明高精度反演比简化式反演结果对储层的纵横向分辨能力更高, 特别是对于薄储层砂体表征能力更强。

5 结论

反射系数偏导方程是进行反演计算的核心, 目前叠前方法中大都应用Zoeppritz方程简化式计算的偏导方程, 由于简化式基于近似假设不能准确代表Zoeppritz方程所描述的反射系数与弹性参数间的关系, 因此造成反演精度低, 适应性差。直接利用Zoeppritz方程精确计算了反射系数偏导方程, 并对其开展了优化, 形成了基于精确反射系数偏导方程的叠前反演方法, 利用模型和实际资料对方法进行了测试, 相比基于简化式的叠前反演方法, 基于精确偏导方程的高精度叠前反演成果精度更高, 对大角度资料以及强变化界面适应性更强。

参考文献

[1] Bortfeld R.Approximations to the reflection andtransmission coefficients of plane longitudinal and transverse wave.Geophys Prosp, 1961; (4) :485—502

[2] Aki K I, Richards P G.Quantitative seismology.W.H.Freeman and Co, 1980:153—154

[3] Shuey, R T.A simplification of the Zoeppritz equations.Geophysics, 1985;50 (5) :609—614

[4] 孟宪军, 刘福平.地震反射系数的相角变化.石油地球物理勘探, 2010;45 (Z) :121—124Men X J, Liu F P.The phase angle change for seismic reflection coefficient.Oil Geophysical Prospecting, 2010;45 (Z) :121—124

[5] 刘福平, 孟宪军.反演纵横波速度的Jacobian矩阵及精确计算方法.中国科学:地球科学, 2010;40 (11) :1608—1616Liu F P, Men X J.Jacobian matrix for the inversion of P-and S-wave velocities and its accurate computation method.Sci China Earhth Sci, 2010;40 (11) :1608—1616

反演计算篇6

关键词：测井资料,低孔低渗储层,含水饱和度,胶结指数,测井评价,南海西部海域

随着南海西部海域大的构造圈闭、岩性圈闭勘探开发的不断深入，勘探开发程度逐渐加大，勘探开发的目标区已经逐渐转移至以低渗透岩性因素起主导作用的砂岩油气藏[1,2,3]。由于特殊的地质特征和复杂的岩石物理特性，油气层与含水层在测井响应特征上区分不明显;并且目前对该类储层的解释方法研究没有系统化，岩石物理解释参数很难确定，从而导致在低孔低渗储层饱和度的精确评价方面存在很大的困难。为了提高此类储层饱和度评价的精度，需要探索研究适应低孔渗储层特点的饱和度计算方法，而评价的关键就是如何准确获取储层的岩电参数。目前对于低孔渗储层岩电参数变化规律的研究，国内外已有测井分析家从岩石物理实验的温压条件、储集空间类型、孔隙结构等方面开展了大量的基础研究工作[4,5,6,7,8,9,10,11,12,13,14]。目前利用Archie公式评价储集层含油气饱和度时，根据岩电实验通常在某一层位取胶结指数m为固定经验值，但由于实际的储层条件与Archie公式假设的地层有一定的差别，特别是在泥质含量高和低孔低渗型储层，如m仍采用固定值，势必会影响饱和度的计算精度，从而影响储层评价的结果。由此，笔者在南海西部海域大量低孔渗储层岩电实验分析研究工作基础上，探索性的建立了利用测井资料优化迭代反演出储层岩电参数m值与含水饱和度的方法。

1 基于测井资料反演m值与饱和度的方法

1.1 方法原理

Archie公式以岩电实验为基础建立了地层因素(F)、孔隙度(Φ)、电阻增大率(RI)、含水饱和度(Sw)之间的关系[15]，由此可得:

对式(1)，令a=1可得到地层胶结指数m值，如式(4)。

参考Archie公式，可得:

同理令a=1，可得:

式中，F为地层因素;m为胶结指数;n为饱和度指数;Φ为地层孔隙度，小数;Fc为拟地层因素;mc为拟胶结指数;Rt为地层电阻率，Ω·m;Rw为地层水电阻率，Ω·m;Φw为含水孔隙度，小数。

1.2 反演mc值的岩电实验基础

为了对低孔低渗储层进行精细评价，近几年有针对性的做了大量地层温压条件下的岩石物理实验。研究数据来源于南海西部海域某气田，储层岩性主要以粉砂岩、泥质粉砂岩为主，选取其中的12口井50块岩心的岩电实验数据，岩心孔隙度分布范围7.9%～21.2%，渗透率分布范围0.023 md～6.3md，属典型的中、低孔低渗储层。

首先参照Archie公式，将岩心的地层因素F与孔隙度Φ在双对数坐标系下作出两者的关系图，如图1;可以看出在相对高孔(18.5%～21.2%)区间，F～Φ具有很好的相关关系;而在低孔(7%～18.5%)区间，F～Φ的相关关系较差，数据点在趋势线两边的分布范围较广，由此可见，这类储层的胶结指数m并不是一个固定值，且变化规律较复杂。

在模拟地层水矿化度条件下，改变岩样的含水饱和度，并通过岩电实验可得到相应电阻率值(Rt)。因此可利用式(1)和式(6)对岩电实验数据进行处理得到每块岩心的胶结指数m值和岩心在不同饱和度下的拟胶结指数mc值。将所有岩心的胶结指数m值和拟胶结指数mc值进行对比分析，结果见图2。图中横坐标代表岩心胶结指数m值，纵坐标为岩心在不同饱和度条件下的拟胶结指数mc值，对于同一块岩心，m值只有一个，而mc值是多个;为了探讨m值与mc值之间的关系，由图2可见:同一块样品，m值与mc值很接近，二者之间的绝对误差小于0.04，即可以认为两者在误差范围内是相等的。这是利用测井资料反演地层胶结指数mc值与含水饱和度的实验基础和前提。

1.3 反演mc值与饱和度的过程

反演拟胶结指数mc值和含水饱和度采用迭代逼近法，资料处理过程有以下四个步骤，详见图3:

(1)首先利用伽玛、中子和密度测井曲线计算得到储层的泥质含量及有效孔隙度参数;

(2)含水饱和度以0.5%为步长(可调整)，从0到100%逐一取值;

(3)储层的初始含水饱和度为Sw，将该值与储层有效孔隙度相乘得到含水孔隙度Φw，然后结合电阻率曲线利用式(5)和式(6)可得到对应于不同含水饱和度的拟胶结指数mc值;

(4)由图2可知，拟胶结指数mc值与岩石胶结指数m值是相等的，因此可把mc值代替胶结指数m值代入式(3)，其中:a=1,b和n值参考气田的区域经验值。由此可以计算出Sw1，此时需要对比分析Sw与Sw1的绝对误差值，如果计算得到的Sw1与假设的Sw值绝对误差大于给定的截止值，则将Sw加上含水饱和度步长值，再代入相应的公式进行循环计算，直到计算的Sw1值与假设的Sw值的绝对误差小于给定的截止值，此时反演计算得到的mc值即为地层胶结指数m值，程序停止迭代循环，输出计算结果。

2 测井解释应用

选取南海西部海域低孔低渗储层测井资料进行处理，图4和图5分别为在A和B井利用本文研究得到的方法进行处理解释得到的成果图。

图4为A井测井解释成果图，图中第5道中蓝色实线为反演得到连续的地层胶结指数mc值，红色点为岩心在地层温压条件下岩电实验得到的m值，反演结果显示在3 068～3 070 m储层反演得到的m值在1.45左右，3 128～3 143 m储层反演得到的m值在1.65左右，两套储层的m值变化特征明显不同，综合分析各测井曲线以及岩屑录井资料可知3 068～3 070 m层段泥质含量高，岩性为泥质粉砂岩;而3 128～3 143 m储层则为细砂岩，储层岩性较纯，泥质含量低;利用岩心资料对测井计算的结果进行标定，二者吻合很好，由此可认为反演结果有效反映储层特征的变化情况。

图5为B井测井解释成果图，图中第5道中蓝色线即为利用上述方法反演得到的含水饱和度曲线，红色实心点为岩心分析得到的束缚水饱和度，可见二者吻合较好;第6道蓝色实线为反演得到连续的地层胶结指数m值，红色点为岩心在地层温压条件下岩电实验得到的m值，可见两者亦吻合很好，计算的结果符合地层实际情况。由此可见，利用本文所用的方法反演得到的m值与含水饱和度是合理的。

3 结论

(1)利用本文的方法可以基于测井资料对低孔低渗储层m值与含水饱和度进行反演计算，提高了低孔渗储层饱和度的解释精度，所用方法在南海西部海域低孔渗储层含水饱和度评价中具有一定的适应性。

(2)对低孔渗地层饱和度进行反演的首要前提是必须先对此类储层的岩电资料进行处理，如果得到的拟胶结指数mc值与岩心胶结指数m值的绝对误差小于0.04时，则这种条件下可以利用本文的方法计算地层含水饱和度。

基于不同算法的温度反演比较篇7

1. 影像准备 :1通过美国地质勘探局网站下载研究区域灞河流域的遥感影像图,下载遥感图为Landsat8影像。2经定位分析得知,屯溪流域位于灞河流域位于陕西省蓝田县灞源乡华岔村西部。查阅该地区TM遥感影像行列号,确定屯溪流域Landsat 8遥感影像所在的行列带号。3按影像行列带号下载影像及头文件。

2. 模板准备 : 准备灞河流域矢量模板,用于对影像进行裁剪获得灞河流域影像图,矢量模板文件为 :灞河裁剪模板 .evf

2 温度反演

Landsa8影像的10波段是热红外波段,接受的是与地表温度高低相对应的地物的热红外辐射,因此,可以通过其来反演出地表温度 , 在温度反演之前需要进行图像拼接与裁剪。

2.1图像拼接和裁剪 :三张原始影像的第10波段影像,进行拼接得到拼接后影像后,打开灞河EVF文件,进行掩膜处理,最后把掩膜的文件放在拼接图像上进行裁剪,得到裁剪结果。

2.2各类算法温度反演

2.2.1 B算法

这种方法简单,易于操作,对影像之外的外来参数依赖性很小。反演效果也较好,是种简单实用的方法。

2.2.2单通道算法

2.2.3单窗算法

2.2.4辐射传输方程法

大气平均作用温度估计误差对地表温度演算精度有较大影响 , 但这种影响直接取决于大气透射率的大小。大气透射率估计误差对数据的地表温度演算精度影响很大,这一误差随地表温度增加而迅速增大 , 温度越高 , 误差越大。

2. 地表比辐射率计算 :水体像元的比辐射率赋值为0.995,自然表面和城镇像元的比辐射率估算则分别根据下式进行计算 :

b1为NDVI值 ;b2为植被覆盖度值,可以通过计算得到地表比辐射率数据。

在NASA官网 (http://atmcorr.gsfc.nasa.gov/) 中输入成影时间以及中心经纬度,可以查询上面公式中所需的参数。利用ENVI主菜单 ->Basic Tools->Band Math,在公式输入栏中输入 :

4. 反演地表温度 :在获取温度为TS的黑体在热红外波段的辐射亮度后,根据普朗克公式的反函数,求得地表真实温度TS :

b1 :公式中的温度为T的黑体在热红外波段的辐射亮度值。由此可以得到真实的地表温度值(摄氏度)。

3 四种反演方法比较

见图1

B算法反演过程较为简便,易操作 ;对影像数据之外的外来参数依赖性较小 ;考虑了地表比辐射率因素的影响 ;单通道算法分类讨论地表比辐射率的影响,对结果反演较好 ;单窗算法既考虑了地表比辐射率的影响,也考虑了大气辐射的影响 ;反演过程所需要的大气参数比传统的辐射传导方程法要少得多,仅需要近地表气温和大气水分含量这两个参数 ; 辐射传输方程法也较为简便,易操作 ;分类讨论了地表比辐射率的影响。

4 论述

单窗算法反演和大气校正法的地表温度最高有54度,误差很大,与其他两者相比,都高出8度左右 ;

1) 总体比较 ;B算法和辐射传输方程法反演出来的温度差异在3度左右 ;两者与单窗算法反演出来的温度相差5-7度左右 ;总体反演温度趋势呈 :单通道算法

2) 各类别比较

见图2

3) 综述 ; 综上可看出,单通道算法反演的温度普遍偏低,但与其他算法最大不超过1度,单窗算法反演温度普遍偏高,比其他算法高出8度左右 ;B算法居中,辐射传输方程法在高温地区反演效果较差,与其他两者相差较多。

东沟尾矿坝渗透系数反演分析篇8

陕西太白黄金矿业有限责任公司东沟尾矿库地处秦岭中部地区,海拔在1240~1945m之间,相对高差在200~800m之间,属中高山侵蚀地貌。矿区地形总体为南西倾向的单面山,北高南低,属太白河左岸的一级支沟,沟长5.2km,沟底宽3~10m,西侧沟坡35°~50°,属“∨”型谷,局部地段沟有高20~40m的第四系冲洪积平台,而且地形切割剧烈。山势陡峻,植被茂密,属秦岭腹地不明显流失造林种草区。尾矿库坝址选在距东沟沟口1400m处的狭窄地段,由两期坝组成。初期坝为碎石土坝,用碎石土逐层碾压筑坝,坝高28m,坝顶宽3.5m,坝上游坡为:1∶1.7,下游坡为:1∶2.2,下游坡设2条宽2m的马道;后期坝为尾矿堆积坝,共设6个平台,平台宽2m,最大堆积高度69m,堆积外坡为1∶5.0,内坡为自然沉积坡1∶100,平均年堆积高度2.76m,最终总坝高97 m,总库容为1279.1×104m3,尾矿库汇水面积:10.125km2,尾矿库类别为:三等库;服务年限25年。

尾矿库的防洪标准:初期洪水校核为50~100年一遇;中、后期为200年一遇校核。

库区下游无人居住,坝底下游600m处有建筑物一处。

尾矿库分为初期坝(28m)和后期堆积坝(69m)形成整个的尾矿库,库内设有5个排水井、一支排洪涵洞、回水泵站一座,坝面设置有平台、排水沟、截洪沟、简易公路(1500m)、马道,在坝体上游设置拦洪坝。

尾矿库采用排水井—排水管—排水涵洞的排水方式,共设有5座(直径:5.5m)框架式排水井,前三座高21m,后两座高24m。其中第5座为2003年10月增补设计。隧洞断面为:高3.5m×宽3.0m的平底圆拱直墙涵洞,最大泻洪量达103.7m3/h,全长1720m ;坝体上游修筑拦洪坝及坝体修筑截洪沟共350m和在坝体上游较大支沟设置长450mФ500mm的截洪引水钢质涵管,同时将地表水引流到坝底而不进入库内。库内溢流水可以通过排水井—排水管—排水涵洞到达回水泵站输送回工艺循环再利用。

为延长该尾矿库使用年限,分析尾矿坝现有运行状况,然后决定是否对该尾矿坝进行筑坝加高扩容改造设计。

尾矿坝的渗流和稳定分析[1,2,3,4,5,6]一直是尾矿坝设计和研究中的最主要问题之一,而尾矿坝加高后的渗流和稳定分析就显得更为重要。

本文基于东沟尾矿坝的实测钻孔水头数据,采用“反问题正算法”,通过三维渗流有限元法以最小二乘法拟合出和实测钻孔水位误差最小的最佳渗透系数[7,8]。该渗透系数综合考虑了尾矿坝各方面影响的是综合渗流参数,用以为尾矿坝运行及加高是否可行方案提供基本依据。

1 东沟尾矿坝渗流参数模型的建立

图1为尾矿坝中包含了所有钻孔实测水位的工程局部图。渗流参数反演分析计算模型建立以工程提供的断面ZK1~ZK8位置和水头为基本依据,钻孔水位资料见表3。计算模型区域为ZK1~ZK8之间实测渗流自由面(浸润面)之下的坝体(含尾矿坝体和初期坝堆石体)和强风化岩体,反演分析计算模型及三维有限元网格分别见图2、图3。

在图2、图3所示计算模型中,整体坐标系的坐标原点取在钻孔ZK1-2处,东西方向为y 轴方向,x轴与y轴正交,z坐标仍为高程。反演分析计算模型共被剖分为51362个结点,八面体尾矿砂单元36000个,风化岩体单元4200个,初期坝单元8000个。

2 渗透系数反演数值分析

考虑到模型渗流参数(渗透系数)和实验室试验渗流参数、现场试验渗流参数的不同,此尾矿坝渗流参数反演分析拟合的模型渗流参数(渗透系数)的取值范围和试算初值如表1所示。

经反演分析和拟合计算,东沟尾矿坝最佳拟合的模型渗流参数见表2,各钻孔水位拟合误差见表3,表3中相对误差为绝对误差和总水头(钻孔1和钻孔8水头差)的比值。最佳拟合模型渗流场分布(各钻孔连线纵剖面)见图5~7。

3 结论

从以上反演分析可以看出:

(1)最佳拟合模型尾矿砂(拟合渗透系数2.0×10-6m/s)渗透系数远大于室内测试平均值9.26×10-8m/s(0.008m/d)。最佳拟合模型渗透系数中已等效反映了排渗措施的作用。

(2)每个剖面中间的钻孔反演误差比两边钻孔的大,这与实测浸润面两边高、中间较低吻合。各钻孔的最佳拟合水位普遍高于实测值,且距离钻孔ZK1和ZK8越远,其误差越大,最大相对误差为7.45%,这说明了虽然反演分析模型可以等效反映排渗措施的作用,但拟合的浸润线位置仍然总体偏高。从另一方面讲,采用等效反映排渗措施作用的最佳拟合模型渗透系数进行计算是偏于安全的。

(3)此尾矿坝加高扩容改造后堆积高程达1359m,渗流分析时尾矿砂、风化岩体和初期坝堆石的模型渗透系数可分别取为2×10-6m/s、4×10-6m/s 和1×10-4m/s。

(4)模型渗透系数中已等效反映了排渗设施的影响,当尾矿坝加高扩容改造堆积高程达1259m时,相应总水头将增加20m,加高方案中应加强排渗设施。

(5)由钻孔实测水位资料可以看出,浸润线位置较高,距坝面较近;加高扩容改造后堆积高程达1359m时浸润线位置将进一步抬高。因此加高方案设计中也应加强排渗设施。

参考文献

[1]赵坚,纪伟,刘志敏.尾矿坝地质剖面概化及其对渗流场计算的影响[J].金属矿山,2003,(12):24~2.

[2]胡明鉴,陈守义,郭爱国,张华.某上游法尾矿坝抗滑稳定性浅析[J].岩土力学,2003,24(s2):254~258.

[3]赵坚,沈振中.尾矿坝复杂排水系统渗流计算方法的改进[J].河海大学学报,1997,25(2):110~113.

[4]蒋卫东,李夕兵.尾矿坝浸润线时空混沌模型及反分析[J].中南工业大学学报自然科学版,2003,34(6):704.

[5]李守义,陈尧隆,胡再强,柴军瑞.米箭沟尾矿坝渗流与稳定分析研究报告[R].西安:西安理工大学,2004.

[6]段蔚平,汪斌.尾矿坝非饱和带滞水曲线模型的建立及应用[J].岩土力学,2003,24(s2):65~68.

[7]柳厚祥,宋军,陈克军.尾矿坝二维固结稳定渗流分析[J].矿冶工程,2002,22(4):8~11.

半干旱地区大气水汽含量反演分析篇9

关键词：半干旱地区,大气水汽含量,反演,降水转化率

水汽是大气中最活跃的成分, 不仅在全球变暖的过程中充当了温室气体的重要作用, 还是天气和气候变化的主要驱动力, 是预测全球气候变化及降雨、中小尺度灾害性天气的一个重要参量。

目前, 大气水汽监测的手段主要有无线电探空仪、微波辐射计、红外卫星遥感、GPS遥感等。除GPS遥感以外的探测手段虽然能获得较高的水汽精度, 但存在运行费用较高、设备原件对环境的要求严格、时空分辨率低等缺陷, GPS虽然弥补了这些观测手段的不足, 但其定位精度有待于进一步提高[1]。使用太阳光度计反演大气水汽含量是地基观测中一种新的方法, 其太阳跟踪精度和时间分辨率较高, 其936 nm通道测得的太阳辐射数据可以反演大气水汽含量。Thome等[2]总结了许多科学家利用水汽吸收率与水汽量的关系反演水汽量的工作。Reagan等[3]采用改进的兰勒法反演了大气水汽含量, 为利用太阳光度计936 nm通道的观测数据反演大气柱水汽总量提供了基本理论和方法。近年来, 国内利用太阳光度计反演大气柱水汽总量的研究有不少成果[4,5,6]。但这些研究主要集中在中东部地区, 而对半干旱地区的研究较少, 同时半干旱地区是气候变化的敏感区, 随着全球气候加剧变暖, 缺水日益严重, 合理开发空中水资源对半干旱地区缺水有一定缓解作用。大气水汽含量是产生大气降水的必要条件, 要合理开发空中水资源, 必须对大气水汽含量进行一定研究。

本文采用兰州大学半干旱气候与环境观测站太阳光度计晴空观测资料, 采用改进的兰勒法反演出大气水汽含量, 开展对半干旱地区大气水汽含量研究, 对于半干旱地区农业生产区及生态区人工增雨 (雪) 作业具有科学的指导作用。

1 资料来源与方法

1.1 资料来源

兰州大学半干旱气候与环境观测站 (SACOL站) 地理坐标为北纬35.946°, 东经104.137°, 属于典型的温带半干旱气候。因此, 利用该站点太阳光度计反演的大气水汽总量能代表方圆数百千米半干旱地区的大气水汽状况。资料为SACOL观测站2006年7—12月和2007年1—4月的CE-318型太阳光度计晴空日观测资料及同期降水资料。

1.2 改进的兰勒法

由于地面测得的太阳辐射信号在936 nm附近水汽吸收带不符合比尔定律, 依据Bruegge等和Halthore等的研究, 水汽透过率用2个参数表达式来模拟:

式 (1) 中, Tw是通道上的水汽平均透过率, w是大气路径水汽总量, a和b是常数, 分别为0.585和0.573。

在936 nm水汽吸收通道, 太阳光度计通过大气到达地面的太阳直射辐射度的响应可表示为:

式 (2) 中, V为太阳光度计的电压输出, V0为大气外界的电压输出, R为测量时刻的日地距离校正量, m为大气质量数, Tw为水汽平均透过率, τ是Rayleigh散射和气溶胶散射光学厚度之和。将 (1) 式代入 (2) 式同时两边取对数:

在稳定和无云的大气条件下, 以mb为X轴, 以ln V+mτ为Y轴画直线, Y截距为ln (V0R-2) , 斜率为-a Pwb, 从而求出大气水汽含量。

2 结果与分析

2.1 半干旱地区大气水汽总量月变化

图1为半干旱地区2006年7月至2007年4月大气水汽含量月变化情况。从图1中可以看出, 2006年7月至2007年4月半干旱地区大气水汽含量月均值分别为2.08、1.90、1.44、0.92、0.73、0.39、0.29、0.46、0.71、0.81 g/cm2, 呈现出季节变化特征, 夏季最为丰富, 秋春季次之, 冬季最少, 这是由于半干旱地区夏季受东亚季风环流和青藏高原的影响, 春季地面温度开始从0℃以下升高至0℃以上, 积雪开始融化, 土壤解冻[6,7]。

2.2 半干旱地区不同季节典型晴空日大气柱水汽总量的日变化

半干旱地区四季典型晴空日大气水汽含量日变化情况见图2。从图2中可以看出, 半干旱地区四季晴空日大气水汽总量的日变化总体为夏季最为剧烈, 秋季次之, 冬春季最小, 主要由于太阳辐射的季节变化造成[8];春季水汽含量的日变化范围为0.47~0.53 g/cm2;夏季水汽含量日变化范围为0.75~1.26 g/cm2, 呈现出中午前后低、早晚高的特征, 原因为中午前后的太阳照射较强;秋季水汽含量日变化范围为0.57~1.00 g/cm2, 呈现出早晨低、晚上高的缓慢渐变特征;冬季水汽总量日变化范围为0.11~0.17 g/cm2, 呈现出中午前后稳定, 早晨和晚上轻微变化的特征, 原因可能为半干旱地区冬季全天逆温层形成的变化和气溶胶的加热效应所引起[7,9]。

注:a为春季, b为夏季, c为秋季, d为冬季。

2.3 半干旱地区降水转化率月变化

大气可降水量是指单位面积内从地表到大气顶层气柱内水汽全部凝结所能形成的降水量, 是评价空中云水资源的一个重要物理量[10,11]。图3为半干旱地区2006年7—12月及2007年1—4月的大气水汽总量月均值、月降水量与月降水转化率, 从图3中可以看出半干旱地区月降水转化率秋季最高, 其中10月最高, 为14.4%;夏季7—8月大气水汽含量虽然较高, 但降水转化率较低, 其中7月最低, 为0.87%;冬季大气柱水汽含量和降水转化率最低;春季降水转化率相对夏、冬季较高, 相对秋季较低。同时, 可以看出半干旱地区大气水汽含量变化与降水量变化不成正比, 说明降水还与水汽、动力抬升和不稳定量有关[4]。以上讨论说明半干旱地区夏、冬季具有人工增雨 (雪) 的潜力, 7月、8月、12月、2月是人工增雨 (雪) 最佳月, 1月虽然降水转化率较低, 但人工增雪的水汽条件很差, 合理开发云水资源对后期及次年的农业生产和生活有很大的影响。

3 结论

(1) 半干旱地区大气水汽含量夏季最丰富, 秋季次之, 冬春季最小。半干旱地区四季典型晴空大气水汽含量日变化情况总体为:夏季最为剧烈, 秋季次之, 冬春季最小。春季、夏季、秋季和冬季水汽含量的日变化范围分别为0.47~0.53、0.75~1.26、0.57~1.00、0.11~0.17 g/cm2。

(2) 半干旱地区夏、冬季具有人工增雨 (雪) 的潜力, 7月、8月、12月、2月是人工增雨 (雪) 最佳时段, 1月虽然降水转化率较低, 但人工增雪的水汽条件差, 合理开发云水资源对后期及次年的农业生产和生活有很大的影响。

参考文献

[1]李国平, 黄丁发.GPS遥感区域大气水汽总量研究回顾与展望[J].气象科技, 2004, 32 (4) :201-205.

[2]THOME K J, HERMAN B M, REAGAN J A.Determination of precipitable water from solar trasmission[J].J Appl Meteor, 1992 (31) :157-165.

[3]REAGAN J, THOM E K J, HERMAN B M.A simple instrument and technique for measuring columnar water vapor via near IR differential solar transmission measurements[J].IEEET rans Geosci Remote Sens, 1992 (30) :825-831.

[4]张文煜, 高润祥, 刘洪韬, 等.利用太阳光度计反演渤海湾西岸大气柱水汽总量[J].南京气象学院学报, 2006, 29 (12) :839-843.

[5]张海鸥, 郑有飞, 蔡子颖, 等.利用太阳光度计反演郑州地区的水汽含量[J].气象科技, 2009, 39 (5) :576-579.

[6]张玉香, 李晓静, 顾行发.利用太阳光度计测值估算北京上空水汽含量[J].遥感学报, 2006, 10 (5) :749-755.

[7]刘世祥, 王遂缠, 刘碧, 等.兰州市空中水汽含量和水汽通量变化研究[J].干旱气象, 2006, 24 (1) :18-22.

[8]梁宏, 刘晶淼.拉萨河谷大气水汽日变化特征[J].水科学进展, 2010, 31 (3) :335-342.

[9]张玉洁.黄土高原半干旱地区气溶胶光学厚度变化特征的初步分析[J].高原气象, 2008, 27 (6) :1416-1422.

[10]李帅, 谢国辉, 何清, 等.阿勒泰地区降水量、可降水量及降水转化率分析[J].冰川冻土, 2008, 30 (4) :675-680.

【反演计算】推荐阅读：

反演算法06-04

三维反演07-11

反演方法09-07

地震反演技术07-22