基于数据分割与主成分分析的 LAI 遥感估算

第30卷第2期2011年4月

文章编号：1001－9014（2011）02－0124－07

红外与毫米波学报

J．Infrared Millim．Waves

Vol．30，No．2April ，2011

基于数据分割与主成分分析的LAI 遥感估算

1，21，22*2222

董莹莹，王纪华，李存军，杨贵军，宋晓宇，顾晓鹤，黄文江

（1．浙江大学农业遥感与信息技术应用研究所，浙江杭州310029；

2．国家农业信息化工程技术研究中心，北京100097）

摘要：针对叶面积指数（LAI ）经典统计反演模型存在估算效果不理想以及反演效率低等问题，提出了一种基于农学

物候的数据分割与主成分分析结合的遥感估算方法．综合了原始光谱和微分（或差分）光谱主成分信息作为自变量，融入了以农学物候为先验的数据分割思想，并引入了多尺度建模方式参与反演过程．以冬小麦为实验对象，进行数值模拟和比较分析．结果显示，该方法既能有效地提高整体估算精度，又能显著地改善数据饱和问题，且在全样本遍历时体现了稳定鲁棒性．

关键词：主成分分析（PCA ）；农学物候；数据分割；多尺度建模；叶面积指数（LAI ）中图分类号：TP79文献标识码：A

Estimating leaf area index from remote sensing data ：based on

data segmentation and principal component analysis

22

DONG Ying-Ying 1，，WANG Ji-Hua 1，，LI Cun-Jun 2*，YANG Gui-Jun 2，SONG Xiao-Yu 2，GU Xiao-He 2，HUANG Wen-Jiang 2

（1．Institute of Agricultural Remote Sensing ＆Information System Application ，Zhejiang University ，Hangzhou 310029，China ；

2．National Engineering Research Center for Information Technology in Agriculture ，Beijing 100097，China ）Abstract ：According to the unsatisfactory and lower efficiency of classical statistical models in leaf area index （LAI ）esti-based data segmentation and principal component analysis was mation ，a new inversion method combined with phenology-proposed in this paper．In the method ，principal components of spectral data and differential （or difference ）spectral data were chosen as independent variables ，and phenology-based data segmentation was integrated into data processing in order

to improve estimation accuracy．In addition ，multi-scale was involved in modeling．Winter wheat was selected as experi-mental object for numerical simulation and comparative analysis．Results not only showed high precision in whole estimation and effectively improved data saturation ，but also manifested stability and robustness under full scan．Key words ：principal component analysis （PCA ）；phenology ；data segmentation ；multi-scale modeling ；leaf area index （LAI ）

PACS ：41．20．-q

引言

LAI ）是一个重要叶面积指数（Leaf Area Index ，

的植物学参数，在农作物长势监测、产量估算等领域应用广泛

．目前遥感反演LAI 主要包括机理模型

［2］

和统计模型．机理模型物理意义明确、普适性好，

［1］

土壤、大气等因素影响，反演精度不高；人工神经网

络对训练样本要求甚高，各层节点数及模型参数选择具有不确定性；主成分分析能够有效压缩数据量、降低数据维数，可靠性及普适性适中．为解决统计模型存在的上述问题，众多研究者在构建新型植被指数、引入多样神经网络、加强主成分分析理论内涵等

［1，3 5］

．方面开展了许多工作，并初步用于LAI 估算

在此，从数据处理模式及建模方式角度出发，提出了

但参数偏多、算法复杂，其应用范围受到局限；统计

模型基于数据的统计特性建模，其经典方法有植被指数、人工神经网络、主成分分析等

［1］

．植被指数受一种基于农学物候的数据分割与主成分分析结合的

Received date ：2010-05-24，revised date ：2010-10-09收稿日期：2010-05-24，修回日期：2010-10-09

基金项目：国家自然科学基金项目（40701120）；国家863计划项目（2006AA120108）；北京市自然科学基金项目（4092016）；北京市科技新星计

划（2008B33）

），E-mail ：cici5201@163．com．作者简介：董莹莹（1984-女，山东淄博人，博士研究生，主要研究方向为农业遥感与信息技术，

*

mail ：licj@nercita．org．cn．通讯作者：E-

2期董莹莹等：基于数据分割与主成分分析的LAI 遥感估算

125

叶面积指数遥感估算方法．

以冬小麦为研究对象，通过在数据处理环节融合

［6，7］

，主成分分析与数据分割，在建模环节引入多尺度

以期实现LAI 的高精度估算和模型的稳定鲁棒性．旨

在通过综合原始光谱及其微分（或差分）光谱增加建

［8］

模信息量，利用主成分分析对数据进行挖掘、融合及重组，提高数据统计的可靠；通过基于以农学物候

以期改善拟为先验的数据分割子集进行多尺度建模，

合精度、抑制数据饱和，并为遥感大面积监测提供参

考依据．其中，单尺度建模是以原始数据或其主成分信息为自变量的反演，多尺度建模是以单尺度建模结果为自变量，并融合数据分割进行反演．

1卫星CCD 相机光谱响应曲线图1HJ-Fig．1Spectral response curve of resources satellite HJ-1

CCD cameras

1

1．1

材料与方法

研究区域与数据采集

实验于2009年在北京通州和顺义开展，两地区

谱而言，虽不存在建模变量多的问题，但由于数据本

身存在自相关性，于是选取主成分变换进行数据融合及重组，以期获得更为有效的光谱信息数据集．

用于数值实验的样本总量为48，其中通州和顺义地区样本数分别为28和20．现将总样本的2/3用作建模，剩余1/3用作检验．表1为样本总体（n =48）、建模样本（n =32）、检验样本（n =16）的统计特征，其中全遍历以及ʃ30% 和ʃ20% 分别代表全遍历样本以及统计参数分布在样本总体ʃ30% 和ʃ20% 之内的样本．

表1LAI 样本总体、建模样本、检验样本统计参数

Table 1Statistical parameters of LAI for all samples ，cali-bration samples and validation samples

样本

统计参数

最大值最小值算术均值几何均值中值标准差方差偏度峰度

样本总体7．491．214．013．633．941．692．860．24－0．90

建模样本7．401．283．993．623．911．682．840．25

7．411．294．013．643．921．702．900．24

7．421．264．003．623．931．712．910．24

7．111．534．043．693．971．672．820．21

检验样本7．311．354．003．633．951．712．930．24

7．311．304．033．653．941．712．910．24

全遍历ʃ 30%ʃ 20%全遍历ʃ 30%ʃ 20%

冬小麦采用正常管理策略，品种包括农大211、烟农19、京9428、京冬8等，在两个地区分别选择了14个和10个大田地块进行数据采集．2009年4月29日（挑旗期）和5月17日（灌浆初期）在田块中选择长进行冠层光谱测定，并同步获取植株势较均匀区域，

样本用于室内叶面积指数测定．

采用美国ASD 公司的ASD FieldSpec Pro FR （350 2500nm ）型光谱仪测定冬小麦冠层光谱，其采样间隔在350 1000nm 内为1．4nm ，在1000 2500nm 内为2nm．观测时探头垂直向下，其高度距地面1．3m ，探头为25ʎ 视场角，每个样点测量20次，每次测量前后用标准参考板校正．在各采集区内获取50cm ˑ 4行的植株样本，并带回实验室采用比叶重法测定叶面积指数．1．21．2．1

数值实验方法数据预处理

对地面高光谱而言，选取受水和二氧化碳强吸收带影响较小的350 1354nm 作为有效数据，进行滤波窗口为5nm 的滑动平均处理．由于未获得环境

1影像，减灾星HJ-故基于光谱响应曲线（见图1）进行数值模拟，将光谱响应曲线及地面高光谱反射率

作为向量，两向量点乘与光谱响应向量模的比值为光谱响应值．

经上述处理后，若直接将地面高光谱各波段作为自变量进行多元线性回归建模，会增加模型复杂度、降低反演效率．为此，引入了既可有效提取波谱信息，又可降低数据维数的主成分分析法（Principal Component Analysis ，PCA ）［1］．对模拟HJ-1

卫星多光

－0．91－1．01－1．01－1．05－0．96－0．96

1．2．2基于农学物候的数据分割与主成分分析结

合的多尺度建模

获取的冬小麦LAI 变化范围是1．21 7．49，较

大变化范围的产生是由于通州和顺义种植区较大，且品种多，因此作物生育期和长势存在差异，此类差

［9］

异表现为覆盖度和植株高度的不同．由于不同覆盖度会导致遥感冠层光谱的不同响应，且植株高度

［10］

也会带来二次散射和多次散射差异，因此引入了以农学物候为先验的数据分割处理．所谓数据分割

126

红外与毫米波学报30卷

（Data Segmentation ），指依据原数据集自身特性，将其合理分割成一组子集，并基于各子集进行曲面拟合及有机重组

［11］

．

结合北京地区冬小麦生长规律可知，一般正常生长的冬小麦在拔节前LAI ＜3，覆盖度较低，植株较矮；从拔节至挑旗3＜LAI ＜5，覆盖度高，植株长高；

［10］

挑旗后LAI ＞5，覆盖度和植株高度达到最大．基于选取3和5为数据分割阈值．之所以融入数据分此，割，是由于其既能削弱数据间的随机扰动，在一定程度上抑制数据饱和，又可提高数值拟合精度．在主成分分析和以农学物候为先验的数据分割

1卫星多基础上，分别以地面高光谱数据和模拟HJ-光谱数据为建模数据源，采用多尺度多元线性回归以原始光谱、一对冬小麦LAI 进行反演估算．首先，

阶微分（或差分）光谱、二阶微分（或差分）光谱主成

分信息作为自变量建模，得到三个LAI 反演模型；然后，基于数据分割子集，对上述反演结果再次建模，并将数据有机拼接，最终得到LAI 估算模型．其中微分（或差分）光谱的计算如下：地面高光谱波段众

1卫星可作为连续函数进行微分运算；模拟HJ-多，

可作为离散函数进行差分运算，多光谱仅四个波段，

其中一阶差分为光谱反射率之差与中心波段波长之

差的比值，二阶差分为一阶差分之差与中心波段波

［12］

长之差的比值．设计的LAI 估算实验方案如图2所示，主要包括：

图2数值实验方案

Fig．2Plan of numerical experiments

精度和抗数据饱和性方面进行文中方法与四类植被

SAVI 、ARVI 、EVI ［1］的分析对比．此四指数即NDVI 、

类植被指数的具体介绍见表3．

表2三类模型检验指标

Table 2Three kinds of assessment indexes

全称均方根误差

缩写RMSE

RMSE =

公式

d i 为模拟值与实测值之差．量次数，

∑i =1n i =1

d i /（n －1），其中n 为测

R 2=

决定系数

R 2

[(∑（x i [(∑（x i

i =1n

－）2（y i －）

2

)]

，

2

（1）以地面高光谱为数据源，基于主成分分析和以农学物候为先验的数据分割进行多尺度多元线

性回归建模，并对模型进行分析检验．

（2）以模拟HJ-1卫星多光谱为数据源，基于主成分变换和以农学物候为先验的数据分割进行多尺度多元线性回归建模，并对模型进行分析检验．（3）所用方法与NDVI 、SAVI 、ARVI 、EVI 在总体模拟精度和抗数据饱和性方面的对比分析．1．2．3

模型检验

选取均方根误差、决定系数和准确度进行模型分析检验，其计算公式见表2．其中均方根误差反映了模拟值偏离实测值的程度，其值越小拟合精度越高；决定系数表示模拟值与实测值的密切程度，其值越接近1，拟合曲面参考价值越高；准确度是截距为0时模拟值和实测值之间线性回归方程的斜率，表征模拟值与实测值的一致程度，其值越接近1，准确度越高．当决定系数与准确度同时趋于1时，模拟值和实测值趋于吻合．

为验证方法的可行性和有效性，

将在整体模拟

－x ）

（y i ∑i =1

n

－y ）

)]

2

x 、y i 、y 分别为实测值、其中x i 、实测值均

n 为测量次数．值、模拟值、模拟值均值、

准确度

Accuracy

y =Accuracy ·x ，y 为实其中x 为模拟值，

测值．

表3四类植被指数

Table 3Four kinds of vegetation indexes

全称

缩写

公式NDVI =

ρNir －ρRed ρNir +ρRed

ρNir －ρRed

ρNir +ρRed +L

入选原由应用最为广泛

归一化

NDVI

植被指数土壤调整

SAVI

植被指数

SAVI =（1+L ）

降低土壤影响

L =0，1ARVI =

ρNir －ρRB

ρNir +ρRB

大气阻抗

ARVI

植被指数改进型土壤大气修正植被指数

ρRB =ρRed －γ（ρBlue －ρRed ），

γ=1EVI =

2（ρNir －ρRed ）

，

+c 1ρRed +c 2ρBlue +L

降低大气影响

EVI ρNir

c 1=6，c 2=7. 5，L =1

降低土壤和

大气影响

ρRed 、ρBlur 分别为近红外波段、红光波段和蓝光波段反射率，上述三个注：ρNir 、

630 690nm 和450 520nm．波段的光谱范围分别是760 900nm 、

2期董莹莹等：基于数据分割与主成分分析的LAI 遥感估算

127

为进行抗饱和性分析，需计算数据饱和点，于是ND-引入了双向建模法．该法分别以LAI 为自变量、NDVI 为自变量进行VI 为因变量和LAI 为因变量、多项式拟合

［12］

，两曲线交点所对应LAI 即为饱和

点，大于此值为饱和数据．之所以使用双向建模，是由于若LAI 与NDVI 之间存在确定的显式函数关系则两曲线应重合，但基于实测数据建模时不会出现此情况，故可采用该法定量地分析两者变化趋势转折的地方，即饱和点．图3为地面高光谱数据饱和点1卫星多光谱亦可采用此法确定对模拟HJ-示意图，饱和点．

图4地面高光谱及其一阶、二阶微分光谱前七个主成分贡

献率及累计贡献率

Fig．4Contribution and cumulative contribution of first seven

order and sec-principal components of ground hyperspectral ，first-ond-order differential Spectral

2

2．1

结果与分析

基于地面高光谱的数值实验

在以地面高光谱为数据源的模拟实验中，将原

始光谱及其一阶、二阶微分光谱主成分信息作为自应用多尺度多元线性回归估算LAI．首先提取变量，

高光谱及其微分数据的主成分；然后统计各主成分见图4．考虑到前贡献率及多个主成分累计贡献率，

64．9%、40．4%的四个主成分分别承载了99．5%、

信息，且再增加主成分对累计贡献率影响不大，故选取前四个主成分为自变量进行基于数据分割子集的多尺度建模

．

模拟结果见表4，其中M1/2/3/4/5为反演模Differential-1、Differ-型，其建模变量依次为Spectral 、ential-2、Spectral +Differential 、Spectral +Differenti-al +Data segmentation （含义见图2）．分析可知，对M2的模拟效果优于M1和三类样本选取策略而言，

M3，这是由于原始光谱混合了土壤信息，

二阶微分

在一定程度放大了噪声，而一阶微分既削弱了土壤

影响又不致过分放大噪声．M4综合了原始光谱与微M2、M3之间．分光谱信息，故其模拟效果介于M1、

2

对M5而言，其RMSE 趋于0的程度以及R 与Accu-racy 同时趋于1的程度均优于其它四个模型．实验M5在一定程度上降低了表明，对地面高光谱而言，

建模数据内部随机扰动的影响，提高了线性拟合精度，体现了较为显著的良好模拟效果．2．2基于模拟HJ-1卫星多光谱的数值实验1卫星多光谱为数据源，以模拟HJ-进行了与

2．1类似的数值实验．首先对多光谱及其差分数据进行主成分变换，然后将主成分信息作为建模变量参与基于数据分割子集的多尺度建模．此外，还给出了多光谱数据未经主成分变换的数值拟合结果进行分析对比．

M1/2/3/4为反演模模拟结果见表5，其中HJ-HJ +Data segmentation 、型，其建模变量依次为HJ 、

HJ +PCA 、HJ +PCA +Data segmentation （含义见HJ-M2图2）．分析可知，对三类样本选取策略而言，

M1，HJ-M4模拟效果好于HJ-M3，模拟效果好于HJ-这是由于基于数据分割子集的模拟可有效提高拟合

M3与HJ-M1拟合结果相当是由于主成分变精度．HJ-换仅对数据进行重组与融合，并未影响建模信息量．

HJ-M4的拟合结果略好于HJ-M2，是由于对本文获取的数据而言，经主成分变换的信息要比未经主成分变

1换的信息对数据分割略敏感．实验表明，对模拟HJ-

图3数据饱和点示意图

Fig．3Saturation point

卫星多光谱而言，文中提出的主成分变换结合数据分割的方法可以显著提高LAI 估算精度．

128

红外与毫米波学报30卷

表4以地面高光谱为数据源的数值实验结果

Table 4Results of numerical experiments based on ground hyperspectral

检验参数

模型编号

M1M2M3M4M5

全遍历0．990．870．890．900．69

均方根误差ʃ 30%1．000．890．900．910．69

ʃ 20%1．000．890．900．930．72

全遍历0．690．750．740．730．84

决定系数ʃ 30%0．710．750．760．740．84

ʃ 20%0．700．750．750．730．83

全遍历0．951．021．191．041．01

准确度ʃ 30%0．971．051．241．061．06

ʃ 20%0．981．071．231．061．08

1卫星多光谱为数据源的数值试验结果表5以模拟HJ-Table 5Results of numerical experiments based on HJ-1multispectral

检验参数

模型编号

HJ-M1HJ-M2HJ-M3HJ-M4

全遍历

1．020．751．020．70

均方根误差ʃ 30%

1．040．711．030．63

ʃ 20%1．050．

731．040．64

全遍历0．700．840．690．86

决定系数ʃ 30%0．710．850．700．87

ʃ 20%0．700．840．700．87

全遍历0．930．900．910．96

准确度ʃ 30%0．950．950．931．02

ʃ 20%0．950．940．931．00

2．3

结果分析与对比

1图5和图6分别是以地面高光谱和模拟HJ-

卫星多光谱为数据源，在全样本遍历时的各方法数

，“△”，“□”值模拟结果，图中符号“○”代表整体，“●”，“▲”，“■”模拟参数代表饱和数据模拟参

NDVI 、SAVI 、ARVI 、EVI 分别代表建模变量为数，

NDVI 、SAVI 、ARVI 、EVI 的反演模型．图7和图8分

其中符号别是实测值与各模型预测值的散点图，“☆”，“★”代表非饱和数据拟合点代表饱和数据拟

合点．

分析图5与图6可知，在整体模拟效果方面，文

R 2与Accuracy 同时趋于1的程中方法RMSE 最小，

度最好，拟合效果最佳．在图5中，文中方法相对于

2

RMSE 降低了其它模型R 提高了13．1 45．8%，

22．1 48．8%．在图6中，文中方法相对于其它模型R 2提高了2．46 50．5%，RMSE 降低了6．53

48．1%．在抗数据饱和性方面，文中方法的RMSE 最R 2与Accuracy 同时趋于1的程度最好，小，拟合效

2

文中方法相对于其它模型R 提果最佳．在图5中，

RMSE 降低了1．0 13．7%．在高了3．69 46．2%，

2

图6中，文中方法相对于其它模型R 提高了13．2

RMSE 降低了16．1 17．6%．55. 7%，

图7和图8以预测值和实测值散点图形式辅助

也直观体现了融合农学物候证实了上述分析结果，

的数据分割与主成分分析反演模型的可行性与有

效性．

图5以地面高光谱为数据源的各方法模拟效果对比图

Fig．5Comparison analysis of ground hyperspectral numerical experiments

2期董莹莹等：基于数据分割与主成分分析的LAI

遥感估算

129

对整体估算效果和抗数据饱和性而言，文

中提出的基于农学物候数据分割和主成分分析结合

的反演模型估算效果最佳，这是由于主成分分析在降低数据维数的同时对数据进行了挖掘、融合及重组，既能有效提取携带原始信息最多的主成分，又能在一定程度上去除背景噪声．此外，将冬小麦农学物候特性融入数据分割，并基于此进行子集划分和拟合重组，此举既有利于增强模型抗噪性，又有助于提高拟合精度．

图5与图6仅给出了全样本遍历时各模型的数值拟合结果，对统计参数分布在样本总体ʃ 30%和ʃ 20%之内的数据而言，文中方法拟合效果较之其

它模型虽有改善和提高，但效果并不明显．这是由于

采用上述样本选取策略时，建模样本与检验样本的方差与标准差均有所增大，虽两者的偏度相对样本总体均无变化，但建模样本的峰度偏离样本总体12．2%，故所建模型不能良好地反映样本总体的情况，拟合效果改善不甚显著．

3结论

在充分考虑主成分分析和以农学物候为先验的数据分割优点的基础上，通过构造结合了原始光谱和微分（或差分）光谱主成分信息的多尺度多元线性回归模型，对叶面积指数经典统计反演模型存在

130

红外与毫米波学报30

卷

1多光谱为数据源的预测值与实测值散点图图8以模拟HJ-Fig．8Results of HJ-1multispectral numerical experiments

的估算效果不理想、反演效率低和数据易饱和等问采用了主成分题进行了针对性处理．在模型构造上，

分析，故能有效进行数据挖掘、重组及融合，融入了

基于农学物候的数据分割，故有助于提高线性回归精度．理论分析和数值实践证实，对遥感多光谱和高光谱而言，提出的数据分割与主成分分析结合的估算方法不仅有效地提高了LAI 的整体估算精度，而

且在全样本遍历时充分且显著地抑制了数据饱和，体现了模型的稳定鲁棒性．

所提方法选用基于指定阈值的农学物候数据分

割来处理样本数据集，由于缺乏分割阈值的自适应选取能力，因此在进一步研究中可考虑采用基于自适应阈值的数据分割方法．此外，在样本选取策略上选用了全遍历和限制样本统计参数变化范围的方法，该法适用于样本量小的数值实验．当样本量较大上述方法会增加计算复杂度，且不利于数值实时，现，因此在后续研究中可对样本选取策略进行深入分析和探讨，以期获得具有普适性的样本选取策略．

REFERENCES

［1］WANG Ji-Hua ，ZHAO Chun-Jiang ，HUANG Wen-Jiang ，et

al．Basis and Application of Agriculture Quantitative Remote Sensing ［M ］．Beijing ：Science Press （王纪华，赵春江，黄文江，等．农业定量遥感基础与应用．北京：科学出版

2008．社），

［2］Verstraete M M ，Pinty B ，Myneni R B．Potential and limi-tations of information extraction on the terrestrial biosphere

from satellite remote sensing ［J ］．Remote Sensing of Envi-ronment ，1996，58（2）：201—214．［3］Rondeaux G ，Steven M ，Baret F．Optimization of soil-ad-justed vegetation indices ［J ］．Remote Sensing of Environ-1996，55（2）：95—107．ment ，

［4］YANG Fei ，ZHANG Bai ，SONG Kai-Shan ，et al．Compari-son of methods for estimating soybean leaf area index ［J ］．

Spectroscopy and Spectral Analysis （杨飞，张柏，宋开山，等．大豆叶面积指数的高光谱估算方法比较．光谱学与光谱

2008，28（12）：2951—2955．分析），

［5］WANG Xiu-Zhen ，HUANG Jing-Feng ，LI Yun-Mei ，et al．

The study on hyperspectral remote sensing estimation models about LAI of rice ［J ］．Journal of Remote Sensing （王秀珍，黄敬峰，李云梅，等．水稻叶面积指数的高光谱遥感估算

2004，8（1）：81—88．模型．遥感学报），

［6］WEN Cheng-Lin ，ZHOU Dong-Hua．Theroy and Application

of Multiscale ［M ］．Beijing ：Tsinghua University Press （文成林，周东华．多尺度估计理论及其应用．北京：清华大

2002．学出版社），

［7］Willsky A S．Multiresolution Markov models for signal and

image processing ［J ］．Proceedings of the IEEE ，2002，90（8）：1396—1458．［8］Somers B ，Delalieux S ，Verstraeten W W ，et al．Magnitude

and shape related feature integration in hyperspectral mix-J ］．IEEE ture analysis to monitor weeds in citrus orchards ［

Transactions on Geoscience and Remote Sensing ，2009，47（11）：3630—3642．［9］Growth Rhythm and Culture Technique of Wheat ［M ］．Bei-jing Academy of Agricultural Sciences．Beijing ：Beijing

Publishing House （北京市农业科学院．小麦生长规律与

1980．栽培技术．北京：北京出版社），

［10］WANG Jin-Di ，ZHANG Li-Xin ，LIU Qin-Huo ，et al．Re-M ］．（王锦地，mote Sensing Spectral Data Base in China ［

张立新，柳钦火，等．中国典型地物波谱知识库．北京：

2009．科学出版社），

［11］Dunham M H．Data Mining ：Introductory and Advanced

M ］．London ：Prentice Hall ，2003．Topics ［

［12］YI Da-Yi ，CHEN Dao-Qi．Introduction of Multiscale Nu-．Hangzhou ：Zhejiang University merical Analysis ［M ］

Press （易大义，陈道琦．数值分析引论．杭州：浙江大学

2003．出版社），