基于模型反演

基于模型反演（精选十篇）

基于模型反演 篇1

随着煤炭工业的发展,高度发展的机械化采煤需要不断提高地震勘探资料解释的精度,以便更加深入细致了解地下开采地质条件。常规地震技术难以实现地震勘探任务对地层纵横分辨率的要求,测井技术虽然分辨率很高,但难以有效实现井间地层参数横向变化对比的地质要求[1]。反演技术就是把地震技术与测井技术有机结合起来,反演结果中既有测井信息亦有地震信息,突破了传统意义上的地震分辨率的限制,理论上可得到与测井资料相同的分辨率[2]。

1 地震反演技术简介

地震反演是利用地表地震观测地震资料,以已知地质规律和钻井、测井资料为约束、对地下岩层空间结构和物理性质进行成像(求解)的过程。钻井资料的特点是纵向精细、横向稀疏,地震资料的特点是纵向粗略、横向密集,包含着丰富的岩性、物性信息。因此,通过地震反演技术把二者的优势有机的结合起来,经过地震反演,把界面型的地震资料转换成岩层性的测井资料,使其能与钻井、测井直接对比,以层岩为单元进行地质解释,研究储层特征的空间变化,可为勘探开发提供重要依据。

2 基于模型地震反演原理

基于模型地震反演的基本原理是建立在地震记录褶积模型基础上,即地震记录S(t)是反射系数R(t)和地震子波w(t)的褶积。

其实质就是从测井资料出发,根据钻井分层数据及时深关系对井进行精细时深标定,合成间隔不足一个采样点的薄层,建立一个初始波阻抗模型,用此模型合成地震剖面与实际地震剖面作比较,然后不断修改模型,使合成剖面最佳地逼近实际剖面,得到最终的地质模型[3],如图1所示。

初始模型的建立是一个人机交互的处理过程,对反演结果的好坏有直接影响。首先要对地震资料进行层位解释;然后通过合成记录,对每口井与井旁地震道做层位标定;最后以层位解释为控制,从井点出发,将测井数据外推内插,在三维空间的每一个点建立初始模型。这个过程实际上是把横向上连续变化的地震界面信息,与垂向具有高分辨率的测井信息相结合的过程。

3 地震岩性反演应用实例

3.1 工区概况

潘东矿西四采区南区三维勘探区面积1.66km2。矿区位于安徽省淮南市西北部,该区8煤层属二叠系上石盒子组煤层且全区分布稳定,从区内的钻孔数据得知8煤层的顶、底板岩性变化较大。

本区构造相对简单,但全区内可用于岩性分析的钻孔柱状及曲线仅有5个,并且分布不均匀。在地质建模时利用合成记录标定好的波阻抗曲线,以解释层位为基础采用了高斯克里金做全区差值,并对所建的地质模型反复进行修改,最终得到一个相对合理的三维地质模型。

3.2 岩性解释的技术步骤

3.2.1 建立准确的地质模型

精细的层位及断层解释,特别是断层的三维空间展布,以及与构造层位的准确接触关系,对于反演过程中的地质建模非常重要。地震资料解释得到的构造层位和断层,是反演处理中建立地质模型的基础资料,也是地震反演处理过程中加入地质认识和地质约束条件的重要依据。本区建立的地质模型如图2所示。

3.2.2 测井曲线校正及分析

为了提高合成记录的制作质量和层位标定精度,需要对本区5个孔的测井资料进行必要的校正。主要内容包括:岩芯资料的深度归位、井曲线的深度对齐、测井资料的环境影响校正以及测井数据的标准化[3]。

3.2.3 波阻抗曲线的重构

淮南潘东矿煤层为石炭二叠纪煤层,8煤属下石合子组,煤层全区分布稳定,平均厚度2.5 m。由于该区视电阻率曲线在砂岩及煤层反映为高电阻率,而泥岩电阻率相对平缓,为了有效的区分潘东矿8煤层顶、底板的岩性,利用视电阻率结合伽码伽码曲线对测井声波波阻抗曲线进行重构,使得重构后的波阻抗曲线能够突出煤层及其顶底板不同岩性的波阻抗特征。

3.2.4 子波估算注意的问题

子波的长度(保证一个主峰,两个旁瓣,最好是零相位),煤田资料一般为128 ms;估算子波的井旁道(最好是走向方向);估算子波的时窗至少是子波长度的三倍以上;边界不要卡在强轴上,避免子波发生畸变;要兼顾子波的波形及频谱(振幅谱,相位谱),特别要注意在地震主频带内,相位接近常相位。

3.2.5 合成记录的制作

根据本区的5个有曲线的孔提取了5个子波进行试验分析,最终确定本区反演所用的子波。根据所提取的子波制作了5个合成记录。7~8孔合成记录及标定,如图3所示。

3.3 潘东矿西四采区南区岩性分析

地质模型建好后,采用了基于模型的地震约束测井反演,得到的南区波阻抗三维数据体,再结合钻孔处分析好的8煤层及上下围岩信息,做全区岩性地质统计属性反演,得到三维岩性数据体。

本区五个井的连井剖面如图4所示,从图4中可以看到8煤底板下的岩性基本都是泥岩,中间大约3 m左右的区域为8煤层,沿8煤层上、下区域为砂岩。图5是将8煤层顶、底上下各20 m范围岩性嵌入在地震剖面中,可在剖面上直观的表现出8煤层上下岩性分布。

8煤层顶板岩性分布图如图6所示。从图6中可以很清楚地反映出砂岩、泥岩、砂质泥岩的分布区域,勘探区内中间部分区域为砂岩分布区,剩余部分全部为泥岩或砂质泥岩分布区域。

4 结 论

本文利用地震反演技术把具有高纵向分辨率的已知测井资料与连续观测的地震资料联系起来,推断地层岩性在平面上的变化,确定煤层的连续性、顶底板的岩性,为煤矿的开发与生产提供了地质依据。要得到好的反演结果,输入参数必须有可靠的保证,解释层位质量、测井数据单位、井数据约束的权值等等,都必须经过严格的落实。总之,反演技术是岩性地震勘探的重要手段之一,应尽快把反演技术应用于煤田地震勘探的工程实践中,实现煤田地震勘探从构造解释阶段向岩性解释阶段的新跨越。

参考文献

[1]张宏兵,杨长春.正则参数控制下的波阻抗约束反演[J].地球物理学报,2003,46(6):827-834.

[2]黄花香,吴战培,熊艳,等.速度反演技术在四川盆地储层预测中的应用[J].石油物探,2002,41(4):456-460.

基于模型反演 篇2

以淮南矿区为例,介绍了基于模型反演的`原理和方法以及利用反演技术对地震资料进行岩性解释的方法.

作 者：陶文朋 董守华 程彦 须振华 TAO Wen-peng DONG Shou-hua CHENG Yan XU Zhen-hua 作者单位：陶文朋,TAO Wen-peng(郑州煤炭工业(集团)有限责任公司,河南,郑州,450000;中国矿业大学资源与地球科学学院,江苏,徐州,221116)

董守华,程彦,须振华,DONG Shou-hua,CHENG Yan,XU Zhen-hua(中国矿业大学资源与地球科学学院,江苏,徐州,221116)

基于模型反演 篇3

关键词：界面反演，剩余重力异常，回归分析

1.引言

密度分界面与区域构造、储油构造、含煤盆地有密切的关系，因此计算密度分界面的起伏和深度的变化在区域构造研究、石油勘探、煤田勘探中具有重要的意义[1][2][3]。

通過分析前人对沉积盆地重震联合反演的研究成果，以及笔者对南华北地区区域地震剖面及构造格架剖面的拟合反演结果，我们发现通常情况下，主要沉积层界面深度与对应的剩余重力异常之间存在一种负相关的关系，即目的层深度越浅，对应异常越大，深度越深，对应异常越小。因此，我们期望运用已有的深度异常信息建立目的层密度界面深度与剩余重力异常之间的回归方程，通过该方程推算出未知区域的深度信息。

2方法原理

2.1线性回归分析

在密度界面起伏平缓的情况下，可以认为重力异常与界面的起伏呈近似线性关系，即

（2-1）

式中： 界面深度， 为界面起伏引起的重力异常； 、 为两个常数，他们与异常起算点处的界面深度和界面上下物质层的密度差有关。

为应用（2-1）式求取深度，至少要知道界面上两个点的深度，以确定 、 两个系数值。若存在n个已知点，它们的深度 ，则根据最小二乘原理，为确定系数 、 ，应使各点的深度 和由（2-1）试计算出的深度 的偏差平方和为最小，即

（2-2）

令 ， 分别等于零，可得：

（2-3）

（2-4）

以上两式联立，解之得：

（2-5）

（2-6）

式中 为 的省略形式。

系数 确定后，就可以由（2-1）式计算出测点下方各界面的深度[1][4]。

2.2抛物线回归分析

与线性回归相比，抛物线回归分析只是给线性回归方程增加了一个二次项，如下式：

（2-7）

应用（2-7）式求取深度，至少要知道界面上三个点的深度，以确定 、 、 三个系数值。对存在n个已知点的情况，同样可以根据最小二乘原理，使各点的深度 和由（2-7）试计算出的深度 的偏差平方和为最小，以确定 、 、 三个系数值[70]。即：

（2-8）

（2-9）

（2-10）

联立以上三式，解之得：

（2-11）

系数 确定后，就可以由（2-7）式计算出测点下方各界面的深度。

2.3算法流程

回归分析的算法流程如图所示，每一个计算环节简单介绍如下：

图2.1 回归分析算法流程图

（1）数据读取

包括剩余重力异常网格数据和已知控制点信息的读取。

（2）搜索控制点

搜索与当前测点距离在指定范围内的已知点，若已知点过少，如对于抛物线回归分析已知点少于4个，则放弃计算该点，若已知点过多，则按距离测点距离远近对已知点排序，取距离最近的指定数目的已知点。

（3）建立回归方程

运用控制半径范围内已知点的深度和异常信息根据前两节所述原理建立界面深度关于剩余异常的回归方程，计算出回归系数。

（4）测点计算

将当前测点的剩余重力异常值代入回归方程，求取其深度值，并对数值的合理性做出判断。

（5）数据输出

若当前深度值求取合理，则输出对应测点的坐标、深度、异常以及相应的回归系数等信息，并进入下一测点的计算，重复1、2、3、4步骤，否则不输出当前测点信息，直接进入下一点的计算。

整个计算流程不是很复杂，在VC6.0中编程实现。计算时需要注意一些细节。首先，对于搜索半径及其范围内制点数目的选取要合适；其次，研究区目的层的深度是有一定范围的，回归分析计算出的深度若超出这个范围应该剔除，而深度范围的确定需要参考地质、钻孔及剖面反演资料。

3约束条件

这里的约束条件包括方法本身的应用条件和对控制点要求。

回归分析的应用前提是密度界面的起伏变化在一定范围内是平缓的，变化越平缓则计算的精度越高。例如当界面起伏最大倾角小于三度，起伏幅度不超过界面最大深度1/10时，由（2-1）式所得的结果的最大相对误差不超过7%；即使界面最大倾角到11度，起伏幅度达到界面最大深度的1/5，带来的最大误差也小于8%[2]，而采用（2-7）式时会更突出一些局部细节，相对来说，误差还会减小。

4剖面回归分析验证

为了验证回归分析的有效性，同时比较线性回归分析和抛物线回归分析的反演效果，以研究区三叠系地层为例，我们把部分地震剖面和区域格架剖面的反演拟合得到的三叠系底界面深度值与对应点剩余重力异常值作为已知控制点，把剖面切割规则剩余重力异常网格得到的异常值作为待求测点，运用两种回归分析方法分别进行了反演计算，部分剖面结果如图4-1所示。图中蓝色十字叉点表示已知控制点，绿色线为线性回归分析反演结果，红色线为抛物线回归分析反演结果。

太康线

EW03线

图4.1 部分剖面深度异常回归分析效果对比图（三叠系底界面）

总体来说，两种回归分析方法求得的深度值都大体反映了剖面下方三叠系底界面深度的变化趋势，因而都具有可行性。在深度变化比较平缓的区域，它们求得的深度值基本没有差异，在深度变化较大的区域，二次回归分析的结果与控制点深度更为接近，更能反映一些深度变化的细节。因此，对于平面的深度回归分析反演，我们优先选取二次回归分析方法。

5区域密度界面反演分析实例

以南华北地区盆地为例。南华北地区（又称华北盆地南部）地处中原和两淮地区，包括河南省和安徽省的大部分以及江苏省的西北部、山东省的西南部。区内诸多盆地是不同构造阶段，多种构造动力体系联合与复合作用的最终產物。对研究区主要密度界面，下古生界底界、上古生界底界深度异常二次回归分析结果如图5.1～5.2所示。

结合研究区区域地层特征分析，反演结果反映了各界面的基本分布格局，即在三门峡—舞钢—信阳—舒城一线以北区域各地层界面起伏变化，该线以北则进入北秦岭逆冲推覆构造带，区内主要分布太古代、元古代区域变质岩及不同时期的侵入岩体，因而主要沉积层深度为零。因此，也从地质角度证明了运用回归分析求取主要目的层残存分布的合理性。

进一步分析结果可以看出，各沉积层底界面深度分布是与区域构造“南北分带、东西分块”的特征是相一致的。在南北方向上，由北侧的济源凹陷、开封坳陷到太康隆起，再到中部的周口坳陷往南经长山隆起进入信阳合肥盆地，各目的层底界面的深度经历了深、浅、较深、浅、深的交替变化，而在东西方向上因为与主要构造单元分布平行，深度变化比较平缓。

图5.1 南华北地区下古生界底界面深度图（单位：m）

图5.2 南华北地区上古生界底界面图（单位：m）

6结论

本文运用基于回归分析的反演模式，建立了盆地、坳陷区剩余重力异常与主要目的层深度之间的回归方程，由此推算未知区域的目的层深度分布情况，并引入实例，计算区域密度界面的分布情况，通过已有地质、地球物理特征认识验证了计算结果的合理性。我们认为这种方法是可行的。

参考文献：

[1]曾华霖，重力场与重力勘探[M]. 地质出版社，2005.6

[2]肖鹏飞，陈生昌，孟令顺. 高精度重力资料的密度界面反演[J]. 物探与化探，2007b，31（1），29-33

[3]韩道范等.利用重力异常反演多层密度分界面的理论和方法[J]. 地球物理学报，1994，37（1），272-281

基于模型反演 篇4

电离层作为地球大气的组成部分之一,对经其传播的电磁波将产生折射、吸收和时延等现象。无论是确保天波雷达系统正常工作[1],还是反演地球大气参数(如利用GPS无线电掩星法)[2],研究电离层的信息都是不可或缺。通常这些信息是利用电离图反演得到的[3],因此,提高反演算法稳健性和自动化程度,降低计算量是实现上述工程应用所必需解决的问题。

电离图是由电离层垂直测高仪扫频得到并利用测量回波和发射脉冲之间的延迟,得到电离层反射高度(称之为虚高)随扫描频率变化的曲线。根据该曲线,以及所建立的电离层模型,求解模型参数,进而反演得到电离层参数信息(主要是电离层各层的临界频率、峰值高度、半层厚度以及一些电离层不规则结构的信息)。如图1所示。

电离图反演的基础是电离层模型。在上个世纪九十年代,CCIR推荐采用修正B/D模型以及多准抛物线(MQP)模型。其中,MQP模型层数可调节,并采用准抛物线层与翻转抛物线层组合保证了层间电子浓度梯度的连续,较好表征了各层电子浓度分布特性,被人们广泛采用[4]。

基于MQP模型的电离图反演方法有切比雪夫多项式法[5]、偏导迭代法[6]等。后者相对前者简洁,参数含义明确而被普遍采用。对某一层参数反演时,偏导迭代法是通过输入至少两点经该层反射的电离图数据进行迭代而求解的。而在使用中发现,一旦数据选择不当,该算法将失效。所以,该方法要求数据录入者具备一定电离层知识,能选择合适的数据点,在算法失效时需人为干预,重新选择数据点,不便进行自动化处理。

本文将从MQP模型虚高表达式入手,在合理近似的基础上,推导出虚高与各电离层参数关系的简约表达式,进而寻求一种电离图数据反射层判断准则。依据该准则,提出一种电离图参数反演的解析算法,减小计算量,以满足反演的实时性要求。

2. 多准抛物线模型的虚高

电离层可以用多层来描述,包括E层、F1层和F2层。每层可用最大等离子频率fci、其所在高度rmi以及rmi与该层底部rbi之间的距离为半层厚度ymi来描述,各层等离子频率随高度变化可用准抛物线来描述[7],其实例如图2所示,表达式见式(1):

式中L表示电离层层数,当i为奇数时,式(1)取负号,代表常规的抛物线层,满足rmi≡rbi+ymi;i为偶数时,式(1)取正号,代表翻转抛物线层,其实现了从低的常规抛物线层顶部,到高的常规抛物线层底部的平滑过渡,其参数可以由上下抛物线层参数计算得到:

其中,

当不考虑地磁场的影响,忽略电子与分子、离子之间的碰撞效应时,电离层媒质的折射率n为[8]:

式中f发射电磁波频率,根据Snell定理[9],经电离层反射的单程传播路径P'的表达式为:

其中,r为电磁波传播前沿到地心的距离;n为r处的电磁波折射率;β为r处电磁波的入射余角(与入射角互余);r0为地球半径(近似为6370km);β0为发射仰角,rt为射线反射顶点到地心的距离。对于电离图中的虚高可表示:

3. 翻转层与常规层反射虚高的判别

由于在电离图中常规抛物线层和与之相连的翻转抛物线层的虚高连在一起,难以区分,易导致偏导迭代法失效。本节将讨论如何在电离图中区分翻转层与常规层。

首先从反射顶点在翻转层时虚高表达式入手。以第一翻转层(第2层)的虚高表达式为例,进行研究,其结果也适用于第二翻转层(第4层)。设电离图读值为(h,f),则由第二层反射的虚高可以表示为:

其中P0'和P1'分别代表自由空间传播路径和第1层(E层)传播路径,这两个值通过第1层参数可直接计算获得。如何反演第1层参数将在下节中介绍。

式中

通常电离层的翻转层很薄,不足百公里(ym2≤100),而rb2>r0=6370km,f>fc2,则有

其中

求解上述一元二次方程,取x的正值解,即ym2与rb2比值为

当读取点(h,f)位于翻转层时,ym2与rb2的比值x不随(h,f)变化而变化;若x随(h,f)变化而变化,则说明该选择点不在翻转层中。由此,可从连续的观测点中找到翻转层的最高频率fb2,其满足

于是,可以利用fb2作为翻转层与常规层的频率分界线,实现了两层反射虚高的区分。

4. 电离层参数反演之解析算法

从式(6)可见,电离图虚高表达式复杂,很难解析法。为此,本节将对表达式进行合理近似处理,在保障计算精度的情况下,提高计算效率。

由于第1层虚高表达式与上层有较大差异,本文将第1层和上层(>1层)反演问题分开讨论。

(1)第1层参数反演

利用式(5),求得第1层反射虚高表达式如下:

其中

由上式可知,化简虚高表达式的关键是对数项的简化。因为rb1>r0=6370km,且F1<1,考虑到通常情况下ym1≤300km,于是,利用rm1≡rb1+ym1,并考虑到rm1/(F1ym1)>rm1/ym1>rb1/ym1>θ>1,则

代入式(10),整理得:

对此二元二次方程,可利用两个观测点数据(h1',f1')和(h2',f2)完成求解。

(2)上层参数反演

借助第三节的确定的fb2和x,可以推导一种反演上层参数的解析法。

根据rh2处电子浓度表达式

可求得rb2为

进而求出ym2=x×rt2。

利用式(2),即可求得rm3、rb3和ym3。同理,若电离层是5层结构,可用类似相应求出第5层参数。

可见上述方法是通过给定一个翻转层观测点(计算x)以及翻转层的最高频率fb2,完成第3层参数的反演。并且,上述方法没有矩阵求逆过程,所以不存在偏导迭代法计算量大、初值选择不当会造成算法失效的情况。输入三组观测点,对比偏导迭代法,其仿真分析结果如表1。仿真参数为:fc1=3.2MHz,rm1=6474km,ym1=14km,fc3=4.61MHz,rm3=6553km,ym3=77km,其电离图如图3所示。

由表1可知:

1)对于第1层(E层)反演,由于E层的高度很低以及层厚度相对很薄,结果可以看出,只要选择的观测点是位于E层(即观测点频率小于fc1),偏导迭代法可得到准确的反演结果,而由于对虚高表达式进行了近似处理,解析算法估值会出现一些偏差,但最大偏差不超过0.2%。由于电离图中fc1很易读出,其为E层反射描迹与高层的描迹分界线,如图3中fc1处所示,所以观测点选取条件很容易满足。

2)由于第一翻转准抛物线层(第2层)与常规准抛物线层第3层(F1层)反射描迹连在一起,如图3所示,无法区分,所以针对F1层两个的观测点位置,做了3组参数反演仿真。第一组中两个观测点分别位于第2层和第3层,第二组两个参考点位于第3层,第三组两个参考点均位于第2层。从仿真结果可以发现,只要有一个观测点位于第3层,偏导迭代法都可得到理想的结果;但是,当两个观测点均位于第2层时,估值失败。研究表明造成此结果的原因是在此时偏导迭代法中迭代矩阵奇异,无法求逆,算法失效。

3)本文提出的解析反演算法对观测点的位置没有限制,所获得的参数估计误差均在3%以内。考虑到经第2层反射时虚高要小于经第3层反射结果,而解析算法是基于近似处理,选择位于第2层的观测点的估值结果要优于选择位于第3层的观测点。

4)由于解析算法采用了解析计算,其计算量要远远小于偏导迭代法,其计算用时通常降为1/15以下,并对求解初值设定宽松,适合于实时性要求较高且专业性不强的场合,如天波雷达。

5)根据本文提出的层判断准则,可以自动判断观测数据所在层,在实时要求不高,可将次准则与偏导迭代法相结合,提高该法的自动计算以及稳定性。

考虑到fb2可能存在读取偏差,可以选择一个第3层观测点作结果校正参考点,即在输入的fb2附近微调fb2,得一组参数估值,用此参数估值计算虚高。当虚高最接近参考点时,此刻的参数估值为最终结果,处理流程如图3。按此思路,做fb2微调处理仿真结果如表3。

如表2所示,当存在fb2读值误差时,只需输入一个第3层的参考点,可以自动纠正fb2偏差,代价是计算量增大,但和偏导迭代法相比,性能改善还是很显著的。

综合本章的第1层参数反演和上层参数反演法,可以得到一种新的电离图反演方法。虽然在估值精度方面稍逊色于偏导迭代法,但是计算性能和稳定性都有显著的改善,适合于实时性电离图参数反演场合。

5. 结论

针对电离图实时自动处理的需求,本文提出了一种电离图解析反演算法。相对现今普遍采用的偏导迭代法,该算法的计算性能和稳健性有了很大改善。同时,本文提出的翻转层和常规层区分的判断准则也可用于偏导迭代法的初值选取中,用以改善其稳定性。

参考文献

[1]D.Bourgeois,C.Morisseau,M.Flecheux. Over-the-Horizon Radar Target Tracking Using Multi-quasi-parabolic Ionospheric Modeling[J].IEE Proc.-Radar Sonar Navig. 2006,153(5):409-416.

[2]蒋虎.GPS无线电掩星技术反演地球大气参数中若干问题的研究[C].中国科学院博士论文.2002.

[3]S.J.Anderson.Generation assimilation of propagation advice for HF skywave radar systems.2008.

[4]Hou Chengyu.A study on dynamic modeling and method to demodulate disturbance of ionosphere in HF band[J].Ph D. dissertation,Harbin Institute of Technology,2008:5-6.

[5]B.W.Reinisch,X.Huang.Automatic Calculation of Electron Density Profiles from Digital 1.Automatic O and X Trace Identification for Topside Ionograms[J]. Radio Sci,1982,17:421-434.

[6]R.J.Norman.An Inversion Technique for Obtaining Quasi-Parabolic Layer Parameters fromⅥIonograms[J].IEEE Radar 2003. 2003:363-367.

[7]D.C.Baker,J.J.Burden.Multisegmented Parabolic Model for Ionospheric Electron Density Distribution[J].IEE Conference Publication.1991,392:183-188.

[8]P.Jiao.Some New Experimental Research of HF Backscatter Propagation in CRIRP[D]. 2004 Asia-Pacific Radio Science Conference-Proceedings. 2004:13-16.

基于模型反演 篇5

为了提高垂直电测深二维反演的收敛速度和反演精度,给出了利用直接反演公式确定初始模型的`方法.在直接反演过程中,通过综合利用三次样条求导法、五点求导法和列维尔插值计算电测深曲线在每个极距处的梯度,并对直接反演结果进行深度校正,获得了更准确的反演结果,为二维反演提供了较好的初始模型.通过对模拟和实测电测深数据进行反演试算,验证了采用直接反演法给定非均匀初始模型的可行性和有效性.

作 者：刘海飞 阮百尧 廖建龙 蒋小红 LIU Hai-fei RUAN Bai-yao LIAO Jian-long JIANG Xiao-hong  作者单位：刘海飞,LIU Hai-fei(中南大学,信息物理工程学院,长沙,410083;桂林理工大学,广西地质工程中心重点实验室,广西,桂林,541004)

阮百尧,RUAN Bai-yao(桂林理工大学,广西地质工程中心重点实验室,广西,桂林,541004)

基于模型反演 篇6

大气水汽含量不仅在人类日常活动中占有重要的角色,而且对全球热带分布具有重要的影响,因此对大气水汽含量进行探测对于水汽循环、水汽输送以及水资源评估等气象研究具有重要意义。常用的水汽探测手段有气象卫星遥感、雷达探测、水汽辐射计和无线电探空探测。20 世纪90 年代以来,人们开始利用GPS理论和技术遥感地球大气,并逐渐形成地基GPS气象学［1］,其基本原理是通过对GPS定位中多路径效应造成的GPS信号延迟进行逆处理,最终利用GPS延迟信号反演出大气水汽含量。地基GPS大气水汽探测具有其它传统水汽探测［2］所不具有的很多优点,如观测数据可以全天候实时获得,反演精度高,并且无须校准［3 － 5］。目前国内GPS观测站网建设正在迅速普及到每个省和地区,每个省和地区的GPS观测站网的站点密度也在增加,从而使GPS网在天气监测和天气预报中发挥的作用越来越大。然而为了满足GPS水汽反演的高精度和实时性要求,加权平均温度起决定性作用,因此本文通过拟合出适合上海地区的加权平均温度曲线的方法来实现水汽反演,并利用上海市2012 年8 月的可降水量反演结果与探空计算结果进行比较,结果证明此方法提高了转换精度并保证了反演的实时性。

1 反演原理与方法

1. 1 反演原理

在高精度GPS定位中,GPS探测信号穿过地球大气层时会产生时间延迟,这是由于信号在传播过程中会因电离层和大气折射造成传播速度减小和路径弯曲［7］。电离层延迟与电磁波的频率平方成反比,可以通过接收双频信号,线性组合后消除。大气延迟一般是指从地面到50km高度的非电离大气对电磁波的折射,又称中性延迟。由于折射的80%发生在对流层,所以通常称为对流层延迟。对流层大气在天顶方向上的总延迟中,其中干空气造成的称为干延迟或静力延迟,占90% 以上,水汽造成的湿延迟占10% 左右。通过对GPS观测数据、GPS卫星星历、GPS测站坐标进行解算,获得GPS信号在对流层传输的天顶总延迟ZTD,同时利用改进的经验模型将GPS信号在对流层传输时受到的干延迟ZHD解算出来,接着将总延迟量ZTD与干延迟量ZHD相减得到湿延迟ZWD,最后利用实时计算的转化因子 Π,将湿延迟量转换为大气可降水量PWV。因此,可得GPS大气可降水量反演步骤如图1 所示。

1. 2 理论模型

假设对流层的折射率是n,那么信号传播中的对流层折射改正为:

式( 1) 中n为大气折射率( 真空中的折射率为1) ,s为信号传播路径。

由于n - 1 的数值很小,为了方便于计算,引进大气折射指数N:

其中大气折射指数N表示大气密度值,它是关于压强、温度和水汽含量的关系式,常用的是1954 年Smith与Weiniranb利用大量实际探测资料拟合出的模型［8］:

由式( 3) 可知,大气折射指数N分为两部分,一部分是与大气压P和绝对温度T相关的静力学折射部分( NDry) ,另一部分则是与水汽压ew和绝对温度T相关的湿气折射部分。

将式( 2) 、式( 3) 带入式( 1) 得到对流层折射改正［9］为:

由式( 4) 可知,要计算出GPS信号传播过程中每点的对流层折射,必须测量出路径上每点的气象因子,这是实际测量中不可能实现的,能够获得的气象因子一般是通过地面测站测得,因此我们必须利用地面测站的气象因子通过数学经验模型计算每点相应的海平面或大气的气象因子,本文采用的是改进的Saastamoinen模型［10］。

1. 2. 1 Saastamoinen经验模型

该经验模型将对流层分为两层进行积分,第一层是地面到12km高空附近的对流层顶,在此大气层中气体温度值随着高度的增加而递减,其递减率为6. 5℃ /km; 第二层是从对流层顶一直到50km附近的平流层顶,此空间大气温度一般设为常数。

Saastamoinen经验模型的天顶总延迟分为干分量和湿分量部分:

将式( 6) 、式( 7) 代入式( 5) 得:

式中,ew为水汽压,Ps为干大气压,f( B,H) 为重力加速度,它随测站纬度  和椭圆地球表面的高度H变化。

1. 2. 2 Saastamoinen改进模型

对Saastamoinen经典模型进行如下改进:

其中,T是测站绝对温度,P代表大气压,ew是指水汽压,z为卫星的天顶距离,B和 δR可以通过插值表获得,式( 10) 中的T、P、ew可以利用标准大气参数获得,并且气压、温度、相对湿度分别和高程的关系可由式( 12) ~ ( 14) 得到:

其中,P0、H0、RH0是标准气压、温度和相对湿度的缺省值,具体数值如下:

P0= 1013. 25mbar,H0= 0m,RH0= 50%

1. 3 反演可降水量及误差分析

1. 3. 1 反演大气可降水量

大气可降水量( PWV) 是指单位底面积空气柱体内所含水汽全部凝结降落的总水量。利用GPS数据资料解算出来的ZTD和ZHD相减得到天顶湿延迟ZWD,ZWD和大气可降水量之间的转换关系为:

式( 15) 中的转换因子 Π 精度取决于大气加权平均温度Tm,Π 的精确计算公式如下:

式( 16) 中 ρw是液态水的密度,Rv为水汽气体常数,k'2和k3是大气物理参数,它们的经验值通常取为:

由式( 16) 可知,天顶湿延迟量转换成可降雨量的过程和加权平均温度Tm有直接关系,Tm的精度直接影响可降水量的精度。目前加权平均温度Tm的计算方法主要分为: 一是利用探空资料中对流层不同高度的各点温度T和水汽压e进行离散化积分,如式( 17) ; 二是1992 年Bevis对美国地区多年的探空资料进行回归分析后得到: Tm= 70. 2+ 0. 72Ts。

其中ei、Ti、Δhi分别为第i层大气的平均水汽压( h Pa) 、平均温度( K) 和厚度( mm) 。

第一种方法的局限在于探空资料的不足,每天只有两次高空探测,因此时间分辨率很低,但是却是目前为止最精确的算法。第二种方法比较方便可取,但是Bevis公式是针对欧美地区的拟合曲线,在中国地区并不适合。所以可以利用第一种方法确定的对流层加权平均温度Tm与地面测站所在点的温度Ts拟合出与上海地区相匹配的关系式。本文就是利用上海市地面与探空资料进行加权平均温度的曲线拟合,基本流程如图2 所示。

1. 3. 2 误差分析

为使式( 16) 中转换因子 Π 的转换精度小于2mm,这里就要讨论加权平均温度的容许误差。通过式( 16) 可知,除了加权平均温度以外,大气物理参数k'2和k3也影响转换因子的精度,取式( 16)对k'2,k3和Tm的微分得:

设k'2、k3和Tm的误差分别为 σk'2、σk3和 σTm,并且这三个参数相互独立。则根据误差传播定律,转换因子 Π 的误差 σΠ为:

由式( 19) 可知,转换因子误差 σΠ是加权平均温度Tm的减函数,所以取Tm的最小值就可以计算出误差 σΠ的最大值,将上海地区Tm最小值270K以及其他物理常数代入式( 19) 得:

由于天顶湿延迟ZWD最大值为500mm,则利用此最大值和式( 20) 代入式( 15) ,得到转换因子误差引起的可降水量误差如表1 所示。由表1 可知,要使可降水量误差在2mm以内,则加权平均温度精度至少要在7K以内。

2 仿真结果

2. 1 GAMIT处理GPS观测数据解算湿延迟

GPS数据分析处理是在Unbuntu下安装的软件GAMIT 10. 4 下进行,主要分为四个部分,首先是包括卫星星历G文件、GPS观测数据o文件、气象观测数据m文件等数据准备及收集工作; 接着用GAMIT软件处理GPS观测数据o文件,获得每个时段对应的解; 然后采用卡尔曼滤波方法对多时段进行综合解算,从而获得网平差结果和测站坐标等参数; 最后一步是在GAMIT中利用前一步平差得到的测站坐标,强约束待求测站,反演天顶对流层总延迟。

2012 年8 月上海市降水较常年偏多,主要是受台风影响,共有5 个台风影响上海。因此,本文采用上海宝山站2012 年8 月共31 天GPS观测数据进行解算,数据准备完成之后,编辑测站相关配置文件,并解算出当天每隔1 小时的湿延迟量,利用命令sh_ nrms检测其均方根残差均在0. 5 以内,从解算结果中提取出每天8 时与20 时的湿延迟量结果,如图3 所示。

2. 2 加权平均温度的曲线拟合

本文选取上海宝山( 站号: 58362) 2012 年1 月~ 2012 年12 月共一年的探测资料,采用回归分析法利用MATLAB编程,拟合出曲线关系式:Tm= 0. 8936Ts+ 19. 1321,如图4 所示。结果表明,加权平均温度的最大误差为6K,那么由表1 所知,可降水量转换误差均在1. 823mm以内。

图 3 2012 年 8 月上海宝山每天 8 时与 20 时的对流层湿延迟量

图 4 上海地区加权平均温度 Tm与地面温度 Ts的关系

2. 3 可降水量反演曲线

利用图3、图4 所示的天顶湿延迟量和加权平均温度曲线,根据式( 15) 计算得到上海市2012年8 月每天8 时与20 时GPS估计的的可降水量与探空计算结果比较,如图5 所示。结果表明,两种方法得到的大气可降水量的误差在3mm以内。

图5 2012年8月GPS估计的可降水量与探空计算的比较

本文在提取参与拟合的观测数据时已经对无效数据进行整行删除,并且由拟合出的曲线方程与实际计算值比较可知其相差最大只有6K,且数量较少,因此本文没有进行拟合方程的系数显著性检验,且通过图5 的最终可降水量计算的应用中也可知结果较准确。

3 结束语

本文分析了GPS可降水量反演的方法与特点,并对反演过程中可能存在的转换误差进行了分析讨论,利用2012 年整年上海市高空和地面探测数据拟合出适应于该地区的Tm与Ts的关系式,并利用其转换出2012 年8 月共31 天的大汽可降水量,最后将每天8 时和20 时的GPS估计的可降水量与对应时间的探空计算结果进行比较,结果表明,GPS估计的可降水量与探空结果一致,误差在3mm以内,该方法解决了反演可降水量的转换精度问题。本方法的下一步工作是对更多年份的数据进行验证分析比较,进而能够进行相关应用研究。

摘要：GPS可降水量反演的本质是利用GPS观测数据o文件解算出垂直方向上传播信号的湿延迟量,并将湿延迟量转换为可降水量,这两个参数的转换需由加权平均温度(T m)来完成。本文利用上海市2012年1月12月的地面和高空探测资料采用MATLAB编程拟合出适合上海市的T m与地面温度(T s)的关系式,并反演了受台风“海葵”严重影响的2012年8月上海地区可降水量。结果证明,该T m和T s的关系式不仅弥补了探空数值积分法的缺点,而且能使GPS可降水量反演的转换误差控制在3mm以内。

基于模型反演 篇7

自20世纪70年代初期开始,遥感技术逐渐应用到陆地水体的研究中,从单纯的水域识别发展到对水质参数进行遥感定量监测。随着遥感技术的发展和对水质参数光谱特征及算法研究的不断深入,水质遥感监测经历了分析方法20世纪(80年代以前,主要针对开阔海洋)一经验方法(80—90年代,多光谱遥感技术的应用)一半分析方法(90年代以后,高光谱遥感技术的应用)的发展过程[1]。经验方法与半分析方法都是通过对遥感数据和与其(准)同步的地面实测数据进行适当的统计分析得到水质参数的估测算法,是目前水质遥感监测最常用的方法。根据研究对象的不同,水质遥感可分为大洋开阔水体遥感、近岸水体遥感和内陆水体遥感[2]。内陆水体遥感监测的主要水质参数包括叶绿素α、悬浮物、黄色物质等。其中叶绿素α是反映内陆水体富营养化程度的一个重要参数,通过遥感定量反演模型提取叶绿素α浓度是内陆水体富营养化监测的重要途径。

内陆水体叶绿素α浓度反演常用的多光谱遥感数据包括Landsat MSS/TM/ETM+、EO—1 ALI、NOAA AVHRR、SPOT HRV、IRS—1数据等[3]。综合空间、时间、光谱分辨率和数据可获得性特征,TM数据是目前叶绿素α浓度反演中使用最广泛的多光谱遥感数据,但缺乏通用的模型。内陆水体叶绿素α浓度反演常用的高光谱遥感数据包括美国的AVIRIS、MODIS和Hyperion、欧空局的Envisat MER-IS、加拿大的CASI、芬兰的AISA、日本的ASTER和GLI以及中国上海光机所的OMIS成像光谱数据[3]。其中MODIS数据具有较高的时间分辨率,适合大范围内陆水体的叶绿素α动态监测[4],但其空间分辨率较低。闻建光等[5]研究表明,Hyperion数据可以满足内陆水体叶绿素α提取的精度要求,但时间分辨率太低,不便于叶绿素α的实时动态监测。环境一号卫星是中国国务院批准的专用于环境和灾害监测的对地观测系统,由两颗光学卫星(HJ—1A、HJ—1B)和一颗雷达卫星(HJ—1C)组成,具有大范围、全天候、全天时、动态的环境和灾害监测能力。该卫星搭载了CCD相机和超光谱成像仪(HSI),可获取高时间分辨率、高光谱分辨率、中等空间分辨率的对地观测数据,数据可获取性强,已经成为我国环境遥感监测的主要数据源之一。近年来,国内开始利用这一数据开展了一些叶绿素α浓度反演模型研究[6,7]。

太湖是我国最大的存在严重蓝藻水华的浅水内陆湖泊,是国务院重点治理的富营养化水域之一,对太湖进行叶绿素α浓度反演模型研究具有重要代表性。本研究选取太湖为实验区,基于环境一号CCD多光谱数据和同步地面实测叶绿素α浓度数据,构建叶绿素α浓度反演的经验模型。并基于环境一号HSI高光谱数据,结合内陆水体中叶绿素α、悬浮物、黄色物质和纯水的固有光学特性以及三波段模型理论,构建三波段模型。对两种定量反演模型对比分析,探讨环境一号卫星影像用于内陆水体叶绿素α遥感定量反演的最优模型。

1 数据获取与预处理

1.1 数据获取

采用的太湖地区环境一号CCD多光谱数据和HSI高光谱数据从中国资源卫星应用中心获取。选取云量覆盖较少的遥感影像,分别为2009年10月6日的HJ—1B CCD数据和2009年8月18日HJ—1A HSI数据,数据均为经过系统几何校正的2级数据。地面同步实测数据包括37个采样点的经纬度数据和叶绿素α浓度数据[8],采样点均匀分布于湖面,见图1。

1.2 遥感数据预处理

环境一号数据预处理前,首先通过ENVI环境小卫星数据处理工具包(可免费下载)读取数据。预处理过程主要包括几何精校正、大气校正和太湖水域的提取。几何精校正选用从国际科学数据服务平台下载的Landsat ETM+数据为参考图像,在ENVI 4.8支持下,均匀选取校正点,图像重采样采用最邻近方法,总误差控制在0.6个像元内,完成校正。大气校正采用ENVI 4.8软件下的FLAASH大气校正模块,输入相关参数,计算得到大气校正后的反射率图像。并在ArcGIS软件中人工数字化提取太湖水域。

2 模型构建与验证

2.1 基于环境一号CCD多光谱数据的叶绿素α浓度遥感定量模型构建与验证

叶绿素α浓度反演主要是利用传感器接收到的辐射值来分析叶绿素α的含量,通过遥感数据单个波段或者多个波段的组合与(准)同步的地面实测数据进行适当的统计分析建立反演模型。

随机选取湖区30个采样点数据用于模型反演,其余7个采样点用于模型验证。利用Excel对CCD四个波段的反射率值和叶绿素α浓度值进行皮尔逊相关性(r)分析,各波段的r依次为-0.616 0,-0.696 1,-0.584 7和0.746 2。可见,仅利用单波段反演叶绿素α浓度的模型相关度不高。有研究表明,采用多波段反射比可以部分消除水表面光滑度和微波随时间和空间变化的干扰,并在一定程度上减小其他污染物的影响[9]。因此,采用最大正相关的近红外波段(CCD4)与其余三个负相关波段(CCD1,CCD2,CCD3)分别进行比值运算,并将运算结果与叶绿素α浓度值进行线性、对数、指数和多项式等回归分析。结果表明,此方法有效的提高了两者间的相关性。其中,CCD4/CCD2与叶绿素α浓度具有最大相关性,r为0.929 1。以CCD4/CCD2为自变量,叶绿素α浓度为因变量建立经验反演模型:

式中,CChla为叶绿素α浓度,单位mg/L;RCCD4,RCCD2分别为环境一号CCD数据第4和第2波段的反射率值。

将建立的叶绿素α经验模型应用于整个湖区,得到叶绿素α浓度的反演值,与7个验证值进行对比分析(如图3)。通过计算,模型估测的最大相对误差为16%,最小相对误差为3%,平均相对误差为6.2%。均方根误差(RMSE)为0.004 2 mg/L,远低于叶绿素α浓度的平均值0.047 4 mg/L。金焰等[10]研究显示,基于环境一号CCD数据,采用RGB:343波段假彩色合成能够对太湖蓝藻水华进行监测,且目视效果较好。因此,采用此方法提取蓝藻水华,并与叶绿素α浓度反演结果对比,两者吻合度较好(如图4)。可见,基于环境一号CCD多光谱数据的经验模型可以应用于太湖叶绿素α浓度的估算。

2.2 基于环境一号卫星HSI高光谱数据的叶绿素α定量反演模型构建与验证

内陆水体大都属于二类水体,其光学特性复杂,主要由叶绿素α、悬浮物、黄色物质和纯水四种因子共同决定,对于浅水水体,还要考虑水体底质的影响,并且各因子间相互干扰。因此,对于内陆水体叶绿素α的遥感定量反演,选择受其他物质干扰小而与叶绿素α显著相关的波段组合算法是关键。Gitelson等[11]曾基于植被、土壤等物质的光学性质提出了三波段模型来估算叶绿素含量。Zmiba[12]、Dall’Olmo[13]等研究了适用于富营养化浑浊内陆水体叶绿素α反演的三波段模型。随着高光谱技术的发展,三波段模型越来越多的被应用于叶绿素α的定量反演,模型综合考虑了内陆水体的固有光学特性,有效的降低了其他物质的干扰,是叶绿素α反演效果较好的一类模型。

基于环境一号HSI高光谱数据,选取太湖富营养化较为严重的南部沿岸湖区和东太湖为研究区,结合内陆水体中叶绿素α、悬浮物、黄色物质和纯水的固有光学特性以及三波段模型理论,选取反演叶绿素α浓度的三个最优波段,构建三波段模型。研究区内共有13个采样点,随机选取9个采样点用于模型反演,其余4个采样点用于模型验证(如图5)。

根据三波段模型的原理,模型的概念模式可表示为:

式中Rrs(λ1),Rrs(λ2),Rrs(λ3)分别表示波段λ1、λ2、λ3的反射率。建模过程中,确定出最优的λ1、λ2、λ3是模型反演的重点和难点。其中λ1要突出叶绿素α的吸收作用,同时满足黄色物质、悬浮物的吸收系数和水体的后向散射系数贡献相对较小,此时叶绿素α的吸收作用对水体的反射率影响最大。叶绿素α存在两个吸收峰,分别为440 nm和670 nm左右,由于黄色物质、悬浮物的吸收系数和水体的后向散射系数随波长的增加而减小,在652 nm以后这三个值都较小[14]。因此,确定λ1的波段范围为650~680 nm;λ2波段将黄质和悬浮物的吸收最小化,但同时要保证叶绿素α的吸收系数相对较小,即λ2应该位于700 nm附近叶绿素α的荧光峰处[15]。故确定λ2的波段范围为680~710 nm;λ3波段要去除后向散射的影响,且叶绿素α、黄色物质、悬浮物的吸收作用对总吸收系数影响甚微,即遥感反射率主要受纯水的吸收特性影响。由于纯水在750—760nm处吸收作用最大[14]。因此,确定λ3的波段范围为710~760 nm。对应环境一号HSI高光谱数据,λ1、λ2、λ3的波段范围分别为HSI66~HSI73、HSI74~HSI79、HSI80~HSI88。

为了确定λ1、λ2、λ3的准确波段,采用逐波段迭代的方法。首先从λ1开始,将λ2、λ3分别固定在HSI77、HSI83,在λ1范围内进行迭代,将迭代结果中R2最大的波段作为λ1的输出波段。将此输出波段加入到λ2的迭代中(仍固定λ3为HSI83),完成λ2的迭代。再将λ1、λ2的迭代结果加入λ3的迭代中,完成λ3的迭代,直到出现相同的迭代结果便停止(如表1)。

由表1可以看出,三波段模型的最优三个波段为HSI73、HSI74、HSI87。根据建模原理,以(1/HSI73-1/HSI74)×HSI87为自变量,叶绿素α浓度为因变量建立线性模型

式中,HSI73、HSI74、HSI87分别为HSI高光谱数据73、74、87波段的遥感反射率,CChla为叶绿素α的浓度值,单位mg/L。

将反演模型应用于研究区的HSI高光谱数据,得到整个研究区内的叶绿素α浓度分布,与4个验证点数据进行对比分析。结果表明,模型估测的最大相对误差为13%,最小相对误差为2%,平均相对误差为7%。RMSE为0.002 9 mg/L,远低于叶绿素α浓度的平均值0.050 3 mg/L。可见,基于环境一号HSI高光谱数据的三波段模型可用于太湖叶绿素α的遥感定量反演。

3 模型对比分析

为了进一步确定以上两类模型的适用性和模拟精度,以图5显示的灰色区域为研究区,将区域内13个实测值与两类模型的模拟值进行对比分析(如图7)。通过计算,环境一号CCD数据叶绿素α反演经验模型模拟的RMSE为0.003 5 mg/L,环境一号HSI数据叶绿素α反演三波段模型模拟的RMSE为0.002 9 mg/L。由此可见,基于HSI高光谱数据的叶绿素α遥感定量反演模型模拟精度更高。

4 结论

选取典型的内陆水体——太湖为实验区,分别基于环境一号CCD多光谱数据和环境一号HSI高光谱数据建立叶绿素α浓度的遥感定量反演模型。通过验证,两类模型均可用于太湖水体的叶绿素α浓度反演,模型拟合度R2分别为0.865 3和0.9366。并通过两类模型对比实验,实验结果表明,基于HSI高光谱数据的三波段模型模拟估测的RMSE更小,精度更高。可见,由于三波段模型综合考虑了水体各组分的光谱特征,模拟精度能够大幅度提高,是内陆水体叶绿素α反演的有效方法。本研究为内陆水体叶绿素α浓度定量反演提供了研究基础。

摘要：叶绿素α是反映内陆水体富营养化程度的一个重要参数,通过遥感定量反演模型提取叶绿素α浓度是内陆水体富营养化监测的重要途径。以典型的内陆水体——太湖为实验区,基于环境一号CCD多光谱数据和同步地面实测叶绿素α浓度数据,建立了以CCD4/CCD2为自变量的叶绿素α浓度反演经验模型。基于环境一号HSI高光谱数据,结合内陆水体中叶绿素α、悬浮物、黄色物质和纯水的固有光学特性以及三波段模型理论,选取了反演叶绿素α浓度的三个最优波段:HSI73、HSI74、HSI87,并建立三波段模型。最后通过两种定量反演模型对比实验。实验结果表明基于HSI高光谱数据的三波段模型反演精度更高,可以有效的用于内陆水体的叶绿素α浓度反演。研究为内陆水体叶绿素α浓度定量反演提供了研究基础。

基于模型反演 篇8

我国多数油田陆相储层非均质性较强。在油田投入开发初期钻井资料少、井距大的条件下,如何在地质建模过程中减少井间储层物性的不确定性成为模型准确程度的关键。应用地震或反演资料对储层地质模型进行井间趋势约束是一种常用的解决办法,这种方法可以充分利用地震或反演资料横向信息丰富、测井资料纵向分辨率高的特点,实现跨学科对储层物性分布进行描述,使得储层物性不仅在井周围纵向上有高分辨率与准确性,同时又能反映出地震或反演资料在平面上大尺度的变化与连续性,确保地质模型能够精准地描述出储层物性的空间变化,为油田的进一步开发提供依据[1]。目前地震或反演资料对地质模型的约束作用主要体现在以下几个方面:

(1)对于构造模型,采用地震解释的深度域构造面对井间构造面起伏进行趋势约束,使构造层面严格通过井上分层点,同时尊重地震资料解释成果,确保构造精确;

(2)对于属性模型,通过数据分析获得地震资料或反演资料与沉积相(微相)或砂泥岩相的相关性,再通过相模型以“相控”的方式实现地震资料或反演资料对属性模型的间接约束[2]。这种方法面临的问题在于:首先,把建立属性模型的过程分为“两步走”,对于最终物性模型的准确度,在实际工作中很可能会出现误差“一加一大于二”的结果;其次,当储层内不同沉积微相之间物性差异小的情况下(如水下分流河道和河口坝),该方法实现效果不理想;

(3)对于属性模型,也可以通过分析反演属性(波阻抗或密度)与储层物性的相关性,用反演资料直接对物性模型进行约束[3]。这种方法面临的问题在于反演属性与储层物性的相关性难以用函数精确地表述,即二者真实的相关关系与定义的相关函数之间不可避免地存在误差。

基于上述原因,本文创新使用遗传反演方法获得孔隙度的反演数据体,并将其直接用来对孔隙度模型进行约束。由于减少了建模过程中部分中间环节,减少了误差出现,使模型能够更加准确地反映地震信息,结果更加精准,并缩短了工作周期,在实际工作中取得了良好的效果。

1遗传反演的基本原理

遗传反演方法是斯伦贝谢公司在其Petrel软件中提出的一种新的反演方法,该方法是基于神经网络与遗传算法合并的一种特殊算法[4],用以得到精确直观的地震反演结果。其具体实现过程分为神经网络训练与权重更新两部分。首先将井数据(如孔隙度)和地震振幅资料作为训练数据,通过神经网络建立二者的非线性相关关系,把通过神经网络迭代后产生的孔隙度计算结果与真实的孔隙度数据通过计算误差函数的方法进行对比。若计算结果不满足收敛条件,则通过基因算法对不同训练数据的权重进行更新后再进入神经网络进行迭代,直到计算结果满足收敛条件,产生出理想的孔隙度数据体为止[5]。

这种方法与传统反演方法的区别与优势有以下几点:

(1)对于遗传反演,需要的输入数据仅限于地震振幅体,及测井曲线数据(如孔隙度)作为训练数据。既不需要子波,也不需要初始属性模型作为运行反演的优先条件。

(2)在后台的遗传算法计算误差以得到用于神经网络的权重系数,遗传算法约束反演收敛性即达到全局最小误差值的机会远远大于以往的基于反演的神经网络算法,因此精确度是显而易见的[6]。

(3)它不限于常规的波阻抗反演或者密度反演,遗传反演可以应用于任何类型的岩石物理属性。更直接地说,所有包含在波动方程中的参数都可以进行反演,例如,速度、密度、孔隙度、反射系数等[7]。本次研究就是基于遗传反演的这种优势,对孔隙度进行了反演,进而对孔隙度模型进行约束。

2实例分析

本次研究以某油田东营组内某一小层为例,对遗传反演及其对储层地质模型的约束作用进行论述。

2.1地质概况

研究区为某油田边部的井区,目的层为东营组三角洲相储层内某一小层。该井区所处部位构造相对简单,无断层发育。小层平均厚度14 m,平均有效厚度9 m,平均有效孔隙度0.28。小层内部发育一条较稳定的夹层。目前该区已钻井37口,建模区域面积2.2 km2,平均井距400 m。

由于目前钻井资料少,井间储层变化难以用数据分析的方式准确获得。在储层非均质性较强的情况下,储层描述工作要强烈地依赖于地震或反演资料开展,因此反演资料的质量对地质模型有巨大影响。

2.2孔隙度遗传反演

应用遗传反演方法,对目标区进行孔隙度反演。为了验证反演准确程度,采用工区范围内各处任意分布的30口井作为源数据,其余不参与运算,用于检验反演效果。从井上效果看,井旁反演道与测井解释孔隙度曲线有明显的对应关系,可以良好地体现出井上单砂体与泥岩夹层的分布(图1),说明遗传反演得到的孔隙度数据体可以很好地表明储层内物性变化,适合用作对孔隙度模型进行趋势约束[8]。相比之下,在目前的大井距条件下,沉积微相图精细程度不足,用沉积相约束储层物性势必造成较大误差。

2.3 孔隙度反演体约束储层地质建模

在遗传反演孔隙度数据体的约束下,进行应用序贯高斯模拟算法对目的小层产生了多个实现(图2)。从图2中可以看出,尽管建立孔隙度模型采用的是序贯高斯随机模拟这种随机建模方法[9],但各个实现之间具有很强的相似性。说明通过反演约束,可以有效降低孔隙度模拟的多解性和井间储层物性的不确定性。在工区南部,孔隙度高值区域沿方位角-45°呈条带状分布,与目前对物源方向与水下分流河道分布的认识大体一致。说明反演约束下的储层孔隙度模型很好地反映出了目前对该井区的地质认识,并且使测井与反演资料在模型中得到了良好的体现。模型同时受到井、震两种数据的约束,既符合井数据的地质统计学特征,同时又能反映出反演数据观测到的大尺度结果和储层连续性,提高了模型的精度。

Por是孔隙度,孔隙度是指岩石中孔隙体积(或岩石中未被固体物质充填的空间体积)与岩石总体积的比值。用百分数或小数表示。图3同。相比之下,在没有反演资料约束的情况下,仅用井数据以相同的序贯高斯模拟算法建立孔隙度模型的几个实现如图3所示。与前者相比,井间孔隙度分布在局部差异较大[10]。这表明在大井距情况下,未经反演约束,尽管用井数据建立的每一个实现都满足井数据的地质统计学特征,但在井间储层的预测还有较大的随机性和不确定性。

3 结论

(1) 应用遗传反演方法得到的孔隙度反演数据体可以在三维空间中较真实地反映出储层物性与已知地质认识,适合对孔隙度模型进行趋势约束;

(2) 采用反演约束的孔隙度模型同时受到井、震两种数据的约束,既符合井数据的地质统计学特征,同时又能反映出反演数据观测到的大尺度结构和储层连续性,有效降低了地质模型的不确定性,能够提高储层孔隙度模型的精度。

(3) 利用孔隙度遗传反演数据体对孔隙度模型进行约束的方法,相比传统地震约束方法,中间步骤减少,准确性更高,工作周期更短,在实际工作中取得了良好的效果。

摘要：针对油田投产初期井距较大、井间储层发育状况难以预测的情况,应用遗传反演方法进行孔隙度反演。在地质建模过程中利用反演数据体约束孔隙度属性。以某油田东营组某小层为例进行地质建模。结果表明:应用遗传反演方法得到的孔隙度反演数据体可以在三维空间中较真实地反映出储层物性。利用孔隙度遗传反演数据体对对孔隙度模型进行约束的方法,相比传统地震约束方法,中间步骤少,准确性高,工作周期短,在实际工作中取得较好效果。

关键词：地质建模,遗传反演,孔隙度模型,反演约束

参考文献

[1]刘文岭.地震约束储层地质建模技术.石油学报,2008;29(1):64—74

[2]刘显太,昌峰,张世明.东海油气田A油藏跟踪研究及效果分析.西南石油学报,2003;25(2):15—18

[3] Seifert D,Jensen J L.Using sequential indicator simulation as a tool inreservoir description:Issues and uncertainties.Mathematical Geology,1999;29(5):527—550

[4]周明辉.储层地质模型的建立及动态实时跟踪研究.东营:中国石油大学(华东),2009

[5]张伟,林承焰,周明辉,等.地质模型动态更新方法在关家堡油田的应用.石油勘探与开发,2010;37(2):220—225

[6]张英波.地震孔隙率反演方法和应用.石油地球物理勘探,1994;29(3):261—273

[7]任殿星,李凡华,李保柱.多条件约束油藏地质建模技术.石油勘探与开发,2008;35(2):205—214

[8] Haldorsen H H,Damsleth E.Stochastic modeling.Journal of Petrole-um Geology,1990;42(4):404—412

[9] Damsleth E,More H,Haldorsen H,et al.A two-stage sto-chasticmodel applied to North Seareservoir.Journal of Petroleum Geology,1992;44(4):402—408

基于模型反演 篇9

岩石的蠕变研究对于合理评价岩体的长期稳定性，即岩石力学行为的时间效应具有重要意义。岩石蠕变研究一般是从试验出发获取有关数据，然后构建岩石蠕变的本构模型并确定相应的模型参数。因此，根据蠕变试验资料，如何正确确定岩石蠕变模型参数是岩石蠕变研究领域的关键课题之一。目前，确定岩石蠕变模型参数的方法主要有最小二乘法[1]、多项式回归法[2]、流变曲线分解法[3]、复合形优化法[4]、模式搜索算法[5]等。但岩石蠕变模型的参数众多，为得到问题的全局最优解，须尝试应用新的优化算法。近些年来，随着优化计算技术的发展，人工神经网络[6]、遗传算法（Genetic Algorithm,GA）[7]、扩展卡尔曼滤波技术[8]、粒子群优化（Particle Swarm Optimization,PSO）算法[9]等逐渐在岩石蠕变模型参数确定中得到了应用，为解决问题提供了一种新的思路。

基于算法简单性和寻优精度方面的考虑，本文给出了岩石蠕变模型非定常参数智能反演的一种新算法——微进化算法（Microevolution Algorithm,MA），并利用页岩的蠕变试验数据进行了参数反演分析。计算结果表明，微进化算法作为蠕变模型参数反演的新方法是可行有效的，且算法计算精度高于混沌粒子群优化算法，具有较高的工程应用价值。

1 微进化算法基本原理

微进化算法，是受人类社会发展进程中人类趋同和趋异学习行为过程的启发构造出来的。它所采用的基本搜索策略是，群体中的每一个个体向群体中的优秀个体进行学习，利用优秀个体积累的经验知识来改变自身，使自身得以处于不断运动之中。具体来说，MA基于实数编码，采用群体社会学习机制，对种群中的每个个体i，以其自身所处历史最优位置为基础，以群体最佳位置与当前个体i的矢量差异为信息参照，进行动态随机搜索，以实现种群的进化。

为便于描述，设Xi=(xi1,xi2,xiD)为个体i=,1,2,…Np的D维矢量，每个Xi代表一个潜在的解；设Xpbesti为个体i迄今为止搜索到的最优状态、Xgbest为整个群体中的所有个体迄今为止搜索到的最优位置状态。这样，在每一次的迭代搜索中，以式（1）更新各个体所处状态：

式中，iter为迭代次数；r为0～1之间的正实数，一般为了使用方便，可令r∈U(0,1)；N(,0σ)为正态分布随机数，其中σ∈(0.5,2.5)。

分析式（1）可知，当Ni(,0σ)0时，群体中的各个体i以矢量差(Xgbest-Xi,iter)为方向趋同；当Ni(,0σ)≺0时，各个体以矢量差(Xgbest-Xi,iter)的反方向趋异。可见，式（1）将趋同与趋异有机地结合了起来，以实现局部的开采与全局的探索进化功能。算法实现步骤如下：

步骤1（初始化）令iter=0，采用实数编码，在可行域空间随机初始化种群Xiter=(X,1iter,X2,iter,…Xi,iter,…XNp,iter)T，若搜索空间为D维，则每个个体Xi,iter中包含了D个变量，

即Xi,iter=(xi,1iter,…xi2,iter,…xij,iter,xi D,iter)，j=,1,2D。

步骤3（种群演化）执行如下操作，产生第iter+1代种群Xiter+1。

步骤4（终止检验）（1）判断算法循环执行次数是否达到最大进化代数ITERMAX；(2）判断最优目标函数值f(Xgbest)是否达到预设精度VTR。若满足终止条件其中之一，则结束算法，输出Xgbest及其目标函数值，否则，转入步骤2直至满足终止条件。

由算法的实现步骤可见，微进化算法采用的基本策略是：（1）采用实数编码；（2）使用固定规模的种群Np；(3）采用贪婪的种群更新策略，即用Np个新个体Xi,iter+1完全替换Np个旧个体Xi,iter；(4）采用精英保留策略，即更新个体自身所处历史最优状态Xpbesti和群体最佳状态Xgbest。

MA需要设置的算法控制参数较少，除算法的运行控制参数（ITERMAX、VTR）外，仅需设置群体数目Np和σ。过大的Np会影响算法的运算速度，经试验，对于一般问题，建议Np取20～40左右；而σ∈(0.5,2.5)，一般取1即可获得较满意解。可见，与GA和PSO算法相比，MA控制参数的选择更为简单，这使得算法易于编程实现和便于用户使用。

2 微进化算法在岩石蠕变模型非定常参数反演中的应用

为探讨本文的微进化算法优化反演岩石蠕变模型参数的可行性，利用文献[10]中的页岩蠕变试验数据，采用标准线性体三元件模型，并引入非定常参数后进行反演分析。

对于标准线性体三元件模型，蠕变试验时，其本构方程为：

在蠕变本构方程中引入非定常参数EK(t)=p1+p2ep3t，并令ηK=p4、EH=p5，式（2）变为[10]：

设式（3）等式左边和右边项分别记为fL和fR，则据式（3），由m对实测蠕变数据{tl,εl}（l=1,2,…m）可建立反演问题的优化准则函数为：

式中，X=(p1,p2,p3,p4,p5)为待反演参数。由于l=式1(4）具有高度的非线性特征，传统优化方法如高斯-牛顿或麦夸特法等，对输入的模型参数初始值有一定的要求，若输入的参数值与参数的真实值偏离较大，可能造成迭代不收敛，而且也很难得到全局优化解。而本文微进化算法对模型是否线性、连续等没有限制，也不受优化变量数目的束缚，直接在优化准则函数的引导下进行全局自适应寻优，该方法简便、通用、适应性强。因此，现利用MA进行反演计算。步骤如下：

步骤1：初始化Np组待反演的蠕变参数Xi=(p 1,p2,p3,p4,p5)，i=1…Np，作为初始种群。

步骤2：将m组蠕变试验数据和初始化得Np组待反演的蠕变参数带入已建立的反演问题的优化准则函数式（4）中，计算得到相应应力水平下的最小二乘误差。

步骤3：评价各组待反演蠕变参数Xi，更新Xgbest。

步骤4：利用式（1）策略进行微进化寻优操作，并重复执行步骤3，直至满足微进化算法的终止条件，然后输出待反演蠕变参数的结果以达到最优化的目的。

根据上述蠕变参数的优化反演步骤，可反演计算得到页岩的蠕变模型非定常参数，如表1所示。为便于与文献[9]提出的混沌粒子群优化（Chaos Particle Swarm Optimization,CPSO）算法进行比较，表1中也列出了CPSO算法反演的参数结果。比较可见，在各应力水平作用下，MA与CPSO算法反演的参数p1、p5结果差别不大，而参数p2、p3、p4差别较大，相差达一个数量级之多。由优化准则函数值min J结果可见，MA计算所得非定常参数的最小二乘误差远小于CPSO计算所得的最小二乘误差，MA可将优化准则函数值min J减少50%左右，这说明MA较CPSO算法在处理式（4）所示这种较为复杂的非线性模型时，具有更高的寻优精度。

为便于直观反映MA与CPSO算法反演不同应力水平各时刻的应变值，图1绘出了两种算法所得蠕变计算值与试验值结果。从图1可以看出，各应力水平下，MA计算所得理论反演变形曲线与试验变形曲线较为吻合，误差较小，说明本文提出的MA是有效可靠的；而在应力水平σc=90MPa、σc=300MPa时，CPSO算法计算所得的后期理论变形曲线偏离试验曲线相对较大。

注：表中CPSO算法计算结果直接来源于文献[10]。

3 结论

(1）以实测蠕变值与理论计算值之间的最小二乘误差为优化准则函数，用微进化算法对岩石蠕变模型非定常参数进行优化反演，能最大限度地利用了所有试验数据，且可避免传统优化算法初始参数选取的困难。算例计算结果表明，微进化算法不受初始选值的影响，计算精度高于混沌粒子群优化算法，具有较高的工程应用价值。

(2）本文虽以标准线性体三元件模型为例，进行了非定常参数的反演分析，但方法同样适用于其他蠕变模型，只需采用与模型对应的蠕变本构方程作为优化准则即可。

(3）微进化算法简单有效，且易于编程实现和便于用户使用。

参考文献

[1]王宏贵,巍丽敏,赫晓光.根据长期单向压缩试验结果确定三维流变模型参数[J].岩土工程学报,2006,28(5):669-673.

[2]许宏发,沉新万.多项式回归间接求解岩石流变力学参数的方法[J].有色金属,1994,46(4):19-22.

[3]BOUKHAROV G N,CHANDA M W,BOUKHAROV N G.The three processes of brittle crystalline rock creep[J].Int.J.Rock Mech.Min.Sci.and Geomech.Abstr.,1995,32(4):325-335.

[4]李云鹏,王芝银,丁秀丽.岩体原位流变荷载试验的力学参数与模型反演[J].实验力学,2005,20(2):297-303.

[5]陈炳瑞,冯夏庭,丁秀丽.基于模式搜索的岩石流变模型参数识别[J].岩石力学与工程学报,2005,24(2):207-211.

[6]陈炳瑞,冯夏庭,丁秀丽.基于模式–遗传–神经网络的流变参数反演[J].岩石力学与工程学报,2005,24(4):553-558.

[7]刘杰,王媛.岩体流变参数反演的加速遗传算法[J].大坝观测与土工测试,2001,256:33-37.

[8]杨成祥,冯夏庭,陈炳瑞.基于扩展卡尔曼滤波的岩石流变模型参数识别[J].岩石力学与工程学报,2007,26(4):754-761.

[9]李志敬,朱珍德,周伟华.基于CPSO算法的岩石蠕变模型非定常参数反演分析[J].河海大学学报(自然科学版),2008,36(3):346-349.

基于模型反演 篇10

传统的水稻穗颈瘟的监测主要基于田间取样、调查,综合其它信息进行预测预报,但对大面积病害的发生,传统的方法不但耗时、费力,而且预报滞后增加了损失程度,从而在一定程度上影响了预报的精确度。与传统病情测报方法相比,遥感监测技术是一种无损测试技术,具有快速、宏观、客观、大面积和无破坏等显著优点,同时可以提取信息制成各种专题报告等特点[2]。由于遥感是在高空尺度下对地面目标宏观状况的一种监测,所以它能在较早的时期发现病害,提供水稻穗颈瘟发生与发展的状况,并能生成病害发生的时空分布图,同时多时相的监测能给决策者提供重要病害蔓延的趋势,从而让决策者有针对性地进行决策、加强重点防治,减少产量损失[3]。

国内外利用遥感技术来监测农作物病害状况已有不少研究。在国外,Malthus等用地物光谱仪研究大豆和蚕豆斑点葡萄孢子感染后的反射光谱,发现其一阶反射率较高,可用来监测病虫害的感染发生情况[4]。Adams等利用大豆黄痿病光谱二阶导数设计的发黄指数对病情评价进行了研究[5]。Rinehart等利用可见/近红外分光镜对匍匐翦股颖的褐斑病和一年生牧草的硬币圆状斑病进行了研究,发现不同病情的冠层光谱反射率一阶导数在700nm、1400nm和1930nm处有着明显的特征[6]。在国内,黄木易等对感染条锈病的小麦光谱进行研究,发现630-687nm、740-890nm及976-1350nm为遥感监测条锈病的敏感波段[7]。蔡成静等研究表明在535nm以后,小麦单叶的病情严重度与反射率达到了显著相关;安虎等研究表明,单叶光谱与小麦病情指数相关性最好的波段为670-690nm[8];黄木易等利用单叶光谱的吸收面积指数AAI(absorption area index)反演小麦条锈病严重度[9];Huang等利用单叶光谱构建的归一化光化学指数NPRI来反演单叶的病情指数,线性方程的决定系数达到极显著水平[10]。上述研究主要基于作物单叶光谱进行的,也有研究者对病害胁迫下的冠层光谱特征进行研究。蔡成静等对感染条锈病的冬小麦冠层光谱特征进行研究,发现在769-938 nm的红外波段,小麦条锈病不同病情的反射率存在很明显而规律的差异,并用930nm处的冠层光谱反射率反演小麦条锈病严重度[11]。蒋金豹等利用微分光谱变量反演条锈病胁迫下小麦冠层叶片色素含量[12]。综上所述,以上研究都取得了较好的效果,但目前利用高光谱反演水稻病害严重度的研究报道较少。

基于以上分析,本文对不同严重度的水稻穗颈瘟进行光谱测定,研究病害严重度与特征光谱之间的定量关系,建立水稻穗颈瘟严重度的高光谱监测模型,研究结果将为水稻穗颈瘟病害的航空、航天遥感监测和预警提供理论依据。

1 材料与方法

1.1 试验材料

水稻品种为2008年浙江省优质早籼新品种28份,对照品种2份,感病对照品种1份,共31份。供试稻瘟病菌为致病力较强的多个菌株的混合孢子,孢子培养采用大麦粒培养基,并将产生的孢子吸附于滤纸上,经干燥后保存在干燥器内备用。

供试材料按常规方法栽培,待水稻稻穗完全抽出时,将穗部剪下,取穗颈5-6cm,放置于培养皿内,并用100mg/kg无毒苯骈咪唑作保绿剂,然后将混有1%粘着剂羧甲基纤维素的稻瘟病菌孢子悬浮液(孢子浓度为1×105个孢子/mL)涂于穗颈处,放入26℃恒温箱内保湿培养,10-12d后进行光谱采集和病情指数调查。

1.2 光谱测定

实验使用美国ASD(analytical spectral device)公司的便携式可见/近红外光谱仪(Handheld Field Spec),其光谱采样间隔为1.5nm,测定范围325-1075nm,扫描次数30次,分辨率3.5nm,探头视场角为25度。光源是配套的14.5V卤素灯。盛样容器用直径是120mm,高10mm的玻璃皿。光谱采集时,先将籽粒从穗上取下,平放在测量的玻璃皿中。对光谱仪进行过标准白板光谱校正之后,进行样本表面扫描。每个样本扫描重复2次,并通过计算机自动保存。

1.3 病情指数(DI)调查

在病症状开始出现后,在每个培养皿中分别调查其发病情况,严重度分为6个梯度,即:0、1.0、3.0、5.0、7.0和9.0,分别记录各严重度的水稻穗颈数,然后利用公式1计算病情指数。

式中n为每个培养皿中穗颈总数,x为各病情梯度的级值,f为各梯度的穗颈数。

1.4 数据分析

1.4.1 光谱数据预处理

光谱数据预处理采用平均平滑法,选用平滑窗口大小为7,进行MSC(multiplicative scatter correction)处理。由于光谱曲线在首端和末端有较大噪音,所以只取400-l000nm波段的光谱用于分析。

1.4.2 光谱反射率一阶导数变换

光谱求导技术被广泛应用于减弱光照强度差异、背景光谱以及仪器噪声对目标物光谱特征的影响。一般导数光谱采用差分法近似计算:

其中,λi表示波段i的波长值,ρλi表示位于波长λi处的光谱值。

1.4.3 模型精度检验

(1)决定系数评价即相关系数R的平方值,用来说明用预测值解释实测值变差的程度。其计算公式为:

其中,y和是病情指数的测定值和平均值,x和是光谱反射率的测定值和平均值。

(2)均方根差(RMSE)评价由单变量和多变量回归模型估计出的参数,其精度可用均方根差(Root Mean Square Error)评价。

式4中,γi和γ'i分别为实测值和理论值,n为样本数。RMSE值越小说明方程的精度越高。

2 结果与分析

2.1 水稻穗颈瘟的光谱特征

在400-700nm和920-1000nm范围内,水稻穗颈瘟病情指数高(DI=7.9)的光谱反射率高于健康穗颈(DI=0)的光谱反射率,而750-850nm范围内穗颈瘟病情指数高的光谱反射率低于健康的穗颈光谱反射率(图1A),这表明,前两个光谱区域和后一光谱区域与病情程度分别呈正和负的相关关系;在720-760nm的微分光谱范围内,穗颈瘟病情指数高的光谱反射率低于健康的穗颈的光谱反射率(图1B),表明光谱反射率与病情程度呈负相关关系。

2.2 光谱特征与水稻穗颈瘟病情指数的相关性分析

由图2可以看出,在400-700nm光谱范围内,病情指数与光谱反射率的相关性达到了极显著正相关(图2A)。对于微分光谱,在520-540nm和710-760nm光谱范围内,病情指数与光谱反射率的相关性达到了极显著负相关,在570-610nm光谱范围内,病情指数与光谱反射率的相关性达到了极显著正相关(图2B)。由于作物在可见-近红外波段呈现特有的光谱反射率特性,即色素吸收决定着可见光波段的光谱反射率,细胞结构决定近红外波段的光谱反射率。在可见光波段区域随DI的增加与特定范围光谱呈现显著的相关性,表明水稻穗颈瘟破坏了细胞的色素特征。

2.3 水稻穗颈瘟病情指数的高光谱模型与精度分析

根据光谱特征与穗颈瘟病情指数的相关性分析,685nm处光谱反射率与711nm处微分光谱与病情指数的相关系数分别达到0.856和-0.867,相关系数的绝对值均为最大值(图2)。因此,利用685nm处光谱反射率与711nm处微分光谱分别与穗颈瘟严重度建立模型,其模拟方程、RMSE值与相关系数如图3和表1所示。由表1可以看出,2个穗颈瘟严重度高光谱模型相关系数分别为0.924和0.928,达到极显著相关水平,RMSE均在0.35左右,说明可以较好预测穗颈瘟的严重度。

3 结论与讨论

水稻受到穗颈瘟侵染时,其穗颈光谱反射率在可见光区有着明显的特征差异,其光谱反射率和微分光谱随病情的加重呈明显的上升或下降趋势,这与叶片的光合色素含量改变有着紧密联系。穗颈瘟菌孢子落在水稻穗颈后形成侵染菌丝并长出吸胞,伸入附近细胞内,从穗颈中吸取养料和水分,叶绿素被破坏。细胞叶绿素与水分含量的下降使得穗颈的光合作用能力随之下降,在光谱曲线上表现为可见光区的反射率上升及近红外反射率下降的趋势,这是穗颈瘟遥感监测的内在机理。

虽然以上结果是对穗颈瘟的严重度进行高光谱预测和模拟,但是对其他种类的稻瘟病(如苗瘟、叶瘟、节瘟等)的严重度反演同样具有借鉴意义。同时,由于本文是在离体条件下进行光谱采集,受植物和环境因素干扰较少;而田间活体条件下的冠层光谱与航天、航空遥感获取的地物混合光谱较为接近,所以本文所建立的穗颈瘟严重度反演模型还需要在田间条件进行检验和完善。当然,为了进一步整体提高水稻穗颈瘟严重度遥感反演的研究和应用水平,还需要开展以下工作:①研究基于冠层高光谱的穗颈瘟严重度反演模型,以增强田间条件下的病害反演能力;②将神经网络、小波分析等方法应用到反演模型中,以提高反演模型的精度[13];③采用多光谱图像(包括卫星影像)来获取水稻穗颈瘟的病情信息,以增强水稻病害宏观监测的能力[11,14]。

摘要：通过人工诱发不同等级水稻穗颈瘟,测定染病稻穗光谱及其病情指数(DI)。分析了病穗反射光谱和一阶微分在绿光区、黄光区和近红外区反射特征差异,将病情指数及光谱数据进行相关分析,建立水稻穗颈瘟严重度的估测模型并进行可靠性检验。结果表明:750-850nm反射光谱及720-760nm一阶微分为水稻穗颈瘟的敏感波段;在400-700nm光谱范围内病情指数与光谱反射率的相关性达到了极显著正相关,对于微分光谱,在520-540nm和710-760nm光谱范围内病情指数与光谱反射率的相关性达到了极显著负相关,在570-610nm光谱范围内病情指数与光谱反射率的相关性达到了极显著正相关,并在此基础上建立了水稻穗颈瘟病情指数的高光谱反演模型,该研究为今后通过航空、航天遥感大面积监测水稻穗颈瘟提供了理论依据。