张玉君
（地质矿产部航空物探遥感中心研究所北京）
摘要：本文介绍一种独特的航放数据图像复原方法。该方法的主要技术关键是：提出航放数据图像复原原理和理论基础；建立航放数据图像复原处理流程；制定重建数据网格文件的途径；进行航放数据图像复原效果及误差评价。
关键词：航放数据，大气本底，图像处理，图像复原技术。
一、引言
自20世纪70年代初期方柱形NaI晶体进入机载综合航空站以来，航放测量的灵敏度和有效性显著提高，地质及地球物理勘探界对于航空放射性测量的需求产生了根本变化。
在近20年的应用实际中，大气氡本底（简称大气本底）的改正很难准确，它始终是困扰该方法应用效果的主要难点；其后果是在图面上造成条带现象，从而严重地影响着图件的可用性及方法的效果。其原因[1]可概括为：空中所测放射性不仅仅来源于地下，而且受飞机硬件环境、宇宙射线、大气中氡及其子体的影响。后者称之为大气本底干扰，它又受气候、风力、风向、温度、季节及一天中何时测量等因素的影响。大气本底干扰的主要表现形式是架次与架次间本底水平不同。受干扰最大的是铀道，钾道次之，钍道和总道虽较小，但也不可忽视，（见彩版附图7中的图3、图4）。由于这种噪声的存在，来源于地质体的信息常常淹没于噪声之中。图3a（彩版附图7）为哈密土墩测区K（红）、Th（绿）、U（蓝）三元素复原图像，图3b为该测区航放原始数据合成图像，图4a（彩版附图7）表示各架次早、晚校准读数，图4b为总道原始数据图像。条带噪声的存在，可以形象地比喻为挂在有用信息图像前面的彩条窗帘，条带的严重性使得该工区原始航放数据无法绘制等值线图。
航放图面条带问题是一个“世界性”的问题[2]。解决得较好的是加拿大，靠星罗棋布的湖泊河流等水域上空测量结果改正本底，取得较好的效果，而且不使用向上探头[1]。美国Geometrics公司及其他航空物探公司则是靠向上探头测量，作为本底改正的依据[3]。1986年Grasty[4]提出当测区内没有湖泊时，可用测线上无异常区的平均值代替本底。
本文介绍的方法与国际上已采用的各种方法全然不同，该方法在数字图像处理学中可称之为航放图像复原技术。图像复原技术的主要目的，是要改善给定的图像。复原是一个过程，它试图利用蜕化现象的某种先验知识，把已经蜕化了的图像加以重建或恢复。因此，复原技术是把蜕化模型化，并运用相反的过程在某种程度上恢复原来的图像。
Cannon博士[5]研究了一种图像复原技术或称图案去除技术，它适用于：从规则图案（如纺织品）上提取指纹图形，改善散焦图像，消除卫片图像探测器与探测器间的噪声，使在曝光过程中相机或物体平移造成模糊的图像清晰化等。Srinivasan也报道了这类研究[6]。张玉君等研究了深海锰结核照片光照不均匀等蜕化现象的图像复原问题[7]。航放数据图像复原处理是数字图像复原技术在地学界成功应用的又一实例，但航放数据图像所存在的蜕化问题与上述各例均不相同。该方法研究成功后，曾在6个测区得到验证。
二、航放数据图像复原技术原理和理论基础
航放所测到的是一幅蜕化了的图像G（x,y）它可视为由真实图像F（x,y）与干扰图像η（x,y）叠加而成，简化了的蜕化过程，见图1。航放图像蜕化现象的先验知识来自对航放测量过程及原始图像的分析。在测量过程中，来自地质体的有用信息是不随时间为转移的。而干扰在本质上是随时间变化的，但在图像上干扰已变为（x,y）的函数，因为：
张玉君地质勘查新方法研究论文集
图1 航放数据图像蜕化示意图
η的变化可分为架次之间的跳变及架次之内的渐变，见图4（彩版附图7），在每一测线上此干扰大致为一个常数，如果将x（即图像上的列）表示垂直测线的方向，则η（x,y）简化为η（x），则有
张玉君地质勘查新方法研究论文集
航放图像复原的目的，就在于设法近似地求出η（x），从而近似地得到F（x,y）。为此，沿测线方向对原始图像进行多次单列
多行窄长窗口褶积：
张玉君地质勘查新方法研究论文集
式中W为褶积模板，是由加权因子组成的矩阵。褶积过程是一种线性运算，其算子H不随空间变化。因为算子为线性的，则两个输入之和的响应等于两个响应之和。
张玉君地质勘查新方法研究论文集
由于假定了η仅与x有关，又由于褶积窗口为单列，则有：
张玉君地质勘查新方法研究论文集
现在分析HF（x,y）的性质，由于沿y方向的多次滑动平均，局部异常“淹没”于近区域特征之中，这种近区域特征表现为沿测线方向的低缓变化；如果用f（x,y）表示局部异常，用L（x,y）表示近区域场，于是：
张玉君地质勘查新方法研究论文集
再经如下处理
张玉君地质勘查新方法研究论文集
由（9）式可见，从原始图像中减去噪声图像后，所得复原图像f（x,y），它从局部异常角度是接近真实图像的，误差取决于所减去的“近区域背景值”在测线方向起伏变化的幅度。
三、航放数据图像复原处理过程
航放数据图像复原技术的研究是以多元统计为理论基础，以图像处理为工具来完成的，并体现了图像处理快速直观的特点，其流程如图2所示。
图2 航放数据图像复原处理流程
此方法假定航放噪声本底沿测线方向不变或呈线性变化。通过沿测线方向的多次滑动平均，使局部异常逐渐淹没在噪声本底之中，得出一幅与噪声本底线性相关的噪声图像。噪声图像尚需进行边缘影响补偿；对去除噪声后的图像，经中值滤波和空间变量反差增强，达到最终复原的效果。这一复原过程集中表示于图2的左半部。
图2的右半部为数据网格文件的重建过程，它是实际应用所必不可少的。经过分类分区，求得各类别在复原前后的均值向量，经最小二乘拟合求出复原图像的元素含量或计数率值，重新建立为在主计算机上绘制等值线图用的网格文件。
本研究曾试验通过沿测线方向取平均值做为噪声水平，结果不及上述方法理想。
四、效果及误差评价
1.航放数据图像复原的效果
（1）图面直观效果的改善。
可以形象地说，航放图像复原好比揭去一层条带窗帘，使原来透过此窗帘隐约可见的图像显示出了真面目，见图3a（彩版附图7）。图面直观效果的改善还表现在由于定位问题所引起的岩体边界上的锯齿状噪声得以消除，见图5（彩版附图7）。图5为总道对比图像，5a为原始数据，5b表示噪声图像，5c为去噪声后的图像，5d为复原图像。
（2）用复原数据所做等值线图真实可信。
以哈密土墩测区为例，原始资料由于条带干扰，在主计算机上，钾、钍、铀道都无法绘等值线图，仅提供了平剖图；只有总道提供了等值线图，但仍可看到条带的影响。
经图像复原、重建网格文件，反馈回主计算机后绘制了TC、K、Th和U等值线图，现以K道经复原后数据等值线图为例示于图6（彩版附图7），与地质图对比，表明异常和地质体对应良好，各类岩性的放射性趋势也都吻合，证实这些等值线图的可靠性。利用复原图像所做分类图也证实了这一点，见图7（彩版附图7），图7中数字分别为：①超基性岩；②基性岩；③花岗岩；④闪长岩；⑤变质岩；⑥混合岩；⑦第四纪沉积；⑧第三、四纪沉积；⑨第三纪沉积。
（3）有用信息增加。
本研究利用多元统计的方法，对航放图像复原的效果给出定量评价。可用一幅图像有用信息构成的变异值的大小来对它做定量评价。为此应计算全图面总变异对于一个象素的平均值，即平均变异值。用C、C´和C＂分别表示原始图像中有用信息平均变异值、原始图像中干扰信息平均变异值和最终复原图像有用信息平均变异值。统计时以G´（x,y）近似代表η（x）；以[G（x,y）-G´（x,y）]近似代表F（x,y）；以P（x,y）表示最终复原图像，并假定它已无干扰存在。
张玉君地质勘查新方法研究论文集
式中，字母上加“—”表示平均值；M、N为图像的行、列数。
表1为哈密土墩测区航放数据图像按上述各式所做定量评价的统计结果。
表1
从表1可见，K、Th、U、TC经图像复原后，有用信息都有十分显著的增长；就此工区而论，TC和K原始图像相对质量较好，Th和U较差。
2.复原图像的准确度及误差评价
复原图像的主要误差来源是“近区域背景值”L（x,y），它是在多次滑动平均时形成的。通过对干扰图像剖面数据的统计，得到以下准确度评价：
K±0.16%（绝对含量）；Th+2.1 pp m；
U±0.15 ppm；TC±869.6计数。
五、结论
（1）本文介绍的方法是在国内外首次提出的方法独特的航放数据图像的复原技术，并在多个工区验证了其可靠性和实用性。
（2）本技术可以基本上消除由于大气本底及阈值变动所造成的图面条带现象，基本复原航放图像的真面目，为进一步图像处理（诸如：求导、增强、分类、逻辑运算等）做了准备，因此本技术也是一种快速预处理方法。
（3）本方法改善了由于飞行往返定位位移所引起的某些地质体边缘呈锯齿状的图像噪声问题。
（4）本研究建立了“有用信息平均变异值”做为定量评价航放数据图像复原效果的尺度。还讨论了图像复原做为一种预处理过程，对于元素含量值可能导入的绝对误差或称为方法的准确度。
参考文献
[1]Grasty, R.L., Gamma ray spectrometric methods in uranium exploration—Thcory and operational procedures, Geophysics and Geochemistry in the Search for Metallic Ores,GSC,Ottawa,147-162,1977.
[2]Creen, A.A.,Leveling airborne gamma-radiation data using between-channel correlation information,Geophysics,52,1557-1562,1987.
[3]Foote, R.S., Improvement in airborne gamma-radiation data analysis by removal of environmental and pedologic radiation changes, in Sympos.on the Use of Nuclear Techniques in Prospecting and Developmcnt of Mineral Resources: Internat.Atomic Energy Agency Mtg.,Buenos Aires,187-196,1968.
[4]Grasty,R.L,Automated system for computing on-line atmospheric backgrounds,GSC paper,1-52,1987.
[5]Cannon,M.,Lehar,A.and Preston,F.,Background pattern removal by power spectral filtering,Applied Optics,22,777-779,1983.
[6]Srinivasan,R.,Software image restoration techniques,Digital Design, 16,4,27-34, 1986.
[7]张玉君，史鉴文.深海多金属结核照片的图像复原和图像处理技术研究.物探与化探，1989,（13）:435～441.
致谢 林振民同志对本文提出了宝贵的意见，史鉴文同志参加了重复工区试验，张志民和谢欣同志分别编制了网格文件转换和最小二乘拟合程序，杨星虹同志拍摄了屏幕图片，水恩海同志搜集了试验工区校准资料，在此一并致谢。
A STUDY ON IMAGE RESTORATION TECHNIQUES FOR AERORADIOMETRIC DATA
Zhang Yu jun
（Research Institute, Center of Aero-Gcophysics and Remote Sensing,Ministry of Geology and Mineral Resources, Beijing）
Abstract
This paper represents a specific methodfor restoration of images of airborne radiometric data.The main technical keys involved in this study are；the advancementof the principles and theory；the establishment of the flow-diagram for processing；the formulation of the means for reestablishment of the gridded data file；the evaluation of the restoration results and the errors, involved by the restoration processing.
Key words Aeroradiometric data, Atmospheric background, Image processing, Image restoration techniques.
原载《地球物理学报》，1990,Vol.33,No.4。
李立
（地质矿产部地球物理地球化学勘查研究所，廊坊 065000）
摘要 利用深大地电磁测深结果研究了地壳、上地幔的地电结构。30km深度的地电结构复杂，由很多大小不等电性各异的块段组成，大多数低阻区对应于高地温区。90km深度的地电结构相对简单，在高阻的背景上分布着一些低阻异常带，其中的松辽、华北地台东部以及扬子地台西部连成一片，构成了北东-南西向的巨大低阻异常带，大兴安岭-太行山重力梯度带、汾渭地堑、龙门山-攀西构造带位于该低阻带西侧，依兰-伊通断裂带、郯庐断裂带位于其东侧，推测该低阻异常带与上述断裂带的形成有关，并有部分熔融的地幔物质存在。150km深处的电性结构则以低阻为背景，在其上零星分布着高阻块体，其轴向线构成的“构造线”将中国大陆分割成若干地质构造单元。在新生代裂谷区和近代活动构造区大多存在壳内低阻层，且深度浅，浅源地震的震源一般分布在壳内低阻层以上。上地幔低阻层的深度变化较大，最浅处仅50～60km，而最深可达200km以上，东部地区平均约100km，西部地区约120km。从东北的松辽地区经华北地台东缘向西南沿汾渭地堑再折向南直到扬子地台西缘有一个北东-南西向的巨大的上地幔低阻层隆起带，大部分上地幔低阻层隆起区与高热流区以及强震震中分布带有较好的对应关系。
关键词 大地电磁测深 地壳 上地幔 电阻率 中国大陆
1 引言
本项研究的数据是基于地质矿产部系统十余年的大地电磁测深（MTS）深部调查结果以及国家地震局、中国科学院、地质院校等单位多年的MTS结果[1～24]。搜集并整理了近1000个测深点的资料得出了30km、90km、150km深度的电阻率值，壳内低阻层以及上地幔低阻层的深度值，并将这些数据按1°×1°的经纬网格分别进行平均，以平均值代表每个网格的电阻率值与低阻层的深度。由于深大地电磁测深的测点分布不均匀，尚有不少空白区，有些空白区的壳内与上地幔内低阻层的深度是用大地热流值估算的[25]，因此只能粗略地给出中国大陆深部的地电结构与格架。
2 不同深度上的地电特征
根据30km、90km以及150km深度的电阻率研究了地壳、上地幔的地电结构，它们基本上反映了地壳内、岩石圈底部以及岩石圈以下的电性结构。
2.1 30km深度
在30km深度上电阻率的变化无规律，由很多大小不等，电阻率各异的块段组成，说明了地壳内构造复杂。根据岩石电阻率的高温高压实验结果表明，地壳内的低阻异常大多是由含流体或石墨化的地层引起[26,27]。多数的低阻异常带与高地温区对应。
2.2 90km深度
图1 中国大陆90km深度上的地电结构及低阻异常
①松辽地块；②华北地台东部；③扬子地台西部与地台东缘；④华南褶皱系东部；⑤腾冲褶皱带；⑥青藏高原；⑦祁连褶皱带；⑧鄂尔多斯西缘的银川地堑；⑨天山-北山褶皱带；⑩阿尔泰褶皱带；⑪班公错-嘎尔地区。图中等值线为松辽—扬子正磁异常[28]
90km深度上的地电结构较为简单，在高阻的背景上分布着一些低阻异常带（图1），它们是松辽地块、华北地台东部、扬子地台西部与该地台东缘、华南褶皱系东部、腾冲褶皱带、青藏高原、祁连褶皱带、鄂尔多斯西缘的银川地堑、天山-北山褶皱带、阿尔泰褶皱带以及西藏西缘的班公错—嘎尔地区。其中的松辽地块、华北地台东部以及扬子地台西部连成一片，构成了一个北东-南西向的巨大低阻带，大兴安岭-太行山重力梯度带、汾渭地堑、龙门山-攀西构造带位于该低阻带西侧；依兰-伊通断裂带、郯庐断裂带位于东侧，推测上述断裂构造带的形成与该低阻异常带有关。这一北东-南西向的低阻带与安振昌编制的卫星磁异常图[28]中的松辽-扬子正异常带基本吻合，推测该低阻带是上地幔的隆起区。
90km深度上的低阻异常大多分布在岩石圈的减薄区，那里可能存在部分熔融的地幔物质。
2.3 150km深度
150km深度在低阻的背景上零星散布着一些高阻块（图2），大部分地区的地幔物质已处于部分熔融状态，只有松辽地块北部、冀鲁陆核东部、扬子地台与华南褶皱系的接合处，右江褶皱带南部、龙门山-攀西构造带、柴达木地块、祁连山褶皱带、阿尔泰-天山地区以及西昆仑山表现为局部的高阻区。推测上述高阻区是由难熔的残余地幔物质形成的“硬块”。这些高阻块体有可能是深部构造带的标志，如果将一些高阻块体的轴向线分别连接，则可得出6条有趣的“构造线”F1—F6。自西昆仑山经柴达木、祁连山及其东延部分直到冀鲁陆核南端构成近东西向的“构造线”F1，它将中国大陆分割成南北两大部分，准噶尔、塔里木、华北地台位于其北侧，青藏高原、扬子地台及华南褶皱系位于南侧。松辽地块北部的轴向线向南延伸则与郯庐断裂带相连构成北东-南西向的“构造线”F2。扬子地台与华南褶皱系间的高阻块体轴向线向西南延伸与右江褶皱带的高阻轴同线相连构成“构造线”F3，它可能是扬子地台与华南褶皱系的深部边界。除此尚有位于龙门山-攀西的“构造线”F4以及天山-阿尔泰“构造线”F5，它们可能分别反映青藏高原与扬子地台的深部边界以及准噶尔地块的深部边界。
图2 中国大陆150km深度上的地电结构及低阻背景上的高阻块体
①松辽地块北部；②冀鲁陆核东部；③扬子地台与华南褶皱系的接合处；④右江褶皱带南部；⑤龙门山-攀西构造带；⑥柴达木地块；⑦祁连褶皱系；⑧阿尔泰-天山地区；⑨西昆仑山；F1—F6为高阻块体轴向线构成的“构造线”
3 壳内低阻层的深度变化
在中国大陆上除华南褶皱系、扬子地台及松辽地块的部分地区外，大多发现了壳内低阻层。壳内低阻层的深度一般为15～30km，厚度几到十几公里，电阻率几到几十欧姆米。大多数的壳内低阻层与活动构造区对应，且地温较高（表1）。此外，壳内低阻层的上隆区还与莫霍面的上隆区（图3）以及上地幔低阻层的上隆区（图4）有一定的对应关系。多数浅源地震的震源位于壳内低阻层以上（图5）。
表1 壳内低阻层上隆区
图3 壳内低阻层上隆区与莫霍面上隆区的分布
①松辽地块北部；②华北地台东缘；③汾渭地堑；④二连-呼和浩特-东胜；⑤南北地震带北段；⑥祁连山；⑦松潘-甘孜-康滇地区；⑧雅鲁藏布江-腾冲地区；⑨天山-准噶尔地区；⑩北山；⑪下扬子地区；⑫泉州－赣州地区；⑬洞庭盆地；⑭建始-恩施盆地；⑮西昆仑北部；⑯班公错-嘎尔地区
图4 汾渭地堑壳内低阻层与上地幔内低阻层的对应关系[5]
图5 拉萨与唐山地区壳内低阻层与浅源地震震源分布[29]
4 上地幔内低阻层的深度变化
中国大陆上地幔低阻层的顶面深度变化很大，从最浅的50～60km到最深的200km以上，平均深度为100～120km，一般情况下，中国的东部地区浅，西部地区深（图6）。从东北的松辽地区经华北地台东缘，向西南沿汾渭地堑再折向南直到扬子地台西缘有一个北东－南西向的巨大的上地幔低阻层隆起带，在90km深度的地电结构图上也清楚地显示了一个巨大异常带。
一般情况下，上地幔低阻层顶界面是岩石圈的底界面，因此根据上地幔低阻层的深度粗略地给出了岩石圈的厚度，并对中国大陆的岩石圈进行了分区，岩石圈厚度小于100km的地区定为岩石圈减薄区，岩石圈厚度大于120km的地区定为岩石圈的增厚区。表2给出了中国大陆岩石圈的厚度。
岩石圈的减薄区大多对应于高热流区（图7），多数的强地震分布在岩石圈减薄区。这些都说明岩石圈减薄区具有活动构造带的特征。有资料[30,31]表明大型内生金属矿多数分布在岩石圈的减薄区或岩石圈厚度的陡变带上。因此，上地幔低阻层的起伏变化对矿产预测具有重要意义。
图6 中国大陆上地幔低阻层的深度（单位：km）
表2 中国大陆岩石圈的厚度
5 结论
在90km深的上地幔内发现一个自松辽盆地直到扬子地台西部的北东-南西向巨大低阻异常带，推测该低阻异常带是上地幔的隆起带。
中国大陆自90km到150km深的上地幔是由塑性的“软体”和致密的“硬块”组成，而不是简单的层状结构。中国大陆150km深度的大部分地区，地幔物质已处于部分熔融状态，反映为低电阻率，只有局部地区为高电阻率，它们构成的一些“构造线”有可能是大地构造单元的深部边界。
图7 岩石圈减薄区与高热流区（金昕）的分布
①松辽-华北地区；②扬子地台东缘-华南褶皱系东部；③鄂尔多斯西缘；④汾渭地堑；⑤康滇隆起区；⑥腾冲褶皱带；⑦金沙江-玉树地区；⑧东昆仑；⑨祁连褶皱带；⑩天山褶皱带；⑪二连-东乌珠穆沁带
在多数的近代活动构造区内均发现了壳内低阻层，且埋藏度浅，一般小于20～25km。
在上地幔内发现一个自松辽至扬子地台西缘的北东-南西向的巨大的低阻层隆起带。
中国大陆岩石圈的平均厚度为100～120km，东部小于西部。岩石圈的减薄区大多位于活动构造区内，并与高热流区以及强地震区对应。
参考文献
[1] 袁学诚，李立，金国元等．西藏洛扎—羊八井地区的磁大地电流地壳测深．地质学报，1985,59（1）:25～32.
[2] 李立，金国元．西藏地区岩石圈的电性与热状态．地质矿产部物化探研究所所刊，（2）,121～128．北京：地质出版社，1987.
[3] 李立，金国元．攀西裂谷及龙门山断裂带的大地电磁测深研究．物探与化探，1987,11（3）:161～168.
[4] 陈沪生．下扬子地区HQ—13线的综合地球物理调查及其地质意义．石油与天然气地质，1988,9（3）:211～222.
[5] 邢集善，姚典群，黎明．试从地球物理资料论山西地堑系的构造特征．山西地质，1989,4（2）:95～108.
[6] 蒋宏堪，战双庆，王宏勋．十堰至洛阳大地电磁观测结果．物探与化探，1990,14（4）:285～291.
[7] 蒋宏堪，战双庆，王宏勋．湖北麻城—九宫山剖面深部地质特征初探．物探与化探，1990,14（5）:357～364.
[8] 詹麒．塔里木盆地东北部的大地电磁测深工作及其初步地质成果．地球科学，1990,15（增刊）：97～106.
[9] 罗志琼．开封—灵壁大地电磁测深研究．地球科学，1990,15（增刊）：87～95．
[10] 张胜业，欧阳昌定．大地电磁测深在鄂中地区区域地电结构探测中的应用．地球科学，1990,15（增刊）：79～85.
[11] 李立，金国元，刘畅往等．可可托海—阿克塞剖面的地电特征．见：中国地球物理学会年刊．北京：地震出版社，1992.
[12] 蒋宏堪，战双庆，王宏勋．四川大足—福建泉州深部地电特征．地球物理学报，1992,35∶214～222.
[13] 顾群，孙洁，史书林等．华北、西北一些地区地壳和上地幔内高导层．地震地质，1980,2（2）:21～29.
[14] 国家地震局《深部物探成果》编写组．中国地壳上地幔地球物理探测成果．北京：地震出版社，1986,271～285.
[15] 吴广跃，曾陆海．大地电磁法在福建漳州地热区的应用．地球科学，1988,13（3）:325～333.
[16] 张云琳，安海静，刘晓玲等．我国西北部分地区地壳—上地幔电性横向变化特征与地震活动的关系．地震地质，1988,10（2）:65～73.
[17] 孙洁，徐常芳，江钊等．滇西地区地壳上地幔电性结构与地壳构造活动的关系．地震地质，1989,11（1）:35～45.
[18] 高文，蒋邦本，白登海．邢台地震区大地电磁观测与研究．地球物理学报，1990,33:291～297.
[19] 邓前辉，刘国栋，刘金汉等．湖北襄樊—福建罗源的大地电磁测量与地壳上地幔电性特征研究．地震地质，1990,12（2）:149～157.
[20] 史书林，徐常芳，王继军等．辽宁义县—内蒙古乐乌珠穆沁旗剖面深部电性研究．地震地质，1991,13（2）:115～125.
[21] 孔祥儒，刘士杰，张建军等．福建东部地区大地电磁测深研究．地球物理学报，1991,34:724～735.
[22] 白登海，张丽，孔祥儒．内蒙古东部古生代块体碰撞区的大地电磁测深研究．地球物理学报，1993,36:773～783.
[23] 秦国卿，陈九辉，刘大建等．昆仑山脉和喀喇昆仑山脉的地壳上地幔电性结构特征．地球物理学报，1994,32:193～198.
[24] Белявский，В．В．．ГЕОЭЛЕКТРИЧЕСКАЯМОдЕЛЬТЯНЬ-ШАНЯ．Физика Земли，1995，（1）：3～13.
[25] Adam,A..Geothermal effects in the formation of electrically conducting zones and temperature distribution in the Earth.Phys.Earth Plan.Int.,1978,17:21～28.
[26] Ваныян，Л．Л．．Мекдународный симиозиум“Природа злектропроводности земнойй коры”．физика Земли，1995,（1）:95～96.
[27] Жамалетдинов，А．А．．Модель Электропроводности Литосферы по Результатам Исследовании с контролируемыми источниками поля．Ленинград：Наука．1990,125～147.
[28] 安振昌，马石庄，谭东海．中国及邻近地区卫星磁异常的球冠谐和分析．地球物理学报，1992,35（增刊）：188～201.
[29] 马杏垣．中国岩石圈动力学图集．北京：中国地图出版社，1989.
[30] Щеглов，А．Д．．Эндогенная металлогения и тектоносфера．Геотектоника，1990,（5）:9～16.
[31] 郭文魁，刘梦庚．中国内生金屑成矿图．北京：中国地图出版社，1989.
[32] 黄汲清，任纪舜，姜春发．中国大地构造及其演化．北京：科学出版社，1981,29～65．
尹祥础
（中国地震局分析预报中心，北京 100036）
王裕仓
（中国科学院力学研究所非线性连续介质力学开放实验室（LNM），北京 100080）
摘要 由于损伤过程的不可逆性，当孕震区介质受到损伤后，其对加载的响应将不同于卸载响应。加载响应与卸载响应的比Y（称之为加卸载响应比，英文为Load/Unload Response Ratio简称LURR）可以度量孕震区介质的损伤程度或接近失稳的程度，因而可以作为一种地震预测的新途径。对数百个地震震例的检验（震级从4级到8.6级）表明：在主震发生前的一段时间里，Y值显著大于1。而对于7个稳定区（指在较长时期内未发生过强震，而小震资料又较丰富的地区），在长达20年的时间内，Y值始终在1附近作轻微的起伏。近年来，利用本方法对发生在中国大陆的十几次中强地震及美国北岭地震（1994.1.17,Mw=6.7）与日本关东地区地震（1996.09.11,Ms=6.6）作出了成功的中期预测。
关键词 地震预报 加卸载响应比理论
1 引言
地震的物理实质是什么？从力学的观点看，它正是震源区介质（岩体）的损伤与快速破坏（失稳）过程，并伴随应力与应变能的快速释放。让我们研究孕震区（含断层或弱化区的岩体）在高温高压下的本构关系，如图1所示。图1中纵坐标为广义载荷P，而横坐标为对于载荷P的响应R。
我们首先定义如下两个参数：响应率X与加卸载响应比Y。
响应率X定义为：
第30届国际地质大会论文集 第5卷 现代岩石圈运动 地震地质
式中：△P表示载荷增量，而△R表示相当于△P的响应增量。
加卸载响应比Y定义为：
第30届国际地质大会论文集 第5卷 现代岩石圈运动 地震地质
式中：X+、X-分别表示在加载及卸载条件下的响应率。
众所周知，材料若处于弹性阶段（图1中的OA段），加载时（△P＞0）的响应率X﹢，将等于卸载时（△P＜0）的响应率，即Y=1。但是若应力超过弹性后，X+＞X-，因而Y＞1。当材料逐步趋向破裂时，Y值也随之逐渐增大。当趋近于图1中的顶点T时，X+趋于无限大；而X-仍保持为有限值，因而Y值也将趋于无穷。所以顶点T可以作为预测失稳的前兆点。
图1 震源区的本构关系
从损伤力学的观点看，地震孕育过程就是孕震区介质的损伤过程。因此有希望采用损伤力学中的损伤参数D来定量地刻划地震的孕育进程。损伤D有多种方法定义。最直接的一种是用弹性模量（4阶张量）的变化来定义D。为简单起见，有时只用弹性模量张量的一个分量来定义D。例如Lamaitre[23]定义D为：
第30届国际地质大会论文集 第5卷 现代岩石圈运动 地震地质
式中：Eo为未损伤材料的杨氏模量，E为已损伤材料的杨氏模量。如果卸载时的模量为Eo，则（3）式可表示为：
第30届国际地质大会论文集 第5卷 现代岩石圈运动 地震地质
式（4）意味着D与Y之间存在着密切的内在关系。也就是说，Y也可以作为孕震区损伤程度的度量。
即使损伤D不按式（3）而另行定义，D与Y之间仍然会存在内在关系。这就明参数Y可以定量地刻划地震孕育进程，因而可以作为地震预报的定量前兆[11～18,24～27]。
要用加卸载响应比理论进行地震预报，必须首先解决几个主要问题。一是如何对地球进行加载与卸载，以及如何判别加载与卸载。其次是怎样选择合适的参数作为响应。以下分别讨论这几个问题。
（1）如何对地壳块体进行加卸载？
孕震区的线性尺度可以达到几百甚至上千千米。对其进行加卸载的方法之一是利用潮汐应力。潮汐力不断地周期性地变化。也就是说它对地球的各部分不停地进行加载与卸载。
（2）用什么准则判定加载与卸载？
加载与卸载问题，在塑性力学中有详细的讨论。对于不同的材料应选择不同的准则。对塑性较好的多种金属（如低中碳钢、铝等）Von Mesis准则比较适用；而对地质材料的破坏，则Coulumb准则[22]更适合。我们在文献[12,13,24,26]中对此作了详细的研究与论述。请参阅上述文献，本文不再赘述。
（3）选择什么参数作为计算Y值的响应
地壳形变、井水水位、地震活动性及其他震源参数以及许多其他地球物理参数都可以作为响应，用于计算Y值。从“八五”期间起，我们与国内外许多地球物理学家合作，开展了多学科的研究，国家科委、国家地震局的“八五”、“九五”攻关项目中，均安排了相应的项目，同时还得到了国家自然基金会及地震科学基金会的支持，取得了比较多的成果[1,4～6,10,19]。在本文中主要介绍以地震能量为响应的加卸载响应YE（在本文中很少涉及其他参数的加卸载响应比，所以仍以Y代替YE）。Y定义为：
第30届国际地质大会论文集 第5卷 现代岩石圈运动 地震地质
式中：E为地震能量，符号“＋”“－”分别表示加载与卸载，N+、N-分别表示在规定的时空范围内发生的加载地震及卸载地震（在加载时段内及卸载时段内发生的地震，在以往的文献中，有时也简称为正地震及负地震）的数目，m为常数，取为0、1/3、1/2、2/3及1。m=1时，Y为所有正地震能量之和与负地震能量之比；m=1/2,Em为贝尼奥夫应变；m=1/3及2/3时，Em分别表示震源的线尺度及面尺度；当m=0时，Y=N﹢/N﹣，即正、负地震数1之比。在本文中恒取m=1/2。
取一定的时、空范围（例如2°×2°，几个月至1年），按式（5），计算出一个Y值。利用该区域Y值随时间的变化，可能预测出该区域内未来发生地震的危险程度。
2 震例检验
我们用国内外数百个已发生的地震，对LURR理论进行了广泛的震例检验，震级范围从Ms=4到Ms=8.6。检验的结果是满意的[13,16,24,26]。以下是1970～1992年期间发生在中国大陆的10个大地震（Ms≥7）的孕震区LURR随时间的变化情况（图2）。
在这段时间内，中国大陆共发生Ms>7级以上地震13个，其中有3个地震（青海地震、通海地震及西藏亦基台错地震）因数据太少无法利用它们计算LURR之外，对其他10个大震全部进行了研究，地震前各孕震区Y随时间的变化示于图2，山图可知，10个震例中，有9个震例，在震前Y明显大于1,Y>1的时间大约持续1～3年。
除了较系统地研究了中国大陆的震例（4≤Ms≤8.6）外，根据我们所能得到的资料还研究了日本、美国和其他国家的若干震例。都取得了好的结果[13,18,26]。
此外，我们选择了中国大陆的7个区域作对比研究，这7个区域在历史上曾发生过强震，但近20多年来，地震活动性较低，没有发生过中强地震，处于低地震活动期。其LURR的变化情况（图3）与图2形成强烈的对比，在所有这7个区域里，在20多年（1970～1992）的时间里，其Y值均在1附近作轻微起伏。
图2 1970～1992年中国大陆所有7级以上大地震前Y随时间的变化曲线
a—1973.02.06四川炉霍地震（Ms=7.6）;b—1974.05.11云南永善地震（Ms=7.1）;c—1974.08.11新疆乌恰地震（Ms=7.3）;d—1975.02.04辽宁海城地震（Ms=7.3）;e—1976.05.29云南龙陵地震（Ms=7.4）;f—1976.07.28河北唐山地震（Ms=7.8）;g—1976.08.16四川松潘地震（Ms=7.2）；h—1985.08.23新疆乌恰地震（Ms=7.1）;i—1988.11.06云南澜沧地震（Ms=7.6）;j—1990.04.26青海共和地震（Ms=7.0）
图3 1970～1992年间，7个平静区的Y-t曲线
a—郯庐断裂带南段（35.5°N士1°，118°E±1°）;b—陕北（40.5°N±1°，119°E±1°）;c—川东（31.0°N±1°，118E士1°）；d—鲁北（37.0°N±1°，119°E±1°）;e—鲁西（37.0°N士1°，118°E±1°）;f—豫北（35.0°N±1°，113°E±1°）;g—鲁南（35.0°N士1°，117°E±1°）
3 地震预报实践
近年来我们尝试用本方法进行地震预报，多次成功对发生在中国大陆的6级以上地震成功地作出了中期预报[5,7,15～18]。此外，还成功地预报了1994年1月17日的美国北岭地震[5.24]及1996年9月11日的日本千叶地震[18]。部分被预测的地震震前的Y-t曲线示于图4。
以下对其中几个典型地震的预测情况略作说明。
1993年夏初，我们得到由USGS所属NEIC（美国国家地震信息中心）供给的加利福尼亚州地震目录，利用此目录研究了该州的加卸载响应比，发现其中有3个地区在较长的时间（长于一年）内Y值显著大于1。经研究后，于1993年10月28日写信给提供我们数据的那位科学家（ISOP项目的负责人），在信中提供了分区的加卸载响应比结果，并且预测：在其中3个区域或其附近，在一年内（1993.10.28～1994.10.28）可能发生中强地震（7>M≥6）。在预测后不到3个月，1994年1月17日发生北岭地震（图4e），发生在预测的一个地区的边缘。再后，1994年9月12日，在另两个地区的中间发生一个Ms=6.0级地震。
图4 用LURR理论成功地预测的某些中外地震的震前Y-t曲线
a—1991.03.26山西大同地震（Ms=6.1）;b—1993.01.27云南普洱地震（Ms=6.3）;c—1993.10.26青海共和地震（Ms=6.0）;d—1994.01.03青海共和地震（Ms=6.0）;e—1994.01.17美国北岭地震（Mw=6.7）;f—1996.09.11日本千叶地震（Ms=6.6）
1996年春，应日本气象厅科学家的要求对日本的关东地区，和歌山地区及兵库地区的加卸载响应比进行了分析研究（资料由对方提供），研究后得到几点结果：①关东地区（按对方提供资料的范围为了35°～36°N,139°～141°E）在一年内发生Ms=6级地震可能性很大；②和歌山地区在近期内不会发生中强以上地震（对方原来预计此地区危险性很高）；③1995.01.17神户地震前，加卸载响应比异常很显著。我们于1996年4月1日将有关结果传真给了对方科学家，同时写成论文[18]于1996年5月间投《中国地震》（季刊）。该文于《中国地震》中文版1996年第三期（1996年9月初出版）及英文版（由美国Allerton出版公司出版）第四期刊出。其后，在1996年9月11日发出了Ms=6.6级千叶地震（35.5°N,140.9°E）。
关于国内的地震预测只讨论一个震例——1994年12月31日的广西北部湾地震（Ms=6.1）。
1993年底我们在“1994年中国大陆地震趋势研究”的报道中将广西沿海地区列为地震危险区[15]。直至于1994年11月分析预报中心召开会商会时，该区未发生过任何中强地震，但当我们研究1992.09.01～1994.08.31时段中国大陆LURR的空间分布时，发现该区域的Y值异常仍非常突出[16]，因此我们认为该区域仍可能在年内发生强烈地震，结果在1994年最后一天发生了北部湾Ms=6.1级地震，并于1995年1月10日再次发生Ms=6.2级地震。
以上震例是成功的例子，但也有些Y值较高的区域，在预测的时段内并未发生强烈地震。这些区域在一定时段内加卸载响应比升高，说明该区域的地震孕育过程正在进行，但随后却可能发生了卸载使地震孕育过程推迟甚至中断，对于这种情况，如何判别是以后要着力研究的课题。
4 加卸载响应比的时空演化特征
大量的震例研究表明，LURR的空间分布图像是很复杂的。当一个地区未来要发生强烈地震前，将出现一系列高Y值区，这些高Y值区往往连成一个环状，形成面包圈图像[7]，大部分未来的地震将会在发生面包圈内或其邻域。图5是1979.03.14云南普洱Ms=7.0级地震震前一年间该区域内Y值的空间分布。由图可见，高Y值区围绕着未来的震中形成一个面积约为5°×5°的面包圈。将LURR的空间分布作成立体图，每一个高Y值区形成一个高峰，很多个高Y值区就形成群峰突起的图像。形成鲜明对比的是，在地震活动性低（指未发生强震）的区域内，Y值在1°上下轻微起伏，所以LURR的空间分布立体图就像平原地区的地形，我们形象地称这种图像为“一马平川”。
图5 1979.03.14云南普洱Ms=7.0级地震震前一年时段内，Y=2.0的等值线图图中符号代表未来的震中
对于同一地区，在地震孕育过程中，不同时段的LURR空间分布图像是不同的，也就是说，空间分布图像随时间发生变化。我们发现一个非常有趣的现象：即强震前多个高Y区向未来的震中迁移，称之为高Y值区的会聚现象[8]。
研究了1970年后中国大陆的12个Ms≥7.0大震[8]。12个震例中，有11个发现了会聚现象，且未来的震中处于面包圈内，只有1992.04.23中缅边界上的Ms=7.0级地震，未来的震中处于面包圈外，距圈的外边界约50km。这可能与该地震发生在两国边境地区，缅甸一侧的数据不好收集有关。
进一步，我们还发现：强震发生前，高Y值区迁移速度在同一地区，近似不变，大致为100km/a的量级。但不同地区的迁移速度有所差异[8]。
Scholz曾经撰文说，华北地区的形变锋（deformation front）传播速度约为150km/a,Press和Allen则观测到美国南加州地区的形变波（deformotion wave）速度为100km/a。这三者在物理上是有关系的，而且其数值在数量级上也是彼此相符的。
5 展望
前已述及，除了地震能量外，其他许多有关的地球物理参数（如地下水位、地壳体应变、地倾斜、地磁参数、尾波Q…）均可选择为响应，进行加卸载响应比的研究[1,2,4,10]。图6为取尾波Q值的例数作为响应，计算LURR—YQ的例子。图6所示为美国加利福尼亚州南部北岭地震（1994.01.17,Mw=6.6）前该区YQ随时间的变化情况，将它与图4e作比较后可以看出，二者在定性上是一致的。
图6 美国加利福尼亚州南部北岭地震的YQ变化图
对同一时空域用众多的参数可以计算出众多的LURR值，然后进行综合预报，必然会提高用加卸载响应比理论进行地震预报的精度。
简而言之，加卸载响应比理论可能为地震预报开辟了一条新途径。现在国内地震界有不少人在研究，应用与改善它[5]。
近来通过研究，表明北京地区的有感地震（指6＞M≥4）前[9]甚至北京地区的矿震（M＞2）前，Y值也有较明显的升高。这说明，除天然地震外，对于诱发地震（矿震、水库地震[4]…）以及某些其他地震灾害（如岩爆、滑坡、火山喷发…），也可能用LURR理论进行预测。
致谢 谨向傅承义、Keiti Aki、秦馨菱、王仁、陈章立、何永年、葛治州、陈鑫连、梅世蓉、罗灼礼、张国民、李宣瑚、张伯民教授及E.A.Bergman、K.Hosono、H.P.Ouyang博士致以诚挚的谢意，感谢他们多方面的支持与帮助。
本项目得到国家自然基金会（批准号 19732006）、国家科委及国家地震局“八五”及“九五”攻关项目、地震学基金会以及中国科学院LNM开放实验室（Lab of Nonlinear Mechanics of Continuous Media,Institute of Mechanics）的资助。
参考文献
[1]陈建民，张昭栋，杨林章，石荣会，张继红.地下水位固体潮响应比的地震异常.地震，1994,10（1）:73～78.
[2]陈学忠，尹祥础.水库地震主震前加卸载响应比的变化特征.中国地震，1994,10（4）:361～367.
[3]陈学忠，尹祥础，K Ake,H Ouyang，宋志平，王裕仓.以介质参数波Q1作为响应的加卸载响应比研究.中国地震，1996,12（3）:243～249.
[4]国家地震局分析预报中心.关于1994年我国地震趋势的研究报告.见：国家地震局分析预报中心编.中国地震趋势预测研究（1994年度）.北京：地震出版社，1994,43.
[5]李宜瑚.“八五”地震预报理论及方法攻关新进展之一：加卸载响应比理论预测洛杉矶地震获得成功.国际地震动态，1994,10（4）:24～25.
[6]刘桂萍，马丽，尹祥础.首都圈地区中等地震前响应比特征的研究.地震，1994,（6）:34～39.
[7]宋志平，尹祥础，陈学忠.加卸载响应比的时空演变特征及其对地震三要素的预测意义.地震学报，1996,18（2）:179～186.
[8]施行觉，许和明，万永中，卢振刚，陈学忠.模拟引潮力作用下的岩石破裂特征：加卸载响应比理论的实验研究之一。地球物理学报，1994,37（5）:633～637.
[9]王丹文.加卸载响应比理论在以磁报震中的应用探索.地震地磁观测与研究，1995,（16）:26～29.
[10]杨林章，何世海，郗钦文，黎凯武，李松阳.用体应变潮汐加卸载响应比确定岩石弹性的变化.中国地震，1994（增刊）：90～94.
[11]尹祥础.地震预测新途径的探索.中国地震，1995,3（1）:1～7.
[12]尹祥础，尹灿.非线性系统的失稳前兆与地震预报.中国科学，1991,（5）:512～518.
[13]尹祥础，陈学忠，宋治平，尹灿.加卸载响应比理论（LURR）：一种新的地震预报方法.地球物理学报，1994.37（6）:767～775.
[14]尹祥础，陈学忠.加卸载响应比理论及其地震预测中的应用研究进展.地球物理学报，1994,37（增刊）：223.
[15]尹祥础，陈学忠，宋志平.加卸载响应比理论的新进展及其在地震趋势研究中的应用.见：国家地震局分析预报中心编.中国地震趋势预测研究（1994年度）.北京：地震出版社.1993.
[16]尹祥础，陈学忠，宋志平.从加卸载响应比的时空分布，研究中国大陆未来的地震大形势.国家地震局分析预报中心编：中国地震趋势预测研究（1995年）.北京：地震出版社，1994.
[17]尹祥础，陈学忠，宋志平，王裕仓.由加卸载响应比的时空变化预测中国大陆地震趋势.见：国家地震局分析预报中心编.中国地震趋势预测研究（1996年）.北京：地震出版社.1995,75～178.
[18]尹祥础，陈学忠，宋治平，王裕仓.关东地区加载响应比的时间变化及其预测意义.中国地震，1996,12（3）:331～335.
[19]张继红.响应比法的地磁异常分析.地震地磁观测与研究，1995,16:61～63.
[20]B.K.Atkinson.The Theory of Subcritical Crack Growth with Application to Minerals and Rocks.In:Fracture Mechamics of Rock（Ed.B.K.Atkinson）,Academic Press,London,1987.
[21]Ding Zhongyi et al..Seismic triggering effect of tidal stress.Tectonophysics,1994,93:319～335.
[22]J.C.Jaeger,N.G.W.Cook.Fundamentals of Rock Mechanics.Chapman and Hall,London,1976,78～99.
[23]J.Lemaitre.How to Use Damage Mechanics.Nuclear Eng.＆.Design,1984,80:233～245.
[24]Yin Xiangchu,Chen Xuezhong,Song Zhiping and Yin Can.The Load-Unload Response Ratio Theory and its Application to Earthquake Prediction.Journal of Earthquake Prediction Research,1994,3:325～333.
[25]Yin Xiangchu,Yin Can and Chen Xuezhong.The Precursor of Instability for Nonlinear Systen and Its Application to Earthquake Prediction——the load-Unload Response Ratio Theory.Non-linear dynamics and predictability of geophysical phenomena（W.I.Newmana,A.M.Gabrelov and D.L.Turcote eds.）.AGU Geophysical Monograph,1994,83:55～60.
[26]Yin Xiangchu,Chen Xuezhong,Song Zhiping and Yin Can.A New Approach to Earthquake Prediction——The Load/Unload Response Ratio（LURR）Theory.PAGEOPH,1995,145（3/4）:701～715.
[27]Иин Ксянчу（Yin Xiangchu）.Новыйнодходк Ирогногу Землетрясеннй.НРИРОД А，1993,（1）:21～27.
