一种基于高光谱技术的致密砂岩储层露头孔隙度表征方法技术领域：
本发明属于高光谱技术在油气精细地质应用领域，涉及一种基于高光谱技术
的致密砂岩储层露头孔隙度表征方法，在油气勘探开发及地质应用中具有较好的
应用前景。
背景技术：
露头是地下储层的真实刻画，在建立露头数字模型的基础上，进行露头储层
物性参数的定量表征，可以为更准确有效的建立地下储层地质模型。孔隙度可以
反映岩石储存流体的能力，是评价储层的一个重要指标。目前，露头表面孔隙度
表征的常用方法是：在露头表面取若干岩心，测得岩心孔隙度，插值得到整个露
头表面的孔隙度数据，这种方法存在岩心采样位置受限，岩心资料少以及插值结
果不准确的问题，不能快速获得储层露头宏观、定量、精确的孔隙度数据。高光
谱技术能够快速准确地获得储层露头的光谱数据，使遥感发生了由宏观到微观探
测，定性解译到定量反演的质的飞跃，在对露头岩石光谱特征的定量分析和理解
的基础上，可被应用于露头矿物精细识别、岩性填图、地质环境信息反演等，为
地质家提供致一种研究密储层露头孔隙度精细与定量表征的新手段。
发明内容：
本发明的目的在于克服现有技术存在的缺点，寻求设计提供一种基于高光谱
技术的致密砂岩储层露头孔隙度表征方法，利用地面激光雷达技术和地面高光谱
成像技术采集野外露头信息，建立高精度数字露头表层模型，在数字露头模型基
础上，利用高光谱数据对孔隙度进行预测，实现致密砂岩储层露头孔隙度的宏观
表征。
为了实现上述目的，本发明实现致密砂岩储层露头孔隙度表征的具体过程
为：
(1)、野外露头数据采集及岩石样品采样：采用地面激光扫描仪和地面高
光谱成像仪分别获取露头表层三维点云和储层露头高光谱图像，并对典型的砂岩
进行采样得到岩石样品，岩石样品的数量根据实际需要确定；
(2)、室内建立孔隙度预测模型：对于步骤(1)采集的岩石样品分别测定
其孔隙度数据和光谱反射率数据，先提取反射率、反射率一阶导数和反射率二阶
导数三个光谱指标，然后计算孔隙度和各波段光谱指标的相关系数，确定出光谱
指标和波段范围，最后利用偏最小二乘方法建立孔隙度对光谱指标的预测模型；
(3)、高光谱图像处理及孔隙度预测：对步骤(1)的野外露头高光谱图像
进行预处理得到，并求取光谱反射率的一阶导数得到处理后的高光谱数据，再根
据步骤(2)中确定的波段范围和预测模型，对处理后的高光谱数据进行孔隙度
预测，得到孔隙度预测结果图；
(4)、三维数字露头表层模型建立及与预测结果图配准：对步骤(1)中地
面激光扫描仪获得的露头表层三维点云数据进行拼接，并采用基于最优趋势面构
建三角网的方法对不规则海量点云进行建模，得到三维数字露头表层模型；并手
动选择同名点，对步骤(3)中预处理过的高光谱图像与三维数字露头表层模型
进行配准，利用配准参数对步骤(3)中的孔隙度预测结果图进行配准，使三维
数字露头表层模型上的任一点都具有预测孔隙度信息，实现致密储层露头孔隙度
的自动定量表征。
本发明与现有技术相比，其工艺简单，操作方便，数据精确，能快速获得储
层露头宏观、定量、精确的孔隙度数据。
附图说明：
图1为本发明实施例所述孔隙度和三个光谱指标的相关系数随着波长的变
化趋势图。
图2为本发明实施例利用反射率建立的预测模型预测孔隙度与实测孔隙度
的相关关系曲线图。
图3为本发明实施例利用反射率一阶导数建立的预测模型预测孔隙度与实
测孔隙度的相关关系曲线图。
具体实施方式：
下面通过实施例并结合附图对本发明作进一步说明。
实施例：
本实施例实现致密砂岩储层露头孔隙度表征的具体过程为：
1.野外露头数据采集及岩石样品采样
先选用奥地利rigel-vz400作为地面激光扫描仪采集露头表层三维点云，扫描
距离为27米，扫描点间距为1mm；选用美国HeadwallPhotonics公司制造的
ExtendedVNIR作为地面高光谱成像仪，其光谱范围是为
600nm-1600nm，波段为213个，数据间隔为4.7nm，图像分辨率为320×256像
元；再对不同层位的典型砂岩进行采样，共采集17块砂岩样品；
2.室内建立孔隙度预测模型
(1)孔隙度数据测定：
将采集的17块砂岩样品取岩芯，岩芯直径为2.509-2.536cm，长度为
1.324-5.963cm，并对岩芯进行常规清洗和烘干预处理后采用岩石孔隙度测定岩芯
孔隙度，岩芯孔隙度数据的分布范围是3.345％-16.329％，平均值为10.1266％；
(2)室内光谱数据测量及预处理
采用美国ASD公司生产的ASDFieldSpec3非成像光谱仪测量室内光谱，所
采用光谱仪的波段范围是350-2500nm，数据间隔为1nm，每次测定前严格按照
操作规范，去除暗电流影响，进行标准白板定标；为使光谱数据具有代表性，对
每一个砂岩样品，测量20次取其算术平均值，得到该砂岩样品的反射光谱曲线；
再采用移动平均法对实测光谱进行去噪处理，即选取测定样本某一点前后反射光
谱曲线上一定范围测定它的平均值作为该点的值，计算公式
为: NEWR i = 1 2 k + 1 ( R i - k + R i - k + 1 + ... + R i + ... + R i + k ) , ]]>i＝k+1,2,3…n-k；式中，Ri-k，Ri-k+1，
Ri，Ri+k为距离波段i最近的2k+1个反射率值，NEWRi为波段i点处的计算反
射率值，n为波段数，本实施例中k＝2；然后确定三个光谱指标反射率NEWRi、
反射率一阶导数NEWRi’＝(NEWRi+Δi-NEWRi)/Δi(式中，NEWRi+Δi为距离波段i
点Δi处的反射率，Δi为光谱微分间隔)、反射率二阶导数
NEWRi”＝(NEWR′i+Δi-NEWRi')/Δi(NEWR′i+Δi为距离波段i点Δi处的反射率一阶导
数，Δi为光谱微分间隔)；
(3)光谱指标提取与波段选择
孔隙度和三个光谱指标的相关系数随着波长的变化趋势如图1所示，从图1
可以看出，孔隙度与光谱反射率在600-1600nm波段范围内呈正相关关系，并且
相关系数从0.48逐渐增加至0.70，说明砂岩样品对近红外光谱具有良好的光谱
响应特征；孔隙度与反射率一阶导数的相关系数曲线存在4个明显的波峰和波
谷，在波峰附近其相关系数大于孔隙度与反射率的相关系数，在波谷附近其相关
系数大于孔隙度与反射率的相关系数；孔隙度与反射率二阶导数的相关系数波动
较大，且大部分处在-0.5到0.5之间，因此舍弃这一光谱指标；以相关系数r＝0.5
为界，将与孔隙度相关系数大于0.5的波长保留，将与孔隙度相关系数小于0.5
的波长舍弃，对于反射率这一指标，选取的波段范围是669-1600nm，其中共198
个波长；对于反射率一阶导数，选取的波段是波峰附近的波长范围，即
739-797nm，895-1042nm，1112-1117nm，1130-1372nm，1413-1423nm，
1437-1600nm，其中共134个波长；
(4)采用偏最小二乘方法建立孔隙度预测模型
采用偏最小二乘方法(PLS)分别建立反射率、反射率一阶导数对孔隙度的
多元线性回归预测模型，PLS能解决自变量之间的多重相关性和样本容量少的问
题，其基本原理如下：
设有p个自变量{x1,x2,…xp}和q个因变量{y1,y2,…yq}，为研究因变量和自变
量的统计关系，观测n个样本点，由此构成自变量与因变量的数据表
X＝{x1,x2,…xp}n×p和Y＝{y1,y2,…yq}n×q，X经标准化处理后的数据矩阵记为
E0＝{E01,E02,…,E0p}n×p,Y经标准化处理后的数据矩阵记为F0＝{F01,F02,…,F0q}n×q，
记t1是E0的第一个成分，t1＝E0w1，w1是E0的第一个轴，它是一个单位向量，
即||w1||＝1，记u1是F0的第一个成分，u1＝F0c1，c1是F0的第一个轴，并且||c1||＝1，
由于本实施例中F0是单因变量，所以c1是一个常数，则c1＝1，即有u1＝F0；
如果要t1，u1能分别很好地代表X与Y中的数据变异信息，根据主成分原
理，应该有Var(t1)→max，Var(u1)→max；由于回归建模的需要，又要求t1与
u1的有最大的解释能力，t1与u1的相关度应该达到最大值，即r(t1,u1)→max。
综上，在偏最小二乘回归中，要求t1与u1的协方差达到最大，即
C o v ( t 1 , u 1 ) = V a r ( t 1 ) V a r ( u 1 ) r ( t 1 , u 1 ) → m a x - - - ( 1 ) ]]>
正规的数学表述应该是求解下列优化问题，即
max w 1 , c 1 < E 0 w 1 , F 0 c 1 > s . t w 1 ′ w 1 = 1 , c 1 ′ c 1 = 1 - - - ( 2 ) ]]>
采用拉格朗日算法，记
s＝w'1E'0F0c1-λ1(w'1w1-1)-λ2(c'1c1-1)(3)
λ1，λ2为拉格朗日乘数，对s分别求关于w1，c1，λ1，λ2的偏导，并令其
为0，有
∂ s ∂ w 1 = E 0 ′ F 0 c 1 - 2 λ 1 w 1 = 0 - - - ( 4 ) ]]>
∂ s ∂ c 1 = F 0 ′ E 0 w 1 - 2 λ 2 c 1 = 0 - - - ( 5 ) ]]>
∂ s ∂ λ 1 = - ( w 1 ′ w 1 - 1 ) = 0 - - - ( 6 ) ]]>
∂ s ∂ λ 2 = - ( c 1 ′ c 1 - 1 ) = 0 - - - ( 7 ) ]]>
由此推出
2λ1＝2λ2＝w'1E'0F0c1＝<E0w1,F0c1>(8)
记θ1＝2λ1＝2λ2＝w'1E'0F0c1，所以θ正是最优化问题的目标函数值，式(4)
和式(5)可以写成
E'0F0c1＝θ1w1(9)
F0'E0w1＝θ1c1(10)
将式(9)代入式(10)有
E 0 ′ F 0 F 0 ′ E 0 w 1 = θ 1 2 w 1 - - - ( 11 ) ]]>
同理可得
F 0 ′ E 0 E 0 ′ F 0 c 1 = θ 1 2 c 1 - - - ( 12 ) ]]>
由此可见w1是矩阵E'0F0F0'E0的特征向量，对应的特征值为θ1是目标函
数值，它要求取最大值，所以w1是对应于E'0F0F0'E0矩阵最大特征值的单位特征
向量，c1是对应于矩阵F0'E0E'0F0最大特征值的单位特征向量；计算出w1和c1
后，可得t1＝E0w1，u1＝F0c1；然后分别求E0和F0对t1的回归方程：
E0＝t1p1'+E1(13)
F0＝t1r1'+F1(14)
式中E1和F1记为残差矩阵，回归系数向量是
p 1 = E 0 ′ t 1 | | t 1 | | 2 - - - ( 15 ) ]]>
r 1 = F 0 ′ t 1 | | t 1 | | 2 - - - ( 16 ) ]]>
E1＝E0-t1p1'(17)
F1＝F0-t1r1'(18)
用残差矩阵E1和F1取代E0和F0，求第二个轴w2和c2及第二个成分t2，u2，
有
t2＝E1w2(19)
u2＝F1c2(20)
实施E1和F1对t2的回归，有
E1＝t2p'2+E2(21)
F1＝t2r2'+F2(22)
式中E2和F2为残差矩阵，
依次类推，最后用交叉有效性确定成分th的提取个数，停止迭代，在得到m
个成分后，有F0关于th(1≤h≤m)的回归模型为
F0＝t1r1'+t2r2'+...+tmr'm+Fm(23)
式中r1，r2，…rm为回归系数向量，Fm为残差矩阵。
由于th均为E0的线性组合，又
t h = E h - 1 w h = E 0 Π j = 1 h - 1 ( I - w j p j ′ ) w h = E 0 w h * - - - ( 24 ) ]]>
记 w h * = Π j = 1 h - 1 ( I - w j p j ′ ) w h ]]>
所以 F 0 = E 0 w 1 * r 1 ′ + E 0 w 2 * r 2 ′ + ... + E 0 w m * r m ′ + F m ]]>
= E 0 ( w 1 * r 1 ′ + w 2 * r 2 ′ + ... + w m * r m ′ ) + F m - - - ( 25 ) ]]>
若记y*＝F0，其中是的第j个分量；
有y*对的回归方程
y * ^ = α 1 x 1 * + α 2 x 2 * + ... α p x p * - - - ( 26 ) ]]>
最后变换成y对原始变量x1,x2,…xp的回归方程；
再利用反射率建立的PLS模型提取了2个主成分，对自变量的累计解释率
为99.9％，对因变量的累计解释率为76.4％，并具有61％的预测能力，利用反射
率建立起的预测模型为：
y＝-4.30361-2.52953X16-2.50922X17-2.4854X18+…+2.58062X213；
式中X16，X17，…，X213分别为选定的第16，17，…213个波段对应的反射
率。
预测孔隙度与实测孔隙度的相关关系如图2所示；利用反射率一阶导数建立
的PLS模型提取了2个主成分，对自变量的累计解释率为79.9％，对因变量的累
计解释率为85.2％，并具有66.8％的预测能力。利用反射率一阶导数建立起的预
测模型为：
y＝-1.22118+1666.15X31+385.667X32+1506.04X33+…+3607.28X212
式中X31，X32，…，X212分别为选定的第31，32，…212个波段对应的反射率一
阶导数；预测孔隙度与实测孔隙度的相关关系如图3所示；由图2和图3可以看
出，利用反射率一阶导数建立的预测模型，主成分更能说明因变量的变化，具有
更高的预测能力，并且预测孔隙度与实测孔隙的相关系数更大，模型效果优于采
用反射率建立的模型，并且一阶导数处理能在一定程度上削弱背景影响，在野外
高光谱数据处理中更具有优势，因此，采用反射率一阶导数建立的预测模型进行
野外储层露头孔隙度预测；
3.高光谱图像处理及孔隙度预测
高光谱测量前先进行白板校正，并对遥感图像中的机内暗电流和零响应偏移
进行校正，在ENVI环境下，所有的图像通过基于地面控制点的正射校正算法的
有理多项式系数进行几何校正，通过暗像素法进行大气校正，并利用移动窗口平
均法进行去噪后对高光谱图像反射率进行一阶导数处理，在2(3)中选择的波段范
围内，利用反射率一阶导数预测模型y＝-1.22118+1666.15X31+385.667
X32+1506.04X33+…+3607.28X212进行孔隙度预测，式中X31，X32，…，X212分别
为选定的第31，32，…212个波段对应的反射率一阶导数；
4.数字露头表层模型建立及配准
对地面激光扫描仪采集的数据进行处理，利用地面激光扫描仪自带处理软
件，将多站扫描的点云数据进行拼接，形成完整的露头点云数据，并采用基于最
优趋势面构建三角网的方法对不规则海量点云进行建模，然后将所有点云向各个
方向进行投影，选择投影面积最大的方向的平面作为最佳趋势面，将所有点云投
影到此最佳趋势面上，在平面上建立Delaunay三角网，再把平面三角网通过高
程值还原到三维空间中，从而形成数字露头表层模型；为在露头模型上显示和定
量分析孔隙度数据，需将预测结果图与数字露头表层模型进行配准，首先将校正
后的高光谱图像与数字露头模型进行配准，选取至少8对同名点，计算高光谱相
机的内外方位元素，利用配准参数对孔隙度预测结果图进行配准，使数字露头模
型上的任一点都具有预测孔隙度信息。