PhoonKK2006_GeoCongress会议_Modeling and Simulation of Stochastic Data

各学科重要国际会议目
录
标准化管理处编码[BBX968T-XBB8968-NNJ668-MM9N]
建筑学院重要国际学术会议一、A类会议
二、B类会议
土木水利学院土木工程系重要国际学术会议一、A类会议
二、B类会议
土木水利学院 建设管理系
重要国际学术会议
一、A 类会议
二、B类会议
土木水利学院水利水电工程系重要国际学术会议
一、A类会议
二、B类会议
环境科学与工程系重要国际学术会议一、A类会议
二、B类会议
机械工程系重要国际学术会议一、A类会议
二、B类会议
精仪系机械工程学科重要国际学术会议一、A类会议
二、B类会议
精仪系仪器科学与技术学科重要国际学术会议一、A类会议
二、B类会议
精仪系光学工程学科重要国际学术会议一、A类会议
二、B类会议
热能工程系重要国际学术会议一、A类会议
二、B类会议
汽车工程系重要国际学术会议一、A类会议
二、B类会议
工业工程系重要国际学术会议一、A类会议
二、B类会议
电机系重要国际学术会议一、A类会议
二、B类会议
电子工程系电子科学与技术一级学科重要国际学术会议汇总一、A类会议
二、B类会议。

第47卷第1期自动化学报Vol.47,No.1 2021年1月ACTA AUTOMATICA SINICA January,2021深度信念网络研究现状与展望王功明1,2乔俊飞1,2关丽娜1贾庆山3摘要深度信念网络(Deep belief network,DBN)是一种基于深度学习的生成模型,克服了传统梯度类学习算法在处理深层结构所面临的梯度消失问题,近几年来已成为深度学习领域的研究热点之一.基于分阶段学习的思想,人们设计了不同结构和学习算法的深度信念网络模型.本文在回顾总结深度信念网络的研究现状基础上,给出了其发展趋势.首先,给出深度信念网络的基本模型结构以及其标准的学习框架,并分析了深度信念网络与其他深度结构的关系与区别;其次,回顾总结深度信念网络研究现状,基于标准模型分析不同深度信念网络结构的性能;第三,给出深度信念网络的不同无监督预训练和有监督调优算法,并分析其性能;最后,给出深度信念网络今后的发展趋势以及未来值得研究的方向.关键词深度信念网络,深度学习,无监督预训练,有监督调优,结构设计引用格式王功明,乔俊飞,关丽娜,贾庆山.深度信念网络研究现状与展望.自动化学报,2021,47(1):35−49DOI10.16383/j.aas.c190102Review and Prospect on Deep Belief NetworkWANG Gong-Ming1,2QIAO Jun-Fei1,2GUAN Li-Na1JIA Qing-Shan3Abstract Deep belief network(DBN)is a generative model based on deep learning and overcomes vanishing gradient problem resulted from traditional gradient-based algorithm when it comes to deep architecture,and it is one of hot issues in thefield of deep learning.Based on the idea of learning in stages,DBN models with different structures and learning algorithms have been proposed.The aim of this paper is to summarize the current research on DBN and gives some views about its development trends in the future.First,the basic structure and standard learning framework of DBN are introduced,the relationship and difference between DBN and other deep structures are analyzed.Second,the current research on DBN is given,the performances of DBN with different structures are analyzed based on standard the model. Thirdly,the different unsupervised pre-training and supervisedfine-tuning of DBN are given,and their performances are also analyzed.Finally,some views about DBN s development trends in the future are presented.Key words Deep belief network(DBN),deep learning,unsupervised pre-training,supervisedfine-tuning,structure designCitation Wang Gong-Ming,Qiao Jun-Fei,Guan Li-Na,Jia Qing-Shan.Review and prospect on deep belief network. Acta Automatica Sinica,2021,47(1):35−49人工神经网络是计算机模拟人类大脑处理信息的一种运算模式,即通过训练输入和输出数据,使网络得到关于输入和输出的非线性映射关系,从而在未来的任务中进行自主计算.因此,人工神经网络是计算机科学、认知科学、脑科学和数学的交叉学科,其在模式识别、智能控制、多种信号处理、优化设计等领域得到较大的发展,并已在信息产业中得收稿日期2019-02-25录用日期2019-05-19Manuscript received February25,2019;accepted May19,2019国家自然科学基金(61533002)资助Supported by National Natural Science Foundation of China (61533002)本文责任编委张敏灵Recommended by Associate Editor ZHANG Min-Ling1.北京工业大学信息学部北京1001242.计算智能与智能系统北京市重点实验室北京1001243.清华大学自动化系智能与网络化系统研究中心北京1000841.Faculty of Information Technology,Beijing University of Technology,Beijing1001242.Beijing Key Laboratory of Computational Intelligence and Intelligent System,Beijing 1001243.Center for Intelligent and Networked Systems,De-partment of Automation,Tsinghua University,Beijing100084到了成功的应用[1−5].20世纪八十年代末期,用于人工神经网络的反向传播算法(Back-propagation, BP)的发明,给机器学习带来了希望,掀起了基于统计模型的机器学习热潮.这个时候的人工神经网络虽然也被称作多层感知器(Multi-layer perceptron, MLP),但实际上是一种只含有一个隐含层的浅层人工神经网络模型.进入21世纪以后,随着互联网的高速发展,对大数据的智能化分析和预测提出了巨大需求.由于浅层网络往往采用梯度类学习算法,人为经验因素较多,缺乏自主学习过程且对初始参数的设定依赖性较强[6−8],这限制了神经网络的特征自动提取能力,使得其在处理大规模不确定性数据时往往误差较大.生物神经系统学研究结果表明,人类的智能主要取决于大脑皮层,而大脑皮层是一个大规模互连的深层生物神经网络[9−11],主要认知方式是无监督自主学习与推理.探求大脑的组织结构和运行机制,从模仿人脑深层学习机制的角度出发,36自动化学报47卷寻求新的信息处理方法是当前人工智能领域发展的优先方向.然而,由于理论分析的难度,加上训练方法需要很多经验和技巧,所以这个时期深层人工神经网络相对较为沉寂.2006年,加拿大多伦多大学教授、机器学习领域泰斗―Geoffrey Hinton和他的学生Ruslan Salakhutdinov在顶尖学术刊物《Science》上发表了一篇文章,开启了深度学习(Deep learning,DL)在学术界和工业界的浪潮[12−14].主要思想是利用“逐层初始化(Layer-wise pre-training)”来完成自主学习与推理过程,从而有效克服深层结构的训练困难.近几年来,深度学习凭借其模拟人脑分层学习和自主推理的认知机理逐渐成为研究热点[15],同时也带动了人工神经网络领域的进一步发展.由于深度学习能够在大量数据任务中快速稳定地计算,这推动了云计算、大数据科学的发展,如今已经在自然语义理解、模式识别问题、机器人学和数据挖掘等方面得到了较好的应用[16−19],甚至在机器情感分析方面也开始被研究,使得该领域朝着图灵机的实现又迈进了一大步.2016年,利用深度学习技术训练过的阿尔法围棋(AlphaGo)击败人类围棋冠军,引起了学术界和科技界的巨大轰动,并激起了人们对深度学习研究的再一次热潮.目前,深度信念网络(Deep belief network, DBN)是深度学习的主要实现方法之一.DBN是具有若干潜变量层的生成模型.潜变量通常是二值的,而可见单元可以是二值或实数[20−21].尽管构造连接比较稀疏的DBN是可能的,但在一般的模型中,每层的每个单元连接到每个相邻层中的每个单元,而层内没有连接.DBN可以通过若干个受限玻尔兹曼机(Restricted Boltzmann machine,RBM)的顺序堆叠来构造,其学习过程分为两个阶段,即首先对RBM进行逐层无监督预训练,再用反向传播算法对整个网络进行有监督的调优.DBN的这种分阶段训练方法使其在学习深层结构上取得了一定的成功,并在图像处理、模式识别、系统建模和预测等任务中得到了关注和研究[20,22−27].近年来,众多学者在现有DBN结构和学习算法的基础上进行了拓展与改进,并提出了多种类型的DBN变种模型.目前,比较常见的DBN变种模型主要有稀疏DBN[28−29]、自组织DBN[26]、增量式DBN[27]、递归DBN[30].与传统的DBN相比,改进型的DBN分别在各自的聚焦点上取得了部分性能上的提升.但是,在结构自主确定方面,DBN仍然存在一些难以解决的瓶颈问题,相关的研究工作还处于刚刚起步状态,在理论、技术以及应用层面上还有很大的提升空间,在未来一段时间内仍将是深度学习研究中比较热门的研究方向之一.1深度信念网络基本模型与概述深度信念网络是为了简化逻辑斯蒂信念网络的推理困难而提出的一种深度模型,也是目前深度学习最主要的实现方式之一.DBN可以通过受限玻尔兹曼机的顺序堆叠来构造,其学习过程分为两个阶段,首先是对RBM进行逐层无监督预训练,然后再用反向传播算法对整个网络进行有监督的调优.本节重点介绍DBN的无监督学习.RBM和DBN的结构分别如图1和图2所示.图1RBM结构图Fig.1Structure of RBM图2DBN结构图Fig.2Structure of DBN给定模型参数θ=(w R,b v,b h),那么可视层和隐含层的联合概率分布P(v,h;θ)用能量函数E(v,h;θ)定义为P(v,h;θ)=1Ze−E(v,h;θ)(1)其中,Z=v,he−E(v,h;θ)是归一化因子,模型关于1期王功明等:深度信念网络研究现状与展望37 v的边缘分布为P(v;θ)=1Zhe−E(v,h;θ)(2)对于一个伯努利(可视层)分布–伯努利(隐含层)分布的RBM,能量函数定义为E(v,h;θ)=−mi=1b vi v i−nj=1b hj h j−m i=1nj=1v i w Rijh j(3)其中,w Rij是RBM的连接权值,b vi和b hj分别表示可视层节点和隐含层节点的偏置.那么条件概率分布可表示为Ph j=1v,θ=σb hj+mi=1v i w Rij(4) Pv i=1h,θ=σb vi+nj=1w Rijh j(5)式中,σ(·)是一个Sigmoid函数.由于可视层和隐含层是伯努利的二值状态,所以判断它们二值概率取值的标准常通过设定一个阈值来实现[31].通过计算对数似然函数log P(v;θ)的梯度,可以得到RBM权值更新公式为w Rij (τ+1)=w Rij(τ)+η∆w Rij(6)∆w Rij=E data(v i h j)−E model(v i h j)(7)式中,τ和η分别表示RBM的迭代次数和学习率, E data(v i h j)和E model(v i h j)分别表示训练集中观测数据的期望和模型所确定分布上的期望[32].特别地, RBM有一个有趣的性质,即当利用基于最大似然的学习规则训练时,连接两个神经元的特定权重的更新仅取决于这两个神经元在不同分布下收集的统计信息:P model(v)和ˆP data(h/v).网络的其余部分参与塑造这些统计信息,但是权值参数可以在完全不知道网络其余部分或这些统计信息如何产生的情况下更新.这意味着学习规则是“局部”的,这使得RBM的学习似乎在某种程度上是符合生物学机理.我们可以设想每个神经元都是RBM中随机变量的情况,那么连接两个随机变量的轴突和树突只能通过观察与它们物理上实际接触细胞的激发模式来学习.特别地,经常发生某种强烈的脉冲激励时的两个神经元之间的连接会被加强,这就是Hebb学习规则的核心思想.Hebb学习规则给出了生理学与心理学之间的内在联系,该规则至今仍被许多神经网络学习算法所使用.作为一种深层网络模型,DBN兼具生成模型和判别模型的双重属性.因为DBN的预训练过程主要用来表达数据的高阶相关性或者描述数据的联合统计分布,具有生成模型的特点;DBN有监督调优过程通常用来分类数据的内在模式或者描述数据的后验分布,具有判别模型的特点.这里的“生成”是指从隐含层到输入数据的的重构过程,而“判别”是指从输入数据到隐含层的归约过程.同时,作为一种生成模型,生成式对抗网络(Generative adversarial network,GAN)近年来同样受到很大的关注并进行了广泛的应用[32−33].GAN实质上属于一种基于深度学习的混合模型,其通过框架中生成模型和判别模型的互相博弈学习产生相当好的输出.从数据生成角度看,GAN的数据生成过程是在有监督信号的反馈作用下完成的.而DBN作为一种生成模型时,其监督信号是数据本身,即通过对原始数据的重构完成网络的训练,从而具有生成能力.具体应用中, DBN常作为GAN的生成模型,与判别模型进行对抗学习[32].DBN学习模型的优点是通过组合许多RBM,把上一层RBM的特征激励作为下一层的训练数据,可以高效地对隐含层进行学习.递归神经网络(Recurrent neural networks,RNN),它的深度甚至可以达到和输入数据序列的长度一致.在无监督学习模式下,RNN被用来根据先前的数据样本预测未来的数据序列,并且学习过程中没有用到类别信息.然而,RNN在近几年才得以广泛使用,部分原因是由于在训练中遇到的梯度弥散或梯度爆炸问题,它很难通过训练来捕捉长时相关性.随着在Hessian-free优化研究方面的进展,在一定程度上解决了这个问题,该方法使用了近似二阶信息或随机曲率估计.另外,RNN没有基于无监督预训练的参数初始化过程,这也是其与DBN在训练原理上的最大区别.卷积神经网络(Convolutional neural net-works,CNN)是另一种具有判别性能的深度学习网络,它的每个模块都是由卷积层(Convolutional layer)和池化层(Pooling layer)组成.卷积层共享权值,池化层对卷积层的输出进行降采样,减少了下一层的数据量.研究发现,CNN的应用主要集中于计算机视觉或者图像识别领域,并且效果较为出色[34].而DBN的应用则广泛分布于计算机视觉和数据建模及预测等领域.另一种与DBN相似的深度结构基本学习模型是自编码器(Auto encoder),自编码器主要用于完成数据转换的学习任务,在本质上是一种无监督学习的非线性特征提取模型.自编38自动化学报47卷码器与DBN也有着重要的区别,这种区别的核心在于:自编码器希望通过非线性变换找到输入数据的特征表示,它是某种确定论性的模型;而DBN的训练则是围绕概率分布进行的,它通过输入数据的概率分布(能量函数)来提取高层表示,是某种概率论性的模型.另外,DBN具有较多的超参数,可分为两类:一类是训练参数(如学习率和动量项);另一类是定义网络结构的参数(如网络层数和每层神经元数).前者的自动调优属于超参数优化(Hyperparame-ter optimization,HO)的范畴,而后者的自动调优一般称为神经网络架构搜索(Neural architecture search,NAS).严格地讲,NAS属于DBN结构设计的方法之一,目前DBN结构设计大多数通过提前赋值来完成,即在网络训练过程中结构不变,只有训练参数在不断调整.本文即将介绍的两种变结构设计策略(自组织结构和增量式结构)对固定结构来讲是一种突破,但是与NAS又存在区别,主要体现在: NAS先定义搜索空间,然后通过搜索策略找出候选网络结构,对它们进行评估,根据反馈进行下一轮的搜索;而变结构策略只要是以某种触发机制或误差导向来实时调整结构规模.2深度信念网络结构分析与性能比较2.1固定结构深度信念网络目前最为常见的DBN应用形式是定结构模型,即在训练过程中DBN结构固定不变.尽管现在与其他无监督或生成学习算法相比,固定结构的DBN 大多已经失去了青睐并很少使用,但它们在深度学习历史中的重要作用仍应该得到承认[20].定结构DBN在处理实际复杂数据时,无监督预训练和反向传播调优算法均具有提升和改进的空间,主要表现在预训练耗时和调优精度两方面.同时,定结构DBN主要是通过足够的经验和充足的数据来确定其结构,且其结构一旦确定将不再调整,这导致定结构DBN无法满足待处理数据的多样性变化要求.众所周知,DBN具有很强的计算和信息处理能力,但是它对于模式识别、感知以及在复杂环境中作决策等问题的处理能力却远不如人.神经生理学研究结果表明,人的智能主要取决于大脑皮层,而大脑皮层是一个大规模互连的生物深度神经网络.在处理不同信息时,生物深度神经网络会启用不同的神经元连接结构,也就是说,其采用的是一种变结构的信息处理机制[35].而在实际过程中,定结构DBN只是通过改变权值参数来适应任务的变化,但如何构造一种DBN使其结构在动态调整的同时不断调整权值参数,是今后DBN发展的趋势,也是一个开放且尚未解决的问题[36].2.2稀疏深度信念网络研究发现,现有的DBN模型在学习过程中内部神经元之间的权值连接均是一种密集表述[37−38].然而,在深度学习算法中,一个主要的目的是独立地表述数据的差异[36],密集表述容易导致网络不稳定,因为任何输入上的扰动都会引起中间隐含层特征表述向量发生变化,甚至是巨变[38].稀疏表述就是用较少的基本信号的线性组合来表述大部分或者全部的原始信号.利用稀疏表述对DBN进行稀疏连接训练,可以有效地降低输入扰动对中间隐含层特征表述向量的影响[39].无监督学习过程中的稀疏表述原理如图3所示.图3稀疏表述原理图Fig.3Sparse representation schemeLee等[40]通过在RBM训练过程中引入一个正则化惩罚项来降低密集表述的程度.具体来讲,首先设置一个隐含层神经元的期望激活强度值,然后惩罚隐含层神经元实际激活强度与期望激活强度之间的偏差.给定m组训练数据集v(1),...,v(m),其实现稀疏表述的优化问题为Maximizeθlog P(v)+λR sparse1(8) R sparse1=−µ−1nnj=1Eh jv2(9)其中,λ是正则化常数,µ是控制着第j个隐含层神经元稀疏度的期望激活强度值,通过这种提前给定期望激活阈值的方法可以实现一定意义上的稀疏表述.为了使所有隐含层神经元能够以一定的概率或者波动性逼近期望激活强度值,Keyvanrad等[41]通过引入正态函数的集中分布思想来控制网络的稀疏度.根据这种思想,对应于稀疏表述优化问题的正则化项可表示为R sparse2=nj=11σ√2πe−(h j−µ)22σ2(10)其中,σ是控制稀疏强度波动性的方差.1期王功明等:深度信念网络研究现状与展望39同时,应该注意到参数的设置对网络学习效果的影响是显著的[38],如果设置不当,要实现较高精度的建模并学习到正确的特征信息往往比较困难.因此上述稀疏表述方法虽然在网络性能的鲁棒性方面取得一定程度的效果,但对无监督学习的迭代次数和神经元数量等有一定依赖.2.3自组织深度信念网络目前DBN在应用中存在一个重要问题,即针对不同的问题,DBN需要提前设置网络深度,然后利用经验法比较各种不同深度的精度和训练效果.这极大地制约了网络解决问题时的效率,使DBN的进一步推广与应用受到很大限制.实际上,著名深度学习专家Bengio在2009年提出了一个与此类似的问题[36],该问题原文描述为:“Is there a depththat is mostly sufficient for the computations nec-essary to approach human-level performance of AItasks?”.意思是,是否存在一个合适深度的DBN,可以用来尽可能像人类解决问题那样去解决大多数的AI问题呢?由于该问题比较笼统,涉及的学科范围太广,很难通过一个有效的数学方法来解决该问题,难以设计出包含较多的特征并具有代表性的实验对其进行验证,因此该问题在短时间内难以得到彻底的解决.目前,针对此问题的初步试探性解决方法有结构自组织策略和凑试法.本节只介绍结构自组织策略.Qiao等[26]提出了一种基于神经元激活强度和误差下降率最小化的结构自组织方法.首先,在无监督预训练阶段将隐含层神经元的激活强度作为神经元的“贡献度”,并根据“贡献度”的大小对神经元进行增加或删减.其次,在有监督调优阶段,将训练误差的下降率作为隐含层的删减标准,当训练误差下降率首次出现递减时删掉一个隐含层,否则增加隐含层.激活强度SI可表示为SI li =αs2i·l1+s2i·l+o2i·l(11)其中,α是正常数,o i·l是第l个隐含层的第i个神经元的输出,i=1,2,3,···,N l,N l是第l个隐含层的神经元个数,s i·l表示第l个隐含层的第i个神经元的输入权值之和,可通过如下公式计算得到s2 i·l =n ij=1w ij r ij(12)其中,r ij是i个神经元的第j个输入量,w ij是第j个输入神经元和第i个神经元之间的连接权值,n i 是第i个神经元的输入神经元个数,s i·l所表示的权值连接过程如图4所示.DBN的结构自组织策略原理如图5所示.在传统浅层神经网络的结构设计方面,研究人员注重结构自组织设计方法[42−43],即根据神经元激活强度的大小来增加或者删减结构.尽管结构自组织设计方法在浅层神经网络中得到了成功的应用并取得了较好的效果,但关于DBN结构自组织方法的研究却非常有限.本节介绍的基于传统自组织方法的变结构DBN模型在学习精度上有所提高,但是在学习效率方面提高不明显,相关研究还需要进一步加强.图4计算激活强度的权值连接过程Fig.4Weights connecting process of computing spikingintensity2.4增量式深度信念网络与传统浅层神经网络的结构自组织相比,DBN 结构自组织策略一直没有得到学术界的广泛关注,主要原因有:1)自组织方法将神经元的激活强度作为增加和删减结构的评判标准,而DBN往往拥有多个隐含层且每个隐含层含有较多的神经元,这导致DBN自组织设计过程复杂且计算量庞大[26];2)预训练好的初始DBN可被视为一种知识源域(Source domain),其中的知识可被视为一种可重复利用的经验[44],但是结构自组织方法未能在知识源域到目标域(Target domain)之间实现知识的转移.因此,在DBN结构自组织过程中需要不间断地对目标域内若干个新增子结构进行参数初始化,从而导致自组织方法在DBN结构设计中应用成本较高,甚至难以实现.通过上述分析可知,DBN结构自组织方法遇到的主要障碍是计算量巨大,而如何在知识源域与目标域之间实现知识的有效转移成为关键.迁移学习(Transfer learning,TL)是一种旨在实现知识转移的学习方法且具有较强的鲁棒性[45−47].常用的迁移学习方法是:首先训练一个模型并将其作为知识源域,然后再利用特定的方法将知识源域中可重复40自动化学报47卷图5结构自组织策略原理图Fig.5Self-organizing structure strategy scheme利用的知识转移到目标域中来加速新结构的学习过程[48−49],从而提高复杂模型的训练效率.近些年来,基于迁移学习的神经网络复合训练方法大批涌现并取得了较好的效果[50−51].为了解决上述问题,Wang等[27]提出了一种基于迁移学习策略的增量式深度信念网络(TL-GDBN)模型.相较于浅层神经网络的结构自组织方法,不同之处在于TL-GDBN没有利用神经元的激活强度作为结构增长或删减的依据.首先,初始化一个单隐含层DBN并对其进行预训练(Pre-training),然后固定预训练好的初始DBN并将其作为知识源域.其次,在初始DBN的基础上不断增加固定规模的隐含层和神经元并将其作为目标域,同时建立基于迁移学习的知识转移规则来加速目标域的训练过程.第三,根据TL-GDBN的预训练的重构误差设置结构增长的停止准则及其阈值,从而获得最优的结构.基于迁移学习的增量式深度信念网络(TL-GDBN)的结构增长过程仅在预训练阶段进行.每一步的结构增长包括神经元和隐含层两部分.数据被分为三部分:训练数据(Training data)、验证数据(Validating data)和测试数据(Testing data).训练数据用来预训练初始DBN并获得知识源域,验证数据用来结合迁移学习实现TL-GDBN结构的增量式变化,测试数据用来测试TL-GDBN.预训练结束后TL-GDBN结构将不再变化.知识在迁移学习规则下持续地被转移到新增结构中,TL-GDBN的一步结构增长过程如下:步骤1.结构初始化和预训练.首先初始化一个单隐含层的DBN结构,然后利用对比散度(Con-1期王功明等:深度信念网络研究现状与展望41trastive divergence,CD)算法和训练数据进行预训练.假设初始化DBN的输入和其隐含层神经元的个数分别为m和n,那么预训练后学习到的知识(权值参数矩阵)w R1∈R m×n将被保存在知识源域中.步骤2.增加神经元.增加两倍于初始DBN隐含层神经元数量的神经元,新的权值参数矩阵变为ˆw R1∈R m×3n.步骤3.增加隐含层.增加与初始DBN具有相同数量神经元的隐含层,对应的新增权值参数矩阵为w R2∈R3n×n.步骤4.计算预训练过程的重构误差,并将重构误差作为预训练过程误差.步骤5.设置结构增长的停止准则.利用验证数据计算重构误差,并将重构误差的连续若干步的减小量作为结构增长的停止准则.同时设置停止准则的阈值,当训练过程中的重构误差满足阈值条件时, TL-GDBN结构停止增长并进入步骤6;否则,跳转到步骤2.步骤6.固定当前TL-GDBN的最优结构,预训练过程结束.TL-GDBN的一步结构增长过程原理如图6所示.结构增长过程一旦结束,TL-GDBN的结构和对应的初始权值参数即被确定.实验结果发现,TL-GDBN的稀疏度随着结构的不断扩大而表现出先增大后稳定的趋势.这种趋势表明在结构增长过程中TL-GDBN的密集表述越来越弱,网络各隐含层提取到的特征向量受输入波动影响的程度也越来越弱,即网络鲁棒性较强.然而,关于如何进行知识迁移仍然是一个难点,究其原因主要在于:在迁移学习中,学习器必须执行两个或更多个不同的任务,但是我们假设能够解释P1变化的许多因素和学习P2需要抓住的变化相关.例如,我们可能在第一种情景中学习了一组数据分布特性,然后在第二种场景中学习了另一组数据分布特性.如果第一种情景P1中具有非常多的数据,那么这有助于学习到能够使得从P2抽取的非常少的样本中快速泛化表示.一般来讲,当不同情景或任务存在有用特征时,并且这些特征对应多个情景出现的潜在因素,迁移学习可以发挥事半功倍的效果.然而,有时不同任务之间共享的不是输入的数据分布特性,而是输出的目标数据分布特征.这种情况下,使用迁移学习往往会得到不尽人意的学习效果.2.5递归深度信念网络从学习策略上看,传统DBN模型是一种前馈网络,堆叠的RBM只能保存暂时的信息(达到能量平衡后的稳态信息),故现有的DBN模型对时间序列的建模与预测精度相对较低[52−55].Ichimura 等[30]提出一种递归深度信念网络(RNN-DBN),其在结构上是由若干个递归受限玻尔兹曼机(RNN-RBM)[56]堆叠组成.而RNN-RBM则是在递归时间RBM(RTRBM)的基础上发展起来的[52],是一种基于能量的时域序列密度估计模型.RTRBM结构如图7所示.图6TL-GDBN的一步增长过程Fig.6Illustration of one-growing step。

and Ni(111) surfaces studied by DFTaqueous environment as function of pH and potential. In the present work we identify the most stable surface structures relevant for oxygen reduction on the (111) surfaces of Ag, Ni and Pt. We refer to the corresponding phase diagrams for the surfaces as surface Pourbaix diagrams. Furthermore, for the three different metals we determine the potential where surface catalysis starts playing a limiting role for ORR. This introduces aself-consistency problem, since the ORR potential will depend on the surface structure and the surface structure depends on the potential. Fortunately only a few surfaces are relevant for ORR. In the context of steady state ORR the cathode surface structure is independent of pH. In other words: the ORR potentials versusthe reversible hydrogen electrode (RHE) are the same in alkaline and acidic solution. Hence, we find that the reason for the activity of Ag and Ni in alkaline solution compared to acidic solution, is not related to a different surface structure, but rather related to higher dissolution potential (vs. RHE) and better stability in alkaline solution than in acidic solution.2. Methods2.1 Surface Pourbaix diagramsWe determine the most stable surfaces without molecular oxygen present, as these are the relevant surfaces for the surface Pourbaix diagrams and normal cyclic voltammograms (CV). To calculate the stability of a surface we use the concept of the theoretical standard hydrogen electrode (SHE),5 which has previously been shown to give reliable results for theoretical cyclic voltammograms for hydrogen on Pt surfaces8 and for predicted ORR activity for Pt alloys.6The method is summarized here. It is assumed the surface is in equilibrium with protons and liquid water at 298 K, so that oxygen and hydroxyl may be exchanged between the surface and a reference electrolyte through the following stepsH2O(l) + * HO* + H+(aq) + e−(1) andHO* O* + H+(aq) + e−,(2) where * denotes a site on the catalyst surface and O* and HO* denote oxygen and hydroxide adsorbed on a site.The oxidation of water to O* and HO* depends on the potential U and pH through the chemical potentials of H+ and e−. The reaction:H+(aq) + e− 1/2H2(g)(3) has G° = 0 at standard conditions (pH = 0, p(H 2) = 1 bar) defining U = 0 V SHE. This allows us to calculate the free energy of adsorbed O* and HO* at U = 0 V and pH = 0 directly from the free energy of the reactionsH2O(l) + * HO* + 1/2H2(g)(4)HO* O* + 1/2H2(g).(5) This has the advantage that H2(g) is much easier to handle computationally than solvated protons. Reactions (4) and (5) define the free energy, G 0, of HO* and O* at zero pH and zero potential vs. SHE.At finite pH and potential, the free energy of reaction (3) is (p(H2) = 1 bar)G = e U SHE + k B T ln 10 pH,(6) Since (4) and (5) are independent of pH and potential, we may calculate the free energy of (1) and (2) at any pH and potential fromG(HO*) = G0(HO*) −e U SHE−k B T ln 10 pH +(7)G field,where G fieldis the change in the adsorption energy due to the electric field in the electrochemical double layer at the cathode. To first order in the electric field G field = −μ·E, where μ is the dipole moment of the adsorbate and E is the electric field at the position of the dipole. G field has been studied in detail for Pt(111),9,10 and the relative change in stability of O* and HO* is less than 0.11 eV when the potential is increased from 0 to 1 V vs. SHE. Since the first term in eqn (6)changes the relative stability of O* and HO* by 1 eV when the potential is increased by 1 V, the trends in adsorption energies are well described by neglecting G field in the construction of the surface Pourbaix diagram.10We calculate G 0 by correcting the DFT energies for zero point energies and entropy viaG0 = E + ZPE −T S,(8) where E is the DFT energy andZPE is the change in zero-point energy of the adsorbates. We have used the zero point energy of the adsorbates on Cu(111) obtained within the harmonic approximation using the PW91 functional.11,12S is approximated from the loss of entropy of the gas phase molecules upon binding them to the surface.13,14 We note that the method described above could equally well have been formulated in terms of OH− instead of H+, without changing the results.The free energy of a given molecular surface structure can be calculated as a function of potential U SHE and pH from eqn (7). The structure with the lowest free energy at a given set of conditions determines the surface structure. An example is shown in Fig. 1, where the stability of different Ag(111) surfaces is plotted against potential at pH = 0. The lowest line determines the surface with the lowest free energy at a given potential. The surface Pourbaix diagram in Fig. 2is then constructed by plotting the most stable surface as function of pH and U SHE.Fig. 1 Stability of O* and HO* on Ag(111) at pH = 0.Dissolution is spontaneous at U RHE > 0.80 V.Neglecting dissolution, a structure with 1/6 ML HOstabilized by water forms at 0.93 V. The hydroxylcoverage increases to 1/3 ML at 1.11 V.Fig. 2 Surface Pourbaix diagram for Ag(111). Thedashed line at U RHE 0.8 V marks the potential of agood ORR catalyst under steady state conditions. It isseen that dissolution is a much larger problem inacidic than in alkaline electrolytes, and that thecoverage of the Ag surface is low at U RHE 0.8 V.For HO* and HOO* we have included the possible stabilization by forming hydrogen bonds to water by incorporating the adsorbates in a hexagonal hydrogen bonded network with water. On Pt(111) this hexagonal water layer is known to be stable and has been observed in UHV experiments.15,16 This approach for including water has previously been used in calculations.5,9We adjust the water layer found on Pt to fit Ni and Ag. The considered structures are shown in the ESI.The theoretical reversible hydrogen electrode (RHE) can be defined by demanding that the equilibrium in reaction (3) defines U RHE = 0 at all conditions. The potential versus the reversible hydrogen electrode becomesU RHE = U SHE + k B T ln10 pH/e.(9) Clearly, at pH = 0, U SHE = U RHE and any fixed potential vs. RHE, e.g. the potential for reversible water oxidation will show up as lines with a slope of −k B T/e ln10 (−59 mV/pH at 298 K) in the Pourbaix diagrams.A special case of this is the formation of HO*, because by rewriting eqn (7) in terms of U RHE and neglecting the effect of the electric field,G(HO*) = G0(HO*) −e U RHE,(10) we can see that hydroxyl adsorbs at a fixed potential vs. the RHE. A similar argument holds for the adsorption of O* formed from water, which means that the coverage and structure of adsorbed O* and HO* at a givenU RHEare independent of pH within this model. As the potential of the cathode relative to the hydrogen electrode determines the fuel cell potential, the coverage and structure of adsorbates at the cathode depend only on the fuel cell potential under steady state conditions and not on pH, provided we disregard dissolution and adsorption of the electrolyte ions.The metals may dissolve in acidic or alkaline solution. The free energy of dissolution has been estimated from standard reduction potentials and thermochemical data, as uniform dissolution of a monolayer of metal atoms on the surface corresponds to dissolution of the bulk. When speaking of a process as spontaneous, it will be in the sense G°(U) < 0, which gives a conservative estimate of the stability vs. dissolution. An example is dissolution of silver in acidic solution:Ag(s) Ag+(aq) + e−.(11) Since this process involves electrons but not protons, the free energy of dissolution of the metal surface depends on the potential vs. SHE but not on pH. Hence the boundary between Ag and Ag+-ions is horizontal in a Pourbaix diagram. Formation of 1/2 ML O on Ag(111) from H2O requires transfer of one electron per Ag atom. The formation of Ag+ions also involves one electron. The difference in free energy of 1/2 ML O on Ag and Ag+ ions therefore do not depend on the potential, so the boundary becomes vertical in a Pourbaix diagram.Adsorption of spectator ions from the solvent, e.g. anions, usually involves only electrons.*+ Cl−(aq) Cl* + e−(12) The free energy of anion adsorption therefore depends on the potential vs. SHE, rather than the potential vs. RHE. If HO* and O* compete with spectator species for the free sites, then this might change the surface coverage of HO* and O* as function of potential versus the reversible hydrogen electrode for different pH. This effect has not been considered in this work.2.2 Calculation detailsAll calculations have been done within density functional theory (DFT),17,18 with the RPBE functional19 for exchange and correlation. The RPBE functional has been shown to provide better19 adsorption energies than the closely related PBE20 functional. The Kohn–Sham equations18 have been solved using a plane wave basis with a cutoff of 350 eV for the eigenstates and the densities have been described using a cutoff corresponding to 1000 eV. The ionic cores are described using Vanderbilt ultrasoft pseudo potentials,21 which give faster convergence with respect to the plane wave basis set than norm conserving pseudo potentials.The adsorption energies on Ag, Pt and Ni(111) have been calculated using (2×2) and (2×3) supercells consisting of 3 close packed layers of metal atoms. The top metal layer is allowed to relax and the two bottom layers are fixed in their bulk positions. The adsorbates are allowed to relax freely on the surface, until the maximum component of the Cartesian forces is below 0.02 eV −1. The periodically repeated slabs have been separated by at least 12.4of vacuum and the dipoles of the images have been decoupled. The Brillouin zone has been sampled using an [8×8×1] and a [10×6×1] Monkhorst-Pack special grid of k-points22 for the (2×2) and the (2×3) cells, respectively. The DFT calculations have been performed using Dacapo with the ASE interface on our local Linux cluster supercomputer.233. Results3.1 Ag(111)To investigate the Ag(111) surface, different structures of hydroxyl and oxygen on the surface have been considered. In acidic solution, dissolution of Ag is spontaneous at potentials above 0.80 V, see Fig. 1, so ORR on this surface is not likely. In alkaline environment the dissolution process is24Ag(s) + H2O(l) AgO−(aq) + 2H+(aq) + e−, U° = 2.2(13)V.In alkaline solution at low potentials, the most stable surface considered is the clean Ag(111) surface, showing that Ag is a noble metal. At pH = 14, the free energy of dissolution in (13) is negative for U SHE > 0.5 V for the clean (111) surface. Dissolution is therefore not as critical as in acidic solution, see Fig. 2. Water oxidation starts at 0.93 V (RHE) forming a structure with 1/6 ML HO* and 1/2 ML water. At 1.11 V (RHE) the HO* coverage increases to 1/3 ML HO*. At potentials above 1.19 V (RHE) hydroxyl is oxidized further to O*. At higher potentials we find that as more oxygen or hydroxyl is adsorbed on the surface, O atoms will spontaneously substitute Ag atoms on the surface when the surface is relaxed. This may be interpreted as an onset of dissolution or oxidation of the Ag(111) surface. We note that the substitution is favored by 0.25 V relative to the alkaline dissolution process at pH = 14. The pH and U range where the observed substitution mechanism is favorable instead of alkaline dissolution is shown in the Pourbaix diagram Fig. 2. The substitution becomes spontaneous at potentials above 1.44 V (RHE).The potential of the onset of dissolution or oxidation is determined by the details of the substituted Ag surface, and it is very likely that surfaces with different atomic configurations turn out to be more stable than the ones we have considered—for example, a surface oxide25,26 or even formation of bulk Ag2O which may form at potentials above 1.18 V (RHE). The kinetics of the oxidation and dissolution mechanisms of theAg(111) surface is without the scope of the present paper.We have added the potentials of the reversible hydrogen and oxygen electrodes to the Pourbaix diagram to mark the stability range of water at U RHE = 0 V and U RHE = 1.23 V, respectively. We have also added a line at U RHE= 0.80 V, because it is a realistic potential for a good ORR catalyst. We emphasize that the surface structure does not vary along the lines of constant U RHE, if the surface does not dissolve. This means that at U RHE = 0.80 V, the same surface structure will be relevant for ORR regardless of pH.CV experiments in alkaline electrolytes show a reversible peak around 0.3 V (RHE) on Ag(111), and at lower potentials on Ag(100) and Ag(110).4,27This peak has been suggested to be caused by hydrogen adsorption, alkali metal deposition or oxidation of water to HO*.28Recent EC-STM experiments have shown that the Ag(111) surface is unmodified up to a potential near 0.7 V (RHE), after which protrusions, which are assigned to 2D Ag2O nuclei, begin to emerge.29 Our calculations do not support the view that the CV peak at 0.3 V (RHE) is due to HO* or H* on the flat Ag(111) surface, as we find water oxidation to start at 0.93 V (RHE). To get HO* at 0.3 V therefore requires a stabilization of HO* by 0.63 eV relative to our calculations. Such stabilization could, in principle, be caused by HO* adsorbing at steps or defects rather than on the Ag(111) terraces. This is however an unlikely explanation for the existence of a peak at 0.3 V (RHE), as the HO* coverage is reported to be 0.2–0.35 ML.4,27 We note that our theoretical HO* formation is in good agreement with the high potential adsorption peak starting at 0.9 V (RHE) also observed in CV experiments.43.2 Pt(111)For Pt the surface structures have been investigated previously.9,10 Here we have redone the calculations with the same setup as for Ag and Ni. At potentials below 0.78 V (RHE) the pure Pt surface without adsorbates is the most stable. Water is oxidized to O* at potentials above 0.78 V (RHE). From 0.78 V (RHE) to 1.0 V (RHE) the O coverage is 1/4 ML. At higher potentials the coverage increases gradually to 1/3 ML O, then 1/2 ML O etc., see the Pourbaix diagram Fig. 3. Potentials near 1.5 V are required before standard dissolution becomes thermodynamically favored.Fig. 3 Surface Pourbaix diagram for Pt(111). IfG field is neglected, the oxygen coverage increasesgradually starting at 0.78 V (RHE). If G field isincluded as found by Karlberg et al.,9 1/3 ML OH isless than 0.05 eV more unstable than O* in thepotential range 0.78 VRHE–0.84 VRHE (black area).We consider O* and HO* to be equally stable in thispotential range within the accuracy of thecalculations.Recent DFT calculations including the electric field in the double layer show that 1/3 ML hydroxyl in a hexagonal layer of water is stable from 0.75 V (RHE) to 0.9 V (RHE), as a result of the electric field stabilizing HO* relative to O*.9,10 This surface has also previously been found to be stable in UHV.15,16 Including the effect of the electric field from ref. 9, we find that hydroxyl in the water layer is unstable byless than 0.05 eV relative to 1/4–1/3 ML O, in the potential range from 0.78 to 0.84 V (RHE). We therefore consider the stability of 1/3 ML hydroxyl and 1/3 ML O to be the same near 0.80 V (RHE) within the accuracy of the calculations.3.3 Ni(111)The Ni(111) surface is quite reactive, and corrosion is a problem in acidic as well as alkaline solutions. We consider the acidic dissolution13 processNi(s) Ni2+(aq) + 2e−, U° = −0.257 V.(14) and the dominant alkaline corrosion process13,30,31Ni(s) + 3H2O(l) 2e− + 3H+(aq) + Ni(OH)3− (aq),(15)U° = 0.63 V.We find that at U RHE= 0 V, 1/4 ML H is adsorbed at the surface up to 0.14 V (RHE). Oxygen starts to adsorb at U RHE > 0.20 V, with 1/3 ML O* stabilized by water being the most stable structure in the potential range from 0.34 V (RHE) to 0.82 V (RHE). The highest resistance against corrosion is obtained at pH = 9.9, where G° < 0 for U RHE > 0.38 V, as seen from Fig. 4 and 5. The free energy of formation for O* is 0.39 eV/O at 1/4 ML O coverage as calculated from (4) and (5), which is in reasonable agreement with the experimental value of 0.49 ± 0.10eV/O*.32Fig. 4 The most stable Ni(111) surfaces at pH = 9.9.The resistance against dissolution is at maximum atthis potential because G of the acidic and alkalinedissolution are equal at this potential. Neglectingdissolution, the oxygen coverage then increasesgradually with 1/3 ML O having a broad stabilityrange from 0.34 V (RHE) to 0.82 V (RHE).Dissolution is spontaneous for U > 0.38 V (RHE).Fig. 5 Surface Pourbaix diagram for Ni(111)including alkaline and acidic dissolution. Dissolutionis a problem in both acidic and alkaline electrolytes.The best resistance against dissolution is obtained atpH = 10, where dissolution becomes spontaneous forU RHE > 0.38 V. Oxygen is adsorbed for U RHE > 0.20V.In contrast to Taylor et al.,33we do not find HO* or co adsorbed O* and HO* to be stable on the Ni(111)-surface. However we note that at 0.18 V the free energies of 1/4 ML H*, 1/4 ML O* and 1/6 ML HO* are equal within 0.1 eV. The difference could therefore be caused by different functionals and atomic structures of the adsorbates.In UHV experiments oxygen on Ni(111) may form an ordered p(2×2) structure with 0.25 ML coverage or a-structure with 0.33 ML coverage.34 Oxidation of the Ni(111) surface initiates at a coverage between 0.33 and 0.5 ML,32,34suggesting that ORR on the metallic Ni surface could occur on a surface with up to 1/2 ML O. Thermodynamically -Ni(OH)2 is however the most stable nickel phase at these conditions.303.4 Oxygen reduction reactionWe now consider the situation in the presence of molecular oxygen. The aim is to determine the intermediates of the oxygen reduction reaction at the relevant surface for the different metals. In an acidic solution we write the electrode reactions as:Anode:H2 2H* 2H+ + 2e−(16) Cathode:O2 + 4H+ + 4e− HOO* + 3H+ + 3e− H2O + O*(17)+ 2H+ + 2e− HIn an alkaline solution the same reactions can be expressed with OH− instead:Anode:2OH− + H2 2H* + 2OH− 2H2O + 2e−(18) Cathode:O2 + 2H2O + 4e− HOO* + H2O + OH− + 3e−(19)O* + H4OH−It is interesting to note that the adsorbed intermediates along the reaction coordinate are independent of the concentration of protons. The only change is the chemical potential of the protons, and that the proton donor changes from being H+ or hydronium ion in acid electrolytes to being H2O in alkaline electrolytes.We now calculate the free energy of the intermediates in eqn (17) or (19). We will consider ORR feasible as long as all reaction steps along the reaction decrease the free energy. We determine the highest potential where this still holds, and refer to this as the ORR potential. The ORR potential depends on the metal and surface structure, which means we have to find the surface structure with the highest ORR potential that is consistent with the stability of the surface structure. We use the fact that the most stable surface at a given potential found above will be one of the intermediates during the ORR. Free energy diagrams of the intermediates on Pt, Ag and Ni are shown in Fig. 6, 7 and 11 at the ORR potential and at 0 V (RHE) and 1.23 V (RHE). The higher the ORR potential, the better the ORR catalyst.Fig. 6 Stability of the intermediates for the oxygenreduction reaction on Pt(111). The stability of theintermediates has been corrected for the effect of theelectric field using the results of Karlberg et al.9 Theformation of HOO* and desorption of HO* becomeexothermic at nearly the same potential, in agreementwith previous findings that Pt is close to the optimalcatalyst for ORR.Fig. 7 Stability of intermediates for the oxygenreduction reaction on Ag(111). The formation ofHOO* determines the overpotential, in agreementwith previous trend studies showing that Ag bindsadsorbates too weakly. The stability of theintermediates has been corrected for the effect of theelectric field using the results of Karlberg et al.9It has previously been shown for Pt and Pt-alloys that the ORR potential is determined by either reduction of HO* (the last step in reaction 17) or the formation of OOH* (the first step in reaction 17).9 For reactive surfaces the first of the two is rate determining and for noble surfaces the latter is rate determining. Both reaction energies scale linearly with the oxygen binding energy, which determines a Sabatier volcano where the lowest overpotential is found for a material with the optimal trade-off between strong and weak binding of oxygen. Pt is close to this optimum, Ni binds too strongly and Ag binds just too weakly. However, previous analysis has been done for the surface relevant for Pt. As seen above, the most stable surface of Ni at a given potential will be different from the most stable Pt surface which, again, is different from the Ag surface.3.4.1 ORR on Pt(111).We have previously studied ORR on Pt in details. The results are summarized here for comparison with Ag and Ni. The most stable surface in the relevant potential region is O* and the mixed HO*-water layer as discussed in details above. We therefore use these two surfaces to represent the two last ORR intermediates in eqn (17). We construct the OOH* intermediate from the HO*-water layer by inserting an oxygen atom between HO* and the surface.The stability of the intermediates is shown in Fig. 6 at potentials 0.0, 0.75 and 1.23 V vs. RHE. The two potential determining steps are O2(g) + H+ + e− + * OOH* and HO* + H+ + e− * + H2O(l).5 That these two steps determine the same ORR potential reflects that Pt has a reactivity close to optimal, as increasingG 0 for one of the steps will decrease G0for the other step due to the linear relations between the adsorption energies. Dissolution of Pt is not a big problem since it requires U SHE > 1.2 V.3.4.2 ORR on Ag(111).We find the Ag(111) surface to be free of adsorbates up to 0.93 V (RHE), where 1/6 ML hydroxyl is adsorbed, followed by additional HO* formation to form a structure with 1/3 ML hydroxyl at 1.11 V. At any reasonable overpotential we would therefore expect the total coverage on the surface to be very low. At 1/6 ML total coverage the rate limiting step is the formation of HOO*, which becomes exothermic at potentials below 0.73 V. This is in agreement with previous estimations based on DFT calculations of the oxygen binding energy on Ag(111).5The adsorption of HOO* being rate limiting reflects the fact that Ag is a noble metal that binds the intermediates too weakly relative to the optimal ORR catalyst. Ag therefore seems like a good alternative to Pt, however the problem in acidic solution is dissolution and not activity. Dissolution involves only electron transfer to Ag, which means that the relevant dissolution potential will not change vs. SHE. Since ORR in alkaline environment runs at lower potentials vs. SHE, Ag-dissolution becomes less of a problem.Assuming the effect of the electric field in the dipole layer is similar on Pt(111) and Ag(111), the electric field destabilizes O* by ca. 0.04 eV and stabilizes HO* by ca. 0.06 eV at the ORR potential. The adsorption energy of HOO*, which determines the overpotential in this model, is destabilized by ca. 0.01 eV by the electric field. The corrections due to the electric field are therefore too small to change the qualitative picture of ORR on Ag(111). The free energy of the intermediates is shown at different potentials in Fig. 7 including the corrections due to the potential.3.4.3 ORR on Ni(111). On the clean Ni(111) surface O 2 dissociates easily, so the rate limiting step on Ni at very low coverage is reduction of O* to HO* or HO* to H2O. With a water layer included to stabilize HO* on the surface, we find G0(HO*) = 0.27 eV and G0(O*) = 0.41 eV. All reduction steps are therefore exothermic at U RHE< 0.14 V, which shows that the clean Ni(111) surface is too reactive to be a good ORR catalyst. From the surface Pourbaix diagram we would however expect some oxygen at the surface at higher potentials. This changes the adsorption energies due to adsorbate interactions.To take adsorbate interactions into account we consider ORR involving the surface, as it is the most stable in the potential range from 0.33 to 0.82 V (RHE). The oxygen reduction reaction might occur by adsorption and subsequent reduction of O2 at vacancies in the structure, as shown in Fig. 8. The free energy of the intermediates is shown in Fig. 9.Fig. 8 Intermediates forming at vacancies on. (a) + (b): The HOO*intermediates dissociates into hydroxyl and oxygen.(c) + (d): Adsorbed oxygen completing thestructure. (e) + (f): Adsorbedhydroxyl. (g) + (h): The vacancy. The Ni atoms andthe molecules adsorbed directly at the surface arerepresented by spacefill whereas the surroundingwater molecules, which are not in direct contact withthe Ni surface, are represented by enlarged CPKspheres. Ni atoms are blue, O atoms red and H atomswhite.Fig. 9 The free energy of the intermediates adsorbedat vacancies in . All stepsdecrease the free energy at potentials below 0.14 V. Awater layer has been used to stabilize all intermediatesas shown in Fig. 8. Without the water layer tostabilize O*, the self-consistent ORR potential is 0.26V.We find that spectator O* makes the Ni(111) surface a bit more noble because the binding energies of O* and HO* are decreased by ca.0.1 eV. The HOO* intermediate dissociates on this surface. All steps are exothermic for U RHE < 0.14 V with the rate limiting step being a reduction of O* to HO*.Another possible reaction mechanism on the surface is that O2 adsorbs and is reduced atfree sites on the surface as shown in Fig. 10. Because of repulsion between adsorbed oxygen and H2O adsorbed at the top sites, a H2O–HO* bonding network directly on the surface is rather unstable. We find HO* to be the most stable with the water layer placed above the surface, with hydrogen bonds between the surface and the water layer. A similar structure is found to be the most stable for O* and HOO*, except that it has been possible to include HOO* into the water layer, as it protrudes from the surface. When HOO* is included in the water layer and allowed to relax, hydrogen moves from HOO* to a nearby water molecule, forming aH 3O molecule. The H–O bond in H3O is 1.05 and the H–OO* bond is 1.51 . The free energy of the intermediates is shown in Fig. 11. The free energies of O* and HO* have increased significantly compared to the clean Ni(111) surface, making the surface much more noble. The significantly increased interaction between adsorbates may be rationalized by the fact that the adsorbates bind to the same Ni atoms on the surface. The rate limiting step is the formation of HOO*, which is exothermic at potentials below 0.40 V (RHE). The binding free energies for oxygen and hydroxyl are 1.64 and 0.80 eV respectively. The O* and HO* binding energy on Ni(111) with 1/3 ML O suggests that this surface is a good ORR catalyst comparable to Pt, but the binding energy of HOO* is weaker than that of Pt. Formation of HOO* therefore becomes the rate limiting step on Ni.Fig. 10 Intermediates adsorbed at free sites on. (a) + (b): HOO*. (c) + (d): O*,(e) + (f): HO*, (g) + (f): Thesurface.Fig. 11 Free energy of the intermediates along thereaction path (17) or (19) adsorbed at the free sites on. Because the formation ofHOO* limits the ORR potential, the overpotentialcould be reduced if the binding energies could beincreased. If O2 dissociates easily on the。

第21卷第4期2006年4月地球科学进展ADVANCES I N E ART H SC I ENCEVol.21　No.4Ap r.,2006文章编号:100128166(2006)0420331207地球圈层相互作用中的深海过程和深海记录(I):研究进展与成果汪品先,翦知　,刘志飞(同济大学海洋地质国家重点实验室,上海　200092)摘　要:由五大系统11个实验室组成的项目组,2000—2005年开展了以“地球圈层相互作用”为主题的深海基础研究。

项目以“热带碳循环”作为核心问题,依靠国际大洋钻探和国内“大洋专项”两大支柱,对西太平洋暖池和南海等海区进行深海过程和深海记录的研究,已圆满完成计划任务。

一方面在南海大洋钻探的基础上,围绕热带海洋在地球系统中的作用向纵深发展,在“热带碳循环”研究中取得了原创性的成果;另一方面依托国内大洋专项和国内外合作航次,在深海研究和圈层相互作用上朝横向发展,取得了一系列国际性成果,在我国形成了与国际接轨的深海研究力量。

对该项目的设计和进展做了简单而又全面的阐述,对于古环境研究中取得的突破性进展将另有续篇介绍。

关　键　词:深海研究;地球系统;圈层相互作用中图分类号:P73 文献标识码:A1　任务与设计水深超过2000m的深海占地球表面积的60%,是地球表面研究最弱的部分,也是研究“地球系统”的关键所在。

由深海发现为先导的重大突破,在近几十年来的地球科学中屡见不鲜。

但是,以大陆为主体的我国地球科学受条件限制,基本上还没有延伸到深海大洋。

1999年春南海大洋钻探ODP184航次的成功,使我国进入国际深海基础研究的前沿阵地;我国的“大洋专项”在太平洋多金属结核调查等多年工作积累的基础上,提出了推动我国地球科学发展的新增目标。

这样,我国系统地开展深海基础研究的时机,终于来到。

正是在这种背景下,国家重点基础研究发展计划学科前沿项目“地球圈层相互作用中的深海过程和深海记录”(简称“深海973”),经过2000年10月至2005年9月的5年工作,已经圆满完成。

GPS HIGH PRECISION ORBIT DETERMIANTION SOFTWARE TOOLS (GHOST)M. Wermuth(1), O. Montenbruck(1), T. van Helleputte(2)(1) German Space Operations Center (DLR/GSOC), Oberpfaffenhofen, 82234 Weßling, Germany, Email:martin.wermuth@dlr.de(2) TU Delft, Email: T.VanHelleputte@tudelft.nlABSTRACTThe GHOST “GPS High precision Orbit determination Software Tools” are specifically designed for GPS based orbit determination of LEO (low Earth orbit) satellites and can furthermore process satellite laser ranging measurements for orbit validation purposes. Orbit solutions are based on a dynamical force model comprising Earth gravity, solid-Earth, polar and ocean tides, luni-solar perturbations, atmospheric drag and solar radiation pressure as well as relativistic effects. Remaining imperfections of the force models are compensated by empirical accelerations, which are adjusted along with other parameters in the orbit determination. Both least-squares estimation and Kalman-filtering are supported by dedicated GHOST programs. In addition purely kinematic solutions can also be computed. The GHOST package comprises tools for the analysis of raw GPS observations as well. This allows a detailed performance analysis of spaceborne GPS receivers. The software tools are demonstrated using examples from the TerraSAR-X and TanDEM-X missions.1.INTRODUCTIONThe GPS High precision Orbit determination Software Tools (GHOST) were developed at the German Space Operations Center (DLR/GSOC) in cooperation with TU Delft. They are a set of tools, which share a common library written in C++. The library contains modules for input and output of all common data formats for GPS observations, auxiliary data and physical model parameters. It also provides modules for the mathematical models used in orbit determination and for modelling all physical forces on a satellite.The tools can be classified in analysis tools and tools for the actual precise orbit determination (POD). Additionally a user can implement new tools based on the library modules.The user interface of the GHOST tools consists of human readable and structured input files. These files contain all the parameters and data filenames which have to be set by the user.Although a software package for orbit determination and prediction existed at GSOC, it was decided with the availability of GPS tracking on missions like CHAMP and GRACE, to implement GHOST as a flexible and modular set of tools, which are dedicated to the processing GPS data for LEO orbit determination. GHOST has been used in the orbit determination of numerous missions like CHAMP, GRACE, TerraSAR-X, GIOVE and Proba-2. The newest tool for relative orbit determination between two spacecrafts (FRNS, see section 5.5) was designed using data from the GRACE mission and will be used in the operation of the TanDEM-X, PRISMA and DEOS missions.Due to its modular design the library offers users a convenient way to create their own tools. As all necessary data formats are supported, tools for handling and organizing data are easily implemented. For example many GPS receivers have their own proprietary data format. Hence for each new satellite mission a tool can be created, which converts the receiver data to the standard RINEX format in order to be compatible with the other tools. This is supported by the library modules.Figure 1: GHOST Software Architecture.2.THE TERRASAR-X AND TANDEM-XMISSIONSGHOST has been used and will be used in the preparation and the data processing of numerous satellite missions including the TerraSAR-X and TanDEM-X missions. The examples shown in this paper are taken from actual TerraSAR-X mission data or from simulations done in preparation of the TanDEM-X mission. Hence these two satellite missions are introduced here.TerraSAR-X is a German synthetic aperture radar (SAR) satellite which was launched in June 2007 from Baikonur on a Dnepr rocket and is operated by GSOC/DLR. Its task is to collect radar images of the Earth’s surface. To support the TerraSAR-X navigation needs, the satellite is equipped with two independent GPS receiver systems. While onboard needs as well as orbit determination accuracies for image processing (typically 1m) can be readily met by the MosaicGNSS single-frequency receiver, a decimeter or better positioning accuracy must be achieved for the analysis of repeat-pass interferometry. To support this goal, a high-end dual-frequency IGOR receiver has been contributed by the German GeoForschungsZentrum (GFZ), Potsdam. Since the launch DLR/GSOC is routinely generating precise rapid and science orbit products using the observations of the IGOR receiver.In mid 2010 the TanDEM-X satellite is scheduled for launch. It is an almost identical twin of the TerraSAR-X satellite. Both satellites will fly in a close formation to acquire a digital elevation model (DEM) of the Earth’s surface by stereo radar data takes. Therefore the baseline vector between the two satellites has to be known with an accuracy of 1 mm. In preparation of the TanDEM-X mission, GHOST has been extended to support high precision baseline determination using single or dual-frequency GPS measurements.Figure 2: The TanDEM-X mission. (Image: EADSAstrium)3.THE GHOST LIBRARYThe GHOST library is written in C++ and fully object-oriented. All data objects are mapped into classes and each module contains one class. Classes are provided for data I/O, mathematical routines, physical force models, coordinate frames and transformations, time frames and plot functions.All data formats necessary for POD are supported by the library. Most important are the format for orbit trajectories SP3-c [1] and the ‘Receiver Independent Exchange Format’ RINEX [2] for GPS observations. The SP3 format is used for the input of GPS ephemerides (as those files are usually provided in SP3-c format) and for the output of the POD results. The raw GPS observations are provided in the RINEX format. At the moment the upgrade from RINEX version 2.20 to version 3.00 is ongoing to allow for the use of multi-constellation and multi-antenna files. Other data formats, which are supported, are the Antenna Exchange Format ANTEX [3] containing antenna phase center offsets and variations and the Consolidated Prediction Format CPF [4] used for the prediction of satellite trajectories to network stations. The Consolidated Laser Ranging Data Format CRD [5] which will replace the old Normal Point Data Format is currently being implemented.The library provides basic mathematical functions needed for orbit determination like a module for matrix and vector operations, statistical functions and quaternion algebra. As numerical integration plays a fundamental role in orbit determination, several numerical integration methods for ordinary differential equations are implemented, like the 4th-order Runge-Kutta method and the variable order variable stepsize multistep method of Shampine & Gordon [6].Forces acting on the satellite are computed by physical models including the Earth's harmonic gravity field, gravitational perturbations of the Sun and Moon, solid Earth and ocean tides, solar radiation pressure, atmospheric drag and relativistic effects.In orbit determination several coordinate frames are used. Most important are the inertial frame, the Earth fixed frame, the orbital frame and the spacecraft frame. The transformation between inertial and Earth fixed frame is quite complex as numerous geophysical terms like the Earth orientation parameters are involved. The orbital frame is defined by the position and velocity vectors of a satellite. The axes are oriented radial, tangential (along-track) and normal to the other axes and often denoted as R-T-N. The spacecraft frame is fixed to the mechanical structure of a satellite and used to express instrument coordinates (e.g. the GPS antenna coordinates). It is connect to the other frames via attitude information. The GHOST library contains transformations between all involved frames.Similar to reference frames, also several time scales like UTC and GPS time are involved in orbit determination.A module provides conversions between different time scales.In order to visualize results of the analysis tools or POD results like orbit differences or residuals, the librarycontains a module dedicated to the generation of post script plots.4.ANALYSIS TOOLSThe GHOST package comprises tools for the analysis of raw GPS observations and POD results. This allows a detailed performance analysis of spaceborne GPS receivers in terms of signal strength, signal quality, statistical distribution of observed satellites and hardware dependent biases. The tools can be used either to characterize the flight hardware already prior to the mission or to analyze the performance of in flight data during the mission. An introduction to the most important tools is given here.4.1.EphCmpOne of the most basic but most versatile tools is the ephemeris comparison tool EphCmp. It simply compares an orbit trajectory with a reference orbit and displays the differences graphically (see Fig. 3). The coordinate frame in which the difference is expressed can be selected. In addition a statistic of the differencesis given. It can be used to visualize orbit differences in various scenarios like the comparison of different orbit solutions, the comparison of overlapping orbit arcs or the evaluation of predicted orbits or navigation solutions against precise orbits.Figure 3: Comparison of two overlapping orbit arcsfrom TerraSAR-X POD.4.2.SkyPlotSkyPlot is a tool to visualize the geometrical distribution of observed GPS satellites in the spacecraft frame. A histogram and a timeline of the number of simultaneously tracked satellites are given as well. Hence the tool can be used to detect outages in the tracking.Fig. 4 shows the output of SkyPlot for the MosaicGNSS single-frequency receiver on TerraSAR-X on 2010/04/05. The antenna of the MosaicGNSS receiver is mounted with a tilt of 33° from the zenith direction. This is very well reflected in the geometrical distribution (upper left) of the observed GPS satellites. It can be seen, that mainly satellites in the left hemisphere have been tracked. The histogram (upper right) shows, that – although the receiver has 8 channels – most of the time only 6 satellites (or less) were tracked. The lower plot in Fig 4. shows the number of observed satellites as timeline. It can be seen, that there was a short outage of GPS tracking around 11h. This is useful information for the evaluation of GPS data and the quality of POD results.Figure 4: Distribution of tracked satellites by the MosaicGNSS receiver on TerraSAR-X.4.3.EvalCN0The tool EvalCN0 is used to analyze the tracking sensitivity of GPS receivers. It plots the carrier-to-noise ratio (C/N0) in dependence of elevation.The example shown in Fig. 5 is taken from a pre-flight test of the IGOR receiver on the TanDEM-X satellite. The test was carried out with a GPS signal simulator connected to the satellite in the assembly hall [7]. It can be seen, that under the given test-setup, the IGOR receiver achieves a peak carrier-to-noise density ratio (C/N0) of about 54 dB-Hz in the direct tracking of the L1 C/A code. The C/N0 decreases gradually at lower elevations but is still above 35 dB-Hz near the cut-off elevation of 10°.For the semi-codeless tracking of the encrypted P-code, the C/N0 values S1 and S2 reported by the IGOR receiver on the L1 and L2 frequency show an even stronger decrease towards the lower elevations. The signal strength of the L2 frequency is about 3dB-Hzlower than that for the L1 frequency. To evaluate the semi-codeless tracking quality, the size of the S1-SA and S2-SA difference is shown in the right plot of Fig.5. Both frequencies show an expected almost linear variation compared to SA. The degradation of the signal due to semi-codeless squaring losses increases with lower elevation.Figure 5: Variation of C/N0 with the elevation of the tracked satellite (left) and semi-codeless tracking losses (right) forthe pre-flight test of the IGOR receiver on TanDEM-X.4.4.SLRRESSatellite Laser Ranging (SLR) is an important tool for the evaluation of the quality of GPS-based precise satellite orbits. It is often the only independent source of observations available, with an accuracy good enough to draw conclusions about the accuracy of the precise GPS orbits.SLRRES computes the residual of the observed distance between satellite and laser station versus the computed distance. The residuals are displayed in a plot (see Fig.6). As output daily statistic, station-wide statistics and an overall RMS residual are given.In order to compute the distance between satellite and laser station, the orbit position of the satellite has to be corrected for the offset of the satellites laser retro reflector (LRR) from the center of mass using attitude information. The coordinates of the laser station are taken from a catalogue and have to be transformed to the epoch of the laser observation. This is done by applying corrections for ocean loading, tides, and plate tectonics. Finally the path of the laser beam has to be modelled considering atmospheric refractions and relativistic effects.Figure 6: SLR Residuals of TerraSAR-X POD forMarch 2010.5.PRECISE ORBIT DETERMINATION TOOLSThe POD tools comprise two different fundamental methods of orbit determination. A reduced-dynamic orbit determination is computed by the RDOD tool, while KIPP produces a kinematic orbit solution. Both tools process the carrier phase GPS observations. They need a reference solution to start with. This is provided by the tools SPPLEO and PosFit.SPPLEO generates a coarse navigation solution by processing the pseudorange GPS observations. The satellite position and receiver clock offset are determined in a least squares adjustment. Next PosFit is run to determine a solution in case of data gaps and to smoothen the SPPLEO solution by taking the satellite dynamics into account.5.1.SPPLEOSPPLEO (Single Point Positioning for LEO satellites) is a kinematic least squares estimator for LEO satellites processing pseudorange GPS observations. The program produces a first navigation solution, with data gaps still present. For each epoch, the satellite position and receiver clock offset are determined in a least-squares adjustment.The tool is able to handle both single and dual-frequency observations. In case of single frequency observations, the C1 code is used without ionosphere correction. In case of dual-frequency observations, the ionosphere free linear combination of P1 and P2 code observations is applied. In case range-rate observations are available, it is also possible to estimate velocities. Before the adjustment, the data is screened and edited. Hence the user can choose hard editing limits for the signal to noise ratio of the observation, the elevation and the difference between code and carrier phase observation. In case of the violation of one limit, the observation is rejected. If the number of observations for one epoch is below the limit set by the user, the whole epoch is rejected. After the adjustment the PDOP is computed, and if it exceeds a limit, the epoch is rejected as well. If the a posteriori RMS of the residuals exceeds the threshold set by the user, the observation yielding the highest value is rejected. This is repeated until the RMS is below the threshold or the number of observations is below the limit.The resulting orbit usually contains data gaps and relatively large errors compared to dynamic orbit solutions. Hence the gaps have to be closed and the orbit has to be smoothed by the dynamic filter tool PosFit.5.2.PosFitPosFit is a dynamic ephemeris filter for processing navigation solutions as those produced by SPPLEO. This is done by an iterated weighted batch least squares estimator with a priori information. The batch filter estimates the initial state vector, drag and solar radiation coefficients. In addition to those model parameters, empirical accelerations are estimated. One empirical parameter is determined for each of the three orthogonal components of the orbital frame (radial, along-track and cross-track) for an interval set by the user. The parameters are assumed to be uncorrelated over those intervals.The positions of the input navigation solution are introduced as pseudo-observations. The filter is fitting an integrated orbit to the positions of the input orbit in an iterative process. In order to obtain initial values for the first iteration, Keplerian elements are computed from the first two positions of the input orbit. All forces that act on the satellite (like atmospheric drag, solar radiation pressure, tides, maneuvers…) are modelled and applied in the integration. Due to imperfections in the force models, the empirical accelerations are introduced, to give the integrated orbit more degrees of freedom, to fit to the observations. The empirical acceleration parameters are estimated in the least squares adjustment together with the initial state vector and model parameters. The partial derivatives of the observations w.r.t. the unknown parameters are obtained by integration of the variational equations along the orbit (for details see [8]). The result of PosFit is a continuous and smooth orbit without data gaps in SP3-c format. It can serve as reference orbit for RDOD and KIPP.Figure 7 displays the graphical output of PosFit. The three upper graphs show the residuals after the adjustment in the three components of the orbital frame. In this example, which is taken from a 30h POD arc of the TerraSAR-X mission, the RMS of the residuals lies between 0.5m and 1.5m. This mainly shows the dispersion of positions of the SPPLEO solution. The three lower graphs show the estimated empirical accelerations in the three components of the orbital frame.Figure 7: Graphical output of PosFit Tool forTerraSAR-X POD.5.3.RDODRDOD is a reduced dynamic orbit determination tool for LEO satellites processing carrier phase GPS observations. This is also done by an iterated weighted batch least squares estimator with a priori information. Similar to PosFit, RDOD estimates the initial state vector, drag and solar radiation coefficients and empirical accelerations. Contrary to PosFit, where the positions of a reference orbit are used as pseude-observations, RDOD directly uses the GPS pseudorange and carrier phase observations. Nevertheless a continuous reference orbit – normally computed by PosFit – is required by RDOD for data editing and for obtaining initial conditions for the first iteration.The tool is able to handle both single and dual-frequency data. In case of single frequency observations, the GRAPHIC (Group and Phase Ionospheric Correction) combination of C/A code and L1 carrier-phase is used as observation. In case of dual-frequency observations, the ionosphere free linear combination of L1 and L2 carrier phase observations is applied.The data editing is crucial to the quality of the results. Hence the data is also screened and edited by RDOD using limits specified by the user. If the signal-to-noise ratio of the observation exceeds the limit, or the elevation is below a cut-off elevation, the observation is rejected. This is done if no GPS ephemerides and clock information is available for that observation. Next outliers in the code and carrier-phase observations are detected. This is done comparing the observations to modelled observations using the reference orbit.The RDOD filter is fitting an integrated orbit to the carrier-phase observations. This is done in a similar way as in PosFit, considering all forces on the satellite and estimating empirical acceleration parameters. But while PosFit uses absolute positions as observations, the carrier-phase observations used in RDOD contain an unknown ambiguity. The ambiguity is considered to be constant over one pass – the time span in which a GPS satellite is tracked continuously. Hence one unknown ambiguity parameter per pass is added to the adjustment.The graphical output of RDOD shows the residuals of the code and carrier phase observations (see Fig. 8). This is an important tool for a fast quality analysis and for detecting systematic errors.Figure 8: Output of RDOD tool for TerraSAR-X POD.Both the antennas of the GPS satellites and the GPS antennas of spaceborne receivers show variations of the phase center dependent on azimuth and elevation of the signal path. It is necessary to model those variations in order to obtain POD results with highest quality. GHOST does not only offer the possibility to use phase center variation maps given in ANTEX format. It also offers the possibility to estimate such phase center variation patterns, and thus can be employed for the in-flight calibration of flight hardware. This was done for the GPS on TerraSAR-X as shown in Fig. 9. The figure shows the phase center variation pattern for the main POD antenna of the IGOR receiver on TerraSAR-X. It was estimated from 30 days of flight data and needs to be applied to carrier phase observations in addition to a pattern which was determined for the antenna type by ground tests (for details cf. [9]).Figure 9: Phase Center Variation Pattern for the MainPOD Antenna on TerraSAR-X.5.4. KIPPKIPP (Kinematic Point Positioning) is a kinematic least squares estimator for LEO satellites. Similar to RDOD carrier phase observations are processed. But in contrast to RDOD no dynamic models are employed, and only the GPS observations are used for orbit determination. For each epoch, the satellite position and receiver clock offset are determined in a weighted least squares adjustment. KIPP also requires a continuous reference orbit, as that computed by PosFit.Like RDOD, the KIPP tool is able to handle both single and dual-frequency data. In case of single frequency observations, the GRAPHIC (Group and Phase Ionospheric Correction) combination of C/A code and L 1 carrier-phase is used as observation. In case of dual-frequency observations, the ionosphere free linear combination of L 1 and L 2 carrier phase observations is applied.5.5. FRNSThe Filter for Relative Navigation of Satellites (FRNS) is designed to perform relative orbit determination between two LEO spacecrafts. This is done using an extended Kalman filter described in [10]. The concept is to achieve a higher accuracy for the relative orbit between two spacecrafts by making use of differenced GPS observations, than by simply differencing two independent POD results. FRNS requires a continuous reference orbit for both spacecrafts, such as computed by RDOD. It then keeps the orbit of one spacecraft fixed, determines the relative orbit between the two spacecrafts and adds it to the positions of the first spacecraft. As result a SP3-c file containing the orbit of both spacecrafts is obtained. The tool is able to process both single and dual-frequency observations.The FRNS tool was developed using data from the GRACE mission and will be applied on a routine basis for the TanDEM-X mission. In contrast to TanDEM-X, GRACE consists of two spacecrafts, which follow each other on a similar orbit with about 200 km distance. The distance between the two spacecrafts is measured by a K-band link, which is considered to be at least one order of magnitude more accurate than GPS observations. Hence the K-band observations can be used to assess the accuracy of the relative navigation results – with the limitation, that the K-band observations only reflect the along-track component, and contain an unknown bias. The differences between a GRACE relative navigation solution and K-band observations are shown in Fig. 10. The standard deviation is about 0.7 mm. As the TanDEM-X mission uses GPS receivers which are follow-on models of those used on GRACE, and the distance between the spacecrafts is less than 1 km, one can expect that the quality of the relative orbit determination will be on the same level of accuracy or even better.Figure 10: Comparison of GRACE relative navigation solution with K-band observations.6.REFERENCES1. Hilla S. (2002). The Extended Standard Product 3Orbit Format (SP3-c, National Geodetic Survey,National Ocean Service, NOAA.2. Gurtner W., Estey L. (2007). The ReceiverIndependent Exchange Format Version 3.00,Astronomical Institute University of Bern.3. Rothacher M., Schmid R. (2006). ANTEX: TheAntenna Exchange Format Version 1.3,Forschungseinrichtung Satellitengeodäsie TUMünchen.4. Rickfels R. L. (2006). Consolidated Laser RangingPrediction Format Version 1.01, The University of Texas at Austin/ Center for Space Research.5. Rickfels R. L. (2009). Consolidated Laser RangingData Format (CRD) Version 1.01, The Universityof Texas at Austin/ Center for Space Research. 6. Shampine G. (1975). Computer solution of OrdinaryDifferential Equations, Freeman and Comp., SanFrancisco.7. Wermuth M. (2009). Integrated GPS Simulator Test,TanDEM-X G/S-S/S Technical Validation Report,Volume 15: Assembly AS-1515, DLROberpfaffenhofen.8. Montenbruck O., Gill E. (2000). Satellite Orbits –Models, Methods and Applications, Springer-Verlag, Berlin, Heidelberg, New York.9. Montenbruck O. Garcia-Fernandez M., Yoon Y.,Schön S., Jäggi A..; Antenna Phase CenterCalibration for Precise Positioning of LEOSatellites; GPS Solutions (2008). DOI10.1007/s10291-008-0094-z.10. Kroes R. (2006). Precise Relative Positioning ofFormation Flying Spacecraft using GPS, PhDThesis, TU Delft.。

67找矿技术P rospecting technology矿区三维地质建模方法研究及深部综合找矿预测王霄霄（河北省地质矿产勘查开发局第一地质大队，河北　邯郸　０５６００１）摘 要：本论文将从矿区三维地质建模方法、三维可视化与分析技术、地质信息集成与分析、模型与算法应用，以及深部矿产资源评价与优选等几个方面进行探讨。

通过对这些关键环节的详细分析和研究，旨在全面展示深部综合找矿预测的理论基础、方法体系以及应用前景，为矿业领域的科学研究和实际应用提供有益的参考和借鉴。

关键词：矿区；三维地质；找矿预测中图分类号：Ｐ６２８ 文献标识码：Ａ 文章编号：１００２－５０６５（２０２３）１７－００６７－３Research on 3D Geological Modeling Methods and Deep Comprehensive Prospecting Prediction in Mining AreasWANG Xiao-xiao（Ｔｈｅ　Ｆｉｒｓｔ　Ｇｅｏｌｏｇｉｃａｌ　Ｂｒｉｇａｄｅ　ｏｆ　ｔｈｅ　Ｇｅｏｌｏｇｉｃａｌ　ａｎｄ　Ｍｉｎｅｒａｌ　Ｅｘｐｌｏｒａｔｉｏｎ　ａｎｄ　Ｄｅｖｅｌｏｐｍｅｎｔ　Ｂｕｒｅａｕ　ｏｆ　Ｈｅｂｅｉ　Ｐｒｏｖｉｎｃｅ，Ｈａｎｄａｎ　０５６００１，Ｃｈｉｎａ）Abstract: Ｔｈｉｓ　ｐａｐｅｒ　ｗｉｌｌ　ｅｘｐｌｏｒｅ　ｓｅｖｅｒａｌ　ａｓｐｅｃｔｓ　ｏｆ　ｍｉｎｉｎｇ　ａｒｅａ　３Ｄ　ｇｅｏｌｏｇｉｃａｌ　ｍｏｄｅｌｉｎｇ　ｍｅｔｈｏｄｓ，　３Ｄ　ｖｉｓｕａｌｉｚａｔｉｏｎ　ａｎｄ　ａｎａｌｙｓｉｓ　ｔｅｃｈｎｉｑｕｅｓ，　ｇｅｏｌｏｇｉｃａｌ　ｉｎｆｏｒｍａｔｉｏｎ　ｉｎｔｅｇｒａｔｉｏｎ　ａｎｄ　ａｎａｌｙｓｉｓ，　ｍｏｄｅｌ　ａｎｄ　ａｌｇｏｒｉｔｈｍ　ａｐｐｌｉｃａｔｉｏｎｓ，　ａｎｄ　ｄｅｅｐ　ｍｉｎｅｒａｌ　ｒｅｓｏｕｒｃｅ　ｅｖａｌｕａｔｉｏｎ　ａｎｄ　ｏｐｔｉｍｉｚａｔｉｏｎ．　Ｔｈｒｏｕｇｈ　ｄｅｔａｉｌｅｄ　ａｎａｌｙｓｉｓ　ａｎｄ　ｒｅｓｅａｒｃｈ　ｏｎ　ｔｈｅｓｅ　ｋｅｙ　ｌｉｎｋｓ，　ｔｈｅ　ａｉｍ　ｉｓ　ｔｏ　ｃｏｍｐｒｅｈｅｎｓｉｖｅｌｙ　ｄｅｍｏｎｓｔｒａｔｅ　ｔｈｅ　ｔｈｅｏｒｅｔｉｃａｌ　ｂａｓｉｓ，　ｍｅｔｈｏｄｏｌｏｇｉｃａｌ　ｓｙｓｔｅｍ，　ａｎｄ　ａｐｐｌｉｃａｔｉｏｎ　ｐｒｏｓｐｅｃｔｓ　ｏｆ　ｄｅｅｐ　ｃｏｍｐｒｅｈｅｎｓｉｖｅ　ｏｒｅ　ｅｘｐｌｏｒａｔｉｏｎ　ｐｒｅｄｉｃｔｉｏｎ，　ｐｒｏｖｉｄｉｎｇ　ｂｅｎｅｆｉｃｉａｌ　ｒｅｆｅｒｅｎｃｅｓ　ａｎｄ　ｒｅｆｅｒｅｎｃｅｓ　ｆｏｒ　ｓｃｉｅｎｔｉｆｉｃ　ｒｅｓｅａｒｃｈ　ａｎｄ　ｐｒａｃｔｉｃａｌ　ａｐｐｌｉｃａｔｉｏｎｓ　ｉｎ　ｔｈｅ　ｍｉｎｉｎｇ　ｆｉｅｌｄ．Keywords: ｍｉｎｉｎｇ　ａｒｅａ；　３Ｄ　ｇｅｏｌｏｇｙ；　Ｐｒｏｓｐｅｃｔｉｎｇ　ｐｒｅｄｉｃｔｉｏｎ收稿日期：２０２３－０６作者简介：王霄霄，女，生于１９９２年，汉族，河北邯郸人，本科，学士学位，矿产地质工程师，研究方向：矿产地质勘查，三维地质建模，地质大数据。

Urban growth and uninsured rural risk:Booming towns in bust timesSteven Poelhekke ⁎,1De Nederlandsche Bank,Research Department,Westeinde 1,1017ZN,Amsterdam,The Netherlandsa b s t r a c ta r t i c l e i n f o Article history:Received 8January 2008Received in revised form 25November 2009Accepted 30July 2010JEL Classi fication:O1R11R23R51D81Keywords:Urbanization RiskNatural resources VolatilityRural –urban migrationRapid urbanization also happens when economic growth and urban job creation are absent,such as in Africa and Latin America during the eighties.Why do some countries urbanize faster while having worse economic growth?This paper finds that higher aggregate agricultural risk induces rural –urban migration,providing an additional channel to explain the urbanization trend.Uninsurable expected risk will lead to rural –urban migration as a form of ex-ante insurance if households are liquidity constrained and cannot overcome adverse shocks.The effect is robust to controlling for the traditional view of urbanization driven by industrialization,and to several alternative explanations such as government spending.©2010Elsevier B.V.All rights reserved.1.IntroductionRapid urbanization even happens when economic growth and urban job creation are absent,such as for example in Africa and Latin America during the eighties.Fig.1shows that growth in GDP per capita slowed signi ficantly or even reversed,while the rate of urbanization continued at a fast pace.Without growth to create jobs or higher wages in cities (such as in East Asia)it seems puzzling that so many rural dwellers choose to become urban inhabitants.Most people end up in slums which do not necessarily offer better living conditions than rural areas for a given income (UN-Habitat,2006).The long time period and crowding should lower the expected income gain from moving to the city.It seems that migration flows are larger and more persistent than the classic Harris-Todaro (1970)model can explain.Big city lights are not always bright.Why do some countries urbanize faster than others while having worse economic growth?If pull-factors are absent,can the observed persistent trend in urbanization be explained by push factors?This paper studies the uninsurable risk involved in dependence on agriculture and natural resource production to explain urbanization occurring even under negative growth.More speci fically,it models rural –urban migration as an insurance mechanism and provides empirical evidence that agricultural risk is an additional explanation for urbanization.Table 1hints at this hypothesis.It ranks regions according to how much faster cities grow with respect to overall population growth (fourth column).On top are Asian regions with fast economic growth,but also regions without much economic growth such as Sub-Saharan Africa.The last three columns show that a third of value added comes from agriculture,exposing a big share of the economy to very high volatility.Such risk became even larger after 1980,the period after which urbanization and economic growth started to diverge.While economic growth is a big driver of urbanization,it may well be that risk offers an additional channel.Unless financial instruments are available to smooth consumption and ‘ride out the bad times ’,the diversity of income sources that cities offer may turn rural –urbanJournal of Development Economics 96(2011)461–475⁎Tel.:+31205243877;fax:+31205242500.E-mail address:steven.poelhekke@dnb.nl .URL:http://www.dnb.nl/onderzoek/onderzoekers/persoonlijke-paginas/auto185449.jsp .1Also af filiated with the Oxford Centre for the Analysis of Resource Rich Economies (OxCarre),Department of Economics,University of Oxford,and CESifo,Munich.Part of this paper was written while the author was af filiated with the European University Institute,Florence.I thank the editor,two anonymous referees,Giancarlo Corsetti,Patrick Eozenou,Marcel Fafchamps,Harry Garretsen,Jan Willem Gunning,Bas Jacobs,Pierre Lafourcade,Rick van der Ploeg,Fabio Schiantarelli,William Strange,Martin Zagler,seminar participants at the European Bank for Reconstruction and Development,and NARSC 2008participants for helpful discussions and suggestions,and INFER 2009for awarding this article the Best Paper Award.All errors are my own.Views expressed are those of the author and do not necessarily re flect of ficial positions of De NederlandscheBank.0304-3878/$–see front matter ©2010Elsevier B.V.All rights reserved.doi:10.1016/j.jdeveco.2010.07.007Contents lists available at ScienceDirectJournal of Development Economicsj o ur n a l h o m e p a ge :ww w.e l s ev i e r.c o m/l o c a t e /d eve cmigration into a crude insurance device for agricultural risk.Such financial services to insure against shocks are relatively absent in rural areas (Collier and Gunning,1999),while natural resources and food prices show very volatile behavior on the world market (Deaton,1999),much more so than (urban)manufactures.Fig.2shows cumulative density functions for four resources,and general price indices for manufacturing products and the world.The x -axis shows the size of yearly standard deviations in %monthly in flation.It is clear that a country with a high dependence on,for example,food products faces much more volatile prices than (OECD)countries which typically trade in manufactures.Food prices are as volatile as the prices of ores and metals,prompting large price stabilization schemes in the 1970s (Newbery and Stiglitz,1981).2Several country examples broadly support this idea.For example,between 1980and 1985Bolivia's urban population grew more than twice as fast as population in total,while manufacturing declined with 5%yearly on average.The same period shows that agricultural risk was 13%:twice as high as urban (manufacturing)risk.Haiti (1985–90)and Paraguay (1995–00)show similar patterns,although less extreme.Also Algeria,Zimbabwe and Mozambique (1990–95)had severe rural risk (resp.11%,18%and 16%standard deviation of yearly growth).This may explain why urban areas increased twice as fast in size as the overall population,even though the urban manufacturing sector declined.The next section presents the main hypotheses in the light of the existing literature and describes possible alternative explanations for continued urban growth.Section 3derives a model of migration as an ex-ante response to risk,leading to the econometric speci fication presented in Section 4.Section 5addresses concerns for endogeneity,and Section 6describes the results.Section 7presents an out-of-sample forecasting exercise to see if the observed general trend can be explained by the model.Section 8concludes.2.Risk,insurance and alternative explanationsWe build on the existing empirical literature on urbanization,which implicitly relies on economic growth of the urban manufac-turing sector to generate a rural –urban income gap and sectoraltransition.Without economic growth,other explanations are needed to explain persistent urbanization,such as risk.For example,if technological progress drives city growth through the industrial sector as simulated in Kelley and Williamson (1984),then urban income may continue to outpace rural income.A resulting positive urban to rural income ratio drives urban population growth in Harris and Todaro (1970)and Brueckner (1990).Moomaw and Shatter (1996)estimate that countries with a higher share of labor in industry are more urbanized,supporting the view that urbanization takes place as a country industrializes.The growth channel lies at the heart of for example China's rapid urbanization (Deng et al.,2008).3Two alternative empirical explanations which do not rely on economic growth are rainfall and government policy.Barrios et al.(2006)use rainfall data to show that low rainfall (low agricultural productivity)is associated with a higher contemporary level of urbanization in Africa.Fay and Opal (2000)and Davis and Henderson (2003)identify government policy resulting in ‘urban bias ’and arti ficially high urban wages as an important cause for high levels of urbanization,in combination with urban poverty.For example,planned economies such as China tend to restrict migration,and policy may affect the sectoral composition through for example import substitution programs that favor cities.Fields (1975)suggests that government involvement may also introduce government jobs or subsidies as winning tickets to the ‘lottery ’for formal urban employment.In that case rural workers may choose to migrate even if living standards are lower in the informal urban sector,as long as it offers the possibility to win a formal job.Becker and Morrison (1988)find empirical evidence for this link,although they are limited to a cross-section of ernment spending may keep cities attractive if it can compensate for lagging job growth under continued rural –urban migration,which would otherwise signi ficantly lower the probability of winning.The possibility that risk provides a third alternative channel is supported by Stark and Levhari (1982)who already noted that,at the micro level,uncertainty and risk can be a motive to migrate (for some family members)but it has not been applied to explain the general trend of urbanization.Daveri and Faini (1999)have estimated this2Crude oil has a different scaling because of the oil crisis.The series ‘OECD Manuf.PPI ’(manufacturing producer price index)starts in 1982.3Cities may also be engines of growth themselves if the bene fits of increasing returns and agglomeration economies outweigh the costs of crowding.See for example Duranton and Puga (2004).10%20%30%40%50%60%70%80%7891011196019802000196019802000196019802000196019802000S u b −S a h a r a n A f r i c aE a s t A s i a & P a c .L . A m e r i c a & C a r i b .W . E u r o p elog GDP/capita (left scale)% Urban (right scale)yearFig.1.Growth and urbanization.462S.Poelhekke /Journal of Development Economics 96(2011)461–475motive for Italian migrants within Italy and internationally and conclude that risk is a significant determinant,driven by risk aversion.4Migration as an ex-ante response to risk seems reasonable, given that the literature suggests that self-insurance(via savings and informal insurance mechanisms)typically only partially succeeds,and against idiosyncratic shocks at best(Besley,1995;Townsend,1994, 1995;Bardhan and Udry,1999).Informal mechanisms require strong information and enforcement institutions within the community (Udry,1990)and transfers across time are limited because of credit constraints(Rosenzweig and Binswanger,1993).Accumulation of buffer stocks(often in the form of bullocks)may also be used to smooth consumption(Deaton,1991),but this can also affect production(Rosenzweig and Wolpin,1993)and is thus a sub-optimal insurance method.Households will have to resort to ex-ante strategies to deal with risk.Elbers et al.(2005)use micro data to quantify the ex-post and ex-ante effects of risk on capital accumu-lation;theyfind that two-thirds of the detrimental effect of risk is due to the ex-ante type which influences households'behavioral deci-sions.5For example,Giles(2006)shows that rural households in China use off-farm labor markets to reduce exposure to ex-ante risk and to increase ex-post smoothing opportunities.This strategy only became possible after(temporary)migration to urban areas became legal in1988,allowing families to diversify income.Similar evidence comes from India where households are more likely to participate in the labor market in regions with higher rainfall risk(Rose,2001). Households obtain additional insurance by letting one or more household members migrate to other areas with less or uncorrelated income risk,expecting remittances to supplement total household income,as happens in Thailand(Paulson,2003)where the destination of choice is the city of Bangkok.The migration choice is essentially a portfolio choice where households decide on the distance(for example home or foreign destination)and on which(and how many)household members to send(Azam and Gubert,2006).Cities,the destination location,may in addition provide relatively better access tofinancial markets.Evidence for this comes from Taylor et al.(1996),who document the failure of local rural credit and risk markets.Conning and Udry(2007)also report the extend of imperfections in rural capital markets,for example pointing out that microfinance has mostly focused on urban or non-farm activities. Commercialfinancial intermediaries are mostly confined to urban areas with more opportunity to diversify their portfolio.Access may not be equally distributed within cities where a dual economy exists, the formal and the informal one(see i.e.Temple,2005).Although slum dwellers live closer to a concentration offinancial services than rural inhabitants they mayfind it harder to smooth consumption than official residents because they cannot provide collateral or a credit history.On the other hand,cities are centers of trade and political power and offer more diverse sources of income than rural areas.At least they offer a chance of improving living conditions.Incomplete markets affect not only households'(ex-post)income, but may also affect their(ex-ante)behavior,possibly resulting in migration to cities.We add to the literature by focusing on house-holds'exposure to aggregate rural risk which is uninsurable by the local informal ex-post methods and hypothesize that periods with more risk induce more migration and urban growth.3.Model of ex-ante risk insuranceThe theoretical reason to look at risk as an explanation for urbanization is derived as follows.We assume that households pool all income from their members.A family may choose to invest in the migration of one or more of its members(workers)who derive income from employment either in the urban or the rural area.6 Household income may then be supplemented with remittances from any members employed outside the home rural area as a return on family investment to migration.The choice of location is directly tied to sectors of the economy.The rural area only offers agriculture and other natural resource production.It is risky because income depends on nature,such as rainfall,and on demand and price shocks for natural products.Moreover,the rural area is a single sector economy with little scope for diversification.The urban area instead consists of formal manufacturing jobs,an informal sector,and the possibility to obtain formal government jobs or ernment income is often spent in cities,closer to the government's supporters(Davis and Henderson,2003).Migration will be influenced by the probability of obtaining formal income in addition to the hypothesized risk channel.7Furthermore,the urban sector typically has better access tofinancial services to insure against shocks in addition to diversification opportunities.It is therefore reasonable to assume that production and employment in the rural area is inherently more risky than employment in the urban area even though its return is not necessarily higher.Some periods are more volatile than others and a country's development over time may change its dependence on natural resources and its ability to cope with external shocks.A time dimension is therefore also important.The goal of this model is to analyze the effect of a risk differential on workers'location choice.A representative household faces a choice to divide its members over two areas which simultaneously requires a choice between sectors.8Household face a liquidity constraint every period because we assume thatfinancial markets are underdeveloped,especially in rural ck of collateral or a credit history(and incomplete4Dustmann(1997)similarly models the duration of international migration as it is determined by risk at home and abroad and the(intertemporal)covariance of labormarket shocks in addition to a wage differential.For example,a temporary migrant may diversify risk if the covariance of shocks is negative.5Ex-ante insurance takes the form of conservative investment decisions,such as postponing adoption of new risky technology,crop diversification,crops of lower yields but faster growth cycles,diversifying family members among different income activities or sharecropping.It is also related to remittances and risk diversification by means of assigning family members to work in a different area,country or sector such as in Stark and Lucas(1988).6Migration as an investment decision goes back to Sjaastad(1962).7Banerjee and Kanbur(1981)use unemployment rates and inequality as measures of risk in a cross-section of Indian regions.Here risk is formulated as shocks to income. 8This section builds heavily on a standard risk and insurance model with a precautionary savings motive(as in Mirrlees,1965;Deaton,1991),see Bardhan and Udry(1999).In these models households make choices between different investments with different risk and return.Here we add the possibility that location is a choice and that each location promises a different stream of income.Table1Economic and urban growth by region(yearly%,1970–2000).Urbanization%Urban–total%pop.growth GDP/cap.growthVolatilityofagriculture%Agri.19702000ΔMean Mean b'80≥'80MeanSub-SaharanAfrica19.036.117.1 2.940.769.09.733.0South Asia12.821.48.6 2.69 2.56 4.8 4.937.4 East Asia&Pacific38.651.312.7 1.40 2.598.87.319.8Middle East&N.Africa60.076.316.3 1.16 2.0812.612.78.6Latin America&Carib.49.763.613.90.97 1.48 5.7 6.112.8WesternEurope65.674.79.10.51 2.54 6.5 6.2 5.9Eastern Europe&47.758.510.80.44 1.77 6.69.919.1 North America74.779.3 4.60.21 2.24 5.6 6.7 2.5 Note:Ordered on column4.Means are calculated as the within-region cross-country-time unweighed average of5-year average growth rates.463S.Poelhekke/Journal of Development Economics96(2011)461–475markets)prohibits borrowing such that in every period the value of assets A t plus expected income y t should be larger than consumption c t .A t +y t −c t ≥0ð1ÞIncome y t depends on the previous period ‘portfolio ’choice of the household migration decision which corresponds to a choice of location and hence of sector.9The household decides to let a share z t −1of its members migrate.If z =1all members will move to the urban area,and if z =0all remain in the rural area.Because each location has a perfectly competitive sector people can always find employment in the city:either formal or informal.10We assume additionally that each location faces a different degree of aggregate (multiplicative)income risk.Income is therefore a function of exogenous shocks taking place in the rural sector t ∼N (1,σ 2)with unit mean and variance σ 2and in the urban sector ηt ∼N (1,ση2)(where the shock ηis a combination of shocks to the manufacturing,the informal and the government sector),with known joint density function f ( ,η).However,migration is costly because it is costly to obtain information about possible destinations (which depends on distance and relates to the cost of searching for a job)and it is more dif ficult to remit income over greater distances.Migrants will have to pay for higher cost of living in urban areas and may have invested in education (which has a higher return outside the rural sector).Azamand Gubert (2006)document that families send their most promising members away.In Lucas (2004),migrant workers may choose to forgo income to allow them to spend more time searching for a better job or to acquire skills.A family with a migrated member will have to spent part of the urban income on (re)paying these household income is therefore a function of the share of family members in the rural area (1−z ),their (stochastic)wages w R ,and (stochastic)income from urban members w U ηnet of costs κ:y t z t −1; t ;ηt ðÞ=1−z t −1ðÞw R ;t t +z t −1w U ηt −κðÞð2ÞThe expected urban wage w U is a combination of formal,informal and government income,weighted by the probabilities of obtaining such income.The marginal income bene fit from letting more household members migrate is increasing in the urban shock:∂2y t /∂ηt ∂z t N 0and decreasing in the rural shock:∂2y t /∂ t ∂z t b 0.It is therefore crucial to form expectations on the relative riskiness of both sectors to make an optimal migration and location choice.Households maximize a discounted flow of expected utility from consumption subject to their budget constraint,where r t is the rate of return on assets:max c t ;z t+1E t ∑Tτ−tβτ−tu c τðÞð3Þs :t :A t +1=1+r t ðÞA t +y t −c t ðÞð4ÞHouseholds are risk averse so u ′N 0,u ″b 0and lim x →0u ′(x )=+∞.Households therefore aim to smooth consumption over time.The corresponding period t value function is given by:V t A t +y t ðÞ=max c tf u c t ðÞ+βE t V t+1½1+r t ðÞA t +y t −c t ðÞ+y z t ; t+1;ηt +1ÀÁ+λt A t +y t −c t ðÞgð5Þwhere λt is the multiplier associated with the liquidity constraint.The income shock is a combination of rural and urban income shocks if.2.4.6.81C u m u l a t i v e P r o b a b i l i t y0246810Agricultural Raw Materials 0.2.4.6.81C u m u l a t i v e P r o b a b i l i t y246810Foods 0.2.4.6.81C u m u l a t i v e P r o b a b i l i t y246810OECD Manuf. PPI0.2.4.6.81C u m u l a t i v e P r o b a b i l i t y246810Ores & Metals 0.2.4.6.81C u m u l a t i v e P r o b a b i l i t y20406080100Crude Petroleum 0.2.4.6.81C u m u l a t i v e P r o b a b i l i t y246810Worldc.d.f NormalFig.2.Densities of yearly standard deviation of monthly price index in flation,1970−2003.9We could also make this choice depend on distance.Fafchamps and Shilpi (2008)show for example how spatial isolation leads to lower subjective welfare,which might be consistent with an increased need for additional income sources.Taking into account distance,the household chooses not only the share of members to migrate but also the distance to the nearest city and thus access to an external market offering more means of diversi fication,but this would not change our main results.See for example Brueckner and Zenou (1999)for an urbanization model with a land market.10We abstract from any unemployment bene fits.464S.Poelhekke /Journal of Development Economics 96(2011)461–4750b z t b 1.The current value of assets and income equals the maximum of current utility from consumption plus the discounted value of future assets and income.Maximization yields:u ′c t ðÞ=βE t V ′t +11+r t ðÞA t +y t −c t ðÞ+y z t ; t+1;ηt +1ÀÁÂÃ+λtð6ÞThe household also chooses the location one period before as a form of ex-ante risk ing the envelop theorem we have:E t −1dV ′t ⋅ðÞdz t −1=E t −1u ′c t ðÞ∂y ∂z t −1=0ð7Þ⇔E t −1β1+r ðÞV ′t +1⋅ðÞ+λt hi ∂y ∂z t −1=0ð8ÞIf the liquidity constraint 1never binds (λt =0),the location is chosen such that there is no incentive to move:E t −1V ′t+1⋅ðÞ∂yt −1=0ð9Þbut if it does bind and λt N 0households chose z t −1such that (rewriting Eq.(8))E t −1β1+r ðÞV ′t+1⋅ðÞ∂y t −1=−E t −1λt ∂yt −10ð10ÞThe last inequality holds only when the liquidity constraint binds,which is when either shock is negative (meaning smaller than 1)but not equal to each other.We look at the short run effects of large shocks rather than the long run effects when shocks are expected to be at their mean of 1.Volatility is then interpreted as a higher chance of receiving a shock that is so large that all savings are wiped out.Households want to avoid being put in that situation.If most family members live and work in the rural area and the shock is suf ficiently bad ( t b b 1),such that the liquidity constraint binds,we have that ∂y /∂z t −1N 0.Households could then improve utility by moving somefamily members to the city:−E t −1λt ∂yt −1b 0.Conversely,if z t −1is closer to one (a higher share of family members in the urban sector)and ηt bb 1(bad urban year)we have that ∂y /∂z t −1b 0and thus that−E t −1λt∂y∂z t −1N 0.In that case the rural sector would be better ifhouseholds expect the constraint to bind.If both shocks are of equal size they cancel,and we are back in the situation of Eq.(9)where there is no incentive to move.Four insights arise:the more likely it is that the liquidity constraint binds,the more likely households will be able to improve consump-tion and utility by letting family members migrate to the area where they expect shocks to be smaller.Secondly,if the variance of shocks to the rural sector is larger than the variance of shocks to the urban sector,then rural households are more likely to suffer an adverse shock that is large enough to hit the liquidity constraint.This increases pressure to let family members migrate to the urban area.Without modern sector job growth this leads to an increase in the informal sector (for given wages)and a lower expected urban wage,which is the balancing force.The government sector can cushion the urban area against shocks but may also provide a direct source of income.Thirdly,if both shocks are equal in size we have that ∂y /∂z t −1=0.Then no improvement can be gained from migrating,even if households hit a liquidity constraint.This is the case if the covariance of both shocks equals stly,higher cost of migration lower theattractiveness of cities because more of the expected flow of urban income has to be spent on migration.114.Estimating the effect of risk on urban growthThe risk and migration model holds for a representative household.Aggregating over all households implies that the growth rate of the rural population R equals natural rural population growth ΔR n minus any rural to urban migration m :ΔR ≡logR t +1−logR t =ΔR n −m .Since growth of the total population P is a weighted function g of the growth rates of R and the urban population U ,it follows that:ΔU =g −1(ΔP )−ΔR n +m .In equilibrium no families have any incentive to move (∂y /∂z t −1=0)and the rate of migration is zero.The national urban population U it for country i and year t (measured with five-year intervals)is then given by the speci fication in levels in Davis and Henderson (2003)based on Brueckner (1990):logU it =δ0logP it +δ1X it +γi +e itð11ÞX it includes measures capturing the country's state of development,rural –urban differences in public service provision,democracy and infrastructure which may affect migration costs,urban cost of living,and ideally the expected urban and rural wages.The γi capture fixed unobserved country characteristics and the e it is the error term.However,we will not assume that countries are in equilibrium every 5years and rather focus on changes when rural –urban migration m ≠0.This happens when the expected rural and urban wages do not equal,and in risky periods when the liquidity constraint binds.In that case ∂y /∂z t −1≠0.The core message is that risk is a strong destabilizing force which in fluences the speed of urbanization rather than the level.Volatile countries are not necessarily more or less urbanized,but volatile periods will induce more migration.Rural –urban migration is therefore a positive function of urban wages and rural risk,and a negative function of rural wages,migration costs and urban risk:m t =m w R t ;w U ηt ;κ;σ ;t ;ση;tð12ÞThe main econometric model is therefore given by:logU i ;t+1−logU i ;t =β0logP i ;t+1−logP it+β1log U it =P it ðÞ+β2X w R ;w U ηðÞit +β3σ ;it +β4ση;it +γi +e itð13Þfor t =1970,1975,...,where X it also contains control variables such asthe growth of the economy.12Overall growth in average GDP per capita captures economic development and the transition process from an agricultural to an industrial economy (including changes to the urban –rural wage gap if growth,as is often assumed,originates in cities).Positive income growth creates more financial leeway and means that people are less vulnerable to shocks.Changes in wage prospects can also be approximated by including separate growth rates of urban and rural sectors of the economy,measured by11Without migration costs,risk averse households would move independently from liquidity constraints.Lower migration costs should induce more migration indepen-dently of shocks if they fall below the utility costs of risk.The premium a household would be willing to pay to get rid of risk depends on the functional form of utility,but is positive for risk averse households.In long-run equilibrium when shocks are of mean size and no one has an incentive to move this means that the cost of migration is equal to the risk premium.12Since we do not observe rural versus urban fertility we can only include overall population growth.465S.Poelhekke /Journal of Development Economics 96(2011)461–475。